Introduction to SimpleITKv4 Registration No description has been provided for this image

First, a reminder of some SimpleITK conventions:

  • The dimensionality of images being registered is required to be the same (i.e. 2D/2D or 3D/3D).
  • The pixel type of images being registered is required to be the same.
  • Supported pixel types for images being registered are sitkFloat32 and sitkFloat64. You can use the SimpleITK Cast function if your image's pixel type is something else.

Registration Components¶

No description has been provided for this image

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.

Here is a list of currently available choices for each group of components:

Optimizers¶

Several optimizer types are supported via the SetOptimizerAsX() methods. These include:

  • Exhaustive
  • Nelder-Mead downhill simplex, a.k.a. Amoeba
  • Powell optimizer
  • 1+1 evolutionary optimizer
  • Variations on gradient descent:
    • GradientDescent
    • GradientDescentLineSearch
    • RegularStepGradientDescent
  • ConjugateGradientLineSearch
  • L-BFGS2 (Limited memory Broyden, Fletcher, Goldfarb, Shannon)
  • L-BFGS-B (Limited memory Broyden, Fletcher, Goldfarb, Shannon - Bound Constrained) - supports the use of simple constraints.

Similarity metrics¶

Several metric types are supported via the SetMetricAsX() methods. These include:

  • MeanSquares
  • Demons
  • Correlation
  • ANTSNeighborhoodCorrelation
  • JointHistogramMutualInformation
  • MattesMutualInformation

Interpolators¶

Several interpolators are supported via the SetInterpolator() method, which can take one of the following interpolators:

  • sitkNearestNeighbor
  • sitkLinear
  • sitkBSpline1, sitkBSpline2, sitkBSpline3, sitkBSpline4, sitkBSpline5, where the number denotes the spline order, with sitkBSpline3 being the go-to option.
  • sitkGaussian
  • sitkHammingWindowedSinc
  • sitkCosineWindowedSinc
  • sitkWelchWindowedSinc
  • sitkLanczosWindowedSinc
  • sitkBlackmanWindowedSinc

See the list of interpolators in the API reference for more information

Data - Retrospective Image Registration Evaluation¶

In this notebook we will be using part of the training data from the Retrospective Image Registration Evaluation (RIRE) project.

In [1]:
import SimpleITK as sitk

# Utility method that either downloads data from the Girder repository or
# if already downloaded returns the file name for reading from disk (cached data).
%run update_path_to_download_script
from downloaddata import fetch_data as fdata

# Always write output to a separate directory, we don't want to pollute the source directory.
import os

OUTPUT_DIR = "Output"

Utility functions¶

First we define a number of utility callback functions for image display and for plotting the similarity metric during registration.

In [2]:
%matplotlib inline
import matplotlib.pyplot as plt

from ipywidgets import interact, fixed
from IPython.display import clear_output


# Callback invoked by the IPython interact method for scrolling through image stacks of
# the two images being registered.
def display_images(fixed_image_z, moving_image_z, fixed_npa, moving_npa):
    # Create a figure with two subplots and the specified size.
    plt.subplots(1, 2, figsize=(10, 8))

    # Draw the fixed image in the first subplot.
    plt.subplot(1, 2, 1)
    plt.imshow(fixed_npa[fixed_image_z, :, :], cmap=plt.cm.Greys_r)
    plt.title("fixed image")
    plt.axis("off")

    # Draw the moving image in the second subplot.
    plt.subplot(1, 2, 2)
    plt.imshow(moving_npa[moving_image_z, :, :], cmap=plt.cm.Greys_r)
    plt.title("moving image")
    plt.axis("off")

    plt.show()


# Callback invoked by the IPython interact method for scrolling and modifying the alpha blending
# of an image stack of two images that occupy the same physical space.
def display_images_with_alpha(image_z, alpha, fixed, moving):
    img = (1.0 - alpha) * fixed[:, :, image_z] + alpha * moving[:, :, image_z]
    plt.imshow(sitk.GetArrayViewFromImage(img), cmap=plt.cm.Greys_r)
    plt.axis("off")
    plt.show()


# Callback invoked when the StartEvent happens, sets up our new data.
def start_plot():
    global metric_values, multires_iterations

    metric_values = []
    multires_iterations = []


# Callback invoked when the EndEvent happens, do cleanup of data and figure.
def end_plot():
    global metric_values, multires_iterations

    del metric_values
    del multires_iterations
    # Close figure, we don't want to get a duplicate of the plot latter on.
    plt.close()


# Callback invoked when the IterationEvent happens, update our data and display new figure.
def plot_values(registration_method):
    global metric_values, multires_iterations

    metric_values.append(registration_method.GetMetricValue())
    # Clear the output area (wait=True, to reduce flickering), and plot current data
    clear_output(wait=True)
    # Plot the similarity metric values
    plt.plot(metric_values, "r")
    plt.plot(
        multires_iterations,
        [metric_values[index] for index in multires_iterations],
        "b*",
    )
    plt.xlabel("Iteration Number", fontsize=12)
    plt.ylabel("Metric Value", fontsize=12)
    plt.show()


# Callback invoked when the sitkMultiResolutionIterationEvent happens, update the index into the
# metric_values list.
def update_multires_iterations():
    global metric_values, multires_iterations
    multires_iterations.append(len(metric_values))

Read images¶

Now we can read the images, casting the pixel type to that required for registration (Float32 or Float64) and have a look at them.

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

interact(
    display_images,
    fixed_image_z=(0, fixed_image.GetSize()[2] - 1),
    moving_image_z=(0, moving_image.GetSize()[2] - 1),
    fixed_npa=fixed(sitk.GetArrayViewFromImage(fixed_image)),
    moving_npa=fixed(sitk.GetArrayViewFromImage(moving_image)),
);

Fetching training_001_ct.mha
Fetching training_001_mr_T1.mha

Initial Alignment¶

To set sensible default values for the center of rotation and initial translation between the two images, the CenteredTransformInitializer is used to align the centers of the two volumes and set the center of rotation to the center of the fixed image.

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

moving_resampled = sitk.Resample(
    moving_image,
    fixed_image,
    initial_transform,
    sitk.sitkLinear,
    0.0,
    moving_image.GetPixelID(),
)

interact(
    display_images_with_alpha,
    image_z=(0, fixed_image.GetSize()[2] - 1),
    alpha=(0.0, 1.0, 0.05),
    fixed=fixed(fixed_image),
    moving=fixed(moving_resampled),
);

Registration¶

The specific registration task at hand estimates a 3D rigid transformation between images of different modalities. There are multiple components from each group (optimizers, similarity metrics, interpolators) that are appropriate for the task. Note that selecting a component from each group requires setting some parameter values. We have made the following choices:

  • 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, maximum 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

In the next code block we perform registration using the settings given above, and take advantage of the built in multi-resolution framework, using a three tier pyramid.

In this example we plot the similarity metric value during registration. Note that the change of scales in the multi-resolution framework is readily visible.

In [5]:
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, start_plot)
registration_method.AddCommand(sitk.sitkEndEvent, end_plot)
registration_method.AddCommand(
    sitk.sitkMultiResolutionIterationEvent, update_multires_iterations
)
registration_method.AddCommand(
    sitk.sitkIterationEvent, lambda: plot_values(registration_method)
)

final_transform = registration_method.Execute(
    sitk.Cast(fixed_image, sitk.sitkFloat32), sitk.Cast(moving_image, sitk.sitkFloat32)
)
No description has been provided for this image

Post registration analysis¶

Having finished registration, we next query the registration method to see the metric value and the reason the optimization terminated.

The metric value allows us to compare multiple registration runs as there is a probabilistic aspect to our registration. This is due to the random sampling used to estimate the similarity metric.

Always remember to query why the optimizer terminated. This will help you understand whether termination is too early, either due to thresholds being too tight, early termination due to small number of iterations - numberOfIterations, or too loose, early termination due to large value for minimal change in similarity measure - convergenceMinimumValue.

In [6]:
print(f"Final metric value: {registration_method.GetMetricValue()}")
print(
    f"Optimizer's stopping condition, {registration_method.GetOptimizerStopConditionDescription()}"
)

Final metric value: -0.648576284144644
Optimizer's stopping condition, GradientDescentOptimizerv4Template: Convergence checker passed at iteration 73.

Now visually inspect the results.

In [7]:
moving_resampled = sitk.Resample(
    moving_image,
    fixed_image,
    final_transform,
    sitk.sitkLinear,
    0.0,
    moving_image.GetPixelID(),
)

interact(
    display_images_with_alpha,
    image_z=(0, fixed_image.GetSize()[2] - 1),
    alpha=(0.0, 1.0, 0.05),
    fixed=fixed(fixed_image),
    moving=fixed(moving_resampled),
);

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

In [8]:
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")
)

In the cell above we separately saved (1) a resampled version of the moving image and (2) the final registration transformation. Another option is to forgo the resampling and combine the transformation into the original moving image, "physically" moving it in space. This operation modifies the image's origin and direction cosine matrix but leaves the intensity information as is, image size and spacing don't change. In Slicer nomenclature this is referred to as "transform hardening".

In [9]:
transformed_moving = sitk.TransformGeometry(moving_image, final_transform)
print(
    f"origin before: {moving_image.GetOrigin()}\norigin after: {transformed_moving.GetOrigin()}"
)
print(
    f"direction cosine before: {moving_image.GetDirection()}\ndirection cosine after: {transformed_moving.GetDirection()}"
)
sitk.WriteImage(
    transformed_moving,
    os.path.join(OUTPUT_DIR, "RIRE_training_001_mr_T1_transformed.mha"),
)

origin before: (0.0, 0.0, 0.0)
origin after: (-8.053812235231259, 19.905474192760018, 23.27020054966579)
direction cosine before: (1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0)
direction cosine after: (0.997386575264692, -0.0722429579617326, 0.0009871710633186331, 0.07222708740007226, 0.9973231436850912, 0.011392757163613307, -0.0018075750250630346, -0.011291682559569386, 0.9999346131510315)