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:
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:
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.
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())
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.
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]);
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).
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);