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:

Do nothing (hope springs eternal) - initialize using the identity transformation.
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.
Use a sampling of the parameter space (useful mostly for low dimensional parameter spaces).
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 [ ]:

    
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 [ ]:

    
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]);

Register using centered transform initializer (assumes orientation is similar)



In [ ]:

    
initial_transform = sitk.CenteredTransformInitializer(fixed_image, 
                                                      moving_image, 
                                                      sitk.Euler3DTransform(), 
                                                      sitk.CenteredTransformInitializerFilter.GEOMETRY)

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

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 [ ]:

    
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);

Register using sampling of the parameter space

As we want to account for significant orientation differences due to erroneous patient position (HFS...) we evaluate the similarity measure at locations corresponding to the various orientation differences. This can be done in two ways which will be illustrated below:

Use the ImageRegistrationMethod.MetricEvaluate() method.
Use the Exhaustive optimizer.

The former approach is more computationally intensive as it constructs and configures a metric object each time it is invoked. It is therefore more appropriate for use if the set of parameter values we want to evaluate are not on a rectilinear grid in the parameter space. The latter approach is appropriate if the set of parameter values are on a rectilinear grid, in which case the approach is more computationally efficient.

In both cases we use the CenteredTransformInitializer to obtain the initial translation.

MetricEvaluate

To use the MetricEvaluate method we create a ImageRegistrationMethod, set its metric and interpolator. We then iterate over all parameter settings, set the initial transform and evaluate the metric. The minimal similarity measure value corresponds to the best parameter settings.



In [ ]:

    
# Dictionary with all the orientations we will try. We omit the identity (x=0, y=0, z=0) as we always use it. This
# set of rotations is arbitrary. For a complete grid coverage we would naively have 64 entries 
# (0, pi/2, pi, 1.5pi for each angle), but we know better, there are only 24 unique rotation matrices defined by
# these parameter value combinations.
all_orientations = {'x=0, y=0, z=180': (0.0,0.0,np.pi),
                    'x=0, y=180, z=0': (0.0,np.pi, 0.0),
                    'x=0, y=180, z=180': (0.0,np.pi, np.pi)}    

# Registration framework setup.
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)

# Evaluate the similarity metric using the rotation parameter space sampling, translation remains the same for all.
initial_transform = sitk.Euler3DTransform(sitk.CenteredTransformInitializer(fixed_image, 
                                                                            moving_image, 
                                                                            sitk.Euler3DTransform(), 
                                                                            sitk.CenteredTransformInitializerFilter.GEOMETRY))
registration_method.SetInitialTransform(initial_transform, inPlace=False)
best_orientation = (0.0,0.0,0.0)
best_similarity_value = registration_method.MetricEvaluate(fixed_image, moving_image)

# Iterate over all other rotation parameter settings. 
for key, orientation in all_orientations.items():
    initial_transform.SetRotation(*orientation)
    registration_method.SetInitialTransform(initial_transform)
    current_similarity_value = registration_method.MetricEvaluate(fixed_image, moving_image)
    if current_similarity_value < best_similarity_value:
        best_similarity_value = current_similarity_value
        best_orientation = orientation
print('best orientation is: ' + str(best_orientation))

Why loop if you can process in parallel

As the metric evaluations are independent of each other, we can easily replace looping with parallelization.



In [ ]:

    
from multiprocessing.pool import ThreadPool
from functools import partial

# This function evaluates the metric value in a thread safe manner
def evaluate_metric(current_rotation, tx, f_image, m_image):
    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)
    current_transform = sitk.Euler3DTransform(tx)
    current_transform.SetRotation(*current_rotation)
    registration_method.SetInitialTransform(current_transform)
    res = registration_method.MetricEvaluate(f_image, m_image)
    return res

p = ThreadPool(len(all_orientations)+1)
orientations_list = [(0,0,0)] + list(all_orientations.values())
all_metric_values = p.map(partial(evaluate_metric, 
                                  tx = initial_transform, 
                                  f_image = fixed_image,
                                  m_image = moving_image),
                          orientations_list)
best_orientation = orientations_list[np.argmin(all_metric_values)]
print('best orientation is: ' + str(best_orientation))



In [ ]:

    
initial_transform.SetRotation(*best_orientation)
final_transform,_ = multires_registration(fixed_image, moving_image, initial_transform)

Visually evaluate the registration results:



In [ ]:

    
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);

Exhaustive optimizer

The exhaustive optimizer evaluates the similarity metric on a grid in parameter space centered on the parameters of the initial transform. This grid is defined using three elements:

numberOfSteps.
stepLength.
optimizer scales.

The similarity metric is evaluated on the resulting parameter grid: initial_parameters ± numberOfSteps × stepLength × optimizerScales

Example:

numberOfSteps=[1,0,2,0,0,0]
stepLength = np.pi
optimizer scales = [1,1,0.5,1,1,1]

Will perform 15 metric evaluations ($\displaystyle\prod_i (2*numberOfSteps[i] + 1)$).

The parameter values for the second parameter and the last three parameters are the initial parameter values. The parameter values for the first parameter are $v_{init}-\pi, v_{init}, v_{init}+\pi$ and the parameter values for the third parameter are $v_{init}-\pi, v_{init}-\pi/2, v_{init}, v_{init}+\pi/2, v_{init}+\pi$.

The transformation corresponding to the lowest similarity metric is returned.

Using this approach we have superfluous evaluations, due to the symmetry of the grid in parameter space. On the other hand this method is often faster than evaluating the metric using the MetricEvaluate method (due to the setup and tear down time).



In [ ]:

    
initial_transform = sitk.CenteredTransformInitializer(fixed_image, 
                                                      moving_image, 
                                                      sitk.Euler3DTransform(), 
                                                      sitk.CenteredTransformInitializerFilter.GEOMETRY)
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)
# The order of parameters for the Euler3DTransform is [angle_x, angle_y, angle_z, t_x, t_y, t_z]. The parameter 
# sampling grid is centered on the initial_transform parameter values, that are all zero for the rotations. Given
# the number of steps and their length and optimizer scales we have:
# angle_x = 0
# angle_y = -pi, 0, pi
# angle_z = -pi, 0, pi
registration_method.SetOptimizerAsExhaustive(numberOfSteps=[0,1,1,0,0,0], stepLength = np.pi)
registration_method.SetOptimizerScales([1,1,1,1,1,1])

#Perform the registration in-place so that the initial_transform is modified.
registration_method.SetInitialTransform(initial_transform, inPlace=True)
registration_method.Execute(fixed_image, moving_image)

print('best initial transformation is: ' + str(initial_transform.GetParameters()))

Run the registration and visually evaluate the results:



In [ ]:

    
final_transform, _ = multires_registration(fixed_image, moving_image, initial_transform)
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);

Exhaustive optimizer - an exploration-exploitation view

In the example above we used the exhaustive optimizer to obtain an initial value for our transformation parameter values. This approach has two limitations:

It assumes we only want to sample the parameter space using a continuous regular grid. This is not appropriate if we want to explore multiple discontinuous regions of the parameter space (e.g. tx = [0.5, 1.0, 1.5] and tx = [ 12.0, 12.5, 13.0] and...).
It assumes that the parameter values corresponding to the best metric value from our sample will enable convergence to the desired optimum. To quote George Gershwin - "It Ain't Necessarily So".

If we consider the exhaustive optimizer in the context of the exploration-exploitation heuristic framework, we first search the parameter space for promising solution(s) and then refine the solution(s), we can readily address these limitations:

Explore multiple discontinuous regions using a single or multiple instances of the ExhaustiveOptimizer or use the MetricEvaluate approach if the samples do not define a regular grid.
Obtain all of the parameter space samples and exploit (run final registration) for each of the k most promising solutions.

Below we implement the latter: (a) We explore the parameter space using the exhaustive optimizer's callback mechanism to obtain all of the parameter values and their corresponding metric values. (b) We exploit the k_most_promising parameter values to obtain the final transformation.

NOTE: This is a heuristic and only increases the probability of convergence to the desired optimum. In some cases this approach may be detrimental when the parameter values we seek correspond to a local optimum and not a global one.



In [ ]:

    
#
# Exploration step.
#
def start_observer():
    global metricvalue_parameters_list 
    metricvalue_parameters_list = []

def iteration_observer(registration_method):    
    metricvalue_parameters_list.append((registration_method.GetMetricValue(), registration_method.GetOptimizerPosition()))
    
initial_transform = sitk.CenteredTransformInitializer(fixed_image, 
                                                      moving_image, 
                                                      sitk.Euler3DTransform(), 
                                                      sitk.CenteredTransformInitializerFilter.GEOMETRY)
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)
# The order of parameters for the Euler3DTransform is [angle_x, angle_y, angle_z, t_x, t_y, t_z]. The parameter 
# sampling grid is centered on the initial_transform parameter values, that are all zero for the rotations. Given
# the number of steps and their length and optimizer scales we have:
# angle_x = 0
# angle_y = -pi, 0, pi
# angle_z = -pi, 0, pi
registration_method.SetOptimizerAsExhaustive(numberOfSteps=[0,1,1,0,0,0], stepLength = np.pi)
registration_method.SetOptimizerScales([1,1,1,1,1,1])

#We don't really care if transformation is modified in place or not, we will select the k 
#best transformations from the parameters_metricvalue_list.
registration_method.SetInitialTransform(initial_transform, inPlace=True)

registration_method.AddCommand(sitk.sitkIterationEvent, lambda: iteration_observer(registration_method))
registration_method.AddCommand(sitk.sitkStartEvent, start_observer )
_ = registration_method.Execute(fixed_image, moving_image)



In [ ]:

    
#
# Exploitation step.
#

#Sort our list from most to least promising solutions (low to high metric values). 
metricvalue_parameters_list.sort(key=lambda x: x[0]) 

# We exploit the k_most_promising parameter value settings.
k_most_promising = min(3, len(metricvalue_parameters_list))
final_results = []
for metricvalue, parameters in metricvalue_parameters_list[0:k_most_promising]:
    initial_transform.SetParameters(parameters)
    final_results.append(multires_registration(fixed_image, moving_image, initial_transform))

final_transform, _ = min(final_results, key=lambda x: x[1]) 
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);

Register using manual initialization

When all else fails, a human in the loop will almost always be able to robustly initialize the registration.

In the example below we identify corresponding points to compute an initial rigid transformation. You will need to click on corresponding points in each of the images, going back and forth between them. The interface will "force" you to create point pairs (you will not be able to add multiple points in one image).

Note:

There is no correspondence between the fiducial markers on the phantom.
After localizing points in the GUI below, comment out the hard-coded point data, two cells below, which is there FOR TESTING.



In [ ]:

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



In [ ]:

    
# Get the manually specified points and compute the transformation.

fixed_image_points, moving_image_points = point_acquisition_interface.get_points()

# FOR TESTING: previously localized points
fixed_image_points = [(24.062587103074605, 14.594981536981521, -58.75), 
                      (6.178716135332678, 53.93949766601378, -58.75), 
                      (74.14383149714774, -69.04462737237648, -76.25), 
                      (109.74899278747029, -14.905272533666817, -76.25)]
moving_image_points = [(4.358707846364581, 60.46357110706131, -71.53120422363281), 
                       (24.09010295252645, 98.21840981673873, -71.53120422363281), 
                       (-52.11888008581127, -26.57984635768439, -58.53120422363281), 
                       (-87.46150681392184, 28.73904765153219, -58.53120422363281)]

fixed_image_points_flat = [c for p in fixed_image_points for c in p]        
moving_image_points_flat = [c for p in moving_image_points for c in p]
initial_transform = sitk.LandmarkBasedTransformInitializer(sitk.VersorRigid3DTransform(), 
                                                                fixed_image_points_flat, 
                                                                moving_image_points_flat)


print('manual initial transformation is: ' + str(initial_transform.GetParameters()))

Run the registration and visually evaluate the results:



In [ ]:

    
final_transform = multires_registration(fixed_image, moving_image, initial_transform)
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);



In [ ]: