python - SciPy空间Delaunay/ConvexHull混淆

我正在尝试生成随机凸多面体。我生成了一组随机的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四面体化，我想我应该不会对此感到惊讶，但不是我所追求的。

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

谢谢

乙

编辑：在下面的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()

现在一切都按我希望的那样工作。我添加了宽高比阈值以确保更好的三角剖分。

乙

最佳答案

一些事情：


您将points[hull.vertices]作为Delaunay的参数，所以tri.simplices中的整数是points[hull.vertices]的索引，而不是points的索引，因此最终会绘制错误的点
四面体有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)