返回3D scipy.spatial.Delaunay表面三角形

Question

我有一个问题。我试图通过scipy.spatial.Delaunay对点云进行三角剖分。我使用了：

tri = Delaunay(points) # points: np.array() of 3d points 
indices = tri.simplices
vertices = points[indices]

但是，这段代码返回的是四面体。如何使其仅返回表面上的三角形?

谢谢。

- s.t.e.a.l.t.h

Answer 1

要使其以代码形式工作，您需要将表面参数化为2D。例如，在球体的情况下（r，theta，psi），半径是恒定的（删除它），并且点由（theta，psi）给出，这是2D。

Scipy Delaunay是N维三角剖分，因此如果您提供3D点，则返回3D对象。给它2D点，它就会返回2D对象。

以下是我用于创建openSCAD多面体的脚本。U和V是我的参数化（x和y），这些是我提供给Delaunay的坐标。请注意，现在“Delaunay三角剖分属性”仅适用于u，v坐标（角度在uv空间而不是xyz空间中最大化等）。

该示例是从http://matplotlib.org/1.3.1/mpl_toolkits/mplot3d/tutorial.html修改而来，该示例最初使用 Triangulation 函数（最终映射到Delaunay？）。

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.tri as mtri
from scipy.spatial import Delaunay

# u, v are parameterisation variables
u = np.array([0,0,0.5,1,1]) 
v = np.array([0,1,0.5,0,1]) 

x = u
y = v
z = np.array([0,0,1,0,0])

# Triangulate parameter space to determine the triangles
#tri = mtri.Triangulation(u, v)
tri = Delaunay(np.array([u,v]).T)

print 'polyhedron(faces = ['
#for vert in tri.triangles:
for vert in tri.simplices:
    print '[%d,%d,%d],' % (vert[0],vert[1],vert[2]),
print '], points = ['
for i in range(x.shape[0]):
    print '[%f,%f,%f],' % (x[i], y[i], z[i]),
print ']);'


fig = plt.figure()
ax = fig.add_subplot(1, 1, 1, projection='3d')

# The triangles in parameter space determine which x, y, z points are
# connected by an edge
#ax.plot_trisurf(x, y, z, triangles=tri.triangles, cmap=plt.cm.Spectral)
ax.plot_trisurf(x, y, z, triangles=tri.simplices, cmap=plt.cm.Spectral)


plt.show()

以下是（稍微更加结构化的）文本输出：

polyhedron(
    faces = [[2,1,0], [3,2,0], [4,2,3], [2,4,1], ], 

    points = [[0.000000,0.000000,0.000000], 
              [0.000000,1.000000,0.000000], 
              [0.500000,0.500000,1.000000], 
              [1.000000,0.000000,0.000000], 
              [1.000000,1.000000,0.000000], ]);

- Juha

---

很好的例子。但它并没有生成一个“封闭表面” - 底部面（2个三角形）缺失了。有没有一种方法可以用三角形完全封闭它？ - Chris Koston

这个问题涉及如何在三维中创建二维表面。这是创建网格时常见的问题（或特征）。您无法将二维表面映射到三维对象，而不创建某种特殊点、点或线。例如，请考虑地球的纬度和经度是如何定义的。在赤道上它们是完全二维的，在极点上完全疯狂。纬度= 90度定义了北极，而在北极，您可以平稳连续地从经度15转移到经度67或经度87。关键是没有标准方法可以做到您所要求的事情。 - Juha

我会将问题分成两个（或更多）部分，并创建两个表面。或者，可以创建3D三角剖分并获取三角形的外表面。如果您找到了一个好的方法，请告诉我。这是一个非常常见的实际问题。 - Juha

我的点云是由自制的3D扫描仪生成的结果。我发现有一些具有相同z值的点云可以连接在一起。然后，我可以尝试在每对平行线之间进行三角剖分。 - Chris Koston

Answer 2

看起来你想计算点云的凸包。我认为这就是你想做的：

from scipy.spatial import ConvexHull

hull = ConvexHull(points)
indices = hull.simplices
vertices = points[indices]

- Jaime

---

在这张图片链接中，您可以看到我想三角化的点云。当我尝试进行凸包时，我会得到像图片中那样的三角剖分。我只想连接最近的点以形成三角形。 - s.t.e.a.l.t.h

@s.t.e.a.l.t.h，您无法使用scipy.spatial函数来完成您所要求的操作，因为它基于qhull，而该库“不支持非凸面的三角剖分...”。您可以将点分解为多个凸面。 - Jonathan A. Gross

是的，谢谢，我通过CGAL泊松重建解决了这个问题。 - s.t.e.a.l.t.h

Answer 3

在Jaime的回答基础上，举例详细说明：

最初的回答：

import matplotlib as mpl
import matplotlib.pyplot as plt
import mpl_toolkits.mplot3d as a3
import numpy as np
import scipy as sp
from scipy import spatial as sp_spatial


def icosahedron():
    h = 0.5*(1+np.sqrt(5))
    p1 = np.array([[0, 1, h], [0, 1, -h], [0, -1, h], [0, -1, -h]])
    p2 = p1[:, [1, 2, 0]]
    p3 = p1[:, [2, 0, 1]]
    return np.vstack((p1, p2, p3))


def cube():
    points = np.array([
        [0, 0, 0], [0, 0, 1], [0, 1, 0], [0, 1, 1],
        [1, 0, 0], [1, 0, 1], [1, 1, 0], [1, 1, 1],
    ])
    return points


points = icosahedron()
# points = cube()

hull = sp_spatial.ConvexHull(points)
indices = hull.simplices
faces = points[indices]

print('area: ', hull.area)
print('volume: ', hull.volume)


fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.dist = 30
ax.azim = -140
ax.set_xlim([0, 2])
ax.set_ylim([0, 2])
ax.set_zlim([0, 2])
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')

for f in faces:
    face = a3.art3d.Poly3DCollection([f])
    face.set_color(mpl.colors.rgb2hex(sp.rand(3)))
    face.set_edgecolor('k')
    face.set_alpha(0.5)
    ax.add_collection3d(face)

plt.show()

这应该描绘以下图示：

最初的回答：

- gmagno

---

这并不能解决原帖作者的问题，因为表面可能是崎岖不平的，所以使用凸包可能会导致一些表面点无法三角化。使用α形状作为替代方案可能会起作用。 - sodiumnitrate