Using ITK in Tomviz

Tomviz is an application for visualizing, processing, and analyzing large datasets from transmission and scanning transmission electron microscopes. It is based on the ParaView platform and includes the Insight Segmentation and Registration Toolkit (ITK). ParaView provides advanced visualization while ITK provides a wealth of image processing and analysis algorithms. Together they provide powerful capabilities for analyzing and understanding images in the materials science domain.

In this post, I will walk through how to use ITK within Tomviz and end by describing one of Tomviz’s ITK-based data operators to provide a particular image processing capability.

ITK

ITK is a C++ library with deep roots in medical image processing and analysis, but its use is far from limited to medical imaging applications. ITK provides a host of general purpose image processing and analysis algorithms for tasks such as image smoothing, morphology operations, image resampling, basic math operations, automated thresholding, segmentation, optimization routines, level sets, edge detection, Fourier analysis, and I/O for a large number of file formats. These image-oriented analysis features complement the visualization and analysis capabilities in the ParaView platform.
ITK’s Python language bindings make it quite accessible from a Python scripting environment like the one in Tomviz. Tomviz offers Python scripting as a first-class way to define data processing operations and includes Python bindings for ITK and VTK. The bindings for each library can be used in a Python script, enabling a mix of VTK and ITK functions to be used in a custom data processing algorithm. NumPy and SciPy are included as well, along with bridges from NumPy to ITK and VTK.

ITK Basics

Let’s take a quick look at some ITK fundamentals. The template class itk::Image is the basic image data unit in ITK. It is templated by image voxel type and image dimension. As an example, itk::Image<unsigned int, 3> specifies an image with unsigned int voxels in three dimensions. An itk::Image also contains metadata such as the number of voxels in each dimension, the starting voxel index, the position of the origin of the image, the voxel spacing, and the orientation of the image.

ITK typically operates on itk::Image instances via filters. A filter takes one or more itk::Image instances  as input and produces one or more itk::Image instances as output. Filters can be chained together in succession to perform complex processing and analysis algorithms. Most ITK filters are templated by input and output image types, e.g., itk::GradientAnisotropicDiffusionImageFilter<InputImageType, OutputImageType> where InputImageType and OutputImageType are convenient typedefs for the fully expanded template instantiation of the input and output image types.

There is a nice mapping between instantiation of an ITK filter in C++ and instantiation in Python. In C++, one writes:

#include <itkImage.h>
#include <itkGradientAnisotropicDiffusionImageFilter.h>

typedef itk::Image<unsigned char, 3> InputImageType;
typedef itk::Image<float, 3> OutputImageType;
typedef itk::GradientAnisotropicDiffusionImageFilter<
  InputImageType, OutputImageType> DiffusionFilterType;

DiffusionFilterType* filter = DiffusionFilterType::New();

This C++ code defines three types: an itk::Image instantiation for an input 3D image with unsigned char voxels, an output 3D image with float voxels, and an ITK filter that performs anistropic diffusion. The last line shows the object instantiation of the filter.

In Python, creating the equivalent filter can be done with

import itk

input_type = itk.Image.UC3
output_type = itk.Image.F3
DiffusionFilterType = itk.GradientAnisotropicDiffusionImageFilter[
  input_type, output_type]
diffusion_filter = DiffusionFilterType.New()

This Python code does the equivalent of the C++ code. It imports the ITK Python module (import itk), defines the input and output image types, the filter type, and instantiates a filter object. Note the convention of the itk.Image type attributes. The first part stands for the voxel type while the number indicates the dimensionality. In this example, the UC stands for unsigned char, the F for float, and the 3 for 3D images.

Python Data Operators in Tomviz

Tomviz provides a variety of data processing operations called data operators. A data operator modifies the selected data set and replaces the input data set with the result of the operator. Data operators can be chained together in Tomviz to form complex image analysis pipelines.

One style of data operator uses a Python script to define its operation. A Python data operator is defined as a single function in a Python script:

def transform_scalars(dataset, conductance=1.0, iterations=100, timestep=0.125):
    # transform dataset here

Here, dataset is the vtkImageData object produced by the prior data operators in the pipeline. The function parameters conductance, iterations, and timestep are defined in a JSON-formatted description of this data operator, as described in this blog post. The body of this function is where the data transformation, including any ITK-based processing, is done.

Converting the Data Set to an itk::Image

To use Tomviz image data in an ITK filter, one must first convert it to an itk::Image. Tomviz has a convenience function in its tomviz.itkutils module for converting vtkImageData to an itk::Image type with a suitable voxel type and image dimension.

from tomviz import itkutils
itk_image = \
  itkutils.convert_vtk_to_itk_image(dataset)

Behind the scenes, a map from VTK image type to ITK image type determines which ITK image type to instantiate. Doing so turns out to be somewhat involved.

ITK’s use of templates presents some challenges when wrapping it in Python. In C++, ITK is nearly a header-only library due to its heavy use of templates. Templated class definitions are created as needed, and only the code needed to instantiate the desired templated types is compiled into binary form. In Python, however, the story is different. instantiating templated classes at the point they are compiled into bytecode is not possible. As a result, a selection of ITK image and image filter types to instantiate into the wrapped libraries must be defined at the time the wrappings are compiled. The Tomviz build process does just this –  various combinations of filter input and output types are also defined. Unfortunately, the combinatorial explosion that occurs when many voxel types are requested can make the size of the binary library files for ITK’s Python wrapping rather large. To mitigate this, a limited selection of voxel types are used in the wrapping, mainly unsigned integral types and single-precision floats.

As a consequence of the limited number of voxel types supported in the ITK Python bindings, the Python script may encounter a vtkImageData voxel type that does not have a direct corresponding ITK voxel type available. The solution in Tomviz is to cast the VTK image data to a supported voxel type that can reasonably represent the VTK image data if possible and then convert the image data to an itk::Image. Usually this involves a type promotion to the next-highest ranking type. For example, if unsigned char is not available but unsigned short is available, then the image data is cast to an unsigned short.

In some cases, it is necessary to specify the desired ITK image data type directly. For example, the itk::GradientAnisotropicDiffusionImageFilter is instantiated only over floats in the ITK Python bindings, so its input must be a float image. In such cases, the target image voxel type may be given explicitly as an optional second argument:

import itkTypes
itk_image = \
  itkutils.convert_vtk_to_itk_image(dataset,
                                    itkTypes.F)

Setting up ITK filters

Much like VTK, ITK filters can be chained together into a pipeline. The ITK image data converted from the Tomviz image data is typically passed to the first filter in the chain and the output of the last filter is typically converted back to a form that can be copied to the dataset object originally passed into the script. In the simplest case, a pipeline consists of a single filter. Continuing with our example from above, we have a pipeline with a single itk.GradientAnisotropicDiffusionImageFilter instance.

diffusion_filter.SetInput(itk_image)

The itk.GradientAnisotropicDiffusionImageFilter has several properties we wish to set: conductance, number of iterations, and timestep value. We simply forward the parameters passed into the transform_scalars function to the ITK filter:

diffusion_filter.SetConductanceParameter(conductance)
diffusion_filter.SetNumberOfIterations(iterations)
diffusion_filter.SetTimeStep(timestep)

Once all the properties of a filter are set up, it can be executed with

diffusion_filter.Update()

This call executes the filter’s algorithm and updates the output image storing the filter results. Note that in a multi-filter pipeline, only the last filter’s Update() method needs to be called to execute the entire ITK processing pipeline.

The pipeline results are stored in an ITK image of type output_type. To get these results into a format Tomviz can use, one can invoke

itkutils.set_array_from_itk_image(dataset, 
                                  diffusion_filter.GetOutput())

This function takes the output from the filter and copies it to the VTK image data structure, completing the script.

The Complete Data Operator

Putting it all together, here is the final Python script for the Perona-Malik Anisotropic Diffusion data operator:

def transform_scalars(dataset, conductance=1.0, iterations=100, timestep=0.125):
    """This filter performs anisotropic diffusion
    on an image using the classic Perona-Malik,
    gradient magnitude-based equation.
    """

    # Make sure we can import needed modules
    try:
        import itk
        import itkTypes
        from tomviz import itkutils
    except Exception as exc:
        print("Could not import necessary module(s)")
        print(exc)

    try:
        # Get the ITK image. The anisotropic diffusion
        # filter is templated over float pixel types
        # only, so explicitly request a float ITK
        # image type.
        itk_image = \
          itkutils.convert_vtk_to_itk_image(
            dataset, itkTypes.F)
        itk_image_type = type(itk_image)

        DiffusionFilterType = \
           itk.GradientAnisotropicDiffusionImageFilter[
             itk_image_type, itk_image_type]
        diffusion_filter = DiffusionFilterType.New()
        diffusion_filter.SetConductanceParameter(conductance)
        diffusion_filter.SetNumberOfIterations(iterations)
        diffusion_filter.SetTimeStep(timestep)
        diffusion_filter.SetInput(itk_image)
        diffusion_filter.Update()
        itkutils.set_array_from_itk_image(
            dataset, diffusion_filter.GetOutput())
    except Exception as exc:
        print("Exception encountered while running"
              "PeronaMalikAnisotropicDiffusion")
        print(exc)
        raise(exc)

Looking Forward

This post covers one example of exposing one ITK filter’s functionality in Tomviz. Half a dozen other ITK-based data operators are now available in Tomviz as of version 0.9.4. In the future, the Tomviz development team will work on expanding the types of ITK filters that can be used in Tomviz. As one example, we will look into adding intuitive 3D placement of seed points for image analysis algorithms that need them.

The Tomviz project is developed as part of a collaboration between Kitware and Cornell University under DOE Office of Science contract DE-SC0011385. This is a community project, and we are very pleased to take input and contributions from all in the community.