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

# ## 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
py.sign_in('username', 'api_key')
py.iplot(fig, filename='isosurface-volume')


# In[2]:


from IPython.core.display import HTML
def  css_styling():
    styles = open("./custom.css", "r").read()
    return HTML(styles)
css_styling()