Starting with a finite set of 3D points, Plotly can generate a Mesh3d object, that depending on a key value can be the convex hull of that set, its Delaunay triangulation or an alpha set.
This notebook is devoted to the presentation of the alpha shape as a computational geometric object, its interpretation, and visualization with Plotly.
Alpha shape of a finite point set $S$ is a polytope whose structure depends only on the set $S$ and a parameter $\alpha$.
Although it is less known in comparison to other computational geometric objects, it has been used in many practical applications in pattern recognition, surface reconstruction, molecular structure modeling, porous media, astrophysics.
In order to understand how the algorithm underlying Mesh3d works, we present shortly a few notions of Computational Geometry.
Let S be a finite set of 2D or 3D points. A point is called $0$-simplex or vertex. The convex hull of:
In [1]:
from IPython.display import IFrame
IFrame('https://plot.ly/~empet/13475/', width=800, height=350)
Out[1]:
If $T$ is the set of points defining a $k$-simplex, then any proper subset of $T$ defines an $\ell$-simplex, $\ell<k$. These $\ell$-simplexes (or $\ell$-simplices) are called faces.
A 2-simplex has three $1$-simplexes, and three 0-simplexes as faces, whereas a tetrahedron has as faces three 2-simplexes, six 1-simplexes and four zero simplexes.
k-simplexes are building blocks for different structures in Computational Geometry, mainly for creating meshes from point clouds.
Let $S$ be a finite set in $\mathbb{R}^d$, $d=2,3$ (i.e. a set of 2D or 3D points). A collection $\mathcal{K}$ of k-simplexes, $0\leq k\leq d$, having as vertices the points of $S$, is a simplicial complex if its simplexes have the following properties:
The next figure illustrates a simplicial complex(left), and a collection of $k$-simplexes (right), $0\leq k\leq 2$ that do not form a simplicial complex because the condition 2 in the definition above is violated.
In [2]:
IFrame('https://plot.ly/~empet/13503/', width=600, height=475)
Out[2]:
Triangular meshes used in computer graphics are examples of simplicial complexes.
The underlying space of a simplicial complex, $\mathcal{K}$, denoted $|\mathcal{K}|$, is the union of its simplexes, i.e. it is a region in plane or in the 3D space, depending on whether d=2 or 3.
A subcomplex of the simplicial complex $\mathcal{K}$ is a collection, $\mathcal{L}$, of simplexes in $\mathcal{K}$ that also form a simplicial complex.
The points of a finite set $S$ in $\mathbb{R}^2$ (respectively $\mathbb{R}^3$) are in general position if no $3$ (resp 4) points are collinear (coplanar), and no 4 (resp 5) points lie on the same circle (sphere).
A particular simplicial complex associated to a finite set of 2D or 3D points, in general position, is the Delaunay triangulation.
A triangulation of a finite point set $S \subset \mathbb{R}^2$ (or $\mathbb{R}^3$) is a collection $\mathcal{T}$ of triangles (tetrahedra), such that:
A Delaunay triangulation of the set $S\subset\mathbb{R}^2$ ($\mathbb{R}^3$) is a triangulation with the property that the open balls bounded by the circumcircles (circumspheres) of the triangulation triangles (tetrahedra) contain no point in $S$. One says that these balls are empty.
If the points of $S$ are in general position, then the Delaunay triangulation of $S$ is unique.
Here is an example of Delaunay triangulation of a set of ten 2D points. It illustrates the emptiness of two balls bounded by circumcircles.
In [4]:
IFrame('https://plot.ly/~empet/13497/', width=550, height=550)
Out[4]:
The notion of Alpha Shape was introduced by Edelsbrunner with the aim to give a mathematical description of the shape of a point set. In this notebook we give a constructive definition of this geometric structure. A more detailed approach of 3D alpha shapes can be found in the original paper.
An intuitive description of the alpha shape was given by Edelsbrunner and his coauthor in a preprint of the last paper mentioned above:
A huge mass of ice-cream fills a region in the 3D space, and the point set $S$ consists in hard chocolate pieces spread in the ice-cream mass. Using a sphere-formed ice-cream spoon we carve out the ice-cream such that to avoid bumping into chocolate pieces. At the end of this operation the region containing the ciocolate pieces and the remaining ice cream is bounded by caps, arcs and points of chocolate. Straightening all round faces to triangles and line segments we get the intuitive image of the alpha shape of the point set $S$.
Now we give the steps of the computational alpha shape construction.
Let $S$ be a finite set of points from $\mathbb{R}^d$, in general position, $\mathcal{D}$ its Delaunay triangulation and $\alpha$ a positive number.
Select the d-simplexes of $\mathcal{D}$ (i.e. triangles in the case d=2, respectively tetrahedra for d=3) whose circumsphere has the radius less than $\alpha$. These simplexes and their faces form a simplicial subcomplex of the Delaunay triangulation, $\mathcal{D}$. It is denoted $\mathcal{C}_\alpha$, and called $\alpha$-complex.
The $\alpha$-shape of the set $S$ is defined by its authors, either as the underlying space of the $\alpha$-complex, i.e. the union of all its simplexes or as the boundary of the $\alpha$-complex.
The boundary of the $\alpha$-complex is the subcomplex consisting in all k-simplexes, $0\leq k<d$, that are faces of a single $d$-simplex (these are called external faces).
In the ice-cream example the alpha shape was defined as the boundary of the alpha-complex.
The underlying space of the $\alpha$-complex is the region where the ice-cream spoon has no access, because its radius ($\alpha$) exceeds the radius of circumscribed spheres to tetrahedra formed by pieces of chocolate.
To get insight into the process of construction of an alpha shape we illustrate it first for a set of 2D points.
The following panel displays the Delaunay triangulation of a set of 2D points, and a sequence of $\alpha$-complexes (and alpha shapes):
In [5]:
IFrame('https://plot.ly/~empet/13479/', width=825, height=950)
Out[5]:
We notice that the Delaunay triangulation has as boundary a convex set (it is a triangulation of the convex hull of the given point set).
Each $\alpha$-complex is obtained from the Delaunay triangulation, removing the triangles whose circumcircle has radius greater or equal to alpha.
In the last subplot the triangles of the $0.115$-complex are filled in with light blue. The filled in region is the underlying space of the $0.115$-complex.
The $0.115$-alpha shape of the given point set can be considered either the filled in region or its boundary.
This example illustrates that the underlying space of an $\alpha$-complex in neither convex nor necessarily connected. It can consist in many connected components (in our illustration above, $|\mathcal{C}_{0.115}|$ has three components).
In a family of alpha shapes, the parameter $\alpha$ controls the level of detail of the associated alpha shape. If $\alpha$ decreases to zero, the corresponding alpha shape degenerates to the point set, $S$, while if it tends to infinity the alpha shape tends to the convex hull of the set $S$.
In order to generate the alpha shape of a given set of 3D points corresponding to a parameter $\alpha$,
the Delaunay triagulation or the convex hull we define an instance of the go.Mesh3d class. The real value
of the key alphahull points out the mesh type to be generated:
alphahull=$1/\alpha$ generates the $\alpha$-shape, -1 corresponds to the Delaunay
triangulation and 0, to the convex hull of the point set.
The other parameters in the definition of a Mesh3d are given here.
Mesh3d generates and displays an $\alpha$-shape as the boundary of the $\alpha$-complex.
An intuitive idea on the topological structure modification, as $\alpha=1/$alphahull varies can be gained from the following three different alpha shapes of the same point set:
In [6]:
IFrame('https://plot.ly/~empet/13481/', width=900, height=950)
Out[6]:
We notice in the subplots above that as alphahull increases, i.e. $\alpha$ decreases, some parts of the alpha shape shrink and
develop enclosed void regions. The last plotted alpha shape points out a polytope that contains faces of tetrahedra,
and patches of triangles.
In some cases as $\alpha$ varies it is also possible to develop components that are strings of edges and even isolated points.
Such experimental results suggested the use of alpha shapes in modeling molecular structure. A search on WEB gives many results related to applications of alpha shapes in structural molecular biology.
Here is an alpha shape illustrating a molecular-like structure associated to a point set of 5000 points.
In [7]:
import numpy as np
import plotly.plotly as py
import plotly.graph_objs as go
from plotly import tools as tls
Load data:
In [8]:
pts = np.loadtxt('Data/data-file.txt')
x, y, z = zip(*pts)
Define two traces: one for plotting the point set and another for the alpha shape:
In [9]:
points = go.Scatter3d(mode='markers',
name='',
x =x,
y= y,
z= z,
marker=dict(size=2, color='#458B00'))
In [10]:
simplexes = go.Mesh3d(alphahull =10.0,
name = '',
x =x,
y= y,
z= z,
color='#90EE90',
opacity=0.15)
In [11]:
axis = dict(showbackground=True,
backgroundcolor="rgb(245, 245, 245)",
gridcolor="rgb(255, 255, 255)",
gridwidth=2,
zerolinecolor="rgb(255, 255, 255)",
tickfont=dict(size=11),
titlefont =dict(size=12))
In [12]:
x_style = dict(axis, range=[-2.85, 4.25], tickvals=np.linspace(-2.85, 4.25, 5)[1:].round(1))
y_style = dict(axis, range=[-2.65, 1.32], tickvals=np.linspace(-2.65, 1.32, 4)[1:].round(1))
z_style = dict(axis, range=[-3.67,1.4], tickvals=np.linspace(-3.67, 1.4, 5).round(1))
In [13]:
layout = go.Layout(title='Alpha shape of a set of 3D points. Alpha=0.1',
width=500,
height=500,
scene=dict(xaxis=x_style,
yaxis=y_style,
zaxis=z_style))
In [14]:
fig = go.FigureWidget(data=[points, simplexes], layout=layout)
#fig
In [16]:
fig = go.FigureWidget(data=[points, simplexes], layout=layout)
#py.plot(fig, filename='3D-AlphaS-ex')
Out[16]:
In [17]:
IFrame('https://plot.ly/~empet/13499/', width=550, height=550)
Out[17]:
We construct the alpha shape of a set of 2D points from the Delaunay triangulation,
defined as a scipy.spatial.Delaunay object.
In [18]:
from scipy.spatial import Delaunay
In [19]:
def sq_norm(v): #squared norm
return np.linalg.norm(v)**2
Compute the circumcenter and circumradius of a triangle (see their definitions here):
In [20]:
def circumcircle(points,simplex):
A = [points[simplex[k]] for k in range(3)]
M = [[1.0]*4]
M += [[sq_norm(A[k]), A[k][0], A[k][1], 1.0 ] for k in range(3)]
M = np.asarray(M, dtype=np.float32)
S = np.array([0.5*np.linalg.det(M[1:, [0,2,3]]), -0.5*np.linalg.det(M[1:, [0,1,3]])])
a = np.linalg.det(M[1:, 1:])
b = np.linalg.det(M[1:, [0,1,2]])
return S/a, np.sqrt(b/a + sq_norm(S)/a**2) #center=S/a, radius=np.sqrt(b/a+sq_norm(S)/a**2)
Filter out the Delaunay triangulation to get the $\alpha$-complex:
In [21]:
def get_alpha_complex(alpha, points, simplexes):
#alpha is the parameter for the alpha shape
#points are given data points
#simplexes is the list of indices in the array of points
#that define 2-simplexes in the Delaunay triangulation
return filter(lambda simplex: circumcircle(points,simplex)[1] < alpha, simplexes)
In [22]:
pts = np.loadtxt('Data/data-ex-2d.txt')
tri = Delaunay(pts)
In [23]:
colors = ['#C0223B', '#404ca0', 'rgba(173,216,230, 0.5)']# colors for vertices, edges and 2-simplexes
Get data for the Plotly plot of a subcomplex of the Delaunay triangulation:
In [24]:
def Plotly_data(points, complex_s):
#points are the given data points,
#complex_s is the list of indices in the array of points defining 2-simplexes(triangles)
#in the simplicial complex to be plotted
X = []
Y = []
for s in complex_s:
X += [points[s[k]][0] for k in [0,1,2,0]] + [None]
Y += [points[s[k]][1] for k in [0,1,2,0]] + [None]
return X, Y
In [25]:
def make_trace(x, y, point_color=colors[0], line_color=colors[1]):# define the trace
#for an alpha complex
return go.Scatter(mode='markers+lines', #vertices and
#edges of the alpha-complex
name='',
x=x,
y=y,
marker=dict(size=6.5, color=point_color),
line=dict(width=1.25, color=line_color))
In [34]:
figure = tls.make_subplots(rows=1, cols=2,
subplot_titles=('Delaunay triangulation', 'Alpha shape, alpha=0.15'),
horizontal_spacing=0.1,
)
In [35]:
title = 'Delaunay triangulation and Alpha Complex/Shape for a Set of 2D Points'
figure.layout.update(title=title,
font=dict(family="Open Sans, sans-serif"),
showlegend=False,
hovermode='closest',
autosize=False,
width=800,
height=460,
margin=dict(l=65,
r=65,
b=85,
t=120));
In [36]:
axis_style = dict(showline=True,
mirror=True,
zeroline=False,
showgrid=False,
showticklabels=True,
range=[-0.1,1.1],
tickvals=[0, 0.2, 0.4, 0.6, 0.8, 1.0],
ticklen=5
)
In [37]:
for s in range(1,3):
figure.layout.update({'xaxis{}'.format(s): axis_style})
figure.layout.update({'yaxis{}'.format(s): axis_style})
In [42]:
alpha_complex = list(get_alpha_complex(0.15, pts, tri.simplices))
In [43]:
X, Y = Plotly_data(pts, tri.simplices)# get data for Delaunay triangulation
figure.append_trace(make_trace(X, Y), 1, 1)
X, Y = Plotly_data(pts, alpha_complex)# data for alpha complex
figure.append_trace(make_trace(X, Y), 1, 2)
In [44]:
shapes = []
for s in alpha_complex: #fill in the triangles of the alpha complex
A = pts[s[0]]
B = pts[s[1]]
C = pts[s[2]]
shapes.append(dict(path=f'M {A[0]}, {A[1]} L {B[0]}, {B[1]} L {C[0]}, {C[1]} Z',
fillcolor='rgba(173, 216, 230, 0.5)',
line=dict(color=colors[1], width=1.25),
xref='x2',
yref='y2'
))
figure.layout.shapes=shapes
In [46]:
py.plot(figure, filename='2D-AlphaS-ex', width=850)
Out[46]:
In [47]:
IFrame('https://plot.ly/~empet/13501', width=800, height=460)
Out[47]:
In [ ]: