我有这个问题。我尝试通过scipy.spatial.Delaunay对点云进行三角测量。我用了:
tri = Delaunay(points) # points: np.array() of 3d points
indices = tri.simplices
vertices = points[indices]
但是,这段代码返回了四面体。怎么可能只返回表面三角形?
谢谢
要使其像代码形式一样工作,您必须将曲面参数化为2D。例如,在球(r,theta,psi)的情况下,半径是常数(将其丢弃)并且点由(θ,psi)给出,即2D。
Scipy Delaunay是N维三角剖分,因此如果你给3D点,它会返回3D对象。给它2D点并返回2D对象。
下面是我用于为openSCAD创建多面体的脚本。 U和V是我的参数化(x和y),这些是我给Delaunay的坐标。请注意,现在“Delaunay三角剖分属性”仅适用于u,v坐标(角度在uv-spatial中最大化而不是xyz -space等)。
该示例是来自http://matplotlib.org/1.3.1/mpl_toolkits/mplot3d/tutorial.html的修改后的副本,该副本最初使用三角测量功能(最终映射到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], ]);
看起来你想要计算点云的convex hull。我想这就是你想要做的:
from scipy.spatial import ConvexHull
hull = ConvexHull(points)
indices = hull.simplices
vertices = points[indices]
按照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()
哪个应该描绘如下图: