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, the setup_shapeworks_env function can found in Examples/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 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

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

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 = True)

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!!!')

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)

In [ ]:

# 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!!!')

Define helper functions¶

Mesh helper functions¶

In [ ]:

# 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¶

In [ ]:

# 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¶

In [ ]:

# 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¶

In [ ]:

# 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¶

In [ ]:

# 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/

In [ ]:

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 to True
Datatype : bool
Defaut value : True
randomize_rotation: randomizes the orientation of the ellipsoid if set to True
Datatype : bool
Defaut value : True
randomize_x_radius: randomizes the radius of the ellipsoid along x-axis if set to True 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 to True 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 to True or else the value is fixed as 10 for all ellipsoids
Datatype : bool
Defaut value : True

In [ ]:

num_samples = 8
meshFiles   = ellipsoid_generator.generate(num_samples)

In [ ]:

# 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.

In [ ]:

# 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 to True
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 to True,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

In [ ]:

segFiles = ellipsoid_generator.generate_segmentations()

Let's visualize the generated segmentations.

In [ ]:

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

In [ ]:

imageFiles = ellipsoid_generator.generate_images()

Let's compare a segmentation to it's corresponding image.

In [ ]:

print("Segmentation:")
seg0 = sw.Image(segFiles[0])
itkw.view(sw2itkImage(seg0))

In [ ]:

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/'

In [ ]:

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 to True
Datatype : bool
Defaut value : True
randomize_rotation: randomizes the orientation of the ellispoids if set to True
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

In [ ]:

num_samples = 8
meshFiles   = ss_generator.generate(num_samples)

In [ ]:

# 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.

In [ ]:

# 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.

In [ ]:

segFiles = ss_generator.generate_segmentations()

Let's visualize the generated segmentations.

In [ ]:

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.

In [ ]:

imageFiles = ss_generator.generate_images()

Let's compare a segmentation to it's corresponding image.

In [ ]:

print("Segmentation:")
seg0 = sw.Image(segFiles[0])
itkw.view(sw2itkImage(seg0))

In [ ]:

print("Image:")
img0 = sw.Image(imageFiles[0])
itkw.view(sw2itkImage(img0))

Getting Started with Shape Cohort Generator¶

Before you start!¶

In this notebook, you will learn:¶

About ShapeCohortGenerator¶

What is a shape cohort ?¶

Why ShapeCohortGenerator?¶

What families of shape can be generated by ShapeCohortGenerator?¶

What you can do with ShapeCohortGenerator?¶

Prerequisites¶

Note about shapeworks APIs¶

Notebook keyboard shortcuts¶

Prerequisites¶

Setting up shapeworks environment¶

Importing shapeworks library¶

Importing ShapeCohortGen library¶

Define helper functions¶

Mesh helper functions¶

Misc helper functions¶

Image helper functions¶

Helper functions for visualization¶

Defining parameters for pyvista plotter¶

Generating an Ellipsoid Cohort¶

Step 1: Initalize Ellipsoid Generator¶

Step 2: Generate Meshes¶

Step 3: Generate Segmentations¶

Step 4: Generate Images - Turning segmentations into non-binary images¶

Generating Supershapes Cohort¶

Step 1: Initialize SuperShapes Generator¶

Step 2: Generate Meshes¶

Step 3: Generate Segmentations¶

Step 4: Generate Images¶

About `ShapeCohortGenerator`¶

Why `ShapeCohortGenerator`?¶

What families of shape can be generated by `ShapeCohortGenerator`?¶

What you can do with `ShapeCohortGenerator`?¶

Importing `shapeworks` library¶

Importing `ShapeCohortGen` library¶

Defining parameters for `pyvista` plotter¶