#!/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()