#!/usr/bin/env python
# coding: utf-8

# ## Plotly Alpha Shapes as Mesh3d instances ##

# Starting with a finite set of 3D points, Plotly can generate a [`Mesh3d`](https://plot.ly/python/reference/#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](https://en.wikipedia.org/wiki/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.

# ### Simplicial complexes and Delaunay triangulation ###

# Let S be a finite set of 2D or 3D points. A point is called $0$-simplex or vertex. The [convex hull](https://en.wikipedia.org/wiki/Convex_hull) of:
#  - two distinct points is  a 1-simplex or edge;
#  - three non-colinear points is a 2-simplex or triangle;
#  - four non-coplanar points in $\mathbb{R}^3$  is a 3-simplex or tetrahedron;

# In[1]:


from IPython.display import IFrame
IFrame('https://plot.ly/~empet/13475/',  width=800, height=350)


# 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:
#   1. If $\sigma$ is a simplex in $\mathcal{K}$, then all its faces are also simplexes in $\mathcal{K}$;
#   2. If $\sigma, \tau$ are two simplexes in $\mathcal{K}$, then their intersection is either empty or a face in both simplexes.
#                                       

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


# 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:
#   1. The union of all triangles (tetrahedra) in $\mathcal{T}$ is the convex hull of $S$.
#   2. The union of all  vertices of triangles (tetrahedra) in  $\mathcal{T}$ is the set  $S$.
#   3. For every distinct pair $\sigma, \tau \in \mathcal{T}$, the intersection $\sigma \cap \tau$ is either empty or  a common  face  of $\sigma$ and $\tau$.

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


# ### Alpha shape of a finite set of points ###

#  The notion of Alpha Shape was introduced by [Edelsbrunner](https://users.cs.duke.edu/~edels/Papers/1983-J-06-ShapeSetPoints.pdf) 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](https://users.cs.duke.edu/~edels/Papers/1994-J-04-3DAlphaShapes.pdf). 

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


# 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$.

# ### Plotly `Mesh3d` ###

# 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](https://plot.ly/python/reference/#mesh3d).
# 
# `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)


# 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](https://plot.ly/~empet/13485/molecular-like-structure-defined-as-an-alpha-shape-alphahull1730/)
# is an alpha shape illustrating a molecular-like structure associated to a point set of 5000 points.

# ### Generating an alpha shape with `Mesh3d` ###

# 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')


# In[17]:


IFrame('https://plot.ly/~empet/13499/',  width=550, height=550) 


# ###  Generating the alpha shape of a set  of 2D points ###

# 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](https://en.wikipedia.org/wiki/Circumscribed_circle#Circumcircle_equations)):

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


# In[47]:


IFrame('https://plot.ly/~empet/13501',  width=800, height=460) 


# In[ ]: