Is there an "alpha shape" function in 3 dimensions in python, other than the CGAL python bindings?
Alternatively, is there a way to extend the example below into 3D?
2D example: draw a smooth polygon around data points in a scatter plot, in matplotlib
I'm currently calculating volume using this ConvexHull example, but for my purposes the volumes are inflated due to the "convex" constraint.
Thanks,
The alpha parameter is defined as the value a , such that an edge of a disk of radius 1/ a can be drawn between any two edge members of a set of points and still contain all the points.
a — Critical alpha radius a is the value of the alpha radius that produces an alpha shape, which either encloses all points (if type is 'all-points' ), or encloses all points within a single region (if type is 'one-region' ).
I wrote some code for finding alpha shape surface. I hope this helps.
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
If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!
Donate Us With