## Isosurface in volumetric data¶

Linear and nonlinear slices in volumetric data, as graphs of functions of two variables, were defined in this Jupyter Notebook http://nbviewer.jupyter.org/github/empet/Plotly-plots/blob/master/Plotly-Slice-in-volumetric-data.ipynb. Here we illustrate how to plot an isosurface in volumetric data.

In [10]:
import plotly.graph_objs as go
import numpy as np
from skimage import measure


We define an isosurface of equation $F(x,y,z)=x^4 + y^4 + z^4 - (x^2+y^2+z^2)^2 + 3(x^2+x^2+z^2) - 3=1.2$ in the volume $[-2, 2]\times[-2, 2]\times [-2, 2]$, on which a scalar field (like density or another physical property) is defined by $\psi(x,y, z)= −xe^−(x^2+y^22+z^2)$.

An isosurface $F(x,y,z)=c$ is plotted as a trisurf surface, having the vertices and faces (triangles) of a triangulation returned by the function skimage.measure.marching_cubes_lewiner(F,c).

The isosurface is colored with a colorscale based on the values of the scalar field, $\psi(x,y,z)$, at its points.

The following function returns the intensities at the vertices of the triangulation of the isosurface:

In [3]:
def intensity_func(x,y,z):
return -x * np.exp(-(x**2 + y**2 + z**2))

In [4]:
def plotly_triangular_mesh(vertices, faces, intensities=intensity_func, colorscale="Viridis",
showscale=False, reversescale=False, plot_edges=False):
# vertices: a numpy array of shape (n_vertices, 3)
# faces:  a numpy array of shape (n_faces, 3)
# intensities can be either a function of (x,y,z) or a list of values

x, y, z = vertices.T
I, J, K = faces.T
if hasattr(intensities, '__call__'):
intensity = intensities(x,y,z) # the intensities are computed here via the passed function,
# that returns a 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 = go.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 = go.Scatter3d(
x=Xe,
y=Ye,
z=Ze,
mode='lines',
name='',
line=dict(color= 'rgb(70,70,70)', width=1)
)

return [mesh, lines]


Define a meshgrid on our volume and the function that for the (iso)surface equation:

In [5]:
X, Y, Z = np.mgrid[-2:2:50j, -2:2:50j, -2:2:50j]
surf_eq = X**4 + Y**4 + Z**4 - (X**2+Y**2+Z**2)**2 + 3*(X**2+Y**2+Z**2) - 3


Although our 3D data is defined in $[-2,2]^3$, the function measure.marching_cubes_lewiner returns verts translated such that they belong to the parallelipiped $[0,4]^3$, provided that the spacing key in this function is the same as the spacing of voxels in the initial parallelipiped,

i.e. 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]):

In [6]:
verts, faces = measure.marching_cubes_lewiner(surf_eq, 1.2,
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]
title = 'Isosurface in volumetric data'


Now we translate the verts back in the original parallelipiped:

In [7]:
verts = verts-2

In [8]:
pl_BrBG=[[0.0, 'rgb(84, 48, 5)'],
[0.1, 'rgb(138, 80, 9)'],
[0.2, 'rgb(191, 129, 45)'],
[0.3, 'rgb(222, 192, 123)'],
[0.4, 'rgb(246, 232, 195)'],
[0.5, 'rgb(244, 244, 244)'],
[0.6, 'rgb(199, 234, 229)'],
[0.7, 'rgb(126, 203, 192)'],
[0.8, 'rgb(53, 151, 143)'],
[0.9, 'rgb(0, 101, 93)'],
[1.0, 'rgb(0, 60, 48)']]

In [11]:
data = plotly_triangular_mesh(verts, faces, colorscale=pl_BrBG,
showscale=True)

In [12]:
axis = dict(showbackground=True,
backgroundcolor="rgb(230, 230,230)",
gridcolor="rgb(255, 255, 255)",
zerolinecolor="rgb(255, 255, 255)")

noaxis = dict(visible=False)

layout = go.Layout(
title=title,
font=dict(family='Balto'),
showlegend=False,
width=800,
height=800,
scene=dict(xaxis=axis,
yaxis=axis,
zaxis=axis,
aspectratio=dict(x=1,
y=1,
z=1)
)
)

In [13]:
fig = go.Figure(data=data, layout=layout)

In [11]:
import plotly.plotly as py

from IPython.core.display import HTML