Getting Started with Shape Cohort Generator¶
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. - See Getting Started with Segmentations to learn how to load and visualize binary segmentations.
- See Getting Started with Exploring Segmentations to learn how to decide the grooming pipeline needed for your dataset.
In this notebook, you will learn:¶
How to use the ShapeCohortGenerator
package to generate meshes and segmentations (binary images) for synthetic shape cohorts, i.e., parameterized families of shapes.
About ShapeCohortGenerator
¶
ShapeCohortGenerator
is a python package that generates synthetic shape cohorts with groundtruth surface correspondences by varying different parameters describing such shape families.
What is a shape cohort ?¶
A shape cohort is a collection of geometric shapes that attain clear differences in shape; however, they share common characteristics that stem from the underlying mechanisms involved in their formation. For real-world shapes, e.g., anatomical structures, such common characteristics (or factor of variations) are not known in advance, hence ShapeWorks discovers such factors of variations directly from surface meshes or binary segmentations of such shapes. ShapeCohortGenerator
uses the true factors of variations known for synthetic shapes that are analytically parameterized.
Why ShapeCohortGenerator
?¶
We require a shape population dataset to run the shape modeling workflow. Each population dataset requires unique grooming steps. Developing and testing complicated grooming pipelines for large-scale datasets can consume a lot of computational resources and time. Hence, having a few toy datasets, which are lightweight and robust in variability can make this development and debugging process easier and simpler.These cohorts can also be used to test the optimization workflow.
What families of shape can be generated by ShapeCohortGenerator
?¶
ShapeCohortGenerator
currently supports two families of synthetic shapes, namely ellipsoids and supershapes.
Ellipsoids An ellipsoid is symmetrical about three mutually perpendicular axes that intersect at the center. If a, b, and c are the principal semiaxes, the general equation of such an ellipsoid is $$\frac{x^2}{a^2} + \frac{y^2}{b^2} + \frac{z^2}{c^2} = 1$$
Supershapes are an extension of superellipses that can exhibit variable symmetry as well as asymmetry. Supershapes can be described through a single equation, the so-called superformula, that parametrizes a wide variety of shapes, including geometric primitives. The superformula is given by :
$$ r(\theta) = \left[ \left| \frac{1}{a} \cos \left( \frac{m\theta}{4} \right) \right|^{n_2} + \left| \frac{1}{b} \sin \left(\frac{m\theta}{4} \right) \right|^{n_3} \right]^{-\frac{1}{n_1}} $$Unlike superellipses, supershapes need not to be symmetric; the parameter $m$ controls the rotational symmetry. The values of $a$ and $b$ control the size, and the exponents $n_1,n_2$ and $n_3$ control the curvature of the sides. The superformula can produce a wide range of shapes,including many shapes found in nature.
The ShapeCohortGenerator
package allows the user to specify the rotational symmetry $m$ and the size. The values of $n_1,n_2$ and $n_3$ are randomly selected to creates shapes with different curvatures. Examples of these supershapes with different $m$ values can be seen below.
What you can do with ShapeCohortGenerator
?¶
The ShapeCohortGenerator
package can be used to generate collections of ellipsoids or supershapes, where the user can control the number of shapes in the cohort and the variability of the members of the cohort.
Each cohort will have mesh data (vtk and ply formats) and segmentation image data (nrrd format). These cohorts generated by the package can be directly run with ShapeWorks
. Generating these cohorts in the Output
folder would be a good way to start.
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. - Helper functions for segmentations. See Getting Started with Segmentations and Getting Started with Exploring Segmentations.
- Helper functions for meshes. See Getting Started with Meshes.
- Helper functions for visualization. See Getting Started with Segmentations, Getting Started with Meshes, and Getting Started with Exploring Segmentations.
- Defining your dataset location. See Getting Started with Exploring Segmentations.
- Loading your dataset. See Getting Started with Exploring Segmentations.
- Defining parameters for
pyvista
plotter. See Getting Started with Exploring Segmentations. - Tentative grooming pipeline. See Getting Started with Exploring Segmentations.
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
Supershapes Examples
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 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 shapeworks_bin_dir
# 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 = True)
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!!!')
Importing ShapeCohortGen
library¶
To use this package, first a generator
is defined, then generate()
is called that generates shapes in mesh format (both .ply and .vtk). Then segmentations (binary image) and images (synthetic intensities that mimic imaging data for real shapes) can be created from those meshes.
Each generator has three functions:
generate()
for mesh generation (function specific to generator type)generate_segmentations()
for segmentation generation based on meshes (general function shared by all generator types)generate_images()
for image generation based on segmentations (general function shared by all generator types)
# let's import ShapeCohortGen library
try:
import ShapeCohortGen
except ImportError:
print('ERROR: ShapeCohortGen library failed to import')
else:
print('SUCCESS: ShapeCohortGen library is successfully imported!!!')
# importing relevant libraries
import os
# 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
Misc helper functions¶
# a helper function to get list of files with specific extensions from a given file list
def get_file_with_ext(file_list,extension):
extList =[]
for file in file_list:
ext = file.split(".")[-1]
if(ext==extension):
extList.append(file)
extList = sorted(extList)
return extList
Image helper functions¶
# importing relevant libraries
import pyvista as pv
import numpy as np
import itk
# a helper function that converts shapeworks Image object to vtk image
def sw2vtkImage(swImg, verbose = False):
# get the numpy array of the shapeworks image
array = swImg.toArray()
# the numpy array needs to be permuted to match the shapeworks image dimensions
array = np.transpose(array,(2,1,0))
# converting a numpy array to a vtk image using pyvista's wrap function
vtkImg = pv.wrap(array)
if verbose:
print('shapeworks image header information: ')
print(swImg)
print('\nvtk image header information: ')
print(vtkImg)
return vtkImg
def sw2itkImage(swImg):
print(swImg)
array = swImg.toArray()
itkImg = itk.GetImageFromArray(array)
return itkImg
Helper functions for visualization¶
# importing itkwidgets to visualize single segmentations
import itkwidgets as itkw
import numpy
# itkwidgets.view returns a Viewer object. And, the IPython Jupyter kernel
# displays the last return value of a cell by default. So we have to use the display function
# to be able to call itkwidgets within a function and if statements
from IPython.display import display
# enable use_ipyvtk by default for interactive plots
pv.rcParams['use_ipyvtk'] = True
# a helper function that addes a vtk image to a pyvista plotter
def add_volume_to_plotter( pvPlotter, # pyvista plotter
vtkImg, # vtk image 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
shade_volumes = True, # use shading when performing volume rendering
color_map = "coolwarm", # color map for volume rendering, e.g., 'bone', 'coolwarm', 'cool', 'viridis', 'magma'
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 volume to
pvPlotter.subplot(rowIdx, colIdx)
# add the volume
pvPlotter.add_volume(vtkImg,
shade = shade_volumes,
cmap = color_map)
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)
# helper functions to define the best grid size for subplots
def postive_factors(num_samples):
factors = []
for whole_number in range(1, num_samples + 1):
if num_samples % whole_number == 0:
factors.append(whole_number)
return factors
def num_subplots(num_samples):
factors = postive_factors(num_samples)
cols = min(int(np.ceil(np.sqrt(num_samples))),max(factors))
rows = int(np.ceil(num_samples/cols))
return rows, cols
# helper function to add and plot a list of volumes
def plot_volumes(volumeList, # list of shapeworks images to be visualized
volumeNames = None, # list of strings of same size as shape list used to add text for each plot window, use None to not show text per window
use_same_window = False, # plot using multiple rendering windows if false
is_interactive = True, # to enable interactive plots
show_borders = True, # show borders for each rendering window
shade_volumes = True, # use shading when performing volume rendering
color_map = "coolwarm", # color map for volume rendering, e.g., 'bone', 'coolwarm', 'cool', 'viridis', 'magma'
show_axes = True, # show a vtk axes widget for each rendering window
show_bounds = True, # 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
):
num_samples = len(volumeList)
if volumeNames is not None:
if use_same_window and (len(volumeNames) > 1):
print('A single title needed when all volumes are to be displayed on the same window')
return
elif (not use_same_window) and (len(volumeNames) != num_samples):
print('volumeNames list is not consistent with number of samples')
return
if use_same_window:
grid_rows, grid_cols = 1, 1
else:
# define grid size for the given number of samples
grid_rows, grid_cols = num_subplots(num_samples)
# define the plotter
plotter = pv.Plotter(shape = (grid_rows, grid_cols),
notebook = is_interactive,
border = show_borders)
# add the given volume list (one at a time) to the plotter
for volumeIdx in range(num_samples):
# which window to add the current volume
if use_same_window:
rowIdx, colIdx = 0, 0
titleIdx = 0
else:
idUnraveled = np.unravel_index(volumeIdx, (grid_rows, grid_cols))
rowIdx, colIdx = idUnraveled[0], idUnraveled[1]
titleIdx = volumeIdx
# which title to use
if volumeNames is not None:
volumeName = volumeNames[titleIdx]
else:
volumeName = None
# convert sw image to vtk image
if type(volumeList[volumeIdx]) == sw.Image:
volume_vtk = sw2vtkImage(volumeList[volumeIdx],
verbose = False)
else:
volume_vtk = volumeList[volumeIdx]
# add the current volume
add_volume_to_plotter( plotter, volume_vtk,
rowIdx = rowIdx, colIdx = colIdx,
title = volumeName,
shade_volumes = shade_volumes,
color_map = color_map,
show_axes = show_axes,
show_bounds = show_bounds,
show_all_edges = show_all_edges,
font_size = font_size)
# link views
if link_views and (not use_same_window):
plotter.link_views()
# now, time to render our volumes
plotter.show(use_ipyvtk = is_interactive)
Defining parameters for pyvista
plotter¶
# define parameters that controls the plotter
# common for volumes and meshes visualization
is_interactive = True # to enable interactive plots
show_borders = True # show borders for each rendering window
show_axes = True # show a vtk axes widget for each rendering window
show_bounds = True # 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
# for volumes
shade_volumes = True # use shading when performing volume rendering
color_map = 'coolwarm' # color map for volume rendering, e.g., 'bone', 'coolwarm', 'cool', 'viridis', 'magma'
# for meshes
meshes_color = 'tan' # color to be used for meshes (can be a list with the same size as meshList if different colors are needed)
mesh_style = 'surface' # visualization style of the mesh. style='surface', style='wireframe', style='points'.
show_mesh_edges = False # show mesh edges
Generating an Ellipsoid Cohort¶
Step 1: Initalize Ellipsoid Generator¶
Here, we will initialize an ellipsoid cohort generator. The output directory needs to be specified, otherwise an output directory will automatically generated.
Arguments:
out_dir
: path where the dataset should be saved
Datatype :string
Default value :current_directory/generated_ellipsoid_cohort/
out_dir = "../Output/Generated_Ellipsoids/"
ellipsoid_generator = ShapeCohortGen.EllipsoidCohortGenerator(out_dir)
Step 2: Generate Meshes¶
For the ellipsoid mesh generation, you can specify the following arguments:
num_samples
: number of samples in the cohort(dataset)
Datatype :int
Default value : 3
randomize_center
: randomizes the centers for ellipsoid mesh generation if set toTrue
Datatype :bool
Defaut value :True
randomize_rotation
: randomizes the orientation of the ellipsoid if set toTrue
Datatype :bool
Defaut value :True
randomize_x_radius
: randomizes the radius of the ellipsoid along x-axis if set toTrue
or else the value is fixed as 20 for all ellipsoids
Datatype :bool
Defaut value :True
randomize_y_radius
: randomizes the radius of the ellipsoid along y-axis if set toTrue
or else the value is fixed as 10 for all ellipsoids
Datatype :bool
Defaut value :True
randomize_z_radius
: randomizes the radius of the ellipsoid along z-axis if set toTrue
or else the value is fixed as 10 for all ellipsoids
Datatype :bool
Defaut value :True
num_samples = 8
meshFiles = ellipsoid_generator.generate(num_samples)
# get all the .vtk files for plotting
VTKmeshFiles = get_file_with_ext(meshFiles,'vtk')
print(VTKmeshFiles)
vtkMeshList = []
for i in range(len(VTKmeshFiles)):
shapeMesh_vtk = sw2vtkMesh(sw.Mesh(VTKmeshFiles[i]), verbose=False)
vtkMeshList.append(shapeMesh_vtk)
Let's visualize the generated meshes using itkwidgets
.
# visualize with axes and auto rotation
itkw.view( geometries = vtkMeshList,
rotate = True, # enable auto rotation
axes = True)
Step 3: Generate Segmentations¶
For segmentation generation, you can specify the following arguments:
randomize_size
: randomize the size of the images to include more background if set toTrue
Datatype :bool
Defaut value :True
spacing
: set the spacing of the segmentation image
Datatype:list
Default value:[1,1,1]
allow_on_boundary
: If set toTrue
,randomly selects 20% samples and ensure that the shapes are touching two random selected axes out of[x,y,z]
Datatype :bool
Defaut value :True
segFiles = ellipsoid_generator.generate_segmentations()
Let's visualize the generated segmentations.
shapeSegList = []
shapeNames = []
for segFile in segFiles:
shapeSegList.append(sw.Image(segFile))
shapeNames.append(segFile.split('/')[-1])
print(shapeNames)
plot_volumes(shapeSegList,
volumeNames = shapeNames,
is_interactive = is_interactive,
show_borders = show_borders,
shade_volumes = shade_volumes,
show_axes = show_axes,
show_bounds = show_bounds,
show_all_edges = show_all_edges,
font_size = font_size,
link_views = True ) #link_views)
Step 4: Generate Images - Turning segmentations into non-binary images¶
For the image generation, a Gaussian distribution is used to define foreground and background pixels values and a blur factor is used to blur the boundary with a Gaussian filter. You can specify the following arguments:
blur_factor
: size of Gaussian filter to use for boundary blurring
Datatype :int
Defaut value :1
foreground_mean
: mean of the foreground pixel value distribution
Datatype:int
Default value:180
foreground_var
: variance of the foreground pixel value distribution
Datatype :int
Defaut value :30
background_mean
: mean of the background pixel value distribution
Datatype:int
Default value:80
background_var
: variance of the foreground pixel value distribution
Datatype :int
Defaut value :30
imageFiles = ellipsoid_generator.generate_images()
Let's compare a segmentation to it's corresponding image.
print("Segmentation:")
seg0 = sw.Image(segFiles[0])
itkw.view(sw2itkImage(seg0))
print("Image:")
img0 = sw.Image(imageFiles[0])
itkw.view(sw2itkImage(img0))
Generating Supershapes Cohort¶
SuperShapes are parameterized shapes that have geometry based on a given number of lobes, $m$.
Step 1: Initialize SuperShapes Generator¶
Here, we will initialize SuperShapes cohort generator. The output directory needs to be specified otherwise an output directory will automatically be generated.
Argument:
out_dir
: path where the dataset should be saved
Datatype :string
Default value : 'current_directory/generated_supershapes_cohort/'
out_dir = "../Output/Generated_Supershapes/"
ss_generator = ShapeCohortGen.SupershapesCohortGenerator(out_dir)
Step 2: Generate Meshes¶
For the supershapes mesh generation, you can specify the following arguments:
num_samples
- number of samples in the cohort(dataset)
Datatype :int
Default value : 3
randomize_center
: randomizes the centers for ellipsoid mesh generation if set toTrue
Datatype :bool
Defaut value :True
randomize_rotation
: randomizes the orientation of the ellispoids if set toTrue
Datatype :bool
Defaut value :True
m
: number of lobes supershapes should have
Datatype :int
Default value:3
size
: size of meshes (won't be more than 'size' away from center in any direction)
Datatype:int
Default value:20
num_samples = 8
meshFiles = ss_generator.generate(num_samples)
# get all the .vtk files for plotting
VTKmeshFiles = get_file_with_ext(meshFiles,'vtk')
print(VTKmeshFiles)
vtkMeshList = []
for i in range(len(VTKmeshFiles)):
shapeMesh_vtk = sw2vtkMesh(sw.Mesh(VTKmeshFiles[i]),verbose=False)
vtkMeshList.append(shapeMesh_vtk)
Let's visualize the generated meshes using itkwidgets
.
# visualize with axes and auto rotation
itkw.view( geometries = vtkMeshList,
rotate = True, # enable auto rotation
axes = True)
Step 3: Generate Segmentations¶
This is data type independent, the options are the same as they were for the ellipsoid.
segFiles = ss_generator.generate_segmentations()
Let's visualize the generated segmentations.
shapeSegList = []
shapeNames = []
for segFile in segFiles:
shapeSegList.append(sw.Image(segFile))
shapeNames.append(segFile.split('/')[-1])
print(shapeNames)
plot_volumes(shapeSegList,
volumeNames = shapeNames,
is_interactive = is_interactive,
show_borders = show_borders,
shade_volumes = shade_volumes,
show_axes = show_axes,
show_bounds = show_bounds,
show_all_edges = show_all_edges,
font_size = font_size,
link_views = True ) #link_views)
Step 4: Generate Images¶
This is also a standard function and has all the same options as listed before.
imageFiles = ss_generator.generate_images()
Let's compare a segmentation to it's corresponding image.
print("Segmentation:")
seg0 = sw.Image(segFiles[0])
itkw.view(sw2itkImage(seg0))
print("Image:")
img0 = sw.Image(imageFiles[0])
itkw.view(sw2itkImage(img0))