3D 中的 Alpha 形状

Posted 2023-02-14

【中文标题】3D 中的 Alpha 形状

【英文标题】：Alpha shapes in 3D

【发布时间】：2014-12-05 21:16:26

【问题描述】：

除了CGAL python绑定之外，python中是否有3维的“alpha shape”函数？

或者，有没有办法将下面的示例扩展到 3D？

2D 示例：draw a smooth polygon around data points in a scatter plot, in matplotlib

我目前正在使用这个 ConvexHull 示例计算交易量，但出于我的目的，由于“凸”约束，交易量被夸大了。

谢谢，

【问题讨论】：

见：Python bindings for Dionysus，Alpha shapes functions

【参考方案1】：

我写了一些代码来寻找 alpha 形状表面。我希望这会有所帮助。

from scipy.spatial import Delaunay
import numpy as np
from collections import defaultdict

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

【讨论】：

谢谢。它适用于我的测试数据。你也算过体积吗？

如何使用此函数的返回值绘制曲面？

【参考方案2】：

您正在寻找“凹形船体”。 marching cube algorithm 可以用来找到这样的船体。你可以找到a full example here。

限制：如果您的数据来自体积数据集，或者如果您有可以轻松转换为体积数据集（类似体素）的点云，则此方法非常有效。这可以使用一组密集的点相对容易地完成，例如，使用像scipy cKDTree 这样的空间索引器，但如果你有一个稀疏的点云，你最终可能会有点挠头以获得好的结果。