Registration Initialization: We Have to Start Somewhere

Initialization is a critical aspect of most registration algorithms, given that most algorithms are formulated as an iterative optimization problem. It effects both the runtime and convergence to the correct minimum. Ideally our transformation is initialized close to the correct solution ensuring convergence in a timely manner. Problem specific initialization will often yield better results than generic initialization approaches.

Rule of thumb: use as much prior information (external to the image content) as you can to initialize your registration.

Common initializations strategies when no prior information is available:

  1. Do nothing (hope springs eternal) - initialize using the identity transformation.
  2. CenteredTransformInitializer (GEOMETRY or MOMENTS) - translation based initialization, align the geometric centers of the images or the intensity based centers of mass of the image contents.
  3. Use a sampling of the parameter space (useful mostly for low dimensional parameter spaces).
  4. Manual initialization - allow an operator to control transformation parameter settings directly using a GUI with visual feedback or identify multiple corresponding points in the two images and compute an initial rigid or affine transformation.

In many cases we perform initialization in an automatic manner by making assumptions with regard to the contents of the image and the imaging protocol. For instance, if we expect that images were acquired with the patient in a known orientation we can align the geometric centers of the two volumes or the center of mass of the image contents if the anatomy is not centered in the image (this is what we previously did in this example).

When the orientation is not known, or is known but incorrect, this approach will not yield a reasonable initial estimate for the registration.

When working with clinical images, the DICOM tags define the orientation and position of the anatomy in the volume. The tags of interest are:

  • (0020|0032) Image Position (Patient) : coordinates of the the first transmitted voxel.
  • (0020|0037) Image Orientation (Patient): directions of first row and column in 3D space.
  • (0018|5100) Patient Position: Patient placement on the table
    • Head First Prone (HFP)
    • Head First Supine (HFS)
    • Head First Decubitus Right (HFDR)
    • Head First Decubitus Left (HFDL)
    • Feet First Prone (FFP)
    • Feet First Supine (FFS)
    • Feet First Decubitus Right (FFDR)
    • Feet First Decubitus Left (FFDL)

The patient position is manually entered by the CT/MR operator and thus can be erroneous (HFP instead of FFP will result in a $180^o$ orientation error). In this notebook we use data acquired using an abdominal phantom which made it hard to identify the "head" and "feet" side, resulting in an incorrect value entered by the technician.

In [1]:
import SimpleITK as sitk

# If the environment variable SIMPLE_ITK_MEMORY_CONSTRAINED_ENVIRONMENT is set, this will override the ReadImage
# function so that it also resamples the image to a smaller size (testing environment is memory constrained).
%run setup_for_testing

import os
import numpy as np

from ipywidgets import interact, fixed
%run update_path_to_download_script
from downloaddata import fetch_data as fdata

%matplotlib notebook
import gui

# This is the registration configuration which we use in all cases. The only parameter that we vary 
# is the initial_transform. 
def multires_registration(fixed_image, moving_image, initial_transform):
    registration_method = sitk.ImageRegistrationMethod()
    registration_method.SetMetricAsMattesMutualInformation(numberOfHistogramBins=50)
    registration_method.SetMetricSamplingStrategy(registration_method.RANDOM)
    registration_method.SetMetricSamplingPercentage(0.01)
    registration_method.SetInterpolator(sitk.sitkLinear)
    registration_method.SetOptimizerAsGradientDescent(learningRate=1.0, numberOfIterations=100, estimateLearningRate=registration_method.Once)
    registration_method.SetOptimizerScalesFromPhysicalShift() 
    registration_method.SetInitialTransform(initial_transform, inPlace=False)
    registration_method.SetShrinkFactorsPerLevel(shrinkFactors = [4,2,1])
    registration_method.SetSmoothingSigmasPerLevel(smoothingSigmas = [2,1,0])
    registration_method.SmoothingSigmasAreSpecifiedInPhysicalUnitsOn()

    final_transform = registration_method.Execute(fixed_image, moving_image)
    print('Final metric value: {0}'.format(registration_method.GetMetricValue()))
    print('Optimizer\'s stopping condition, {0}'.format(registration_method.GetOptimizerStopConditionDescription()))
    return (final_transform, registration_method.GetMetricValue())

Loading Data

Note: While the images are of the same phantom, they were acquired at different times and the fiducial markers visible on the phantom are not in the same locations.

Scroll through the data to gain an understanding of the spatial relationship along the viewing (z) axis.

In [2]:
data_directory = os.path.dirname(fdata("CIRS057A_MR_CT_DICOM/readme.txt"))

fixed_series_ID = "1.2.840.113619.2.290.3.3233817346.783.1399004564.515"
moving_series_ID = "1.3.12.2.1107.5.2.18.41548.30000014030519285935000000933"

reader = sitk.ImageSeriesReader()
fixed_image = sitk.ReadImage(reader.GetGDCMSeriesFileNames(data_directory, fixed_series_ID), sitk.sitkFloat32)
moving_image = sitk.ReadImage(reader.GetGDCMSeriesFileNames(data_directory, moving_series_ID), sitk.sitkFloat32)

# To provide a reasonable display we need to window/level the images. By default we could have used the intensity
# ranges found in the images [SimpleITK's StatisticsImageFilter], but these are not the best values for viewing.
# Using an external viewer we identified the following settings.
ct_window_level = [1727,-320]
mr_window_level = [355,178]

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

Fetching CIRS057A_MR_CT_DICOM/readme.txt

Register using centered transform initializer (assumes orientation is similar)

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

final_transform,_ = multires_registration(fixed_image, moving_image, initial_transform)

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

Visually evaluate the results using a linked cursor approach, a mouse click in one image will create the "corresponding" point in the other image. Don't be fooled by clicking on the "ribs" (symmetry is the bane of registration).

In [4]:
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);