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 tovtk
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
anditkwidgets
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, thesetup_shapeworks_env
function can found inExamples/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 shortcutsEsc + A
: insert new cell above the current cellEsc + B
: insert new cell below the current cellEsc + D + D
: delete current cellEsc + Z
: undoShift + enter
: run current cell and move to next- To show a function's argument list (i.e., signature), use
(
thenshift-tab
- Use
shift-tab-tab
to show more help for a function - To show the help of a function, use
help(function)
orfunction?
- 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
.
# 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¶
# 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.
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).
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 .
# 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
# 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.
# 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.
# 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
.
# 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.
# 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.
# use the plot function for the vtk mesh ... notice the static view!
shapeMesh_vtk.plot()
# 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
# 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.
# 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.
# 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.
# 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!
# 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.
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.
# 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.
# 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
¶
# 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)