Basic Registration

Summary:

  1. Creating an instance of the registration framework requires selection of the following components:
    • Optimizer.
    • Similarity metric.
    • Interpolator.
  2. The registration framework only supports images with sitkFloat32 and sitkFloat64 pixel types (use the SimpleITK Cast() function if your image's pixel type is something else).

  3. Successful registration is highly dependent on initialization. In general you can:

    • Use auxiliary information or user interaction to obtain an initial transformation (avoid resampling).
    • Center the images using the CenteredTransformInitializer.
    • Coarsely sample the parameter space using the Exhaustive Optimizer to obtain one or more initial transformation estimates.
    • Manually initialize, via direct manipulation of transformation parameters and visualization or localization of corresponding points in the two images and then use the LandmarkBasedTransformInitializer.

Registration Components



There are many options for creating an instance of the registration framework, all of which are configured in SimpleITK via methods of the ImageRegistrationMethod class. This class encapsulates many of the components available in ITK for constructing a registration instance.

Currently, the available choices from the following groups of ITK components are:

Optimizers

The SimpleITK registration framework supports several optimizer types via the SetOptimizerAsX() methods, these include:

Similarity metrics

The SimpleITK registration framework supports several metric types via the SetMetricAsX() methods, these include:

Interpolators

The SimpleITK registration framework supports several interpolators via the SetInterpolator() method, which receives one of the following enumerations:

  • sitkNearestNeighbor
  • sitkLinear
  • sitkBSpline
  • sitkGaussian
  • sitkHammingWindowedSinc
  • sitkCosineWindowedSinc
  • sitkWelchWindowedSinc
  • sitkLanczosWindowedSinc
  • sitkBlackmanWindowedSinc
In [1]:
import SimpleITK as sitk
from downloaddata import fetch_data as fdata
%matplotlib nbagg
import gui
import registration_gui as rgui

import numpy as np
import os
OUTPUT_DIR = 'output'

Read images

We first read the images, specifying the pixel type that is required for registration (Float32 or Float64) and look at them. In this notebook we use a CT and MR image from the same patient. These are part of the training data from the Retrospective Image Registration Evaluation (RIRE) project.

In [2]:
fixed_image =  sitk.ReadImage(fdata("training_001_ct.mha"), sitk.sitkFloat32)
moving_image = sitk.ReadImage(fdata("training_001_mr_T1.mha"), sitk.sitkFloat32)

ct_window_level = [835,162]
mr_window_level = [1036,520]

gui.MultiImageDisplay(image_list = [fixed_image, moving_image],                   
                      title_list = ['fixed', 'moving'], figure_size=(8,4), window_level_list=[ct_window_level, mr_window_level]);
Fetching training_001_ct.mha
Fetching training_001_mr_T1.mha

Classic Registration

Estimate a 3D rigid transformation between images of different modalities.

We have made the following choices with respect to initialization and registration component settings:

  • Similarity metric, mutual information (Mattes MI):
    • Number of histogram bins, 50.
    • Sampling strategy, random.
    • Sampling percentage, 1%.
  • Interpolator, sitkLinear.
  • Optimizer, gradient descent:
    • Learning rate, step size along traversal direction in parameter space, 1.0 .
    • Number of iterations, maximal number of iterations, 100.
    • Convergence minimum value, value used for convergence checking in conjunction with the energy profile of the similarity metric that is estimated in the given window size, 1e-6.
    • Convergence window size, number of values of the similarity metric which are used to estimate the energy profile of the similarity metric, 10.

We initialize registration by aligning the centers of the two volumes. To qualitatively evaluate the result we use a linked cursor approach, click on one image and the corresponding point is added to the other image.

In [3]:
initial_transform = sitk.CenteredTransformInitializer(fixed_image, 
                                                      moving_image, 
                                                      sitk.Euler3DTransform(), 
                                                      sitk.CenteredTransformInitializerFilter.GEOMETRY)

gui.RegistrationPointDataAquisition(fixed_image, moving_image, figure_size=(8,4), known_transformation=initial_transform, fixed_window_level=ct_window_level, moving_window_level=mr_window_level);

Run the next cell three times:

  1. As is.
  2. Uncomment the automated optimizer scale setting so that a change in rotation (radians) has a similar effect to a change in translation (mm).
  3. Uncomment the multi-resolution settings.
In [4]:
registration_method = sitk.ImageRegistrationMethod()

# Similarity metric settings.
registration_method.SetMetricAsMattesMutualInformation(numberOfHistogramBins=50)
registration_method.SetMetricSamplingStrategy(registration_method.RANDOM)
registration_method.SetMetricSamplingPercentage(0.01)

registration_method.SetInterpolator(sitk.sitkLinear)

# Optimizer settings.
registration_method.SetOptimizerAsGradientDescent(learningRate=1.0, numberOfIterations=100, convergenceMinimumValue=1e-6, convergenceWindowSize=10)
#registration_method.SetOptimizerScalesFromPhysicalShift()

# Setup for the multi-resolution framework.            
#registration_method.SetShrinkFactorsPerLevel(shrinkFactors = [4,2,1])
#registration_method.SetSmoothingSigmasPerLevel(smoothingSigmas=[2,1,0])
#registration_method.SmoothingSigmasAreSpecifiedInPhysicalUnitsOn()

# Don't optimize in-place, we would possibly like to run this cell multiple times.
registration_method.SetInitialTransform(initial_transform, inPlace=False)

# Connect all of the observers so that we can perform plotting during registration.
registration_method.AddCommand(sitk.sitkStartEvent, rgui.start_plot)
registration_method.AddCommand(sitk.sitkEndEvent, rgui.end_plot)
registration_method.AddCommand(sitk.sitkMultiResolutionIterationEvent, rgui.update_multires_iterations) 
registration_method.AddCommand(sitk.sitkIterationEvent, lambda: rgui.plot_values(registration_method))

final_transform = registration_method.Execute(fixed_image, moving_image)

# Always check the reason optimization terminated.
print('Final metric value: {0}'.format(registration_method.GetMetricValue()))
print('Optimizer\'s stopping condition, {0}'.format(registration_method.GetOptimizerStopConditionDescription()))
Final metric value: -0.35954723072306366
Optimizer's stopping condition, GradientDescentOptimizerv4Template: Convergence checker passed at iteration 15.

Qualitatively evaluate the result using a linked cursor approach (visual evaluation):

In [5]:
gui.RegistrationPointDataAquisition(fixed_image, moving_image, figure_size=(8,4), known_transformation=final_transform,fixed_window_level=ct_window_level, moving_window_level=mr_window_level);

If we are satisfied with the results, save them to file.

In [6]:
moving_resampled = sitk.Resample(moving_image, fixed_image, final_transform, sitk.sitkLinear, 0.0, moving_image.GetPixelID())
sitk.WriteImage(moving_resampled, os.path.join(OUTPUT_DIR, 'RIRE_training_001_mr_T1_resampled.mha'))
sitk.WriteTransform(final_transform, os.path.join(OUTPUT_DIR, 'RIRE_training_001_CT_2_mr_T1.tfm'))

ITKv4 Coordinate Systems

Unlike the classical registration approach where the fixed and moving images are treated differently, the ITKv4 registration framework allows you to treat both images in the same manner. This is achieved by introducing a third coordinate system, the virtual image domain.



Thus, the ITK v4 registration framework deals with three transformations:

  • SetInitialTransform, $T_{opt}$ - composed with the moving initial transform, maps points from the virtual image domain to the moving image domain, modified during optimization.
  • SetFixedInitialTransform $T_f$- maps points from the virtual image domain to the fixed image domain, never modified.
  • SetMovingInitialTransform $T_m$- maps points from the virtual image domain to the moving image domain, never modified.

The transformation that maps points from the fixed to moving image domains is thus: $^M\mathbf{p} = T_{opt}(T_m(T_f^{-1}(^F\mathbf{p})))$

We now modify the previous example to use $T_{opt}$ and $T_m$.

In [7]:
registration_method = sitk.ImageRegistrationMethod()

# Similarity metric settings.
registration_method.SetMetricAsMattesMutualInformation(numberOfHistogramBins=50)
registration_method.SetMetricSamplingStrategy(registration_method.RANDOM)
registration_method.SetMetricSamplingPercentage(0.01)

registration_method.SetInterpolator(sitk.sitkLinear)

# Optimizer settings.
registration_method.SetOptimizerAsGradientDescent(learningRate=1.0, numberOfIterations=100, convergenceMinimumValue=1e-6, convergenceWindowSize=10)
registration_method.SetOptimizerScalesFromPhysicalShift()

# Setup for the multi-resolution framework.            
registration_method.SetShrinkFactorsPerLevel(shrinkFactors = [4,2,1])
registration_method.SetSmoothingSigmasPerLevel(smoothingSigmas=[2,1,0])
registration_method.SmoothingSigmasAreSpecifiedInPhysicalUnitsOn()

# Set the initial moving and optimized transforms.
optimized_transform = sitk.Euler3DTransform()    
registration_method.SetMovingInitialTransform(initial_transform)
registration_method.SetInitialTransform(optimized_transform, inPlace=False)

# Connect all of the observers so that we can perform plotting during registration.
registration_method.AddCommand(sitk.sitkStartEvent, rgui.start_plot)
registration_method.AddCommand(sitk.sitkEndEvent, rgui.end_plot)
registration_method.AddCommand(sitk.sitkMultiResolutionIterationEvent, rgui.update_multires_iterations) 
registration_method.AddCommand(sitk.sitkIterationEvent, lambda: rgui.plot_values(registration_method))

# Need to compose the transformations after registration.
final_transform_v4 = registration_method.Execute(fixed_image, moving_image)
#final_transform_v4.AddTransform(initial_transform)

# Always check the reason optimization terminated.
print('Final metric value: {0}'.format(registration_method.GetMetricValue()))
print('Optimizer\'s stopping condition, {0}'.format(registration_method.GetOptimizerStopConditionDescription()))
Final metric value: -0.6187114209525496
Optimizer's stopping condition, GradientDescentOptimizerv4Template: Convergence checker passed at iteration 11.
In [8]:
gui.RegistrationPointDataAquisition(fixed_image, moving_image, figure_size=(8,4), known_transformation=final_transform_v4, fixed_window_level=ct_window_level, moving_window_level=mr_window_level);