问题
我正在尝试为 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：
DT计算点云的 Delaunay 三角剖分- 查找 alpha 复形：搜索 Delaunay 三角剖分中的所有单纯形，并且 (a) 如果单纯形周围的任何球为空且半径小于
alpha（称为“alpha 测试”），然后将其添加到 alpha 复形- 阿尔法复合体的边界是阿尔法形状
代码
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 文档）