Getting Started with Meshes¶

Before you start!¶

This notebook assumes that shapeworks conda environment has been activated using conda activate shapeworks on the terminal.
See Setting Up ShapeWorks Environment to learn how to set up your environment to start using shapeworks library. Please note, the prerequisite steps will use the same code to setup the environment for this notebook and import shapeworks library.

In this notebook, you will learn:¶

How to define your dataset location and explore what is available in it
How to load a single mesh
How to convert shapeworks mesh to vtk mesh for visualization
How to visualize a surface mesh using itkwidgets
How to visualize a surface mesh using pyvista
How to visualize two meshes side-by-side using pyvista
How to visualize two meshes in the same rendering window using pyvista and itkwidgets

We will also define modular/generic helper functions as we walk through these items to reuse functionalities without duplicating code.

Prerequisites¶

Setting up shapeworks environment. See Setting Up ShapeWorks Environment. To avoid code clutter, the setup_shapeworks_env function can found in Examples/Python/setupenv.py module.
Importing shapeworks library. See Setting Up ShapeWorks Environment.

Note about `shapeworks` APIs¶

shapeworks functions are inplace, i.e., <swObject>.<function>() applies that function to the swObject data. To keep the original data unchanged, you have first to copy it to another variable before applying the function.

Notebook keyboard shortcuts¶

Esc + H: displays a complete list of keyboard shortcuts
Esc + A: insert new cell above the current cell
Esc + B: insert new cell below the current cell
Esc + D + D: delete current cell
Esc + Z: undo
Shift + enter: run current cell and move to next
To show a function's argument list (i.e., signature), use ( then shift-tab
Use shift-tab-tab to show more help for a function
To show the help of a function, use help(function) or function?
To show all functions supported by an object, use dot-tab after the variable name

Prerequisites¶

Setting up `shapeworks` environment¶

Here, we will append both your PYTHONPATH and your system PATH to setup shapeworks environment for this notebook. See Setting Up ShapeWorks Environment for more details.

In this notebook, we assume the following.

This notebook is located in Examples/Python/notebooks/tutorials
You have built shapeworks from source in build directory within the shapeworks code directory
You have built shapeworks dependencies (using build_dependencies.sh) in the same parent directory of shapeworks code

Note: If you run from a ShapeWorks installation, you don't need to set the dependencies path and the shapeworks_bin_dir would be set as ../../../../bin.

In [ ]:

# import relevant libraries 
import sys 

# add parent-parent directory (where setupenv.py is) to python path
sys.path.insert(0,'../..')

# importing setupenv from Examples/Python
import setupenv

# indicate the bin directories for shapeworks and its dependencies
shapeworks_bin_dir   = "../../../../build/bin"
dependencies_bin_dir = "../../../../../shapeworks-dependencies/bin"

# set up shapeworks environment
setupenv.setup_shapeworks_env(shapeworks_bin_dir,  
                              dependencies_bin_dir, 
                              verbose = False)

Importing `shapeworks` library¶

In [ ]:

# let's import shapeworks library to test whether shapeworks is now set
try:
    import shapeworks as sw
except ImportError:
    print('ERROR: shapeworks library failed to import')
else:
    print('SUCCESS: shapeworks library is successfully imported!!!')

1. Defining and exploring your dataset¶

Defining dataset location¶

You can download exemplar datasets from ShapeWorks data portal after you login. For new users, you can register an account for free. Please do not use an important password.

After you login, click Collections on the left panel and then use-case-data-v2. Select the dataset you would like to download by clicking on the checkbox on the left of the dataset name. See the video below.

This notebook assumes that you have downloaded ellipsoid-v1 in Examples/Python/Data. Feel free to use your own dataset.

In [ ]:

import os # for paths and mkdir

# dataset name is the folder name for your dataset
datasetName  = 'ellipsoid-v2'

# path to the dataset where we can find shape data 
# here we assume shape data are given as surface meshes
shapeDir      = '../../Data/' + datasetName + '/meshes/'
    
print('Dataset Name:     ' + datasetName)
print('Shape Directory:  ' + shapeDir)

What is available in the dataset?¶

First let's see how many shapes we have in the dataset.

File formats: For surface meshes, all vtk-supported mesh formats can be used (e.g., vtk, ply, and stl).

In [ ]:

import glob # for paths and file-directory search

# file extension for the shape data
shapeExtention = '.vtk'

# let's get a list of files for available meshes in this dataset
# * here is a wild character used to retrieve all filenames 
# in the shape directory with the file extensnion
shapeFilenames = sorted(glob.glob(shapeDir + '*' + shapeExtention)) 

print ('Number of shapes: ' + str(len(shapeFilenames)))
print('Shape files found:')
for shapeFilename in shapeFilenames:
    print('\t' + shapeFilename)

2. Loading a single mesh¶

We will select one mesh to explore for now. We will then use shapeworks Mesh class to load this surface mesh and print out its header information that includes .

In [ ]:

# select a shape by setting the shape index (in the filenames list)
shapeIdx       = 0

# the filename for the selected shape
shapeFilename  = shapeFilenames[shapeIdx]

# use shapeworks Mesh class to load it
print('Loading: ' + shapeFilename)
shapeMesh = sw.Mesh(shapeFilename)

# let's print out header information of this mesh - TODO: #828 
print('Header information: ')
print(shapeMesh)

3. Converting `shapeworks` mesh to `vtk` mesh for visualization¶

We can use python libraries such as itkwidgets and pyvista for interactive 3D visualization. These libraries support, among others, vtk data structures for images and meshes. Hence, to visualize our shapeworks mesh, we need first to convert it to a vtk data structure.

This conversion can be performed by first extracting the mesh vertices and faces arrays from the shapeworks mesh, then constructing a vtk mesh from these arrays and pyvista's wrap function to construct a vtk mesh. (see issue #825)

For now, we will save and re-read mesh for visualization, issue #825

In [ ]:

# importing pyvista
import pyvista as pv

# save mesh
shapeMesh.write('temp.vtk')

# read mesh into an itk mesh data
shapeMesh_vtk = pv.read('temp.vtk')

# remove the temp mesh file
os.remove('temp.vtk')

Defining a helper function¶

As converting between shapeworks Mesh object and vtk mesh is a step that we will need frequently, let's add a helper function for this purpose.

In [ ]:

# a helper function that converts shapeworks Mesh object to vtk mesh 
# TODO: to be modifed when #825 is addressed
def sw2vtkMesh(swMesh, verbose = False):
    
    if verbose:
        print('Header information: ')
        print(swMesh)

    # save mesh
    swMesh.write('temp.vtk')

    # read mesh into an itk mesh data
    vtkMesh = pv.read('temp.vtk')
    
    # remove the temp mesh file
    os.remove('temp.vtk')
    
    return vtkMesh

4. Visualizing surface mesh using `itkwidgets`¶

itkwidgets is a python library that supports interactive Jupyter widgets to visualize images, point sets, and meshes.

itkwidgets supports itk, vtk, and pyvista data structures. Hence, to visualize a shapeworks mesh, we need first to convert it to a vtk mesh.

In [ ]:

# convert shapeworks mesh to a vtk mesh
shapeMesh_vtk = sw2vtkMesh(shapeMesh)

Now we have an vtk mesh, we can simply visualize and interact with it using the view function of itkwidgets.

In [ ]:

# importing itkwidgets
import itkwidgets as itkw

# visualize - this is a volume rendering of the shape mesh
itkw.view(geometries = shapeMesh_vtk)

Let's enable axes and auto rotation.

In [ ]:

# visualize with axes and auto rotation
itkw.view(  geometries     = shapeMesh_vtk, 
            rotate         = True, # enable auto rotation
            axes           = True)

5. Visualizing surface mesh using `pyvista`¶

pyvista is a python library for 3D visualization and analysis. It is built on top of vtk and brings a paraview-like visualizations to notebooks. It also supports multiple rendering windows that can be linked. This feature is very useful when visualizing multiple samples from your dataset side-by-side and making them share the same camera view.

pyvista supports vtk data structures. Hence, to visualize a shapeworks mesh, we need first to convert it to a vtk mesh.

In [ ]:

# use the plot function for the vtk mesh ... notice the static view!
shapeMesh_vtk.plot()

In [ ]:

# to have an interactive visualization, 
# we need enable use_ipyvtk for this plot 
# click r to reset the view after zooming
shapeMesh_vtk.plot(use_ipyvtk = True) # enable interactive plots

In [ ]:

# Or, we can enable use_ipyvtk by default for interactive plots
pv.rcParams['use_ipyvtk'] = True 

# click r to reset the view after zooming
# click w to show wireframe and s to return back to surface
shapeMesh_vtk.plot()

6. Visualizing two meshes side-by-side using `pyvista`¶

When exploring datasets and results of different grooming (data preprocessing) steps, it is important to simultaneously visualize multiple shape samples. Here, we will learn how to visualize two meshes side-by-side and link their views using pyvista. This linking is useful to make all rendering windows share the same camera view.

Loading the second mesh and convert it to `vtk` mesh¶

First, let's select another mesh and load it.

In [ ]:

# select a shape by setting the shape index (in the filenames list)
shapeIdx2       = 1

# the filename for the selected shape
shapeFilename2  = shapeFilenames[shapeIdx2]

# use shapeworks Mesh class to load it
print('Loading: ' + shapeFilename2)
shapeMesh2 = sw.Mesh(shapeFilename2)

# let's print out header information of this mesh 
print('Header information: ')
print(shapeMesh2)

Then, let's convert this shapeworks mesh to a vtk mesh for visualization.

In [ ]:

# sw to vtk
shapeMesh2_vtk = sw2vtkMesh(shapeMesh2)

Defining `pyvista` plotter¶

Next, we will define a pyvista plotter to render multiple windows, each with a single mesh. The multiple rendering windows will be visualized as a grid of plots. Since, we have only two meshes, the grid size will be one row and two columns.

In [ ]:

# define grid size for two meshes
grid_rows  = 1
grid_cols  = 2

# define parameters that controls the plotter
is_interactive  = True   # to enable interactive plots
show_borders    = True   # show borders for each rendering window
mesh_color      = "tan"  # string or 3 item list
mesh_style      = "surface" # visualization style of the mesh. style='surface', style='wireframe', style='points'. 
show_mesh_edges = False  # show mesh edges
show_axes       = True   # show a vtk axes widget for each rendering window
show_bounds     = False  # show volume bounding box
show_all_edges  = True   # add an unlabeled and unticked box at the boundaries of plot. 
font_size       = 10     # text font size for windows
link_views      = True   # link all rendering windows so that they share same camera and axes boundaries

# define the plotter
plotter = pv.Plotter(shape    = (grid_rows, grid_cols),
                     notebook = is_interactive, 
                     border   = show_borders)

Adding meshes to the plotter and start rendering¶

Let's add the two meshes to the plotter and start the viz fun!

In [ ]:

# add the first mesh
plotter.subplot(0, 0)
plotter.add_mesh(shapeMesh_vtk, 
                 color      = mesh_color, 
                 style      = mesh_style,
                 show_edges = show_mesh_edges)

if show_axes:
    plotter.show_axes()
    
if show_bounds:
    plotter.show_bounds(all_edges = show_all_edges)

# add a text to this subplot to indicate which mesh is being visualized
meshFilename = shapeFilenames[shapeIdx].split('/')[-1] 
shapeName    = meshFilename[:-len(shapeExtention)]
plotter.add_text(shapeName, font_size = font_size)

# now, add the second mesh, 
# note that we repeat the same exact code but with a different mesh 
# ---> perfect scenario to define a helper function 
# to reuse this code without having to duplicate the code    
plotter.subplot(0, 1)
plotter.add_mesh(shapeMesh2_vtk, 
                 color      = mesh_color, 
                 style      = mesh_style,
                 show_edges = show_mesh_edges)

if show_axes:
    plotter.show_axes()
    
if show_bounds:
    plotter.show_bounds(all_edges = show_all_edges)

# add a text to this subplot to indicate which mesh is being visualized
meshFilename2 = shapeFilenames[shapeIdx2].split('/')[-1] 
shapeName2    = meshFilename2[:-len(shapeExtention)]
plotter.add_text(shapeName2, font_size = font_size)

# link views
if link_views:
    plotter.link_views()  

# now, time to render our meshes
# for any rendering window, click w to show wireframe and s to return back to surface
plotter.show() #use_ipyvtk=True) #is already enabled by default using pv.rcParams

Defining a helper function¶

Let's define a helper function that adds a mesh to a pyvista plotter.

In [ ]:

def add_mesh_to_plotter( pvPlotter,      # pyvista plotter
                         vtkMesh,         # vtk mesh to be added
                         rowIdx, colIdx, # subplot row and column index
                         title = None,    # text to be added to the subplot, use None to not show text 
                         mesh_color      = "tan",  # string or 3 item list
                         mesh_style      = "surface", # visualization style of the mesh. style='surface', style='wireframe', style='points'. 
                         show_mesh_edges = False, # show mesh edges
                         opacity         = 1,
                         show_axes       = True,  # show a vtk axes widget for each rendering window
                         show_bounds     = False, # show volume bounding box
                         show_all_edges  = True,  # add an unlabeled and unticked box at the boundaries of plot. 
                         font_size       = 10     # text font size for windows
                         ):
    
    # which subplot to add the mesh to
    pvPlotter.subplot(rowIdx, colIdx)

    # add the surface mesh
    pvPlotter.add_mesh(vtkMesh, 
                       color      = mesh_color, 
                       style      = mesh_style,
                       show_edges = show_mesh_edges,
                       opacity    = opacity)

    if show_axes:
        pvPlotter.show_axes()

    if show_bounds:
        pvPlotter.show_bounds(all_edges = show_all_edges)

    # add a text to this subplot to indicate which volume is being visualized
    if title is not None:
        pvPlotter.add_text(title, font_size = font_size)

Let's test the helper functions by adding both meshes then render.

In [ ]:

# clear the plotter first
plotter.clear()  # clears all actors

# add the first mesh
add_mesh_to_plotter( plotter, shapeMesh_vtk,   
                     rowIdx = 0, colIdx = 0, 
                     title           = shapeName,
                     mesh_color      = mesh_color,
                     mesh_style      = mesh_style,
                     show_mesh_edges = show_mesh_edges, 
                     show_axes       = show_axes, 
                     show_bounds     = show_bounds, 
                     show_all_edges  = show_all_edges, 
                     font_size       = font_size)

# add the second mesh - note that we could define a loop to avoid code repetition
add_mesh_to_plotter( plotter, shapeMesh2_vtk,   
                     rowIdx = 0, colIdx = 1, 
                     title          = shapeName2,
                     mesh_color      = mesh_color,
                     mesh_style      = mesh_style,
                     show_mesh_edges = show_mesh_edges, 
                     show_axes       = show_axes, 
                     show_bounds     = show_bounds, 
                     show_all_edges  = show_all_edges, 
                     font_size       = font_size)

# link views
if link_views:
    plotter.link_views()  

# now, time to render our meshes
plotter.show() #use_ipyvtk=True is already enabled by default using pv.rcParams

7. Visualizing two meshes in the same rendering window¶

This type of visualization is useful when exploring differences between more than one mesh, e.g., when inspecting the impact of a grooming/preprocessing step or the spatial relation of multiple samples. This is also useful if your shape data contains multiple domains (or compartments) such as anatomical joints.

Using `pyvista`¶

Note that, since we have a single rendering window (view), linking views is not necessary. But, if this multi-surface visualization is used in conjuction with multiple rendering windows, linking views should be considered.

In [ ]:

# since we want to visualize the two meshes in the same rendering window
# we define grid size for two meshes as 1 x 1 grid
grid_rows  = 1
grid_cols  = 1

# define parameters that controls the plotter
is_interactive  = True   # to enable interactive plots
show_borders    = True   # show borders for each rendering window
mesh1_color     = "blue" # string or 3 item list
mesh2_color     = "tan"  # string or 3 item list
mesh_style      = "surface" # visualization style of the mesh. style='surface', style='wireframe', style='points'. 
show_mesh_edges = False  # show mesh edges
show_axes       = True   # show a vtk axes widget for each rendering window
show_bounds     = False  # show volume bounding box
show_all_edges  = True   # add an unlabeled and unticked box at the boundaries of plot. 
font_size       = 10     # text font size for windows
link_views      = True   # link all rendering windows so that they share same camera and axes boundaries

one_title       = "%s (%s), %s (%s)" % (shapeName, mesh1_color, shapeName2, mesh2_color)

# define the plotter
plotter = pv.Plotter(shape    = (grid_rows, grid_cols),
                     notebook = is_interactive, 
                     border   = show_borders) 

# add the first mesh
add_mesh_to_plotter( plotter, shapeMesh_vtk,   
                     rowIdx = 0, colIdx = 0, 
                     title           = None,
                     mesh_color      = mesh1_color,
                     mesh_style      = mesh_style,
                     show_mesh_edges = show_mesh_edges, 
                     show_axes       = show_axes, 
                     show_bounds     = show_bounds, 
                     show_all_edges  = show_all_edges, 
                     font_size       = font_size)

# add the second mesh - note that we could define a loop to avoid code repetition
add_mesh_to_plotter( plotter, shapeMesh2_vtk,   
                     rowIdx = 0, colIdx = 0, 
                     title           = one_title,
                     mesh_color      = mesh2_color,
                     mesh_style      = mesh_style,
                     show_mesh_edges = show_mesh_edges, 
                     show_axes       = show_axes, 
                     show_bounds     = show_bounds, 
                     show_all_edges  = show_all_edges, 
                     font_size       = font_size)

# now, time to render our meshes
plotter.show() #use_ipyvtk=True is already enabled by default using pv.rcParams

Using `itkwidgets`¶

In [ ]:

# visualize with axes and auto rotation
itkw.view(  geometries       = [shapeMesh_vtk, shapeMesh2_vtk], 
            geometry_colors  = ['blue', 'tan'], # this feature is broken somehow
            rotate           = True, # enable auto rotation
            axes             = True)

Getting Started with Meshes¶

Before you start!¶

In this notebook, you will learn:¶

Prerequisites¶

Note about shapeworks APIs¶

Notebook keyboard shortcuts¶

Prerequisites¶

Setting up shapeworks environment¶

Importing shapeworks library¶

1. Defining and exploring your dataset¶

Defining dataset location¶

What is available in the dataset?¶

2. Loading a single mesh¶

3. Converting shapeworks mesh to vtk mesh for visualization¶

Defining a helper function¶

4. Visualizing surface mesh using itkwidgets¶

5. Visualizing surface mesh using pyvista¶

6. Visualizing two meshes side-by-side using pyvista¶

Loading the second mesh and convert it to vtk mesh¶

Defining pyvista plotter¶

Adding meshes to the plotter and start rendering¶

Defining a helper function¶

7. Visualizing two meshes in the same rendering window¶

Using pyvista¶

Using itkwidgets¶

Note about `shapeworks` APIs¶

Setting up `shapeworks` environment¶

Importing `shapeworks` library¶

3. Converting `shapeworks` mesh to `vtk` mesh for visualization¶

4. Visualizing surface mesh using `itkwidgets`¶

5. Visualizing surface mesh using `pyvista`¶

6. Visualizing two meshes side-by-side using `pyvista`¶

Loading the second mesh and convert it to `vtk` mesh¶

Defining `pyvista` plotter¶

Using `pyvista`¶

Using `itkwidgets`¶