### Author:  [@Adam\_Hendry](https://discourse.vtk.org/u/Adam_Hendry)
#### Post date: [September 9, 2021, 1:44am UTC](https://discourse.vtk.org/t/plot-3d-alpha-shape/6584/1 "2021-09-09T01:44:47Z")
### 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](https://stackoverflow.com/questions/26303878/alpha-shapes-in-3d). 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](https://graphics.stanford.edu/courses/cs268-11-spring/handouts/AlphaShapes/as_fisher.pdf) and [Discrete Differential Geometry](https://www.youtube.com/watch?v=mas-PUA3OvA)):
- **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
* * *
```python
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
### Author:  [@Adam\_Hendry](https://discourse.vtk.org/u/Adam_Hendry)
#### Post date: [September 9, 2021, 3:15am UTC](https://discourse.vtk.org/t/plot-3d-alpha-shape/6584/2 "2021-09-09T03:15:05Z")
@MatthewFlamm @mwestphal @ben.boeckel Any ideas?
---
### Author:  [@Adam\_Hendry](https://discourse.vtk.org/u/Adam_Hendry)
#### Post date: [September 9, 2021, 1:39pm UTC](https://discourse.vtk.org/t/plot-3d-alpha-shape/6584/3 "2021-09-09T13:39:56Z")
### Update 09-SEP-2021
Per @akaszynski, the problem does indeed appear to be a combination of unique and negative points. He fixed this issue with the following:
```python
pts = np.stack((x.ravel(), y.ravel(), z.ravel()), axis=1) + 10
return np.unique(pts, axis=0) - 10
```
However, if someone could perform a deeper investigation to determine why this is the problem, that would help.
---
### Author:  [@cory.quammen](https://discourse.vtk.org/u/cory.quammen)
#### Post date: [September 9, 2021, 1:44pm UTC](https://discourse.vtk.org/t/plot-3d-alpha-shape/6584/4 "2021-09-09T13:44:09Z")
As far as this question relates to VTK, please not that `vtkDelaunay3D` incorporates some support for alpha shapes. Please see the [documentation](https://vtk.org/doc/nightly/html/classvtkDelaunay3D.html) for more details.
As for debugging your code, you can try a couple things:
- make all coordinates positive to rule out negative coordinates being an issue (I doubt that is an issue)
- build a subset of the final PolyData from a subset of the faces and examine point ids to see where things are going wrong.
That is the extent of the suggestions I am able to provide.
---
### Author:  [@Adam\_Hendry](https://discourse.vtk.org/u/Adam_Hendry)
#### Post date: [September 9, 2021, 4:57pm UTC](https://discourse.vtk.org/t/plot-3d-alpha-shape/6584/5 "2021-09-09T16:57:40Z")
@cory.quammen Thank you for the suggestion regarding `vtkDelaunay3D`. I am noticing on my system that this will take roughly 5 hours to complete for a dataset I am using that has ~12 million points.
Is there a way to make this run faster? I am using the `pip install` of `vtk-9.0.3` on Windows 10 x64.
---
### Author:  [@Adam\_Hendry](https://discourse.vtk.org/u/Adam_Hendry)
#### Post date: [September 9, 2021, 5:09pm UTC](https://discourse.vtk.org/t/plot-3d-alpha-shape/6584/6 "2021-09-09T17:09:41Z")
A quick Google search seems to indicate that parallel versions of Delaunay 3D exist. For example,
[Paper 01](https://www.cs.cmu.edu/~./blelloch/papers/BBK06.pdf)
[Paper 02](https://cse.buffalo.edu/faculty/miller/Courses/CSE633/adarsh-prakash-Spring-2017-CSE633.pdf)
is there a parallel vtkDelaunay3D?
---
### Author:  [@cory.quammen](https://discourse.vtk.org/u/cory.quammen)
#### Post date: [September 9, 2021, 5:12pm UTC](https://discourse.vtk.org/t/plot-3d-alpha-shape/6584/7 "2021-09-09T17:12:29Z")
> [@Adam\_Hendry](#):
>
> is there a parallel vtkDelaunay3D?
Alas, no.