python - SciPy 空间 Delaunay/ConvexHull 混淆

Question

我正在尝试生成随机凸多面体。我生成一组随机 3D 坐标，然后找到它们的凸包(到目前为止一切顺利)。

然后我想我会使用 Delaunay 三角剖分来给我凸包的三角剖分。这是我的基本理解开始显现的地方!

这是代码

import numpy as np
from scipy.spatial import ConvexHull
from scipy.spatial import Delaunay
import matplotlib as mpl
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# Generate random points & convex hull
points = np.random.rand(20,3)
hull = ConvexHull(points)

fig = plt.figure()
ax = fig.gca(projection = '3d')

# Plot hull's vertices
for vert in hull.vertices:    
    ax.scatter(points[vert,0], points[vert,1], zs=points[vert,2])#, 'ro')

# Calculate Delaunay triangulation & plot
tri = Delaunay(points[hull.vertices])
for simplex in tri.simplices:
    vert1 = [points[simplex[0],0], points[simplex[0],1], points[simplex[0],2]]
    vert2 = [points[simplex[1],0], points[simplex[1],1], points[simplex[1],2]]
    vert3 = [points[simplex[2],0], points[simplex[2],1], points[simplex[2],2]]
    vert4 = [points[simplex[3],0], points[simplex[3],1], points[simplex[3],2]]
    ax.plot([vert1[0], vert2[0]], [vert1[1], vert2[1]], zs = [vert1[2], vert2[2]])
    ax.plot([vert2[0], vert3[0]], [vert2[1], vert3[1]], zs = [vert2[2], vert3[2]])
    ax.plot([vert3[0], vert4[0]], [vert3[1], vert4[1]], zs = [vert3[2], vert4[2]])
    ax.plot([vert4[0], vert1[0]], [vert4[1], vert1[1]], zs = [vert4[2], vert1[2]])

plt.show()

有几件事让我担心，情节有时会漏掉船体上的一些点，这似乎是 Delaunay 四面体化，我想我不应该对此感到惊讶，但不是我想要的。

我只想对船体表面进行三角剖分，所以我想是一个包含表面刻面的单纯形？这可能吗？

谢谢

B

编辑:在下面 pv 的启示性帖子之后，我修改了代码如下；

import numpy as np
import pylab as pl
import scipy as sp
from scipy.spatial import ConvexHull
from scipy.spatial.distance import euclidean
import matplotlib.pyplot as plt
import mpl_toolkits.mplot3d as a3

aspect = 0
while aspect == 0:

    # Generate random points & convex hull
    points = np.random.rand(20,3)
    hull = ConvexHull(points)

    # Check aspect ratios of surface facets
    aspectRatio = []
    for simplex in hull.simplices:
        a = euclidean(points[simplex[0],:], points[simplex[1],:])
        b = euclidean(points[simplex[1],:], points[simplex[2],:])
        c = euclidean(points[simplex[2],:], points[simplex[0],:])
        circRad = (a*b*c)/(np.sqrt((a+b+c)*(b+c-a)*(c+a-b)*(a+b-c)))
        inRad = 0.5*np.sqrt(((b+c-a)*(c+a-b)*(a+b-c))/(a+b+c))
        aspectRatio.append(inRad/circRad)

    # Threshold for minium allowable aspect raio of surface facets
    if np.amin(aspectRatio) > 0.3:
        aspect = 1

ax = a3.Axes3D(pl.figure())
facetCol = sp.rand(3) #[0.0, 1.0, 0.0]

# Plot hull's vertices
#for vert in hull.vertices:    
#    ax.scatter(points[vert,0], points[vert,1], zs=points[vert,2])

# Plot surface traingulation
for simplex in hull.simplices:
    vtx = [points[simplex[0],:], points[simplex[1],:], points[simplex[2],:]]
    tri = a3.art3d.Poly3DCollection([vtx], linewidths = 2, alpha = 0.8)
    tri.set_color(facetCol)
    tri.set_edgecolor('k')
    ax.add_collection3d(tri)

plt.axis('off')
plt.show()

现在一切都如我所愿。我添加了纵横比阈值以确保更好的三角测量。

B

Answer 1

一些事情:

• 您将 points[hull.vertices] 作为 Delaunay 的参数，因此 tri.simplices 中的整数是 points[hull.vertices] 的索引，而不是点，这样你最终会画错点
• 四面体有 6 个脊，但您只绘制 4 个
• 如果您只需要对凸包曲面进行三角剖分，可以使用 hull.simplices

即，

for simplex in hull.simplices:
    xs, ys, zs = points[simplex].T
    xs = np.r_[xs, xs[0]] # close polygons
    ys = np.r_[ys, ys[0]]
    zs = np.r_[zs, zs[0]]
    ax.plot(xs, ys, zs)

或者只是:

ax.plot_trisurf(points[:,0], points[:,1], points[:,2],
                triangles=hull.simplices)