diff --git a/Examples/Python/notebooks/tutorials/getting-started-with-grooming-segmentations.ipynb b/Examples/Python/notebooks/tutorials/getting-started-with-grooming-segmentations.ipynb deleted file mode 100644 index 3e5391b98d..0000000000 --- a/Examples/Python/notebooks/tutorials/getting-started-with-grooming-segmentations.ipynb +++ /dev/null @@ -1,2923 +0,0 @@ -{ - "cells": [ - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# Getting Started with Grooming Segmentations " - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Before you start!\n", - "\n", - "- This [notebook](getting-started-with-grooming-segmentations.ipynb) assumes that shapeworks conda environment has been activated using `conda activate shapeworks` on the terminal.\n", - "- See [Setting Up ShapeWorks Environment](setting-up-shapeworks-environment.ipynb) 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.\n", - "- See [Getting Started with Segmentations](getting-started-with-segmentations.ipynb) to learn how to load and visualize binary segmentations.\n", - "- See [Getting Started with Exploring Segmentations](getting-started-with-exploring-segmentations.ipynb) to learn how to decide the grooming pipeline needed for your dataset.\n", - "\n", - "\n", - "## In this notebook, you will learn:\n", - "\n", - "1. How to resample segmentations to have smaller and isotropic voxel spacing \n", - "2. How to rigidly align segmentations to factor out global transformations \n", - "3. How to crop segmentations to remove unnecessary background voxels that might increase the memory footprint when optimizing the shape model and pad segmentations to create more room along each dimension for correspondences optimization \n", - "4. How to convert segmentations to smooth signed distance transforms as numerically stable inputs for correspondences optimization \n", - "5. How to run `ShapeWorksStudio` and use the groomed segmentations to optimize its shape model\n", - "\n", - "\n", - "We will also define modular/generic helper functions as we walk through these items to reuse functionalities without duplicating code.\n", - "\n", - "## Prerequisites\n", - "\n", - "- Setting up `shapeworks` environment. See [Setting Up ShapeWorks Environment](setting-up-shapeworks-environment.ipynb). To avoid code clutter, the `setup_shapeworks_env` function can found in `Examples/Python/setupenv.py` module.\n", - "- Importing `shapeworks` library. See [Setting Up ShapeWorks Environment](setting-up-shapeworks-environment.ipynb).\n", - "- Helper functions for segmentations. See [Getting Started with Segmentations](getting-started-with-segmentations.ipynb) and [Getting Started with Exploring Segmentations](getting-started-with-exploring-segmentations.ipynb).\n", - "- Helper functions for meshes. See [Getting Started with Meshes](getting-started-with-meshes.ipynb).\n", - "- Helper functions for visualization. See [Getting Started with Segmentations](getting-started-with-segmentations.ipynb), [Getting Started with Meshes](getting-started-with-meshes.ipynb), and [Getting Started with Exploring Segmentations](getting-started-with-exploring-segmentations.ipynb).\n", - "- Defining your dataset location. See [Getting Started with Exploring Segmentations](getting-started-with-exploring-segmentations.ipynb).\n", - "- Loading your dataset. See [Getting Started with Exploring Segmentations](getting-started-with-exploring-segmentations.ipynb).\n", - "- Defining parameters for `pyvista` plotter. See [Getting Started with Exploring Segmentations](getting-started-with-exploring-segmentations.ipynb).\n", - "- Tentative grooming pipeline. See [Getting Started with Exploring Segmentations](getting-started-with-exploring-segmentations.ipynb). \n", - "\n", - "## Note about `shapeworks` APIs\n", - "\n", - "shapeworks functions are inplace, i.e., `.()` 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.\n", - "\n", - "## Notebook keyboard shortcuts\n", - "\n", - "- `Esc + H`: displays a complete list of keyboard shortcuts\n", - "- `Esc + A`: insert new cell above the current cell\n", - "- `Esc + B`: insert new cell below the current cell\n", - "- `Esc + D + D`: delete current cell\n", - "- `Esc + Z`: undo\n", - "- `Shift + enter`: run current cell and move to next\n", - "- To show a function's argument list (i.e., signature), use `(` then `shift-tab`\n", - "- Use `shift-tab-tab` to show more help for a function\n", - "- To show the help of a function, use `help(function)` or `function?`\n", - "- To show all functions supported by an object, use `dot-tab` after the variable name\n", - "\n" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Prerequisites\n", - "\n", - "### Setting up `shapeworks` environment \n", - "\n", - "Here, we will append both your `PYTHONPATH` and your system `PATH` to setup shapeworks environment for this notebook. See [Setting Up ShapeWorks Environment](setting-up-shapeworks-environment.ipynb) for more details.\n", - "\n", - "In this notebook, we assume the following.\n", - "\n", - "- This notebook is located in `Examples/Python/notebooks/tutorials`\n", - "- You have built shapeworks from source in `build` directory within the shapeworks code directory\n", - "\n", - "**Note:** If you run from a ShapeWorks installation, you don't need to set the `shapeworks_bin_dir`" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# import relevant libraries \n", - "import sys \n", - "\n", - "# add parent-parent directory (where setupenv.py is) to python path\n", - "sys.path.insert(0,'../..')\n", - "\n", - "# importing setupenv from Examples/Python\n", - "import setupenv\n", - "\n", - "# indicate the bin directories for shapeworks and its dependencies\n", - "shapeworks_bin_dir = None # default\n", - "#shapeworks_bin_dir = \"../../../../build/bin\"\n", - "\n", - "# set up shapeworks environment\n", - "setupenv.setup_shapeworks_env(shapeworks_bin_dir, \n", - " verbose = False)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Importing `shapeworks` library" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# let's import shapeworks library to test whether shapeworks is now set\n", - "try:\n", - " import shapeworks as sw\n", - "except ImportError:\n", - " print('ERROR: shapeworks library failed to import')\n", - "else:\n", - " print('SUCCESS: shapeworks library is successfully imported!!!')" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Helper functions for IO" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# a helper function that saves a list of shapeworks images in a directory\n", - "# this could be used to save final and intermediate results (if needed)\n", - "def save_images(outDir, # path to the directory where we want to save the images\n", - " swImageList, # list of shapeworks images to be saved\n", - " swImageNames, # list of image names to be used as filenames\n", - " extension = 'nrrd',\n", - " compressed = False, # use false to load in paraview \n", - " verbose = True):\n", - "\n", - " if (len(swImageList) != len(swImageNames)):\n", - " print('swImageNames list is not consistent with number of images in swImageList')\n", - " return\n", - "\n", - " # create the output directory in case it does not exist\n", - " if not os.path.exists(outDir):\n", - " os.makedirs(outDir)\n", - "\n", - " for curImg, curName in zip(swImageList, swImageNames):\n", - " filename = Path(outDir + curName + '.' + extension)\n", - " if verbose:\n", - " print('Writing: ' + str(filename))\n", - " if not os.path.exists(os.path.dirname(filename)):\n", - " os.makedirs(os.path.dirname(filename))\n", - " curImg.write(str(filename), compressed=compressed) " - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Helper functions for segmentations" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# importing relevant libraries\n", - "import pyvista as pv\n", - "import numpy as np\n", - "\n", - "# a helper function that converts shapeworks Image object to vtk image\n", - "def sw2vtkImage(swImg, verbose = False):\n", - " \n", - " # get the numpy array of the shapeworks image\n", - " array = swImg.toArray()\n", - " \n", - " # the numpy array needs to be permuted to match the shapeworks image dimensions\n", - " array = np.transpose(array,(2,1,0))\n", - " \n", - " # converting a numpy array to a vtk image using pyvista's wrap function\n", - " vtkImg = pv.wrap(array)\n", - " \n", - " if verbose:\n", - " print('shapeworks image header information: ')\n", - " print(swImg)\n", - "\n", - " print('\\nvtk image header information: ')\n", - " print(vtkImg) \n", - " \n", - " return vtkImg" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Helper functions for meshes" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# importing relevant libraries\n", - "import os\n", - "\n", - "# a helper function that converts shapeworks Mesh object to vtk mesh \n", - "# TODO: to be modifed when #825 is addressed\n", - "def sw2vtkMesh(swMesh, verbose = False):\n", - " \n", - " if verbose:\n", - " print('Header information: ')\n", - " print(swMesh)\n", - "\n", - " # save mesh\n", - " swMesh.write('temp.vtk')\n", - "\n", - " # read mesh into an itk mesh data\n", - " vtkMesh = pv.read('temp.vtk')\n", - " \n", - " # remove the temp mesh file\n", - " os.remove('temp.vtk')\n", - " \n", - " return vtkMesh" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Helper functions for visualization" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# importing itkwidgets to visualize single segmentations\n", - "import itkwidgets as itkw\n", - "\n", - "# itkwidgets.view returns a Viewer object. And, the IPython Jupyter kernel \n", - "# displays the last return value of a cell by default. So we have to use the display function\n", - "# to be able to call itkwidgets within a function and if statements\n", - "from IPython.display import display\n", - "\n", - "# enable use_ipyvtk by default for interactive plots\n", - "pv.rcParams['use_ipyvtk'] = True \n", - " \n", - "# a helper function that addes a vtk image to a pyvista plotter\n", - "def add_volume_to_plotter( pvPlotter, # pyvista plotter\n", - " vtkImg, # vtk image to be added\n", - " rowIdx, colIdx, # subplot row and column index\n", - " title = None, # text to be added to the subplot, use None to not show text \n", - " shade_volumes = True, # use shading when performing volume rendering\n", - " color_map = \"coolwarm\", # color map for volume rendering, e.g., 'bone', 'coolwarm', 'cool', 'viridis', 'magma'\n", - " show_axes = True, # show a vtk axes widget for each rendering window\n", - " show_bounds = False, # show volume bounding box\n", - " show_all_edges = True, # add an unlabeled and unticked box at the boundaries of plot. \n", - " font_size = 10 # text font size for windows\n", - " ):\n", - " \n", - " # which subplot to add the volume to\n", - " pvPlotter.subplot(rowIdx, colIdx)\n", - " \n", - " # add the volume\n", - " pvPlotter.add_volume(vtkImg, \n", - " shade = shade_volumes, \n", - " cmap = color_map)\n", - "\n", - " if show_axes:\n", - " pvPlotter.show_axes()\n", - "\n", - " if show_bounds:\n", - " pvPlotter.show_bounds(all_edges = show_all_edges)\n", - "\n", - " # add a text to this subplot to indicate which volume is being visualized\n", - " if title is not None:\n", - " pvPlotter.add_text(title, font_size = font_size)\n", - " \n", - "# a helper function that adds a mesh to a `pyvista` plotter.\n", - "def add_mesh_to_plotter( pvPlotter, # pyvista plotter\n", - " vtkMesh, # vtk mesh to be added\n", - " rowIdx, colIdx, # subplot row and column index\n", - " title = None, # text to be added to the subplot, use None to not show text \n", - " mesh_color = \"tan\", # string or 3 item list\n", - " mesh_style = \"surface\", # visualization style of the mesh. style='surface', style='wireframe', style='points'. \n", - " show_mesh_edges = False, # show mesh edges\n", - " opacity = 1,\n", - " show_axes = True, # show a vtk axes widget for each rendering window\n", - " show_bounds = False, # show volume bounding box\n", - " show_all_edges = True, # add an unlabeled and unticked box at the boundaries of plot. \n", - " font_size = 10 # text font size for windows\n", - " ):\n", - " \n", - " # which subplot to add the mesh to\n", - " pvPlotter.subplot(rowIdx, colIdx)\n", - "\n", - " # add the surface mesh\n", - " pvPlotter.add_mesh(vtkMesh,\n", - " color = mesh_color, \n", - " style = mesh_style,\n", - " show_edges = show_mesh_edges,\n", - " opacity = opacity)\n", - "\n", - " if show_axes:\n", - " pvPlotter.show_axes()\n", - "\n", - " if show_bounds:\n", - " pvPlotter.show_bounds(all_edges = show_all_edges)\n", - "\n", - " # add a text to this subplot to indicate which volume is being visualized\n", - " if title is not None:\n", - " pvPlotter.add_text(title, font_size = font_size)\n", - " \n", - "# helper functions to define the best grid size for subplots\n", - "def postive_factors(num_samples):\n", - " factors = []\n", - " \n", - " for whole_number in range(1, num_samples + 1):\n", - " if num_samples % whole_number == 0:\n", - " factors.append(whole_number)\n", - " \n", - " return factors\n", - "\n", - "def num_subplots(num_samples):\n", - " factors = postive_factors(num_samples)\n", - " cols = min(int(np.ceil(np.sqrt(num_samples))),max(factors))\n", - " rows = int(np.ceil(num_samples/cols))\n", - " \n", - " return rows, cols\n", - "\n", - "# helper function to add and plot a list of volumes\n", - "def plot_volumes(volumeList, # list of shapeworks images to be visualized\n", - " 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 \n", - " use_same_window = False, # plot using multiple rendering windows if false\n", - " is_interactive = True, # to enable interactive plots\n", - " show_borders = True, # show borders for each rendering window\n", - " shade_volumes = True, # use shading when performing volume rendering\n", - " color_map = \"coolwarm\", # color map for volume rendering, e.g., 'bone', 'coolwarm', 'cool', 'viridis', 'magma'\n", - " show_axes = True, # show a vtk axes widget for each rendering window\n", - " show_bounds = True, # show volume bounding box\n", - " show_all_edges = True, # add an unlabeled and unticked box at the boundaries of plot. \n", - " font_size = 10, # text font size for windows\n", - " link_views = True # link all rendering windows so that they share same camera and axes boundaries\n", - " ):\n", - " \n", - " num_samples = len(volumeList)\n", - " \n", - " if volumeNames is not None:\n", - " if use_same_window and (len(volumeNames) > 1):\n", - " print('A single title needed when all volumes are to be displayed on the same window')\n", - " return\n", - " elif (not use_same_window) and (len(volumeNames) != num_samples):\n", - " print('volumeNames list is not consistent with number of samples')\n", - " return\n", - " \n", - " if use_same_window:\n", - " grid_rows, grid_cols = 1, 1\n", - " else:\n", - " # define grid size for the given number of samples\n", - " grid_rows, grid_cols = num_subplots(num_samples)\n", - "\n", - " # define the plotter\n", - " plotter = pv.Plotter(shape = (grid_rows, grid_cols),\n", - " notebook = is_interactive,\n", - " border = show_borders) \n", - " \n", - " # add the given volume list (one at a time) to the plotter\n", - " for volumeIdx in range(num_samples):\n", - " \n", - " # which window to add the current volume\n", - " if use_same_window:\n", - " rowIdx, colIdx = 0, 0\n", - " titleIdx = 0\n", - " else:\n", - " idUnraveled = np.unravel_index(volumeIdx, (grid_rows, grid_cols))\n", - " rowIdx, colIdx = idUnraveled[0], idUnraveled[1]\n", - " titleIdx = volumeIdx\n", - " \n", - " # which title to use\n", - " if volumeNames is not None:\n", - " volumeName = volumeNames[titleIdx]\n", - " else:\n", - " volumeName = None\n", - "\n", - " # convert sw image to vtk image\n", - " if type(volumeList[volumeIdx]) == sw.Image:\n", - " volume_vtk = sw2vtkImage(volumeList[volumeIdx], \n", - " verbose = False)\n", - " else:\n", - " volume_vtk = volumeList[volumeIdx]\n", - "\n", - " # add the current volume\n", - " add_volume_to_plotter( plotter, volume_vtk, \n", - " rowIdx = rowIdx, colIdx = colIdx,\n", - " title = volumeName,\n", - " shade_volumes = shade_volumes,\n", - " color_map = color_map,\n", - " show_axes = show_axes,\n", - " show_bounds = show_bounds,\n", - " show_all_edges = show_all_edges,\n", - " font_size = font_size)\n", - " # link views\n", - " if link_views and (not use_same_window):\n", - " plotter.link_views() \n", - "\n", - " # now, time to render our volumes\n", - " plotter.show(use_ipyvtk = is_interactive)\n", - " \n", - "# helper function to add and plot a list of meshes\n", - "def plot_meshes(meshList, # list of shapeworks meshes to be visualized\n", - " meshNames = 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 \n", - " use_same_window = False, # plot using multiple rendering windows if false\n", - " is_interactive = True, # to enable interactive plots\n", - " show_borders = True, # show borders for each rendering window\n", - " meshes_color = 'tan', # color to be used for meshes (can be a list with the same size as meshList if different colors are needed)\n", - " mesh_style = \"surface\", # visualization style of the mesh. style='surface', style='wireframe', style='points'. \n", - " show_mesh_edges = False, # show mesh edges\n", - " opacities = 1, # opacity to be used for meshes (can be a list with the same size as meshList if different opacities are needed) \n", - " show_axes = True, # show a vtk axes widget for each rendering window\n", - " show_bounds = True, # show volume bounding box\n", - " show_all_edges = True, # add an unlabeled and unticked box at the boundaries of plot. \n", - " font_size = 10, # text font size for windows\n", - " link_views = True # link all rendering windows so that they share same camera and axes boundaries\n", - " ):\n", - " \n", - " num_samples = len(meshList)\n", - " \n", - " if meshNames is not None:\n", - " if use_same_window and (len(meshNames) > 1):\n", - " print('A single title needed when all meshes are to be displayed on the same window')\n", - " return\n", - " elif (not use_same_window) and (len(meshNames) != num_samples):\n", - " print('meshNames list is not consistent with number of samples')\n", - " return\n", - " \n", - " if type(meshes_color) is not list: # single color given\n", - " meshes_color = [meshes_color] * num_samples\n", - " \n", - " if type(opacities) is not list: # single opacity given\n", - " opacities = [opacities] * num_samples\n", - " \n", - " if use_same_window:\n", - " grid_rows, grid_cols = 1, 1\n", - " else:\n", - " # define grid size for the given number of samples\n", - " grid_rows, grid_cols = num_subplots(num_samples)\n", - "\n", - " # define the plotter\n", - " plotter = pv.Plotter(shape = (grid_rows, grid_cols),\n", - " notebook = is_interactive,\n", - " border = show_borders) \n", - " \n", - " # add the given volume list (one at a time) to the plotter\n", - " for meshIdx in range(num_samples):\n", - " \n", - " # which window to add the current mesh\n", - " if use_same_window:\n", - " rowIdx, colIdx = 0, 0\n", - " titleIdx = 0\n", - " else:\n", - " idUnraveled = np.unravel_index(meshIdx, (grid_rows, grid_cols))\n", - " rowIdx, colIdx = idUnraveled[0], idUnraveled[1]\n", - " titleIdx = meshIdx\n", - " \n", - " # which title to use\n", - " if meshNames is not None:\n", - " meshName = meshNames[titleIdx]\n", - " else:\n", - " meshName = None\n", - "\n", - " # convert sw mesh to vtk mesh\n", - " if type(meshList[meshIdx]) == sw.Mesh:\n", - " mesh_vtk = sw2vtkMesh(meshList[meshIdx], \n", - " verbose = False)\n", - " else:\n", - " mesh_vtk = meshList[meshIdx]\n", - "\n", - " # add the current mesh\n", - " add_mesh_to_plotter( plotter, mesh_vtk,\n", - " rowIdx = rowIdx, colIdx = colIdx,\n", - " title = meshName,\n", - " mesh_color = meshes_color[meshIdx],\n", - " mesh_style = mesh_style,\n", - " show_mesh_edges = show_mesh_edges,\n", - " opacity = opacities[meshIdx],\n", - " show_axes = show_axes,\n", - " show_bounds = show_bounds,\n", - " show_all_edges = show_all_edges,\n", - " font_size = font_size)\n", - " \n", - " # link views\n", - " if link_views and (not use_same_window):\n", - " plotter.link_views() \n", - "\n", - " # now, time to render our meshes\n", - " plotter.show(use_ipyvtk = is_interactive)\n", - " \n", - " \n", - "# helper function to add and plot a list of meshes/volumes mix\n", - "def plot_meshes_volumes_mix(objectList, # list of shapeworks meshes to be visualized\n", - " objectsType, # list of 'vol', 'mesh' of same size as objectList\n", - " objectNames = 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 \n", - " use_same_window = False, # plot using multiple rendering windows if false\n", - " is_interactive = True, # to enable interactive plots\n", - " show_borders = True, # show borders for each rendering window\n", - " meshes_color = 'tan', # color to be used for meshes (can be a list with the same size as meshList if different colors are needed)\n", - " mesh_style = \"surface\", # visualization style of the mesh. style='surface', style='wireframe', style='points'. \n", - " shade_volumes = True, # use shading when performing volume rendering\n", - " color_map = \"coolwarm\", # color map for volume rendering, e.g., 'bone', 'coolwarm', 'cool', 'viridis', 'magma' \n", - " show_mesh_edges = False, # show mesh edges\n", - " opacities = 1, # opacity to be used for meshes (can be a list with the same size as meshList if different opacities are needed) \n", - " show_axes = True, # show a vtk axes widget for each rendering window\n", - " show_bounds = True, # show volume bounding box\n", - " show_all_edges = True, # add an unlabeled and unticked box at the boundaries of plot. \n", - " font_size = 10, # text font size for windows\n", - " link_views = True # link all rendering windows so that they share same camera and axes boundaries\n", - " ):\n", - " \n", - " num_samples = len(objectList)\n", - " \n", - " if objectNames is not None:\n", - " if use_same_window and (len(objectNames) > 1):\n", - " print('A single title needed when all objects are to be displayed on the same window')\n", - " return\n", - " elif (not use_same_window) and (len(objectNames) != num_samples):\n", - " print('objectNames list is not consistent with number of samples')\n", - " return\n", - " \n", - " if type(meshes_color) is not list: # single color given\n", - " meshes_color = [meshes_color] * num_samples\n", - " \n", - " if type(opacities) is not list: # single opacity given\n", - " opacities = [opacities] * num_samples\n", - " \n", - " if use_same_window:\n", - " grid_rows, grid_cols = 1, 1\n", - " else:\n", - " # define grid size for the given number of samples\n", - " grid_rows, grid_cols = num_subplots(num_samples)\n", - "\n", - " # define the plotter\n", - " plotter = pv.Plotter(shape = (grid_rows, grid_cols),\n", - " notebook = is_interactive,\n", - " border = show_borders) \n", - " \n", - " # add the given volume list (one at a time) to the plotter\n", - " for objectIdx in range(num_samples):\n", - " \n", - " # which window to add the current mesh\n", - " if use_same_window:\n", - " rowIdx, colIdx = 0, 0\n", - " titleIdx = 0\n", - " else:\n", - " idUnraveled = np.unravel_index(objectIdx, (grid_rows, grid_cols))\n", - " rowIdx, colIdx = idUnraveled[0], idUnraveled[1]\n", - " titleIdx = objectIdx\n", - " \n", - " # which title to use\n", - " if objectNames is not None:\n", - " objectName = objectNames[titleIdx]\n", - " else:\n", - " objectName = None\n", - "\n", - " if objectsType[objectIdx] == 'vol':\n", - " \n", - " # convert sw image to vtk image\n", - " if type(objectList[objectIdx]) == sw.Image:\n", - " object_vtk = sw2vtkImage(objectList[objectIdx],\n", - " verbose = False)\n", - " else:\n", - " object_vtk = objectList[objectIdx]\n", - " \n", - " # add the current volume\n", - " add_volume_to_plotter( plotter, object_vtk,\n", - " rowIdx = rowIdx, colIdx = colIdx,\n", - " title = objectName,\n", - " shade_volumes = shade_volumes,\n", - " color_map = color_map,\n", - " show_axes = show_axes,\n", - " show_bounds = show_bounds,\n", - " show_all_edges = show_all_edges,\n", - " font_size = font_size)\n", - "\n", - " else: # 'mesh'\n", - " # convert sw mesh to vtk image\n", - " if type(objectList[objectIdx]) == sw.Mesh:\n", - " object_vtk = sw2vtkMesh(objectList[objectIdx],\n", - " verbose = False)\n", - " else:\n", - " object_vtk = objectList[objectIdx]\n", - "\n", - " # add the current mesh\n", - " add_mesh_to_plotter( plotter, object_vtk,\n", - " rowIdx = rowIdx, colIdx = colIdx,\n", - " title = objectName,\n", - " mesh_color = meshes_color[objectIdx],\n", - " mesh_style = mesh_style,\n", - " show_mesh_edges = show_mesh_edges,\n", - " opacity = opacities[objectIdx],\n", - " show_axes = show_axes,\n", - " show_bounds = show_bounds,\n", - " show_all_edges = show_all_edges,\n", - " font_size = font_size)\n", - " \n", - " # link views\n", - " if link_views and (not use_same_window):\n", - " plotter.link_views() \n", - "\n", - " # now, time to render our mesh/volume mix\n", - " plotter.show(use_ipyvtk = is_interactive)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "\n", - "\n", - "### Defining dataset location\n", - "\n", - "You can download exemplar datasets from [ShapeWorks data portal](http://cibc1.sci.utah.edu:8080) after you login. For new users, you can [register](http://cibc1.sci.utah.edu:8080/#?dialog=register) an account for free. Please do not use an important password.\n", - "\n", - "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.\n", - "After you download the dataset zip file, make sure you unzip/extract the contents in the appropriate location.\n", - "\n", - "**This notebook assumes that you have downloaded `ellipsoid_1mode` and you have placed the unzipped folder `ellipsoid_1mode` in `Examples/Python/Data`.** Feel free to use your own dataset. \n", - "\n", - "\n", - "

" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "import glob # for paths and file-directory search\n", - "from pathlib import Path # for generating robust paths irrespective of the platform:Win/Linux/Mac\n", - "# dataset name is the folder name for your dataset\n", - "datasetName = 'ellipsoid_1mode'\n", - "\n", - "# file extension for the shape data\n", - "shapeExtention = '.nrrd'\n", - "\n", - "# path to the dataset where we can find shape data \n", - "# here we assume shape data are given as binary segmentations\n", - "shapeDir = '../../Data/' + datasetName + '/segmentations/'\n", - "\n", - "# name for this notebook to use for output directory\n", - "notebookName = 'getting-started-groom-segmentations'\n", - "\n", - "# path to the directory where we want to save results from this notebook\n", - "outDir = '../../Output/Notebooks/' + notebookName + '/' + datasetName + '/'\n", - "\n", - "# create the output directory in case it does not exist\n", - "if not os.path.exists(outDir):\n", - " os.makedirs(outDir)\n", - "\n", - "# let's get a list of files for available segmentations in this dataset\n", - "# * here is a wild character used to retrieve all filenames \n", - "# in the shape directory with the file extensnion\n", - "shapeFilenames = sorted(glob.glob(shapeDir + '*' + shapeExtention)) \n", - "\n", - "print('Dataset Name: ' + datasetName)\n", - "print('Shape Directory: ' + shapeDir)\n", - "print('Output Directory: ' + outDir)\n", - "print('Number of shapes: ' + str(len(shapeFilenames)))\n", - "print('Shape files found:')\n", - "for shapeFilename in shapeFilenames:\n", - " shapeFilename = Path(shapeFilename)\n", - " print(shapeFilename)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Loading your dataset" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# list of shape segmentations\n", - "shapeSegList = []\n", - "\n", - "# list of shape names (shape files prefixes) to be used \n", - "# for saving outputs and visualizations\n", - "shapeNames = [] \n", - "\n", - "# loop over all shape files and load individual segmentations\n", - "for shapeFilename in shapeFilenames:\n", - " print('Loading: ' + shapeFilename)\n", - " \n", - " # current shape name\n", - " shapeFilename = shapeFilename.replace(\"\\\\\",\"/\")\n", - " segFilename = shapeFilename.split('/')[-1] \n", - " shapeName = segFilename[:-len(shapeExtention)]\n", - " shapeNames.append(shapeName)\n", - " \n", - " # load segmentation\n", - " shapeSeg = sw.Image(shapeFilename)\n", - " \n", - " # append to the shape list\n", - " shapeSegList.append(shapeSeg)\n", - "\n", - "num_samples = len(shapeSegList)\n", - "print('\\n' + str(num_samples) + \n", - " ' segmentations are loaded for the ' + datasetName + ' dataset ...')" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Defining parameters for `pyvista` plotter" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# define parameters that controls the plotter\n", - "\n", - "# common for volumes and meshes visualization\n", - "is_interactive = True # to enable interactive plots\n", - "show_borders = True # show borders for each rendering window\n", - "show_axes = True # show a vtk axes widget for each rendering window\n", - "show_bounds = True # show volume bounding box\n", - "show_all_edges = True # add an unlabeled and unticked box at the boundaries of plot. \n", - "font_size = 10 # text font size for windows\n", - "link_views = True # link all rendering windows so that they share same camera and axes boundaries\n", - "\n", - "# for volumes\n", - "shade_volumes = True # use shading when performing volume rendering\n", - "color_map = 'coolwarm' # color map for volume rendering, e.g., 'bone', 'coolwarm', 'cool', 'viridis', 'magma'\n", - "\n", - "# for meshes\n", - "meshes_color = 'tan' # color to be used for meshes (can be a list with the same size as meshList if different colors are needed)\n", - "mesh_style = 'surface' # visualization style of the mesh. style='surface', style='wireframe', style='points'. \n", - "show_mesh_edges = False # show mesh edges \n", - " " - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Tentative grooming\n", - "\n", - "Hence, a tentative grooming pipeline entails the following steps: \n", - "\n", - "1. Resampling segmentations to have smaller and isotropic voxel spacing \n", - "2. Rigidly aligning shapes \n", - "3. Cropping and padding segmentations \n", - "4. Converting segmentations to smooth signed distance transforms \n", - "\n", - "\n", - "Let the fun begins!!!" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## 1. Resampling segmentations\n", - "\n", - "\n", - "This grooming step resamples all the binary volumes which in a raw setting could be in different physical spaces (different dimensions and voxel spacing). This brings all segmentations to the same voxel spacing. \n", - "\n", - "If a smaller voxel spacing is used, this improves the resolution of the segmentations and reduce the staircase effect seen in the volume rendering.\n", - "\n", - "Since image resampling entails interpolation, directly resampling binary segmentations will not result in a binary segmentation, but rather an interpolated version that does not have two distinct labels (i.e., foreground and background).\n", - "\n", - "To mitigate this behavior, we need first to convert the binary segmentations (with zero-one voxels) to a continuous-valued (gray-scale) image. This can be done by antialiasing the segmentations, which smoothes the foreground-background interface. \n", - "\n", - "Hence, the resampling pipeline for a binary segmentation includes the following steps:\n", - "\n", - "- Antialiasing the binary segmentation to convert it to a smooth continuous-valued image\n", - "- Resampling the antialiased image using the same (and possible smaller) voxel spacing for all dimensions\n", - "- Binarizing (aka thresholding) the resampled image to results in a binary segmentation with an isotropic voxel spacing\n", - "\n" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Antialiasing: which `shapeworks` API to use?" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# let's see if there's a function that antialias a shapeworks image\n", - "# use dot-tap to get a list of functions/apis available for shapeSeg\n", - "\n", - "shapeIdx = 6\n", - "shapeSeg = shapeSegList[shapeIdx]\n", - "\n", - "# found it - resample, let's see its help\n", - "help(shapeSeg.antialias)\n", - "#shapeSeg.antialias?\n" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Exploring antialiasing for a single segmentation" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# parameters for antialiasing\n", - "antialiasIterations = 50 # number of iterations to perform antialiasing\n", - "\n", - "# note that shapeworks APIs use inplace computations, so if we do shapeSeg.antialias, \n", - "# shapeSeg will have the antialiased segmentation\n", - "\n", - "# to visualize the effect of padding, let's keep the original segmentation unchanged\n", - "shapeSegAntialiased = shapeSeg.copy()\n", - "\n", - "# let's antialias this segmentation\n", - "shapeSegAntialiased.antialias(antialiasIterations)\n", - "\n", - "# note antialiasing does not change header information\n", - "print('\\nHeader information before antialiasing: ')\n", - "print(shapeSeg)\n", - "\n", - "print('\\nHeader information after antialiasing: ')\n", - "print(shapeSegAntialiased)\n", - " \n", - "plot_volumes([shapeSeg, shapeSegAntialiased], \n", - " ['before antialiasing','after antialiasing'],\n", - " use_same_window = False,\n", - " is_interactive = is_interactive,\n", - " show_borders = show_borders,\n", - " shade_volumes = shade_volumes, \n", - " color_map = color_map,\n", - " show_axes = show_axes,\n", - " show_bounds = show_bounds,\n", - " show_all_edges = show_all_edges,\n", - " font_size = font_size,\n", - " link_views = link_views)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "

" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Note that the antialiasing operation has converted a 0-1 binary volume to a continous-valued volume with voxel values ranging from -4.0 to 4.0. The zero level set (i.e., voxels with zero value) is the foreground-background interface.\n", - "\n", - "Let's extract the isosurface of this continous-valued volume and compare it with the segmentation isosurface to make sure that the shape information is not lost during antialiasing." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# extract the segmentation isosurface\n", - "shapeSegIso = shapeSeg.toMesh(isovalue = 0.5)\n", - "\n", - "# extract the antialiased isosurface\n", - "shapeSegAntialiasedIso = shapeSegAntialiased.toMesh(isovalue = 0.0)\n", - "\n", - "\n", - "plot_meshes([shapeSegIso, shapeSegAntialiasedIso],\n", - " ['before (red) and after (tan) antialiasing'],\n", - " use_same_window = True,\n", - " is_interactive = is_interactive,\n", - " show_borders = show_borders,\n", - " meshes_color = ['red', 'tan'],\n", - " opacities = [0.5, 1],\n", - " mesh_style = mesh_style,\n", - " show_mesh_edges = show_mesh_edges,\n", - " show_axes = show_axes,\n", - " show_bounds = show_bounds,\n", - " show_all_edges = show_all_edges,\n", - " font_size = font_size,\n", - " link_views = link_views)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "

" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Resampling: which `shapeworks` API to use?" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# let's see if there's a function that resamples a shapeworks image\n", - "# use dot-tap to get a list of functions/apis available for shapeSeg\n", - "\n", - "# found it - resample, let's see its help\n", - "help(shapeSegAntialiased.resample) # issue #827" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Exploring resampling for a single segmentation\n", - "\n", - "There are three options to perform resampling:\n", - "- Isotropic resampling with all dimensions having the same voxel size (set isSpacing parameter)\n", - "- Anisotropic resampling by setting spacex, spacey, and spacez parameters to indicate a different voxel size for each dimension\n", - "- Resampling given the output image dimensions by setting sizex, sizey, and sizez parameters. Voxel spacing along each dimension is thus computed based on the output image size.\n", - "\n", - "In this case, we will perform isotropic resampling. " - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# define the voxel spacing for resampling\n", - "\n", - "# let's checkout the spacing of the original segmentation\n", - "spacing = shapeSegAntialiased.spacing()\n", - "\n", - "print(spacing) # voxel spacing in the z-dimension is double that of x- and y-dimension" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# let's define the isospacing to be the minimum of the spacing of all dimensions\n", - "# this will improve the segmentation resolution along \n", - "# the z-dimension (i.e., reduce slice thickness)\n", - "# #830 - add toArray to swVectors\n", - "spacing_array = [spacing[d] for d in range(3)]\n", - "isoSpacing = min(spacing_array) # The isotropic spacing in all dimensions.\n", - "\n", - "print('isoSpacing = ' + str(isoSpacing))" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# to visualize the effect of resampling, let's keep the antialiased segmentation unchanged\n", - "shapeSegResampled = shapeSegAntialiased.copy()\n", - "\n", - "# define the interpolation type\n", - "interp = sw.InterpolationType.Linear \n", - "\n", - "# perform image resampling\n", - "#shapeSegResampled.resample(voxelSpacing, interp)\n", - "shapeSegResampled.resample([isoSpacing,isoSpacing,isoSpacing], interp)\n", - "\n", - "print('\\nHeader information before resampling: ')\n", - "print(shapeSegAntialiased)\n", - "\n", - "print('\\nHeader information after resampling: ')\n", - "print(shapeSegResampled)\n", - " \n", - "# note pyvista's plots don't reflect physical coordinates, you are seeing a larger\n", - "# 3D image (matrix)\n", - "plot_volumes([shapeSegAntialiased, shapeSegResampled],\n", - " ['before resampling','after resampling'],\n", - " is_interactive = is_interactive,\n", - " show_borders = show_borders,\n", - " shade_volumes = shade_volumes,\n", - " show_axes = show_axes,\n", - " show_bounds = show_bounds,\n", - " show_all_edges = show_all_edges,\n", - " font_size = font_size, \n", - " link_views = link_views)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "

" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# visualizing isosurfaces rather respects the physical coordinates\n", - "\n", - "# extract the resampled isosurface \n", - "shapeSegResampledIso = shapeSegResampled.toMesh(isovalue = 1e-20) # Issue #833\n", - "\n", - "plot_meshes([shapeSegIso, shapeSegResampledIso],\n", - " ['segmentation (red) and resampled (tan)'],\n", - " use_same_window = True,\n", - " is_interactive = is_interactive,\n", - " show_borders = show_borders,\n", - " meshes_color = ['red', 'tan'],\n", - " opacities = [0.5, 1],\n", - " mesh_style = mesh_style,\n", - " show_mesh_edges = show_mesh_edges,\n", - " show_axes = show_axes,\n", - " show_bounds = show_bounds,\n", - " show_all_edges = show_all_edges,\n", - " font_size = font_size,\n", - " link_views = link_views)\n" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "

" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Binarizing: which shapeworks API to use?" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# let's see if there's a function that resamples a shapeworks image\n", - "# use dot-tap to get a list of functions/apis available for shapeSeg\n", - "\n", - "# found it - resample, let's see its help\n", - "help(shapeSegResampled.binarize)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Exploring binarization for a single segmentation" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# parameters for binarization\n", - "# all voxels between minVal and maxVal are set to innerVal\n", - "minVal = 0\n", - "maxVal = 5\n", - "innerVal = 1\n", - "outerVal = 0\n", - "\n", - "# to visualize the effect of padding, let's keep the original segmentation unchanged\n", - "shapeSegResampledBin = shapeSegResampled.copy()\n", - "\n", - "shapeSegResampledBin.binarize(minVal, maxVal, innerVal, outerVal)\n", - " \n", - "plot_volumes([shapeSegResampled, shapeSegResampledBin],\n", - " ['before binarization','after binarization'],\n", - " is_interactive = is_interactive,\n", - " show_borders = show_borders,\n", - " shade_volumes = shade_volumes,\n", - " show_axes = show_axes,\n", - " show_bounds = show_bounds,\n", - " show_all_edges = show_all_edges,\n", - " font_size = font_size, \n", - " link_views = link_views)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "

" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Comparing the original segmentation and the resampled one" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "plot_meshes([shapeSegIso, shapeSegResampledBin.toMesh(0.5)],\n", - " ['segmentation (red) and resampled (tan)'],\n", - " use_same_window = True,\n", - " is_interactive = is_interactive,\n", - " show_borders = show_borders,\n", - " meshes_color = ['red', 'tan'],\n", - " opacities = [0.5, 1],\n", - " mesh_style = mesh_style,\n", - " show_mesh_edges = show_mesh_edges,\n", - " show_axes = show_axes,\n", - " show_bounds = show_bounds,\n", - " show_all_edges = show_all_edges,\n", - " font_size = font_size,\n", - " link_views = link_views)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "

" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Resampling all segmentations" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# let's now resample all segmentations, here we will use inplace update for shapeworks image\n", - "# this step will bring all segmentation to the same voxel spacing\n", - "\n", - "# parameters for resampling\n", - "iso_spacing = 1 # use 1 for memory constraints on notebooks\n", - "antialias_iterations = 50\n", - "\n", - "# all voxels between minVal and maxVal are set to innerVal\n", - "minVal = 0\n", - "maxVal = 5\n", - "innerVal = 1\n", - "outerVal = 0\n", - "\n", - "for shapeSeg, shapeName in zip(shapeSegList, shapeNames):\n", - " print('Resampling segmentation: ' + shapeName)\n", - " \n", - " # we can chain aliasing then resample then binarize\n", - " \n", - " # antialiasing binary segmentation\n", - " shapeSeg.antialias(antialias_iterations)\n", - " \n", - " # perform image resampling\n", - " shapeSeg.resample([iso_spacing, iso_spacing, iso_spacing], sw.InterpolationType.Linear)\n", - "\n", - " # binarize resampled image\n", - " # all voxels between minVal and maxVal are set to innerVal\n", - " shapeSeg.binarize(minVal, maxVal, innerVal, outerVal)\n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "shapeSeg = shapeSegList[10]\n", - "itkw.view( image = sw2vtkImage(shapeSeg),\n", - " slicing_planes = True,\n", - " axes = True,\n", - " rotate = True,\n", - " interpolation = True)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "

" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Visualizing resampled segmentations" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# plot all segmentations in the shape list as surfaces\n", - "shapeSegIsoList = []\n", - "for shapeSeg, shapeName in zip(shapeSegList, shapeNames):\n", - " print('Isosurfacing segmentation: ' + shapeName)\n", - " shapeSegIsoList.append(shapeSeg.toMesh(1e-10))" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "plot_meshes(shapeSegIsoList,\n", - " shapeNames,\n", - " use_same_window = False,\n", - " is_interactive = is_interactive, \n", - " show_borders = show_borders,\n", - " meshes_color = meshes_color,\n", - " mesh_style = mesh_style,\n", - " show_mesh_edges = show_mesh_edges,\n", - " show_axes = show_axes,\n", - " show_bounds = False,\n", - " show_all_edges = show_all_edges,\n", - " font_size = font_size, \n", - " link_views = False) # try out True to link views\n" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "

" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Saving resampled segmentations" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "save_images(outDir + 'resampled/', \n", - " shapeSegList,\n", - " shapeNames,\n", - " extension = 'nrrd',\n", - " compressed = False, # use false to load in paraview \n", - " verbose = True)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## 2. Aligning shapes \n", - "\n", - "Rigidly aligning a cohort of shapes entails removing differences across these shapes pertaining to global transformations, i.e., translation and rotation. This step requires a reference coordinate frame to align all shapes to, where one of the shapes can be selected as a reference. \n", - "\n", - "Rigid alignment (aka registration) is an optimization process that might get stuck in a bad local minima if shapes are significantly out of alignment. To bring shapes closer, we can remove translation differences using center-of-mass alignment. This factors out translations to reduce the risk of misalignment and allow for a medoid sample to be automatically selected as the reference for subsequent rigid alignment.\n", - "\n", - "Hence, the shapes alignment pipeline includes the following steps:\n", - "- Center-of-mass alignment\n", - "- Reference shape selection\n", - "- Rigid alignment\n", - "\n" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Center-of-mass alignment\n", - "\n", - "This step takes in a binary volume and translates the center of mass of the shape to the center of the 3D volume space." - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "#### Exploring center-of-mass alignment for a single segmentation" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "shapeIdx = 3\n", - "shapeSeg = shapeSegList[shapeIdx].copy()\n", - "\n", - "# computing the translation vector for this registration step, we need to\n", - "# get the image's center and the shape's center of mass.\n", - "# shapeworks image class has apis to provide both (returning shapeworks point)\n", - "help(shapeSeg.centerOfMass)\n", - "help(shapeSeg.center)\n", - "\n", - "itkw.view( image = sw2vtkImage(shapeSeg),\n", - " slicing_planes = True,\n", - " axes = True,\n", - " rotate = True,\n", - " interpolation = True)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "

" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# compute the center of mass of this segmentation\n", - "shapeCenter = shapeSeg.centerOfMass() # retruns shapeworks point\n", - "print('shapeCenter: ' + str(shapeCenter))\n", - "\n", - "# get the center of the image domain - physical coordinates of the \"image center\" (middle voxel) \n", - "imageCenter = shapeSeg.center() # retruns shapeworks point\n", - "print('imageCenter: ' + str(imageCenter))\n", - "\n", - "# image origin - physical coordinates of voxel(0,0,0)\n", - "imageOrigin = shapeSeg.origin() # retruns shapeworks point\n", - "print('imageOrigin: ' + str(imageOrigin))\n", - "\n", - "# physical coordinates of voxel (row-1,cols-1,slices-1) - end point of the image domain diagonal\n", - "imageEnd = imageOrigin + shapeSeg.size()\n", - "print('imageEnd: ' + str(imageEnd))" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# let' visualize those centers in relation to the image domain and shape/segmentation\n", - "# to better understand this relation, we need to visualize meshes and volumes in one view\n", - "shapeCenter_vtk = pv.Sphere(radius = 2, \n", - " center = (shapeCenter[0], shapeCenter[1], shapeCenter[2]))\n", - "imageCenter_vtk = pv.Sphere(radius = 2, \n", - " center = (imageCenter[0], imageCenter[1], imageCenter[2]))\n", - "\n", - "origin_vtk = pv.Sphere(radius = 2, \n", - " center = (imageOrigin[0], imageOrigin[1], imageOrigin[2]))\n", - "\n", - "end_vtk = pv.Sphere(radius = 2, \n", - " center = (imageEnd[0], imageEnd[1], imageEnd[2]))\n", - "\n", - "plot_meshes([shapeSeg.toMesh(0.5), shapeCenter_vtk, imageCenter_vtk, origin_vtk, end_vtk],\n", - " ['seg (tan), shape center (red), image center (green), image origin (blue), image end (cyan)'],\n", - " use_same_window = True,\n", - " is_interactive = is_interactive,\n", - " show_borders = show_borders,\n", - " meshes_color = ['tan', 'red', 'green', 'blue', 'cyan'],\n", - " opacities = [0.5, 1, 1, 1, 1],\n", - " mesh_style = mesh_style,\n", - " show_mesh_edges = show_mesh_edges,\n", - " show_axes = show_axes,\n", - " show_bounds = show_bounds,\n", - " show_all_edges = show_all_edges,\n", - " font_size = font_size,\n", - " link_views = link_views)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "

" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# now define the translation to move the shape center to the image center\n", - "translationVector = imageCenter - shapeCenter\n", - "print('translationVector: ' + str(translationVector))" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Applying this translation to segmentations entails interpolation due to image resampling in the new coordinate frame. Directly resampling binary segmentations will not result in a binary segmentation, but rather an interpolated version that does not have two distinct labels (i.e., foreground and background).\n", - "\n", - "To mitigate this behavior, similar to the resampling workflow, we will first antialias the segmentation to convert it to a continuous-valued image with a smooth foreground-background interface, then apply the translation, and finally binarize the translated image." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# let's first inspect the api for image translation\n", - "help(shapeSeg.translate)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# let's keep shapeSeg unchanged\n", - "shapeSegCentered = shapeSeg.copy()\n", - "\n", - "translationVector\n", - "# here we will perform antialias-translate-binarize by chaining the apis\n", - "#shapeSegCentered.antialias(antialias_iterations).translate(translationVector).binarize()\n", - "shapeSegCentered.antialias(antialias_iterations).translate(translationVector).binarize()\n", - "\n", - "# notice no change in header information\n", - "print('\\nHeader information before center of mass alignment: ')\n", - "print(shapeSeg)\n", - "\n", - "print('\\nHeader information after center of mass alignment: ')\n", - "print(shapeSegCentered)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# let's visualize the segmentation before and after center of mass alignment\n", - "# notice how the shape is moved (i.e., translated) to the center of the image\n", - "\n", - "print('shape segmentation before center of mass alignment')\n", - "display(itkw.view( image = sw2vtkImage(shapeSeg), # for orthoginal image plane\n", - " label_image = sw2vtkImage(shapeSeg), # for volume rendering segmentation\n", - " slicing_planes = True,\n", - " axes = True,\n", - " rotate = True, # enable auto rotation\n", - " interpolation = True)\n", - " )\n", - "\n", - "print('shape segmentation after center of mass alignment')\n", - "display(itkw.view( image = sw2vtkImage(shapeSegCentered), # for orthoginal image plane\n", - " label_image = sw2vtkImage(shapeSegCentered), # for volume rendering segmentation\n", - " slicing_planes = True,\n", - " axes = True,\n", - " rotate = True, # enable auto rotation\n", - " interpolation = True)\n", - " )" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "

" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "#### Center-of-mass aligning all shapes" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# let's apply center of mass alignment to all shapes\n", - "\n", - "for shapeSeg, shapeName in zip(shapeSegList, shapeNames):\n", - " print('Center of mass alignment: ' + shapeName) \n", - " \n", - " # compute the center of mass of this segmentation\n", - " shapeCenter = shapeSeg.centerOfMass() # retruns shapeworks point\n", - " #print('\\tshapeCenter: ' + str(shapeCenter))\n", - "\n", - " # get the center of the image domain\n", - " imageCenter = shapeSeg.center() # retruns shapeworks point\n", - " #print('\\timageCenter: ' + str(imageCenter))\n", - "\n", - " # now define the translation to move the shape to its center\n", - " translationVector = imageCenter - shapeCenter\n", - " #print('\\ttranslationVector: ' + str(translationVector))\n", - " \n", - " # perform antialias-translate-binarize by chaining the apis\n", - " shapeSeg.antialias(antialias_iterations).translate(translationVector).binarize()" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "#### Visualizing center-of-mass aligned segmentations" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# .. and visualize center of mass aligned dataset\n", - "plot_volumes(shapeSegList, \n", - " volumeNames = shapeNames,\n", - " is_interactive = is_interactive,\n", - " show_borders = show_borders,\n", - " shade_volumes = shade_volumes,\n", - " show_axes = show_axes,\n", - " show_bounds = show_bounds,\n", - " show_all_edges = show_all_edges,\n", - " font_size = font_size,\n", - " link_views = True ) #link_views)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "

" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "#### Saving center-of-mass aligned segmentations" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "save_images(outDir + 'com/', \n", - " shapeSegList, \n", - " shapeNames,\n", - " extension = 'nrrd',\n", - " compressed = False, # use false to load in paraview \n", - " verbose = True)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Reference shape selection\n", - "\n", - "One option for a reference is to select the shape that is closest to the mean shape, i.e., the medoid shape. Note that computing the mean shape and distance to the mean shape require that all segmentations to have the same image dimensions.\n", - "\n", - "Select the medoid shape entails the following steps:\n", - "\n", - "- Bring all segmentations to the same image size if needed \n", - "- Compute mean shape\n", - "- Compute distance to mean shape\n", - "- Select the shape sample that is closest to the mean shape" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "#### Bring all segmentations to the same image size if needed " - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# compute maximum image size in each dimension\n", - "x_dims, y_dims, z_dims = [], [], []\n", - "for shapeSeg in shapeSegList:\n", - " dims = shapeSeg.dims()\n", - " x_dims.append(dims[0])\n", - " y_dims.append(dims[1])\n", - " z_dims.append(dims[2])\n", - "\n", - "print('x_dims: ' + str(x_dims))\n", - "print('y_dims: ' + str(y_dims))\n", - "print('z_dims: ' + str(z_dims))\n", - "\n", - "# define the maximum image size in the given dataset - notice different images have different sizes\n", - "x_max_dim = max(x_dims)\n", - "y_max_dim = max(y_dims)\n", - "z_max_dim = max(z_dims) \n", - "\n", - "# this is the common size that we need to pad individual segmentations \n", - "# to compute the mean shape and distance to the mean shape\n", - "print('x_max_dim = ' + str(x_max_dim))\n", - "print('y_max_dim = ' + str(y_max_dim))\n", - "print('z_max_dim = ' + str(z_max_dim))" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "#### Compute mean shape" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# now let's compute the mean shape\n", - "meanShape = np.zeros((x_max_dim, y_max_dim, z_max_dim))\n", - "padValue = 0\n", - "\n", - "# let's keep a list for the padded ones to compute distances to the mean shape next\n", - "paddedList = [] \n", - "for shapeSeg in shapeSegList:\n", - " dims = shapeSeg.dims()\n", - "\n", - " # convert to numpy array\n", - " shapeSeg_array = shapeSeg.toArray()\n", - " \n", - " # the numpy array needs to be permuted to match the shapeworks image dimensions\n", - " shapeSeg_array = np.transpose(shapeSeg_array,(2,1,0))\n", - " \n", - " # now let's pad and still maintain the shape in the center of the image domain\n", - " x_diff = x_max_dim - dims[0]\n", - " before_x = x_diff // 2\n", - " after_x = x_max_dim - (dims[0] + before_x)\n", - " \n", - " y_diff = y_max_dim - dims[1]\n", - " before_y = y_diff // 2\n", - " after_y = y_max_dim - (dims[1] + before_y)\n", - " \n", - " z_diff = z_max_dim - dims[2]\n", - " before_z = z_diff // 2\n", - " after_z = z_max_dim - (dims[2] + before_z)\n", - " \n", - " paddedList.append( np.pad(shapeSeg_array, ( (before_x, after_x), \n", - " (before_y, after_y), \n", - " (before_z, after_z)\n", - " )\n", - " )\n", - " )\n", - " meanShape += paddedList[-1]\n", - "meanShape /= len(shapeSegList)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# ... and visualize the mean shape\n", - "meanShape_vtk = pv.wrap(meanShape)\n", - "\n", - "print('mean shape')\n", - "itkw.view( #image = meanShape_vtk, # for orthoginal image plane\n", - " label_image = meanShape_vtk, # for volume rendering segmentation\n", - " slicing_planes = True,\n", - " axes = True,\n", - " rotate = True, # enable auto rotation\n", - " interpolation = True\n", - " )" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "

" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "#### Compute distance to mean shape" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "distances = np.zeros((num_samples,))\n", - "for ii, segPadded in enumerate(paddedList):\n", - " distances[ii] = np.linalg.norm(meanShape - segPadded)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "#### Select the shape sample that is closest to the mean shape" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "import matplotlib.pyplot as plt\n", - "\n", - "plt.figure(dpi=150)\n", - "plt.bar(range(num_samples), distances)\n", - "plt.xlabel('shape #')\n", - "plt.ylabel('distance to mean shape')\n", - "\n", - "referenceIdx = distances.argmin()\n", - "plt.bar(referenceIdx, distances[referenceIdx]);" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "

" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Rigid alignment\n", - "\n", - "This step rigidly aligns each shape to the selected references. Rigid alignment involves interpolation, hence we need to convert binary segmentations to continuous-valued images." - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "#### Exploring alignment for a single segmentation\n", - "\n", - "Aligning to a reference shape entails two steps:\n", - "- computing the rigid transformation parameters that would align a segmentation to the reference shape\n", - "- applying the rigid transformation to the segmentation" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# this API takes in \n", - "# (1) target/reference image \n", - "# (2) type of transform (we want to use IterativeClosestPoint transformation for alignment)\n", - "# (3) isovalue that defines the shape's surface in the given images\n", - "# (4) number of iterations for ICP alignment\n", - "help(sw.Image.createTransform)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# first step is to compute rigid transformation parameters\n", - "\n", - "shapeIdx = 5 #0 #referenceIdx #13\n", - "shapeSeg = shapeSegList[shapeIdx].copy() # need copy to avoid inplace antialiasing\n", - "refSeg = shapeSegList[referenceIdx].copy() # need copy to avoid inplace antialiasing \n", - "\n", - "# notice the differences in origin, dims and size\n", - "print('\\nHeader information before rigid alignment: ')\n", - "print(shapeSeg)\n", - "\n", - "print('\\nHeader information of reference: ')\n", - "print(refSeg)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# here we will antialias segmentations to convert to continuous-valued images\n", - "isoValue = 1e-20\n", - "icp_iterations = 200\n", - "\n", - "# compute rigid transformation\n", - "rigidTransform = shapeSeg.antialias(antialias_iterations).createTransform(refSeg.antialias(antialias_iterations),\n", - " sw.TransformType.IterativeClosestPoint,\n", - " isoValue,\n", - " icp_iterations\n", - " )\n", - "\n", - "print(rigidTransform)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# second we apply the computed transformation\n", - "# to avoid segmentations going out of image bounds when applying the transformation\n", - "# we need the aligned shape to have the same image information \n", - "# of the reference shape for subsequent steps\n", - "shapeSegAligned = shapeSeg.copy()\n", - "shapeSegAligned.antialias(antialias_iterations)\n", - "shapeSegAligned.applyTransform(rigidTransform, \n", - " refSeg.origin(), refSeg.dims(), \n", - " refSeg.spacing(), refSeg.coordsys(), \n", - " sw.InterpolationType.Linear)\n", - "shapeSegAligned.binarize()\n", - "\n", - "# notice the change in header information\n", - "print('\\nHeader information before rigid alignment: ')\n", - "print(shapeSeg)\n", - "\n", - "print('\\nHeader information after rigid alignment: ')\n", - "print(shapeSegAligned)\n", - "\n", - "print('\\nHeader information of reference: ')\n", - "print(refSeg) " - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "plot_meshes([shapeSeg.toMesh(0.5), shapeSegAligned.toMesh(0.5), refSeg.toMesh(0.5)],\n", - " ['seg (tan), aligned (red), reference (green)'],\n", - " use_same_window = True,\n", - " is_interactive = is_interactive,\n", - " show_borders = show_borders,\n", - " meshes_color = ['tan', 'red', 'green'],\n", - " opacities = [0.5, 0.5, 0.5],\n", - " mesh_style = mesh_style,\n", - " show_mesh_edges = show_mesh_edges,\n", - " show_axes = show_axes,\n", - " show_bounds = show_bounds,\n", - " show_all_edges = show_all_edges,\n", - " font_size = font_size,\n", - " link_views = link_views)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "

" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "#### Aligning all segmentations to the reference shape\n", - "\n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# get a copy for the reference segmentation to avoid inplace antialiasing \n", - "refSeg = shapeSegList[referenceIdx].copy()\n", - "\n", - "# antialias reference segmentation onnce\n", - "refSeg.antialias(antialias_iterations)\n", - "\n", - "# alignment parameters\n", - "isoValue = 1e-20\n", - "icp_iterations = 200\n", - "\n", - "for shapeSeg, shapeName in zip(shapeSegList, shapeNames):\n", - " print('Aligning ' + shapeName + ' to ' + shapeNames[referenceIdx]) \n", - "\n", - " # compute rigid transformation\n", - " rigidTransform = shapeSeg.antialias(antialias_iterations).createTransform(refSeg,\n", - " sw.TransformType.IterativeClosestPoint,\n", - " isoValue,\n", - " icp_iterations\n", - " )\n", - " \n", - " # second we apply the computed transformation, note that shapeSeg has \n", - " # already been antialiased, so we can directly apply the transformation \n", - " shapeSeg.applyTransform(rigidTransform, \n", - " refSeg.origin(), refSeg.dims(), \n", - " refSeg.spacing(), refSeg.coordsys(), \n", - " sw.InterpolationType.Linear)\n", - " \n", - " # then turn antialized-tranformed segmentation to a binary segmentation\n", - " shapeSeg.binarize()" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "#### Visualizing aligned shapes" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# .. and visualize aligned dataset\n", - "plot_volumes(shapeSegList, \n", - " volumeNames = shapeNames,\n", - " is_interactive = is_interactive,\n", - " show_borders = show_borders,\n", - " shade_volumes = shade_volumes,\n", - " show_axes = show_axes,\n", - " show_bounds = show_bounds,\n", - " show_all_edges = show_all_edges,\n", - " font_size = font_size, \n", - " link_views = link_views)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "

" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "#### Comparing mean shapes before & after alignment" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# let's compute the mean shape - note that all images have the same reference\n", - "# dimensions so no need to worry about different size images and no need for padding\n", - "\n", - "dims = refSeg.dims()\n", - "meanShapeAfterAlignment = shapeSegList[0].copy() # start with a copy of the first shape\n", - "\n", - "# add all the other shapes\n", - "for i in range(1, num_samples):\n", - " meanShapeAfterAlignment += shapeSegList[i]\n", - " \n", - "# compute their average\n", - "meanShapeAfterAlignment /= num_samples" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# let's visualize the mean shape (segmentation) before and after alignment\n", - "# notice how the mean shape after alignment reflects a better ellipsoidal shape\n", - "meanShapeAfterAlignment_vtk = sw2vtkImage(meanShapeAfterAlignment)\n", - "\n", - "print('mean shape before alignment')\n", - "display(\n", - " itkw.view( #image = meanShape_vtk, # for orthoginal image plane\n", - " label_image = meanShape_vtk, # for volume rendering segmentation\n", - " slicing_planes = True,\n", - " axes = True,\n", - " rotate = True, # enable auto rotation\n", - " interpolation = True\n", - " )\n", - " )\n", - "\n", - "print('mean shape after alignment')\n", - "display(\n", - " itkw.view( #image = meanShapeAfterAlignment_vtk, # for orthoginal image plane\n", - " label_image = meanShapeAfterAlignment_vtk, # for volume rendering segmentation\n", - " slicing_planes = True,\n", - " axes = True,\n", - " rotate = True, # enable auto rotation\n", - " interpolation = True\n", - " )\n", - " )" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "

" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "#### Saving aligned segmentations" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "save_images(outDir + 'aligned/', \n", - " shapeSegList,\n", - " shapeNames,\n", - " extension = 'nrrd',\n", - " compressed = False, # use false to load in paraview \n", - " verbose = True)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## 3. Cropping and padding segmentations \n", - "\n", - "\n", - "As you can observed from the mean shape (after alignment), image boundaries are not tight around shapes. This leaves too much irrelevant background voxels that might increase the memory footprint when optimizing the shape model. Hence, let's learn how can we remove this irrelevant background while keeping our segmentations intact and avoid cropped segmentations to touch image boundaries.\n", - "\n", - "This step entails:\n", - "- Finding the largest bounding box in which all segmentations are inscribed\n", - "- Cropping the segmentations using the computing bounding box\n", - "- Padding cropped segmentations" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "itkw.view( image = meanShapeAfterAlignment_vtk, # for orthoginal image plane\n", - " #label_image = meanShapeAfterAlignment_vtk, # for volume rendering segmentation\n", - " slicing_planes = True,\n", - " axes = True,\n", - " rotate = True, # enable auto rotation\n", - " interpolation = True\n", - " )" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Finding the largest bounding box\n", - "\n", - "As finding the bounding box is a process that needs to take into account all the segmentations in our dataset, the `ImageUtils` library in `shapeworks` is the right place to look for an API for this process. " - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# finding the bounding box takes in a list of images and the isovalue that \n", - "# defines the level set of shapes implicitly defined in these images \n", - "help(sw.ImageUtils.boundingBox)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# note that the aligned segmentations are still binary images, so an good isovalue \n", - "# that reflect the foreground-background interface would be 0.5\n", - "isoValue = 0.5\n", - "segsBoundingBox = sw.ImageUtils.boundingBox(shapeSegList, isoValue)\n", - "\n", - "print('Computing bounding box:')\n", - "print(segsBoundingBox)\n", - "print(type(segsBoundingBox))" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Cropping segmentations\n", - "\n", - "Cropping a segmentation is a process that can be applied for a single segmentation. So, let's find out the API from the `Image` library." - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "#### Exploring cropping for a single segmentation" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# let's try this out first on a single segmentation\n", - "shapeIdx = 0\n", - "\n", - "# keep a copy to avoid inplace update for now\n", - "shapeSeg = shapeSegList[shapeIdx].copy() \n", - "\n", - "# dot-tap magic\n", - "help(shapeSeg.crop)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# to see the impact of cropping, let's keep a copy\n", - "shapeSegCropped = shapeSeg.copy()\n", - "\n", - "shapeSegCropped.crop(segsBoundingBox)\n", - "\n", - "print('Header information before cropping')\n", - "print(shapeSeg)\n", - "\n", - "print('Header information after cropping')\n", - "print(shapeSegCropped)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# let's visualize them\n", - "plot_volumes([shapeSeg, shapeSegCropped],\n", - " ['before cropping','after cropping'],\n", - " is_interactive = is_interactive,\n", - " show_borders = show_borders,\n", - " shade_volumes = shade_volumes,\n", - " show_axes = show_axes,\n", - " show_bounds = show_bounds,\n", - " show_all_edges = show_all_edges,\n", - " font_size = font_size, \n", - " link_views = link_views)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "

" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "#### Inspecting how tight the cropping is\n", - "\n", - "It is worth noting here that cropping here is a tight cropping without any extra padding. To see this effect, let's find the segmentation that has the largest volume and inspect the effect of cropping on it." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# compute the volume (number of voxels = 1) of each segmentation\n", - "volumes = np.zeros((num_samples,))\n", - "for ii, shapeSeg in enumerate(shapeSegList):\n", - " volumes[ii] = np.sum(shapeSeg.toArray())\n", - " \n", - "plt.figure(dpi=150)\n", - "plt.bar(range(num_samples), volumes)\n", - "plt.xlabel('shape #')\n", - "plt.ylabel('volume')\n", - "\n", - "largestIdx = volumes.argmax()\n", - "plt.bar(largestIdx, volumes[largestIdx]); " - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "

" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# let's inspect the cropping of the segmentation with the largest volume\n", - "# keep a copy to avoid inplace update for now\n", - "shapeSeg = shapeSegList[largestIdx].copy()\n", - "\n", - "# keep a copy to see the impact\n", - "shapeSegCropped = shapeSeg.copy()\n", - "\n", - "# apply cropping\n", - "shapeSegCropped.crop(segsBoundingBox)\n", - "\n", - "shapeSeg_vtk = sw2vtkImage(shapeSeg)\n", - "shapeSegCropped_vtk = sw2vtkImage(shapeSegCropped)\n", - "\n", - "# to see how the segmentation relate to image boundaries before and after cropping\n", - "# let's visualize the volume with orthogonal image slices\n", - "# note how tight the cropping is\n", - "print('before cropping')\n", - "display(\n", - " itkw.view( image = shapeSeg_vtk, # for orthoginal image plane\n", - " label_image = shapeSeg_vtk, # for volume rendering segmentation\n", - " slicing_planes = True,\n", - " axes = True,\n", - " rotate = True, # enable auto rotation\n", - " interpolation = True\n", - " )\n", - " )\n", - "\n", - "print('after cropping')\n", - "display(\n", - " itkw.view( image = shapeSegCropped_vtk, # for orthoginal image plane\n", - " label_image = shapeSegCropped_vtk, # for volume rendering segmentation\n", - " slicing_planes = True,\n", - " axes = True,\n", - " rotate = True, # enable auto rotation\n", - " interpolation = True\n", - " )\n", - " )" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "

" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Hence, to avoid cropped segmentations to touch the image boundaries, we will crop then pad the segmentations. " - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Padding segmentations\n", - "\n", - "Given a segmentation volume, we would like to pad the volume in the three dimensions with a constant value. In this case, we will use zero padding, which will add more background voxels to all dimensions.\n", - "\n", - "Let's start with one segmentation and figure out which `shapeworks` API to use and what parameters are expected." - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "#### Which `shapeworks` API to use?" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# let's see if there's a function that pads a shapeworks image\n", - "# use dot-tap to get a list of functions/apis available for shapeSeg\n", - "\n", - "# found it - pad, let's see its help\n", - "help(shapeSeg.pad)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "#### Exploring padding for a single segmentation" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# parameters for padding \n", - "# The padding size and padding value are the two parameters for this step.\n", - "paddingSize = 10 # number of voxels to pad for each dimension\n", - "paddingValue = 0 # the constant value used to pad the segmentations\n", - "\n", - "# to visualize the effect of padding, let's keep the original segmentation unchanged\n", - "shapeSegPadded = shapeSegCropped.copy()\n", - "\n", - "shapeSegPadded.pad(paddingSize, paddingValue)\n", - "\n", - "print('\\nHeader information before padding: ')\n", - "print(shapeSegCropped)\n", - "\n", - "print('\\nHeader information after padding: ')\n", - "print(shapeSegPadded)\n", - " \n", - "\n", - "plot_volumes([shapeSegCropped, shapeSegPadded],\n", - " ['before padding','after padding'],\n", - " is_interactive = is_interactive,\n", - " show_borders = show_borders,\n", - " shade_volumes = shade_volumes,\n", - " show_axes = show_axes,\n", - " show_bounds = show_bounds,\n", - " show_all_edges = show_all_edges,\n", - " font_size = font_size,\n", - " link_views = link_views)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "

" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Cropping and padding all segmentations" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# let's now crop-pad all segmentations, here we will use inplace update for shapeworks image\n", - "\n", - "# parameters for padding \n", - "paddingSize = 10 # number of voxels to pad for each dimension\n", - "paddingValue = 0 # the constant value used to pad the segmentations\n", - "\n", - "for shapeSeg, shapeName in zip(shapeSegList, shapeNames):\n", - " print('Cropping & padding segmentation: ' + shapeName) \n", - "\n", - " # crop-pad segmentation by chaining apis\n", - " shapeSeg.crop(segsBoundingBox).pad(paddingSize, paddingValue)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Visualizing cropped and padded segmentations" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# .. and visualize cropped dataset\n", - "plot_volumes(shapeSegList,\n", - " volumeNames = shapeNames,\n", - " is_interactive = is_interactive,\n", - " show_borders = show_borders,\n", - " shade_volumes = shade_volumes,\n", - " show_axes = show_axes, \n", - " show_bounds = show_bounds,\n", - " show_all_edges = show_all_edges,\n", - " font_size = font_size,\n", - " link_views = link_views)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "

" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Saving cropped segmentations" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "save_images(outDir + 'cropped/',\n", - " shapeSegList,\n", - " shapeNames,\n", - " extension = 'nrrd',\n", - " compressed = False, # use false to load in paraview \n", - " verbose = True)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## 4. Converting segmentations to smooth signed distance transforms\n", - "\n", - "We are now one-step away from optimizing our shape model. Stay tuned! \n", - "\n", - "For numerical computations for correspondences optimization, we need to convert binary segmentations to a continuous-valued image that satisfies the following requirements.\n", - "- smooth for gradient updates stability\n", - "- reflect the shape's surface (i.e., foreground-background) interface\n", - "- provide a signal for the particle to snap (move back) to the surface in case particles gets off the surface during optimization, which is a typical scenario when using gradient descent based optimization\n", - "\n", - "So far, we have been antialiasing segmentations every time we need to convert binary segmentations to a continuous-valued image (e.g., resampling and alignment). An antialiased segmentation satisfies the first two requirements. However, if a particles leaves the surface (i.e., the zero-level set), it would be challenging to snap it back to the surface.\n", - "\n", - "Another representation that satisfies all the requirements is the *signed distance transform*.\n", - "\n", - "- A signed distance transform assigns to each voxel the physical distance to the closest point on the surface (i.e., the minimum distance from that voxel to nearest voxel on the foreground-background interface). \n", - "- The sign is used to indicate whether that voxel is inside or outside the foreground object. \n", - "- The zero-level set (zero-distance to the surface) indicates the foreground-background interface (i.e., the shape's surface).\n", - "- The gradient of a signed distance transform at a voxels indicats what direction to move in from that voxels to most rapidly increase the value of this distance. Hence, we can use the negative of this gradient as a signal to move a particle back to the surface." - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Exploring distance transforms for a single segmentation" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "shapeIdx = 0\n", - "\n", - "# let's keep a copy to avoid inplace update\n", - "shapeSeg = shapeSegList[shapeIdx].copy()\n", - "\n", - "# computing distance transforms is a single image operation - let's use the dot-tap magic\n", - "help(shapeSeg.computeDT)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "The `computeDT` API needs an isovalue that defines the foreground-background interface. For binary segmentations, this interface will not be smooth (notice the staircase effect) due to the aliasing effect of binarization. Hence, a smoother interface can be defined by first antialiasing the segmentation then compute the distance transform at the zero-level set." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "isoValue = 0\n", - "shapeDT = shapeSeg.copy()\n", - "shapeDT.antialias(antialias_iterations).computeDT(isoValue)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# let's visualize them\n", - "plot_volumes([shapeSeg, shapeDT],\n", - " ['segmentation','signed distance transform'],\n", - " is_interactive = is_interactive,\n", - " show_borders = show_borders, \n", - " shade_volumes = shade_volumes,\n", - " show_axes = show_axes,\n", - " show_bounds = show_bounds,\n", - " show_all_edges = show_all_edges,\n", - " font_size = font_size,\n", - " link_views = link_views)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "

" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# let's compare isosurfaces\n", - "plot_meshes([shapeSeg.toMesh(0.5), shapeDT.toMesh(0)],\n", - " ['seg (tan), dt (red)'],\n", - " use_same_window = True,\n", - " is_interactive = is_interactive,\n", - " show_borders = show_borders,\n", - " meshes_color = ['yellow', 'red'],\n", - " opacities = [0.7, 0.7],\n", - " mesh_style = mesh_style,\n", - " show_mesh_edges = show_mesh_edges,\n", - " show_axes = show_axes,\n", - " show_bounds = show_bounds,\n", - " show_all_edges = show_all_edges,\n", - " font_size = font_size,\n", - " link_views = link_views)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "

" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# let's have a closer look, is it smooth?\n", - "plot_meshes([shapeDT.toMesh(0)],\n", - " ['isosurface of dt'],\n", - " use_same_window = True,\n", - " is_interactive = is_interactive,\n", - " show_borders = show_borders,\n", - " meshes_color = ['tan'],\n", - " opacities = [1],\n", - " mesh_style = mesh_style,\n", - " show_mesh_edges = show_mesh_edges,\n", - " show_axes = show_axes,\n", - " show_bounds = show_bounds,\n", - " show_all_edges = show_all_edges,\n", - " font_size = font_size,\n", - " link_views = link_views)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "

" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Exploring smoothing the distance transform of a single segmentation\n", - "\n", - "We can still see some leftovers from the aliasing effect of binarization. Let's try to smooth this out without significantly impacting the shape features that we want to model.\n", - "\n", - "`shapeworks` library supports two methods of smoothing, gaussian-blur and topology-preserving smoothing. The gaussian blur method filters/convolves the image with a 3D gaussian filter with a given sigma (in physical coordinates). This method could be use for blobby-like structures. However, for shapes with thin features and high curvature regions, the gaussian blurring method could impact the underlying geometrical features. For these shapes, topology-preserving smoothing is recommended. \n", - "\n", - "\n", - "**Note:** topology-preserving smoothing is currently being debugged and tested. Issue #850" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# let's try gaussian blur first\n", - "shapeDTgauss = shapeDT.copy()\n", - "help(shapeDTgauss.gaussianBlur)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# gaussian blur needs the convolution fiter size (i.e., sigma)\n", - "blur_sigma = 2 # physcial coordinates \n", - "shapeDTgauss.gaussianBlur(blur_sigma)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# let's have a look\n", - "plot_meshes([shapeDT.toMesh(0), shapeDTgauss.toMesh(0)],\n", - " ['distance tranform', 'gaussian blur'],\n", - " use_same_window = False,\n", - " is_interactive = is_interactive,\n", - " show_borders = show_borders,\n", - " meshes_color = 'tan',\n", - " opacities = 1,\n", - " mesh_style = mesh_style,\n", - " show_mesh_edges = show_mesh_edges,\n", - " show_axes = show_axes,\n", - " show_bounds = show_bounds,\n", - " show_all_edges = show_all_edges,\n", - " font_size = font_size,\n", - " link_views = link_views)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "

" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# this seems to be a nicely smooth shape, let's compare isosurfaces\n", - "plot_meshes([shapeDT.toMesh(0), shapeDTgauss.toMesh(0)],\n", - " ['dt (tan), gaussian-blur (red)'],\n", - " use_same_window = True,\n", - " is_interactive = is_interactive,\n", - " show_borders = show_borders,\n", - " meshes_color = ['yellow', 'red'],\n", - " opacities = [0.7, 0.7],\n", - " mesh_style = mesh_style,\n", - " show_mesh_edges = show_mesh_edges,\n", - " show_axes = show_axes,\n", - " show_bounds = show_bounds,\n", - " show_all_edges = show_all_edges,\n", - " font_size = font_size,\n", - " link_views = link_views)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "

" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# this smoothing seems to impact more the major radius of the ellispoid\n", - "# a shape feature that we would like to model\n", - "# let's try to reduce the gaussian's sigma\n", - "\n", - "blur_sigma = 1.3 # physcial coordinates \n", - "shapeDTgauss = shapeDT.copy()\n", - "shapeDTgauss.gaussianBlur(blur_sigma)\n", - "\n", - "plot_meshes([shapeDT.toMesh(0), shapeDTgauss.toMesh(0)],\n", - " ['dt (tan), gaussian-blur (red)'],\n", - " use_same_window = True,\n", - " is_interactive = is_interactive,\n", - " show_borders = show_borders,\n", - " meshes_color = ['yellow', 'red'],\n", - " opacities = [0.7, 0.7],\n", - " mesh_style = mesh_style,\n", - " show_mesh_edges = show_mesh_edges,\n", - " show_axes = show_axes,\n", - " show_bounds = show_bounds,\n", - " show_all_edges = show_all_edges,\n", - " font_size = font_size, \n", - " link_views = link_views)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "

" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# not that bad\n", - "# let's have a look on how smooth it is\n", - "plot_meshes([shapeDT.toMesh(0), shapeDTgauss.toMesh(0)],\n", - " ['distance tranform', 'gaussian blur'],\n", - " use_same_window = False,\n", - " is_interactive = is_interactive,\n", - " show_borders = show_borders,\n", - " meshes_color = 'tan',\n", - " opacities = 1,\n", - " mesh_style = mesh_style,\n", - " show_mesh_edges = show_mesh_edges,\n", - " show_axes = show_axes,\n", - " show_bounds = show_bounds,\n", - " show_all_edges = show_all_edges,\n", - " font_size = font_size,\n", - " link_views = link_views)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "

" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Computing smooth distance transforms for all segmentations" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "isoValue = 0\n", - "sigma = 1.3\n", - "for shapeSeg, shapeName in zip(shapeSegList, shapeNames):\n", - " print('Compute DT for segmentation: ' + shapeName) \n", - "\n", - " # antialias-dt computing\n", - " shapeSeg.antialias(antialias_iterations).computeDT(isoValue).gaussianBlur(sigma)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Saving distance transforms" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "save_images(outDir + 'dts/',\n", - " shapeSegList, \n", - " shapeNames,\n", - " extension = 'nrrd',\n", - " compressed = False, # use false to load in paraview \n", - " verbose = True)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## 5. Optimizing correspondence: just scratching the surface" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# compose the list of groomed segmentations\n", - "filenames = ''\n", - "for curName in shapeNames:\n", - " filename = Path(outDir + 'dts/' + curName + '.nrrd')\n", - " filenames += ' ' + str(filename)\n", - "print(filenames)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# let's launch studio to optimize our shape model \n", - "# see the video in the next cell to an illustration\n", - "\n", - "!ShapeWorksStudio {str(filenames)}" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "

" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Congrats! You have been able to groom your dataset ..." - ] - } - ], - "metadata": { - "kernelspec": { - "display_name": "Python 3", - "language": "python", - "name": "python3" - }, - "language_info": { - "codemirror_mode": { - "name": "ipython", - "version": 3 - }, - "file_extension": ".py", - "mimetype": "text/x-python", - "name": "python", - "nbconvert_exporter": "python", - "pygments_lexer": "ipython3", - "version": "3.7.8" - }, - "toc": { - "base_numbering": 1, - "nav_menu": {}, - "number_sections": false, - "sideBar": true, - "skip_h1_title": true, - "title_cell": "Table of Contents", - "title_sidebar": "Contents", - "toc_cell": false, - "toc_position": { - "height": "calc(100% - 180px)", - "left": "10px", - "top": "150px", - "width": "358.390625px" - }, - "toc_section_display": true, - "toc_window_display": false - } - }, - "nbformat": 4, - "nbformat_minor": 4 -}