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]:
library(SimpleITK)
# 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).
source("setup_for_testing.R")
# Utility method that either downloads data from the Girder repository or
# if already downloaded returns the file name for reading from disk (cached data).
source("downloaddata.R")
# Always write output to a separate directory, we don't want to pollute the source directory.
OUTPUT_DIR = "Output"
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]:
#
# Write the given transformation to file, resample the moving_image onto the fixed_images grid and save the
# result to file.
# @param transform (SimpleITK Transform): transform that maps points from the fixed image coordinate system to the moving.
# @param fixed_image (SimpleITK Image): resample onto the spatial grid defined by this image.
# @param moving_image (SimpleITK Image): resample this image.
# @param outputfile_prefix (string): transform is written to outputfile_prefix.tfm and resampled image is written to
# outputfile_prefix.mha.
#
save_transform_and_image <- function(transform, fixed_image, moving_image, outputfile_prefix)
{
resample <- ResampleImageFilter()
resample$SetReferenceImage(fixed_image)
# SimpleITK supports several interpolation options, we go with the simplest that gives reasonable results.
resample$SetInterpolator("sitkLinear")
resample$SetTransform(transform)
WriteImage(resample$Execute(moving_image), paste0(outputfile_prefix,".mha"))
WriteTransform(transform, paste0(outputfile_prefix,".tfm"))
}
#
# Get the DICOM tag values for a given file.
# @param file_name (string): Name of the DICOM file.
# @param tags (list(string)): List of strings containing the DICOM tags of interest. These are not the human readable
# values (i.e. They are of the form 0010|0010).
# @return (list(string)): List of strings containing the value of each of the DICOM tags.
#
get_DICOM_tag_values <- function(file_name, tags)
{
img <- ReadImage(file_name)
lapply(tags, function(tag) img$GetMetaData(tag))
}
#
# Get a dataframe containing DICOM tag information for the given set of series.
#
display_DICOM_information <- function(series_file_lists)
{
tags_to_print <- list(c('0010|0010', 'Patient name'),
c('0008|0060', 'Modality'),
c('0008|0021', 'Series date'),
c('0008|0031', 'Series time'),
c('0008|0070', 'Manufacturer'))
tags_to_print <- c('0010|0010', '0008|0060', '0008|0021', '0008|0031', '0008|0070')
names(tags_to_print) <- c('Patient name', 'Modality', 'Series date', 'Series time', 'Manufacturer')
data <- lapply(series_file_lists,
function(x, tag_list) get_DICOM_tag_values(x[1], tag_list),
tag_list=tags_to_print)
df <- as.data.frame(do.call(rbind, data))
df <- cbind(SeriesID=rownames(df), df)
rownames(df) <- NULL
return(df)
}
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 <- dirname(fetch_data("CIRS057A_MR_CT_DICOM/readme.txt"))
# Directory contains multiple DICOM studies/series, store the file names in a list
reader <- ImageSeriesReader()
series_IDs <- ImageSeriesReader_GetGDCMSeriesIDs(data_directory) #vector with strings/series IDs
names(series_IDs) <- series_IDs
series_file_names <- lapply(series_IDs,
function(sid, data_dir) ImageSeriesReader_GetGDCMSeriesFileNames(data_dir, sid),
data_dir=data_directory)
display_DICOM_information(series_file_names)
In [4]:
# Based on the DICOM image meta-information, select the fixed and moving images to register.
series_index_fixed_image <- "1.2.840.113619.2.290.3.3233817346.783.1399004564.515"
series_index_moving_image <- "1.3.12.2.1107.5.2.18.41548.30000014030519285935000000933"
# Actually read the data based on the user's selection.
reader$SetFileNames(series_file_names[[series_index_fixed_image]])
fixed_image <- reader$Execute()
reader$SetFileNames(series_file_names[[series_index_moving_image]])
moving_image <- reader$Execute()
# Save images to file and view overlap using external viewer.
WriteImage(fixed_image, file.path(OUTPUT_DIR, "fixedImage.mha"))
WriteImage(moving_image, file.path(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 <- CenteredTransformInitializer(Cast(fixed_image,moving_image$GetPixelID()),
moving_image,
Euler3DTransform(),
"GEOMETRY")
# Save moving image after initial transform and view overlap using external viewer.
save_transform_and_image(initial_transform, fixed_image, moving_image, file.path(OUTPUT_DIR, "initialAlignment"))
Look at the transformation, what type is it?
In [6]:
print(initial_transform)
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 <- ImageRegistrationMethod()
registration_method$SetMetricAsMattesMutualInformation(numberOfHistogramBins=50)
registration_method$SetMetricSamplingStrategy("RANDOM")
registration_method$SetMetricSamplingPercentage(0.01)
registration_method$SetInterpolator("sitkLinear")
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)
final_transform_v1 = registration_method$Execute(Cast(fixed_image, "sitkFloat32"),
Cast(moving_image, "sitkFloat32"))
In [8]:
cat(paste0("Optimizer\'s stopping condition, ",registration_method$GetOptimizerStopConditionDescription(),"\n"))
cat(paste0("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, file.path(OUTPUT_DIR, "finalAlignment-v1"))
Look at the final transformation, what type is it?
In [9]:
print(final_transform_v1)
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_{opt}(T_m(T_f^{-1}(^F\mathbf{p})))$
We now modify the previous example to use $T_{opt}$ and $T_m$.
In [10]:
registration_method <- ImageRegistrationMethod()
registration_method$SetMetricAsMattesMutualInformation(numberOfHistogramBins=50)
registration_method$SetMetricSamplingStrategy("RANDOM")
registration_method$SetMetricSamplingPercentage(0.01)
registration_method$SetInterpolator("sitkLinear")
registration_method$SetOptimizerAsGradientDescent(learningRate=1.0, numberOfIterations=100)
registration_method$SetOptimizerScalesFromPhysicalShift()
# Set the initial moving and optimized transforms.
optimized_transform <- Euler3DTransform()
registration_method$SetMovingInitialTransform(initial_transform)
registration_method$SetInitialTransform(optimized_transform)
dev_null <- registration_method$Execute(Cast(fixed_image, "sitkFloat32"),
Cast(moving_image, "sitkFloat32"))
# Need to compose the transformations after registration.
final_transform_v11 <- Transform(optimized_transform)
dev_null <- final_transform_v11$AddTransform(initial_transform)
In [11]:
cat(paste0("Optimizer\'s stopping condition, ",registration_method$GetOptimizerStopConditionDescription(),"\n"))
cat(paste0("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, file.path(OUTPUT_DIR, "finalAlignment-v1.1"))
Look at the final transformation, what type is it? Why is it different from the previous example?
In [12]:
print(final_transform_v11)
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 <- ImageRegistrationMethod()
registration_method$SetMetricAsMattesMutualInformation(numberOfHistogramBins=50)
registration_method$SetMetricSamplingStrategy("RANDOM")
registration_method$SetMetricSamplingPercentage(0.01)
registration_method$SetInterpolator("sitkLinear")
registration_method$SetOptimizerAsGradientDescent(learningRate=1.0, numberOfIterations=100) #, estimateLearningRate=registration_method.EachIteration)
registration_method$SetOptimizerScalesFromPhysicalShift()
final_transform <- Euler3DTransform(initial_transform)
registration_method$SetInitialTransform(final_transform)
registration_method$SetShrinkFactorsPerLevel(shrinkFactors = c(4,2,1))
registration_method$SetSmoothingSigmasPerLevel(smoothingSigmas = c(2,1,0))
registration_method$SmoothingSigmasAreSpecifiedInPhysicalUnitsOn()
dev_null <- registration_method$Execute(Cast(fixed_image, "sitkFloat32"),
Cast(moving_image, "sitkFloat32"))
In [14]:
cat(paste0("Optimizer\'s stopping condition, ",registration_method$GetOptimizerStopConditionDescription(),"\n"))
cat(paste0("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, file.path(OUTPUT_DIR, "finalAlignment-v2"))
Look at the final transformation, what type is it?
In [15]:
print(final_transform)