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()

有些事情与我有关,情节有时会忽略船体上的点,这似乎是德劳内四面体化,我想我不应该对此感到惊讶,但不是我所追求的.

我想只对船体表面进行三角测量,所以我猜一个单面包含表面刻面？这可能吗？

谢谢

乙

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

现在一切都按照我的希望运作了.我添加了宽高比阈值以确保更好的三角测量.

乙

Answer 1

有些事情:

• 你给points[hull.vertices]Delaunay作为一个参数,所以整数tri.simplices是索引points[hull.vertices],而不是points,所以你最终绘制错误的点
• Tetrahedra有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)