Plot 3D Alpha Shape

[adam-grant-hendry]

Problem

---

I am trying to implement Edelsbrunner's Algorithm for the alpha shape of a 3D point cloud in python as presented in this SO post. However, I'm having trouble plotting results. Half my sphere looks good, and the other half garbled.

I suspect it may have something to do with the fact that I have negative coordinates, but I'm not sure.

---

Definitions

---

I'm adding these so programmers can contribute without being bogged down by math. These are simplified definitions and not meant to be precise (Feel free to skip this part; for more, see see Introduction to Alpha Shapes and Discrete Differential Geometry):

• Delaunay triangulation: for a 2D point set, a tesselation of the points into triangles (i.e. "triangulation") where the circumscribed circle (i.e. "circumcircle") about every triangle contains no other points in the set. For 3D points, replace "triangle" with "tetrahedron" and "circumcircle" with "circumsphere".

• Affinely independent: a collection of points p0, ..., pk such that all vectors vi := pi-p0 are linearly independent (i.e. in 2D not collinear, in 3D not coplanar); also called "points in general position"

• k-simplex: the convex hull of k+1 affinely-independent points; we call the points its vertices.

0-simplex = point (consists of 0+1 = 1 points)
1-simplex = line (consists of 1+1 = 2 points)
2-simplex = triangle (consists of 2+1 = 3 points)
3-simplex = tetrahedron (consists of 3+1 = 4 points)

• face: any simplex whose vertices are a subset of the vertices of another simplex; i.e. "a part of a simplex"

• (geometric) simplicial complex: a collection of simplices where (1) the intersection of two simplices is a simplex, and (2) every face of a simplex is in the complex; i.e. "a bunch of simplices"

• alpha-exposed: a simplex within a point set where the circle (2D) or ball (3D) of radius alpha through its vertices doesn't contain any other point in the point set

• alpha shape: the boundary of all alpha-exposed simplices of a point set

Algorithm

---

Edelsbrunner's Algorithm is as follows:

Given a point cloud pts:

1. Compute the Delaunay triangulation DT of the point cloud
2. Find the alpha-complex: search all simplices in the Delaunay triangulation and (a) if any ball around a simplex is empty and has a radius less than alpha (called the "alpha test"), then add it to the alpha complex
3. The boundary of the alpha complex is the alpha shape

Code

---

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

Output

Point Cloud:

Alpha Shape:

[adam-grant-hendry]

@banesullivan Would you be able to help with this?

[adam-grant-hendry]

@akaszynski Any ideas?

[akaszynski]

I don't have enough spare time conduct a through review of this, but I fixed this with

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

The problem appears to be a combination of unique and negative points.

[adam-grant-hendry]

Thank you @akaszynski I confirm!

I'll post on VTK discourse and SO to see if anyone has an idea as to why.

[adam-grant-hendry]

Closing as solution was found.