{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Isosurface in volumetric data ##"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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."
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [],
"source": [
"import plotly.graph_objs as go\n",
"import numpy as np\n",
"from skimage import measure"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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)$. "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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)`. \n",
"\n",
"The isosurface is colored with a colorscale based on the values of the scalar field, $\\psi(x,y,z)$, at its points."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The following function returns the intensities at the vertices of the triangulation of the isosurface:"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"def intensity_func(x,y,z):\n",
" return -x * np.exp(-(x**2 + y**2 + z**2))"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"def plotly_triangular_mesh(vertices, faces, intensities=intensity_func, colorscale=\"Viridis\", \n",
" showscale=False, reversescale=False, plot_edges=False):\n",
" # vertices: a numpy array of shape (n_vertices, 3)\n",
" # faces: a numpy array of shape (n_faces, 3)\n",
" # intensities can be either a function of (x,y,z) or a list of values\n",
" \n",
" x, y, z = vertices.T\n",
" I, J, K = faces.T\n",
" if hasattr(intensities, '__call__'):\n",
" intensity = intensities(x,y,z) # the intensities are computed here via the passed function, \n",
" # that returns a list of vertices intensities\n",
" \n",
" elif isinstance(intensities, (list, np.ndarray)):\n",
" intensity = intensities #intensities are given in a list\n",
" else:\n",
" raise ValueError(\"intensities can be either a function or a list, np.array\")\n",
" \n",
" mesh = go.Mesh3d(x=x,\n",
" y=y,\n",
" z=z,\n",
" colorscale=colorscale, \n",
" reversescale=reversescale,\n",
" intensity= intensity,\n",
" i=I,\n",
" j=J,\n",
" k=K,\n",
" name='',\n",
" showscale=showscale\n",
" )\n",
" \n",
" \n",
" if showscale is True:\n",
" mesh.update(colorbar=dict(thickness=20, ticklen=4, len=0.75))\n",
" \n",
" if plot_edges is False: # the triangle sides are not plotted \n",
" return [mesh]\n",
" else: #plot edges\n",
" #define the lists Xe, Ye, Ze, of x, y, resp z coordinates of edge end points for each triangle\n",
" #None separates data corresponding to two consecutive triangles\n",
" tri_vertices= vertices[faces]\n",
" Xe=[]\n",
" Ye=[]\n",
" Ze=[]\n",
" for T in tri_vertices:\n",
" Xe += [T[k%3][0] for k in range(4)]+[ None]\n",
" Ye += [T[k%3][1] for k in range(4)]+[ None]\n",
" Ze += [T[k%3][2] for k in range(4)]+[ None]\n",
" \n",
" #define the lines to be plotted\n",
" lines = go.Scatter3d(\n",
" x=Xe,\n",
" y=Ye,\n",
" z=Ze,\n",
" mode='lines',\n",
" name='',\n",
" line=dict(color= 'rgb(70,70,70)', width=1) \n",
" )\n",
" \n",
" return [mesh, lines]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Define a meshgrid on our volume and the function that for the (iso)surface equation:"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"X, Y, Z = np.mgrid[-2:2:50j, -2:2:50j, -2:2:50j]\n",
"surf_eq = X**4 + Y**4 + Z**4 - (X**2+Y**2+Z**2)**2 + 3*(X**2+Y**2+Z**2) - 3 "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Although our 3D data is defined in $[-2,2]^3$, the function `measure.marching_cubes_lewiner` returns `verts` translated such that they belong\n",
"to the parallelipiped $[0,4]^3$, provided that the spacing key in this function is the same as the spacing\n",
"of voxels in the initial parallelipiped, \n",
"\n",
"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])`:"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"verts, faces = measure.marching_cubes_lewiner(surf_eq, 1.2, \n",
" spacing=(X[1,0, 0]-X[0,0,0], Y[0,1, 0]-Y[0,0,0], \n",
" Z[0,0, 1]-Z[0,0,0]))[:2]\n",
"title = 'Isosurface in volumetric data' "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now we translate the verts back in the original parallelipiped:"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"verts = verts-2"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [],
"source": [
"pl_BrBG=[[0.0, 'rgb(84, 48, 5)'],\n",
" [0.1, 'rgb(138, 80, 9)'],\n",
" [0.2, 'rgb(191, 129, 45)'],\n",
" [0.3, 'rgb(222, 192, 123)'],\n",
" [0.4, 'rgb(246, 232, 195)'],\n",
" [0.5, 'rgb(244, 244, 244)'],\n",
" [0.6, 'rgb(199, 234, 229)'],\n",
" [0.7, 'rgb(126, 203, 192)'],\n",
" [0.8, 'rgb(53, 151, 143)'],\n",
" [0.9, 'rgb(0, 101, 93)'],\n",
" [1.0, 'rgb(0, 60, 48)']]"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [],
"source": [
"data = plotly_triangular_mesh(verts, faces, colorscale=pl_BrBG, \n",
" showscale=True)\n"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [],
"source": [
"axis = dict(showbackground=True, \n",
" backgroundcolor=\"rgb(230, 230,230)\",\n",
" gridcolor=\"rgb(255, 255, 255)\", \n",
" zerolinecolor=\"rgb(255, 255, 255)\")\n",
"\n",
"noaxis = dict(visible=False)\n",
"\n",
"layout = go.Layout(\n",
" title=title, \n",
" font=dict(family='Balto'),\n",
" showlegend=False,\n",
" width=800,\n",
" height=800,\n",
" scene=dict(xaxis=axis,\n",
" yaxis=axis, \n",
" zaxis=axis, \n",
" aspectratio=dict(x=1,\n",
" y=1, \n",
" z=1)\n",
" )\n",
" )\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [],
"source": [
"fig = go.Figure(data=data, layout=layout)"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
""
]
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"import plotly.plotly as py\n",
"py.sign_in('username', 'api_key')\n",
"py.iplot(fig, filename='isosurface-volume')"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n"
],
"text/plain": [
""
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"from IPython.core.display import HTML\n",
"def css_styling():\n",
" styles = open(\"./custom.css\", \"r\").read()\n",
" return HTML(styles)\n",
"css_styling()"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.7.1"
}
},
"nbformat": 4,
"nbformat_minor": 1
}