python - 3D 中的 Alpha 形状 - 跟进

Question

问题

---

我正在尝试为 python 中的 3D 点云的 alpha 形状实现 Edelsbrunner 算法，如这篇 SO post中所述。但是，我无法绘制结果。我的球体有一半看起来不错，而另一半则乱码。

我怀疑这可能与我有负坐标的事实有关，但我不确定。

---

定义

---

我正在添加这些，以便程序员可以做出贡献而不会被数学所困。这些是简化的定义，并不意味着精确（请随意跳过这部分；有关更多信息，请参阅Alpha 形状和离散微分几何简介）：

• Delaunay 三角剖分： 对于 2D 点集，将点细分为三角形（即“三角剖分”），其中每个三角形的外接圆（即“外接圆”）不包含集合中的其他点。对于 3D 点，将“triangle”替换为“tetrahedron”，将“circumcircle”替换为“circumsphere”。

• 仿射独立： 点的集合，p0, ..., pk使得所有向量vi := pi-p0都是线性独立的（即在 2D 中不共线，在 3D 中不共面）；也称为“一般位置点”

• k-单纯形： k+1 个仿射独立点的凸包；我们称这些点为顶点。

0-单纯形 = 点（由 0+1 = 1 个点组成）
1-单纯形 = 线（由 1+1 = 2 个点组成）
2-单纯形 = 三角形（由 2+1 = 3 个点组成）
3-单纯形 = 四面体（由 3+1 = 4 分组成）

• 面： 任何单纯形，其顶点是另一个单纯形的顶点的子集；即“单纯形的一部分”

• （几何）单纯复形： 单纯形的集合，其中（1）两个单纯形的交点是单纯形，（2）单纯形的每个面都在复合形中；即“一堆单纯形”

• alpha-exposed： 点集中的单纯形，其中半径为 alpha 通过其顶点的圆 (2D) 或球 (3D) 不包含点集中的任何其他点

• alpha 形状： 点集的所有 alpha 暴露单纯形的边界

算法

---

Edelsbrunner 算法如下：

给定一个点云pts：

1. DT计算点云的 Delaunay 三角剖分
2. 查找 alpha 复形：搜索 Delaunay 三角剖分中的所有单纯形，并且 (a) 如果单纯形周围的任何球为空且半径小于alpha（称为“alpha 测试”），然后将其添加到 alpha 复形
3. 阿尔法复合体的边界是阿尔法形状

代码

---

from scipy.spatial import Delaunay
import numpy as np
from collections import defaultdict
from matplotlib import pyplot as plt
import pyvista as pv

fig = plt.figure()
ax = plt.axes(projection="3d")

plotter = pv.Plotter()

def alpha_shape_3D(pos, alpha):
    """
    Compute the alpha shape (concave hull) of a set of 3D points.
    Parameters:
        pos - np.array of shape (n,3) points.
        alpha - alpha value.
    return
        outer surface vertex indices, edge indices, and triangle indices
    """

    tetra = Delaunay(pos)
    # Find radius of the circumsphere.
    # By definition, radius of the sphere fitting inside the tetrahedral needs 
    # to be smaller than alpha value
    # http://mathworld.wolfram.com/Circumsphere.html
    tetrapos = np.take(pos,tetra.vertices,axis=0)
    normsq = np.sum(tetrapos**2,axis=2)[:,:,None]
    ones = np.ones((tetrapos.shape[0],tetrapos.shape[1],1))
    a = np.linalg.det(np.concatenate((tetrapos,ones),axis=2))
    Dx = np.linalg.det(np.concatenate((normsq,tetrapos[:,:,[1,2]],ones),axis=2))
    Dy = -np.linalg.det(np.concatenate((normsq,tetrapos[:,:,[0,2]],ones),axis=2))
    Dz = np.linalg.det(np.concatenate((normsq,tetrapos[:,:,[0,1]],ones),axis=2))
    c = np.linalg.det(np.concatenate((normsq,tetrapos),axis=2))
    r = np.sqrt(Dx**2+Dy**2+Dz**2-4*a*c)/(2*np.abs(a))

    # Find tetrahedrals
    tetras = tetra.vertices[r<alpha,:]
    # triangles
    TriComb = np.array([(0, 1, 2), (0, 1, 3), (0, 2, 3), (1, 2, 3)])
    Triangles = tetras[:,TriComb].reshape(-1,3)
    Triangles = np.sort(Triangles,axis=1)
    # Remove triangles that occurs twice, because they are within shapes
    TrianglesDict = defaultdict(int)
    for tri in Triangles:TrianglesDict[tuple(tri)] += 1
    Triangles=np.array([tri for tri in TrianglesDict if TrianglesDict[tri] ==1])
    #edges
    EdgeComb=np.array([(0, 1), (0, 2), (1, 2)])
    Edges=Triangles[:,EdgeComb].reshape(-1,2)
    Edges=np.sort(Edges,axis=1)
    Edges=np.unique(Edges,axis=0)

    Vertices = np.unique(Edges)
    return Vertices,Edges,Triangles

def ptcloud_sphere():
    r = 3
    phi = np.linspace(0, np.pi, 18)
    theta = np.linspace(0, 2 * np.pi, 36)

    PHI, THETA = np.meshgrid(phi, theta)

    x = r * np.sin(PHI) * np.cos(THETA)
    y = r * np.sin(PHI) * np.sin(THETA)
    z = r * np.cos(PHI)

    ax.scatter(x, y, z)
    plt.show()
    
    pts = np.stack((x.ravel(), y.ravel(), z.ravel()), axis=1)

    return np.unique(pts, axis=0)


if __name__ == "__main__":
    pts = ptcloud_sphere()

    verts, edges, faces = alpha_shape_3D(pts, alpha=10)

    faces_conn_list = np.insert(faces, 0, 3, axis=1)
    num_faces = faces.shape[0]

    mesh = pv.PolyData(pts[verts], faces_conn_list, n_faces=num_faces)

    plotter.add_mesh(mesh, reset_camera=True)
    plotter.show()

输出

点云：

阿尔法形状：

2021 年 9 月 9 日更新

根据@akaszynski，问题确实似乎是独特点和负面点的结合。他通过以下方式解决了这个问题：

pts = np.stack((x.ravel(), y.ravel(), z.ravel()), axis=1) + 10

return np.unique(pts, axis=0) - 10

但是，如果有人可以进行更深入的调查以确定问题的原因，那将有所帮助。

2021 年 9 月 9 日更新 #2

根据@AndrasDeak，两者都pyvista支持VTK创建 2d 和 3d alpha 形状。该pyvista函数在底层delaunay_3d使用vtkDelaunay3D，它接受一个alpha参数。

（参见vtkDelaunay3D 文档）

标签： pythonmatplotlibconcave-hullpyvistaalpha-shape