用于在python中查找点云边界的3D alpha形状

问题描述

我正在尝试使用3D alpha形状算法查找点云的表面。当我计算外接半径时，行列式a等于零，从而导致错误“在double_scalars中遇到的被零除”。我该怎么办？非常感谢！这是代码：

def alpha_shape_3D(points,alpha):
    tetra = delaunay(points)
    edge = []
    for i1,i2,i3,i4 in tetra.vertices:
        pa = points[i1]
        pb = points[i2]
        pc = points[i3]
        pd = points[i4]
        pa2 = np.dot(pa,pa)
        pb2 = np.dot(pb,pb)
        pc2 = np.dot(pc,pc)
        pd2 = np.dot(pd,pd)
        a = np.linalg.det(np.array([np.append(pa,1),np.append(pb,np.append(pc,np.append(pd,1)]))

        Dx = np.linalg.det(np.array([np.array([pa2,pa[1],pa[2],1]),np.array([pb2,pb[1],pb[2],np.array([pc2,pc[1],pc[2],np.array([pd2,pd[1],pd[2],1])]))

        Dy = - np.linalg.det(np.array([np.array([pa2,pa[0],pb[0],pc[0],pd[0],1])]))

        Dz = np.linalg.det(np.array([np.array([pa2,1])]))

        c = np.linalg.det(np.array([np.array([pa2,pa[2]]),pb[2]]),pc[2]]),pd[2]])]))
       
        r = math.sqrt(math.pow(Dx,2) + math.pow(Dy,2) + math.pow(Dz,2) - 4 * a * c) / (2 * abs(a))  # error occurs here,since some value of a is zero.
        if r<alpha:
            edge.append([pa,pb,pc,pd])
    tetras = np.array(edge)
    # triangles
    TriComb = np.array([(0,1,2),(0,3),2,(1,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,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

解决方法

如果没有您的 points 数据，很难确定，但看起来您的数据已退化。例如，如果您的所有点都在同一平面上，您就可以得到这个结果。