Introduction to SimpleITKv4 Registration - Continued¶
ITK v4 Registration Components¶
Before starting with this notebook, please go over the first introductory notebook found here.
In this notebook we will visually assess registration by viewing the overlap between images using external viewers. The two viewers we recommend for this task are ITK-SNAP and 3D Slicer. ITK-SNAP supports concurrent linked viewing between multiple instances of the program. 3D Slicer supports concurrent viewing of multiple volumes via alpha blending.
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
# Utility method that either downloads data from the network 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"
# GUI components (sliders, dropdown...).
from ipywidgets import interact, fixed
# Enable display of HTML.
from IPython.display import display, HTML
# Plots will be inlined.
%matplotlib inline
# Callbacks for plotting registration progress.
import registration_callbacks
Utility functions¶
A number of utility functions, saving a transform and corresponding resampled image, callback for selecting a DICOM series from several series found in the same directory.
In [2]:
def save_transform_and_image(transform, fixed_image, moving_image, outputfile_prefix):
"""
Write the given transformation to file, resample the moving_image onto the fixed_images grid and save the
result to file.
Args:
transform (SimpleITK Transform): transform that maps points from the fixed image coordinate system to the moving.
fixed_image (SimpleITK Image): resample onto the spatial grid defined by this image.
moving_image (SimpleITK Image): resample this image.
outputfile_prefix (string): transform is written to outputfile_prefix.tfm and resampled image is written to
outputfile_prefix.mha.
"""
resample = sitk.ResampleImageFilter()
resample.SetReferenceImage(fixed_image)
# SimpleITK supports several interpolation options, we go with the simplest that gives reasonable results.
resample.SetInterpolator(sitk.sitkLinear)
resample.SetTransform(transform)
sitk.WriteImage(resample.Execute(moving_image), outputfile_prefix + ".mha")
sitk.WriteTransform(transform, outputfile_prefix + ".tfm")
def DICOM_series_dropdown_callback(fixed_image, moving_image, series_dictionary):
"""
Callback from dropbox which selects the two series which will be used for registration.
The callback prints out some information about each of the series from the meta-data dictionary.
For a list of all meta-dictionary tags and their human readable names see DICOM standard part 6,
Data Dictionary (http://medical.nema.org/medical/dicom/current/output/pdf/part06.pdf)
"""
# The callback will update these global variables with the user selection.
global selected_series_fixed
global selected_series_moving
img_fixed = sitk.ReadImage(series_dictionary[fixed_image][0])
img_moving = sitk.ReadImage(series_dictionary[moving_image][0])
# There are many interesting tags in the DICOM data dictionary, display a selected few.
tags_to_print = {
"0010|0010": "Patient name: ",
"0008|0060": "Modality: ",
"0008|0021": "Series date: ",
"0008|0031": "Series time:",
"0008|0070": "Manufacturer: ",
}
html_table = []
html_table.append(
"<table><tr><td><b>Tag</b></td><td><b>Fixed Image</b></td><td><b>Moving Image</b></td></tr>"
)
for tag in tags_to_print:
fixed_tag = ""
moving_tag = ""
try:
fixed_tag = img_fixed.GetMetaData(tag)
except: # ignore if the tag isn't in the dictionary
pass
try:
moving_tag = img_moving.GetMetaData(tag)
except: # ignore if the tag isn't in the dictionary
pass
html_table.append(
"<tr><td>"
+ tags_to_print[tag]
+ "</td><td>"
+ fixed_tag
+ "</td><td>"
+ moving_tag
+ "</td></tr>"
)
html_table.append("</table>")
display(HTML("".join(html_table)))
selected_series_fixed = fixed_image
selected_series_moving = moving_image
Loading Data¶
In this notebook we will work with CT and MR scans of the CIRS 057A multi-modality abdominal phantom. The scans are multi-slice DICOM images. The data is stored in a zip archive which is automatically retrieved and extracted when we request a file which is part of the archive.
In [3]:
data_directory = os.path.dirname(fdata("CIRS057A_MR_CT_DICOM/readme.txt"))
# 'selected_series_moving/fixed' will be updated by the interact function.
selected_series_fixed = ""
selected_series_moving = ""
# Directory contains multiple DICOM studies/series, store the file names
# in dictionary with the key being the series ID.
reader = sitk.ImageSeriesReader()
series_file_names = {}
series_IDs = list(reader.GetGDCMSeriesIDs(data_directory)) # list of all series
if series_IDs: # check that we have at least one series
for series in series_IDs:
series_file_names[series] = reader.GetGDCMSeriesFileNames(
data_directory, series
)
interact(
DICOM_series_dropdown_callback,
fixed_image=series_IDs,
moving_image=series_IDs,
series_dictionary=fixed(series_file_names),
)
else:
print("This is surprising, data directory does not contain any DICOM series.")
Fetching CIRS057A_MR_CT_DICOM/readme.txt
In [4]:
# Actually read the data based on the user's selection.
fixed_image = sitk.ReadImage(series_file_names[selected_series_fixed])
moving_image = sitk.ReadImage(series_file_names[selected_series_moving])
# Save images to file and view overlap using external viewer.
sitk.WriteImage(fixed_image, os.path.join(OUTPUT_DIR, "fixedImage.mha"))
sitk.WriteImage(moving_image, os.path.join(OUTPUT_DIR, "preAlignment.mha"))
Initial Alignment¶
A reasonable guesstimate for the initial translational alignment can be obtained by using the CenteredTransformInitializer (functional interface to the CenteredTransformInitializerFilter).
The resulting transformation is centered with respect to the fixed image and the translation aligns the centers of the two images. There are two options for defining the centers of the images, either the physical centers of the two data sets (GEOMETRY), or the centers defined by the intensity moments (MOMENTS).
Two things to note about this filter, it requires the fixed and moving image have the same type even though it is not algorithmically required, and its return type is the generic SimpleITK.Transform.
In [5]:
initial_transform = sitk.CenteredTransformInitializer(
sitk.Cast(fixed_image, moving_image.GetPixelID()),
moving_image,
sitk.Euler3DTransform(),
sitk.CenteredTransformInitializerFilter.GEOMETRY,
)
# Save moving image after initial transform and view overlap using external viewer.
save_transform_and_image(
initial_transform,
fixed_image,
moving_image,
os.path.join(OUTPUT_DIR, "initialAlignment"),
)
Look at the transformation, what type is it?
In [6]:
print(initial_transform)
itk::simple::Euler3DTransform
Euler3DTransform (0x7f87e8f0d8f0)
RTTI typeinfo: itk::Euler3DTransform<double>
Reference Count: 1
Modified Time: 23371
Debug: Off
Object Name:
Observers:
none
Matrix:
1 0 0
0 1 0
0 0 1
Offset: [0, 0, 0]
Center: [-0.328125, -0.328125, -106.875]
Translation: [0, 0, 0]
Inverse:
1 0 0
0 1 0
0 0 1
Singular: 0
AngleX: 0
AngleY: 0
AngleZ: 0
ComputeZYX: Off
Final registration¶
Version 1¶
- Single scale (not using image pyramid).
- Initial transformation is not modified in place.
-
Illustrate the need for scaling the step size differently for each parameter:
- SetOptimizerScalesFromIndexShift - estimated from maximum shift of voxel indexes (only use if data is isotropic).
- SetOptimizerScalesFromPhysicalShift - estimated from maximum shift of physical locations of voxels.
- SetOptimizerScalesFromJacobian - estimated from the averaged squared norm of the Jacobian w.r.t. parameters.
- Look at the optimizer's stopping condition to ensure we have not terminated prematurely.
In [7]:
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 learningRate parameter is always required. Using the default
# configuration this parameter is ignored because it is overridden
# by the default setting of the estimateLearningRate parameter which
# is sitk.ImageRegistrationMethod.Once. For the user selected
# learningRate to take effect you need to also set the
# estimateLearningRate parameter to sitk.ImageRegistrationMethod.Never
registration_method.SetOptimizerAsGradientDescent(
learningRate=1.0, numberOfIterations=100
)
# Scale the step size differently for each parameter, this is critical!!!
registration_method.SetOptimizerScalesFromPhysicalShift()
registration_method.SetInitialTransform(initial_transform, inPlace=False)
registration_method.AddCommand(
sitk.sitkStartEvent, registration_callbacks.metric_start_plot
)
registration_method.AddCommand(
sitk.sitkEndEvent, registration_callbacks.metric_end_plot
)
registration_method.AddCommand(
sitk.sitkIterationEvent,
lambda: registration_callbacks.metric_plot_values(registration_method),
)
final_transform_v1 = registration_method.Execute(
sitk.Cast(fixed_image, sitk.sitkFloat32), sitk.Cast(moving_image, sitk.sitkFloat32)
)
In [8]:
print(
f"Optimizer's stopping condition, {registration_method.GetOptimizerStopConditionDescription()}"
)
print(f"Final metric value: {registration_method.GetMetricValue()}")
# Save moving image after registration and view overlap using external viewer.
save_transform_and_image(
final_transform_v1,
fixed_image,
moving_image,
os.path.join(OUTPUT_DIR, "finalAlignment-v1"),
)
Optimizer's stopping condition, GradientDescentOptimizerv4Template: Convergence checker passed at iteration 9. Final metric value: -0.9062210635475781
Look at the final transformation, what type is it?
In [9]:
print(final_transform_v1)
itk::simple::CompositeTransform
CompositeTransform (0x7f87d8fa56a0)
RTTI typeinfo: itk::CompositeTransform<double, 3u>
Reference Count: 1
Modified Time: 1508806
Debug: Off
Object Name:
Observers:
none
TransformQueue:
>>>>>>>>>
Euler3DTransform (0x7f87a8f2bee0)
RTTI typeinfo: itk::Euler3DTransform<double>
Reference Count: 1
Modified Time: 1508652
Debug: Off
Object Name:
Observers:
none
Matrix:
1 0.000194693 -1.88627e-05
-0.000194703 1 -0.000517573
1.87619e-05 0.000517576 1
Offset: [0.012608, -4.07142, 1.035]
Center: [-0.328125, -0.328125, -106.875]
Translation: [0.0145601, -4.01604, 1.03484]
Inverse:
1 -0.000194703 1.87619e-05
0.000194693 1 0.000517576
-1.88627e-05 -0.000517573 1
Singular: 0
AngleX: 0.000517576
AngleY: -1.87619e-05
AngleZ: -0.000194693
ComputeZYX: Off
TransformsToOptimizeFlags:
1
TransformsToOptimizeQueue:
PreviousTransformsToOptimizeUpdateTime: 0
Version 1.1¶
The previous example illustrated the use of the ITK v4 registration framework in an ITK v3 manner. We only referred to a single transformation which was what we optimized.
In ITK v4 the registration method accepts three transformations (if you look at the diagram above you will only see two transformations, Moving transform represents $T_{opt} \circ T_m$):
- 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_m(T_{opt}(T_f^{-1}(^F\mathbf{p})))$
We now modify the previous example to use $T_{opt}$ and $T_m$.
In [10]:
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
)
registration_method.SetOptimizerScalesFromPhysicalShift()
# Set the initial moving and optimized transforms.
optimized_transform = sitk.Euler3DTransform()
registration_method.SetMovingInitialTransform(initial_transform)
registration_method.SetInitialTransform(optimized_transform)
registration_method.AddCommand(
sitk.sitkStartEvent, registration_callbacks.metric_start_plot
)
registration_method.AddCommand(
sitk.sitkEndEvent, registration_callbacks.metric_end_plot
)
registration_method.AddCommand(
sitk.sitkIterationEvent,
lambda: registration_callbacks.metric_plot_values(registration_method),
)
registration_method.Execute(
sitk.Cast(fixed_image, sitk.sitkFloat32), sitk.Cast(moving_image, sitk.sitkFloat32)
)
# Need to compose the transformations after registration.
final_transform_v11 = sitk.CompositeTransform(optimized_transform)
final_transform_v11.AddTransform(initial_transform)
In [11]:
print(
f"Optimizer's stopping condition, {registration_method.GetOptimizerStopConditionDescription()}"
)
print(f"Final metric value: {registration_method.GetMetricValue()}")
# Save moving image after registration and view overlap using external viewer.
save_transform_and_image(
final_transform_v11,
fixed_image,
moving_image,
os.path.join(OUTPUT_DIR, "finalAlignment-v1.1"),
)
Optimizer's stopping condition, GradientDescentOptimizerv4Template: Convergence checker passed at iteration 9. Final metric value: -1.0332024329467742
Look at the final transformation, what type is it? Why is it different from the previous example?
In [12]:
print(final_transform_v11)
itk::simple::CompositeTransform
CompositeTransform (0x7f87d8f1e800)
RTTI typeinfo: itk::CompositeTransform<double, 3u>
Reference Count: 1
Modified Time: 2994265
Debug: Off
Object Name:
Observers:
none
TransformQueue:
>>>>>>>>>
Euler3DTransform (0x7f8779587390)
RTTI typeinfo: itk::Euler3DTransform<double>
Reference Count: 3
Modified Time: 2994097
Debug: Off
Object Name:
Observers:
none
Matrix:
1 -0.000275379 6.40485e-05
0.000275418 1 -0.000605101
-6.38819e-05 0.000605119 1
Offset: [0.651138, -2.12368, 0.492245]
Center: [-0.328125, -0.328125, -106.875]
Translation: [0.644383, -2.0591, 0.492087]
Inverse:
1 0.000275418 -6.38819e-05
-0.000275379 1 0.000605119
6.40485e-05 -0.000605101 1
Singular: 0
AngleX: 0.000605119
AngleY: 6.38819e-05
AngleZ: 0.000275379
ComputeZYX: Off
>>>>>>>>>
Euler3DTransform (0x7f87d8f46ed0)
RTTI typeinfo: itk::Euler3DTransform<double>
Reference Count: 1
Modified Time: 2994261
Debug: Off
Object Name:
Observers:
none
Matrix:
1 0 0
0 1 0
0 0 1
Offset: [0, 0, 0]
Center: [-0.328125, -0.328125, -106.875]
Translation: [0, 0, 0]
Inverse:
1 0 0
0 1 0
0 0 1
Singular: 0
AngleX: 0
AngleY: 0
AngleZ: 0
ComputeZYX: Off
TransformsToOptimizeFlags:
0 1
TransformsToOptimizeQueue:
PreviousTransformsToOptimizeUpdateTime: 0
Version 2¶
- Multi scale - specify both scale, and how much to smooth with respect to original image.
- Initial transformation modified in place, so in the end we have the same type of transformation in hand.
In [13]:
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.EachIteration)
registration_method.SetOptimizerScalesFromPhysicalShift()
final_transform = sitk.Euler3DTransform(initial_transform)
registration_method.SetInitialTransform(final_transform)
registration_method.SetShrinkFactorsPerLevel(shrinkFactors=[4, 2, 1])
registration_method.SetSmoothingSigmasPerLevel(smoothingSigmas=[2, 1, 0])
registration_method.SmoothingSigmasAreSpecifiedInPhysicalUnitsOn()
registration_method.AddCommand(
sitk.sitkStartEvent, registration_callbacks.metric_start_plot
)
registration_method.AddCommand(
sitk.sitkEndEvent, registration_callbacks.metric_end_plot
)
registration_method.AddCommand(
sitk.sitkMultiResolutionIterationEvent,
registration_callbacks.metric_update_multires_iterations,
)
registration_method.AddCommand(
sitk.sitkIterationEvent,
lambda: registration_callbacks.metric_plot_values(registration_method),
)
registration_method.Execute(
sitk.Cast(fixed_image, sitk.sitkFloat32), sitk.Cast(moving_image, sitk.sitkFloat32)
)
Out[13]:
<SimpleITK.SimpleITK.Euler3DTransform; proxy of <Swig Object of type 'itk::simple::Euler3DTransform *' at 0x116007540> >
In [14]:
print(
f"Optimizer's stopping condition, {registration_method.GetOptimizerStopConditionDescription()}"
)
print(f"Final metric value: {registration_method.GetMetricValue()}")
# Save moving image after registration and view overlap using external viewer.
save_transform_and_image(
final_transform,
fixed_image,
moving_image,
os.path.join(OUTPUT_DIR, "finalAlignment-v2"),
)
Optimizer's stopping condition, GradientDescentOptimizerv4Template: Convergence checker passed at iteration 9. Final metric value: -0.7964225864936663
Look at the final transformation, what type is it?
In [15]:
print(final_transform)
itk::simple::Euler3DTransform
Euler3DTransform (0x7f877954c600)
RTTI typeinfo: itk::Euler3DTransform<double>
Reference Count: 3
Modified Time: 4689431
Debug: Off
Object Name:
Observers:
none
Matrix:
0.999997 -0.00241544 0.000157833
0.0024157 0.999996 -0.00171882
-0.000153681 0.0017192 0.999999
Offset: [4.81257, -3.66813, -3.8047]
Center: [-0.328125, -0.328125, -106.875]
Translation: [4.79649, -3.48522, -3.80505]
Inverse:
0.999997 0.0024157 -0.000153681
-0.00241544 0.999996 0.0017192
0.000157833 -0.00171882 0.999999
Singular: 0
AngleX: 0.0017192
AngleY: 0.000153681
AngleZ: 0.00241544
ComputeZYX: Off