#!/usr/bin/env python
# coding: utf-8
# ## Plotting trisurfs with Plotly ##
# A triangulation of a compact surface is a finite collection of triangles that cover the surface in such a way that every point on the surface is in a triangle, and the intersection of any two triangles is either void, a common edge or a common vertex. A triangulated surface is called tri-surface.
#
# The triangulation of a surface defined as the graph of a continuous function, $z=f(x,y), (x,y)\in D\subset\mathbb{R}^2$ or in a parametric form:
# $$x=x(u,v), y=y(u,v), z=z(u,v), (u,v)\in U\subset\mathbb{R}^2,$$
# is the image through $f$, respectively through the parameterization, of the Delaunay triangulation or an user defined triangulation of the planar domain $D$, respectively $U$.
#
# The Delaunay triangulation of a planar region is defined and illustrated in this [Jupyter Notebook](http://nbviewer.jupyter.org/github/empet/Plotly-plots/blob/master/Plotly-Mesh3d.ipynb).
#
# If the planar region $D$ ($U$) is rectangular, then one defines a meshgrid on it, and the points
# of the grid are the input points for the `scipy.spatial.Delaunay` function that defines the triangulation of $D$, respectively $U$.
#
# ### Triangulation of the Moebius band ###
# The Moebius band is parameterized by:
#
# $$\begin{align*}
# x(u,v)&=(1+0.5 v\cos(u/2))\cos(u)\\
# y(u,v)&=(1+0.5 v\cos(u/2))\sin(u)\quad\quad u\in[0,2\pi],\: v\in[-1,1]\\
# z(u,v)&=0.5 v\sin(u/2)
# \end{align*}
# $$
# Define a meshgrid on the rectangle $U=[0,2\pi]\times[-1,1]$:
# In[1]:
import numpy as np
from scipy.spatial import Delaunay
# In[2]:
import plotly.plotly as py
py.sign_in('empet','api_key')
# In[3]:
u=np.linspace(0,2*np.pi, 24)
v=np.linspace(-1,1, 8)
u,v=np.meshgrid(u,v)
u=u.flatten()
v=v.flatten()
#evaluate the parameterization at the flattened u and v
tp=1+0.5*v*np.cos(u/2.)
x=tp*np.cos(u)
y=tp*np.sin(u)
z=0.5*v*np.sin(u/2.)
#define 2D points, as input data for the Delaunay triangulation of U
points2D=np.vstack([u,v]).T
tri = Delaunay(points2D)#triangulate the rectangle U
points3D=np.vstack((x,y,z)).T
# `tri.simplices` is a `np.array` of integers, of shape (`ntri`,3), where `ntri` is the number of triangles generated by `scipy.spatial.Delaunay`.
# Each row in this array contains three indices, i, j, k, such that points2D[i,:], points2D[j,:], points2D[k,:] are vertices of a triangle in the Delaunay triangulation of the rectangle $U$.
#
# The images of the `points2D` through the surface parameterization are 3D points. The same simplices define the triangles on the surface.
# Setting a combination of keys in `Mesh3d` leads to generating and plotting a tri-surface, in the same way as `plot_trisurf` in matplotlib or `trisurf` in Matlab does.
#
# We note that `Mesh3d` with different combination of keys can generate [alpha-shapes](http://nbviewer.jupyter.org/github/empet/Plotly-plots/blob/master/Plotly-Mesh3d.ipynb).
# To plot the triangles on a surface, we set in Plotly `Mesh3d` the lists of x, y, respectively z- coordinates of the vertices, and the lists of indices, i, j, k, for x, y, z coordinates of all vertices:
# Now we define a function that returns data for a Plotly plot of a tri-surface:
# In[4]:
from plotly.graph_objs import *
# In[5]:
def standard_intensity(x,y,z):
return z
# In[6]:
def plotly_triangular_mesh(x,y,z, faces, intensities=standard_intensity, colorscale="Viridis",
showscale=False, reversescale=False, plot_edges=False):
#x,y,z lists or np.arrays of vertex coordinates
#faces = a numpy array of shape (n_faces, 3)
#intensities can be either a function of (x,y,z) or a list of values
vertices=np.vstack((x,y,z)).T
I,J,K=faces.T
if hasattr(intensities, '__call__'):
intensity=intensities(x,y,z)#the intensities are computed here via the set function,
#that returns the list of vertices intensities
elif isinstance(intensities, (list, np.ndarray)):
intensity=intensities#intensities are given in a list
else:
raise ValueError("intensities can be either a function or a list, np.array")
mesh=dict(type='mesh3d',
x=x,
y=y,
z=z,
colorscale=colorscale,
reversescale=reversescale,
intensity= intensity,
i=I,
j=J,
k=K,
name='',
showscale=showscale
)
if showscale is True:
mesh.update(colorbar=dict(thickness=20, ticklen=4, len=0.75))
if plot_edges is False: # the triangle sides are not plotted
return [mesh]
else:#plot edges
#define the lists Xe, Ye, Ze, of x, y, resp z coordinates of edge end points for each triangle
#None separates data corresponding to two consecutive triangles
tri_vertices= vertices[faces]
Xe=[]
Ye=[]
Ze=[]
for T in tri_vertices:
Xe+=[T[k%3][0] for k in range(4)]+[ None]
Ye+=[T[k%3][1] for k in range(4)]+[ None]
Ze+=[T[k%3][2] for k in range(4)]+[ None]
#define the lines to be plotted
lines=dict(type='scatter3d',
x=Xe,
y=Ye,
z=Ze,
mode='lines',
name='',
line=dict(color= 'rgb(50,50,50)', width=1.5)
)
return [mesh, lines]
# Call this function for data associated to Moebius band:
# In[7]:
pl_RdBu=[[0.0, 'rgb(103, 0, 31)'],
[0.1, 'rgb(176, 23, 42)'],
[0.2, 'rgb(214, 96, 77)'],
[0.3, 'rgb(243, 163, 128)'],
[0.4, 'rgb(253, 219, 199)'],
[0.5, 'rgb(246, 246, 246)'],
[0.6, 'rgb(209, 229, 240)'],
[0.7, 'rgb(144, 196, 221)'],
[0.8, 'rgb(67, 147, 195)'],
[0.9, 'rgb(32, 100, 170)'],
[1.0, 'rgb(5, 48, 97)']]
# In[8]:
data1=plotly_triangular_mesh(x,y,z, tri.simplices, intensities=standard_intensity, colorscale=pl_RdBu,
showscale=True, plot_edges=True)
# Set the layout of the plot:
# In[9]:
axis = dict(
showbackground=True,
backgroundcolor="rgb(230, 230,230)",
gridcolor="rgb(255, 255, 255)",
zerolinecolor="rgb(255, 255, 255)",
)
layout = Layout(
title='Moebius band triangulation',
width=800,
height=800,
showlegend=False,
scene=Scene(xaxis=XAxis(axis),
yaxis=YAxis(axis),
zaxis=ZAxis(axis),
aspectratio=dict(x=1,
y=1,
z=0.5
),
)
)
fig1 = Figure(data=data1, layout=layout)
# In[10]:
py.iplot(fig1, filename='Mobius-band-trisurf')
# ### Triangulation of the surface $z=\sin(-xy)$, defined over a disk ###
# We consider polar coordinates on the disk, $D(0, 1)$, centered at origin and of radius 1, and define
# a meshgrid on the set of points $(r, \theta)$, with $r\in[0,1]$ and $\theta\in[0,2\pi]$:
# In[11]:
n=12# number of radii
h=1.0/(n-1)
r = np.linspace(h, 1.0, n)
theta= np.linspace(0, 2*np.pi, 36)
r,theta=np.meshgrid(r,theta)
r=r.flatten()
theta=theta.flatten()
#Convert polar coordinates to cartesian coordinates (x,y)
x=r*np.cos(theta)
y=r*np.sin(theta)
x=np.append(x, 0)# a trick to include the center of the disk in the set of points. It was avoided
# initially when we defined r=np.linspace(h, 1.0, n)
y=np.append(y,0)
z = np.sin(-x*y)
points2D=np.vstack([x,y]).T
tri=Delaunay(points2D)
# Plot the surface with a modified layout:
# In[12]:
pl_cubehelix=[[0.0, 'rgb(0, 0, 0)'],
[0.1, 'rgb(25, 20, 47)'],
[0.2, 'rgb(21, 60, 77)'],
[0.3, 'rgb(30, 101, 66)'],
[0.4, 'rgb(83, 121, 46)'],
[0.5, 'rgb(161, 121, 74)'],
[0.6, 'rgb(207, 126, 146)'],
[0.7, 'rgb(207, 157, 218)'],
[0.8, 'rgb(193, 202, 243)'],
[0.9, 'rgb(210, 238, 238)'],
[1.0, 'rgb(255, 255, 255)']]
# In[13]:
data2=plotly_triangular_mesh(x,y,z, tri.simplices, intensities=standard_intensity, colorscale=pl_cubehelix,
showscale=True, reversescale=False, plot_edges=False)
fig2 = Figure(data=data2, layout=layout)
fig2['layout'].update(dict(title='Triangulated surface',
scene=dict(camera=dict(eye=dict(x=1.75,
y=-0.7,
z= 0.75)
)
)))
# In[14]:
py.iplot(fig2, filename='cubehexn')
# ### Plotting tri-surfaces from data stored in ply-files ###
# A PLY (Polygon File Format or Stanford Triangle Format) format is a format for storing graphical objects
# that are represented by a triangulation of an object, resulted usually from scanning that object. A Ply file contains the coordinates of vertices, the codes for faces (triangles) and other elements, as well as the color for faces or the normal direction to faces.
#
# In the following we show how we can read a ply file via the Python package, `plyfile`. This package can be installed with `pip`.
# In[15]:
from plyfile import PlyData, PlyElement
# Define a function that extract from plydata the vertices and the faces of a triangulated 3D object:
# In[16]:
def extract_data(plydata):
vertices=list(plydata['vertex'])
vertices=np.asarray(map(list, vertices))
nr_faces=plydata.elements[1].count
faces=np.array([plydata['face'][k][0] for k in range(nr_faces)])
return vertices, faces
# We choose a ply file from a list provided [here](http://people.sc.fsu.edu/~jburkardt/data/ply/ply.html).
# In[17]:
import urllib2
req = urllib2.Request('http://people.sc.fsu.edu/~jburkardt/data/ply/chopper.ply')
opener = urllib2.build_opener()
f = opener.open(req)
plydata = PlyData.read(f)
# In[18]:
vertices, faces=extract_data(plydata)
x, y, z=vertices.T
# Get data for a Plotly plot of the graphical object read from the ply file:
# In[19]:
data3=plotly_triangular_mesh(x,y,z, faces, intensities=standard_intensity, colorscale=pl_RdBu,
showscale=True, reversescale=False, plot_edges=True)
# Update the layout for this new plot:
# In[20]:
title="Trisurf from a PLY file
"+\
"Data Source: [1]"
# In[21]:
noaxis=dict(showbackground=False,
showline=False,
zeroline=False,
showgrid=False,
showticklabels=False,
title=''
)
fig3 = Figure(data=data3, layout=layout)
fig3['layout'].update(dict(title=title,
width=1000,
height=1000,
scene=dict(xaxis=noaxis,
yaxis=noaxis,
zaxis=noaxis,
aspectratio=dict(x=1, y=1, z=0.4),
camera=dict(eye=dict(x=1.25, y=1.25, z= 1.25)
)
)
))
# In[23]:
py.iplot(fig3, filename='Chopper-Ply-lines')
# ### Plotly plot of an isosurface ###
# An isosurface, F(x,y,z) = c, is discretized by a triangular mesh, extracted by the [Marching cubes algorithm](http://scikit-image.org/docs/dev/auto_examples/edges/plot_marching_cubes.html) from a volume given as a (M, N, P) array of doubles.
#
# The scikit image function, `measure.marching_cubes_lewiner(F, c)`
# returns the vertices and simplices of the triangulated surface.
#
# In[24]:
from skimage import measure
X,Y,Z = np.mgrid[-2:2:40j, -2:2:40j, -2:2:40j]
F = X**4 + Y**4 + Z**4 - (X**2+Y**2+Z**2)**2 + 3*(X**2+Y**2+Z**2) - 3
vertices, simplices = measure.marching_cubes_lewiner(F, 0,
spacing=(X[1,0, 0]-X[0,0,0], Y[0,1, 0]-Y[0,0,0], Z[0,0, 1]-Z[0,0,0]))[:2]
x,y,z = zip(*vertices)
# In[25]:
pl_amp=[[0.0, 'rgb(241, 236, 236)'],
[0.1, 'rgb(229, 207, 200)'],
[0.2, 'rgb(219, 177, 163)'],
[0.3, 'rgb(211, 148, 127)'],
[0.4, 'rgb(201, 119, 91)'],
[0.5, 'rgb(191, 88, 58)'],
[0.6, 'rgb(179, 55, 38)'],
[0.7, 'rgb(157, 24, 38)'],
[0.8, 'rgb(126, 13, 41)'],
[0.9, 'rgb(92, 14, 32)'],
[1.0, 'rgb(60, 9, 17)']]
# In[27]:
data4=plotly_triangular_mesh(x,y,z,simplices, intensities=standard_intensity, colorscale=pl_amp,
showscale=True, reversescale=True, plot_edges=False)
fig4 = Figure(data=data4, layout=layout)
fig4['layout'].update(dict(title='Isosurface',
width=600,
height=600,
scene=dict(camera=dict(eye=dict(x=1,
y=1,
z=1)
),
aspectratio=dict(x=1, y=1, z=1)
)))
py.iplot(fig4,filename='Isosurface')
# In[28]:
from IPython.core.display import HTML
def css_styling():
styles = open("./custom.css", "r").read()
return HTML(styles)
css_styling()