Notebook

2. Image Registration with Multiple and Custom Parameter Maps¶

Registration¶

In [1]:

# First two imports are currently necessary to run ITKElastix on MacOs
from itk import itkElastixRegistrationMethodPython
import itk

We can run the registration with multiple stages, using the spatial transformation result from the current stage to initialize registration at the next stage.

Typically, we start with a simple spatial transformation and advance to a more complex spatial transformation.

The default, conservative, registration parameters progress from a translation -> affine -> bspline transformation, with 4 resolutions for each transformation type.

When multiple resolutions are used, a multi-scale image pyramid is generated with downscaled versions of the image. Registration results at a lower resolution is used to initialize registration at a higher resolution. This improves speed and robustness.

The number of resolutions can be manually specified aswell.

A default set of registration parameters are available for

translation
rigid
affine
bspline
spline
groupwise

transformations. More information on these transformation and on all possible parameter settings can be found in the elastix manual.

In [2]:

# Import Images
fixed_image = itk.imread('data/CT_2D_head_fixed.mha', itk.F)
moving_image = itk.imread('data/CT_2D_head_moving.mha', itk.F)

# Import Default Parameter Map
parameter_object = itk.ParameterObject.New()
resolutions = 3
parameter_map_rigid = parameter_object.GetDefaultParameterMap('rigid',3)
parameter_object.AddParameterMap(parameter_map_rigid)

# For the bspline default parameter map, an extra argument can be specified that define the final bspline grid spacing in physical space. 
parameter_map_bspline = parameter_object.GetDefaultParameterMap("bspline", resolutions, 20.0)
parameter_object.AddParameterMap(parameter_map_bspline)


# .. and/or load custom parameter maps from .txt file
parameter_object.AddParameterFile('data/parameters_BSpline.txt')

# ... and/or load and customize more
parameter_map_custom = parameter_object.GetDefaultParameterMap('rigid')
parameter_map_custom['Transform'] = ['BSplineTransform']
parameter_map_custom['Metric'] = ['AdvancedMattesMutualInformation']
parameter_object.AddParameterMap(parameter_map_custom)

# ... or customize parameter maps afterwards
# here the 'NumberOfResolutions' parameter of the 2nd parameter map of the parameter_object is set to 1.
parameter_object.SetParameter(2, "NumberOfResolutions", "1")

# here the 'DefaultPixelValue' of all parameter maps is in the parameter_object set to 0
parameter_object.SetParameter("DefaultPixelValue", "0")

# Remove a parameter
parameter_object.RemoveParameter("ResultImageFormat")

print(parameter_object)

ParameterObject (00000293D3DFE780)
  RTTI typeinfo:   class elastix::ParameterObject
  Reference Count: 1
  Modified Time: 60
  Debug: Off
  Object Name: 
  Observers: 
    none
ParameterMap 0: 
  (AutomaticParameterEstimation "true")
  (AutomaticScalesEstimation "true")
  (CheckNumberOfSamples "true")
  (DefaultPixelValue 0)
  (FinalBSplineInterpolationOrder 3)
  (FixedImagePyramid "FixedSmoothingImagePyramid")
  (ImageSampler "RandomCoordinate")
  (Interpolator "LinearInterpolator")
  (MaximumNumberOfIterations 256)
  (MaximumNumberOfSamplingAttempts 8)
  (Metric "AdvancedMattesMutualInformation")
  (MovingImagePyramid "MovingSmoothingImagePyramid")
  (NewSamplesEveryIteration "true")
  (NumberOfResolutions 3)
  (NumberOfSamplesForExactGradient 4096)
  (NumberOfSpatialSamples 2048)
  (Optimizer "AdaptiveStochasticGradientDescent")
  (Registration "MultiResolutionRegistration")
  (ResampleInterpolator "FinalBSplineInterpolator")
  (Resampler "DefaultResampler")
  (Transform "EulerTransform")
  (WriteIterationInfo "false")
  (WriteResultImage "true")
ParameterMap 1: 
  (AutomaticParameterEstimation "true")
  (CheckNumberOfSamples "true")
  (DefaultPixelValue 0)
  (FinalBSplineInterpolationOrder 3)
  (FinalGridSpacingInPhysicalUnits 20)
  (FixedImagePyramid "FixedSmoothingImagePyramid")
  (GridSpacingSchedule 1.9881 1.41 1)
  (ImageSampler "RandomCoordinate")
  (Interpolator "LinearInterpolator")
  (MaximumNumberOfIterations 256)
  (MaximumNumberOfSamplingAttempts 8)
  (Metric "AdvancedMattesMutualInformation" "TransformBendingEnergyPenalty")
  (Metric0Weight 1)
  (Metric1Weight 1)
  (MovingImagePyramid "MovingSmoothingImagePyramid")
  (NewSamplesEveryIteration "true")
  (NumberOfResolutions 3)
  (NumberOfSamplesForExactGradient 4096)
  (NumberOfSpatialSamples 2048)
  (Optimizer "AdaptiveStochasticGradientDescent")
  (Registration "MultiMetricMultiResolutionRegistration")
  (ResampleInterpolator "FinalBSplineInterpolator")
  (Resampler "DefaultResampler")
  (Transform "BSplineTransform")
  (WriteIterationInfo "false")
  (WriteResultImage "true")
ParameterMap 2: 
  (BSplineInterpolationOrder 1)
  (DefaultPixelValue 0)
  (ErodeMask "false")
  (FinalBSplineInterpolationOrder 3)
  (FinalGridSpacingInPhysicalUnits 16)
  (FixedImagePyramid "FixedRecursiveImagePyramid")
  (FixedInternalImagePixelType "float")
  (HowToCombineTransforms "Compose")
  (ImageSampler "Random")
  (Interpolator "BSplineInterpolator")
  (MaximumNumberOfIterations 500)
  (Metric "AdvancedMattesMutualInformation")
  (MovingImagePyramid "MovingRecursiveImagePyramid")
  (MovingInternalImagePixelType "float")
  (NewSamplesEveryIteration "true")
  (NumberOfHistogramBins 32)
  (NumberOfResolutions 1)
  (NumberOfSpatialSamples 2048)
  (Optimizer "AdaptiveStochasticGradientDescent")
  (Registration "MultiResolutionRegistration")
  (ResampleInterpolator "FinalBSplineInterpolator")
  (Resampler "DefaultResampler")
  (ResultImagePixelType "short")
  (Transform "BSplineTransform")
  (UseDirectionCosines "true")
  (WriteResultImage "true")
ParameterMap 3: 
  (AutomaticParameterEstimation "true")
  (AutomaticScalesEstimation "true")
  (CheckNumberOfSamples "true")
  (DefaultPixelValue 0)
  (FinalBSplineInterpolationOrder 3)
  (FixedImagePyramid "FixedSmoothingImagePyramid")
  (ImageSampler "RandomCoordinate")
  (Interpolator "LinearInterpolator")
  (MaximumNumberOfIterations 256)
  (MaximumNumberOfSamplingAttempts 8)
  (Metric "AdvancedMattesMutualInformation")
  (MovingImagePyramid "MovingSmoothingImagePyramid")
  (NewSamplesEveryIteration "true")
  (NumberOfResolutions 4)
  (NumberOfSamplesForExactGradient 4096)
  (NumberOfSpatialSamples 2048)
  (Optimizer "AdaptiveStochasticGradientDescent")
  (Registration "MultiResolutionRegistration")
  (ResampleInterpolator "FinalBSplineInterpolator")
  (Resampler "DefaultResampler")
  (Transform "BSplineTransform")
  (WriteIterationInfo "false")
  (WriteResultImage "true")

The parameter object or the seperate parameter maps can be saved or serialized in several ways

In [10]:

# Save custom parameter map
parameter_object.WriteParameterFile(
    parameter_map_custom, 'exampleoutput/parameters_custom.txt')

# Or serialize each parameter map to a file.
for index in range(parameter_object.GetNumberOfParameterMaps()):
    parameter_map = parameter_object.GetParameterMap(index)
    parameter_object.WriteParameterFile(parameter_map, "exampleoutput/Parameters.{0}.txt".format(index))

In [11]:

# To deserialize the parameters
parameter_files = ["exampleoutput/Parameters.{0}.txt".format(i) for i in range(4)]
parameter_object = itk.ParameterObject.New()
parameter_object.ReadParameterFile(parameter_files)
print(parameter_object)

ParameterObject (00000293EF258E70)
  RTTI typeinfo:   class elastix::ParameterObject
  Reference Count: 1
  Modified Time: 229588
  Debug: Off
  Object Name: 
  Observers: 
    none
ParameterMap 0: 
  (AutomaticParameterEstimation "true")
  (AutomaticScalesEstimation "true")
  (CheckNumberOfSamples "true")
  (DefaultPixelValue 0)
  (FinalBSplineInterpolationOrder 3)
  (FixedImagePyramid "FixedSmoothingImagePyramid")
  (ImageSampler "RandomCoordinate")
  (Interpolator "LinearInterpolator")
  (MaximumNumberOfIterations 256)
  (MaximumNumberOfSamplingAttempts 8)
  (Metric "AdvancedMattesMutualInformation")
  (MovingImagePyramid "MovingSmoothingImagePyramid")
  (NewSamplesEveryIteration "true")
  (NumberOfResolutions 3)
  (NumberOfSamplesForExactGradient 4096)
  (NumberOfSpatialSamples 2048)
  (Optimizer "AdaptiveStochasticGradientDescent")
  (Registration "MultiResolutionRegistration")
  (ResampleInterpolator "FinalBSplineInterpolator")
  (Resampler "DefaultResampler")
  (Transform "EulerTransform")
  (WriteIterationInfo "false")
  (WriteResultImage "true")
ParameterMap 1: 
  (AutomaticParameterEstimation "true")
  (CheckNumberOfSamples "true")
  (DefaultPixelValue 0)
  (FinalBSplineInterpolationOrder 3)
  (FinalGridSpacingInPhysicalUnits 20)
  (FixedImagePyramid "FixedSmoothingImagePyramid")
  (GridSpacingSchedule 1.9881 1.41 1)
  (ImageSampler "RandomCoordinate")
  (Interpolator "LinearInterpolator")
  (MaximumNumberOfIterations 256)
  (MaximumNumberOfSamplingAttempts 8)
  (Metric "AdvancedMattesMutualInformation" "TransformBendingEnergyPenalty")
  (Metric0Weight 1)
  (Metric1Weight 1)
  (MovingImagePyramid "MovingSmoothingImagePyramid")
  (NewSamplesEveryIteration "true")
  (NumberOfResolutions 3)
  (NumberOfSamplesForExactGradient 4096)
  (NumberOfSpatialSamples 2048)
  (Optimizer "AdaptiveStochasticGradientDescent")
  (Registration "MultiMetricMultiResolutionRegistration")
  (ResampleInterpolator "FinalBSplineInterpolator")
  (Resampler "DefaultResampler")
  (Transform "BSplineTransform")
  (WriteIterationInfo "false")
  (WriteResultImage "true")
ParameterMap 2: 
  (BSplineInterpolationOrder 1)
  (DefaultPixelValue 0)
  (ErodeMask "false")
  (FinalBSplineInterpolationOrder 3)
  (FinalGridSpacingInPhysicalUnits 16)
  (FixedImagePyramid "FixedRecursiveImagePyramid")
  (FixedInternalImagePixelType "float")
  (HowToCombineTransforms "Compose")
  (ImageSampler "Random")
  (Interpolator "BSplineInterpolator")
  (MaximumNumberOfIterations 500)
  (Metric "AdvancedMattesMutualInformation")
  (MovingImagePyramid "MovingRecursiveImagePyramid")
  (MovingInternalImagePixelType "float")
  (NewSamplesEveryIteration "true")
  (NumberOfHistogramBins 32)
  (NumberOfResolutions 1)
  (NumberOfSpatialSamples 2048)
  (Optimizer "AdaptiveStochasticGradientDescent")
  (Registration "MultiResolutionRegistration")
  (ResampleInterpolator "FinalBSplineInterpolator")
  (Resampler "DefaultResampler")
  (ResultImagePixelType "short")
  (Transform "BSplineTransform")
  (UseDirectionCosines "true")
  (WriteResultImage "true")
ParameterMap 3: 
  (AutomaticParameterEstimation "true")
  (AutomaticScalesEstimation "true")
  (CheckNumberOfSamples "true")
  (DefaultPixelValue 0)
  (FinalBSplineInterpolationOrder 3)
  (FixedImagePyramid "FixedSmoothingImagePyramid")
  (ImageSampler "RandomCoordinate")
  (Interpolator "LinearInterpolator")
  (MaximumNumberOfIterations 256)
  (MaximumNumberOfSamplingAttempts 8)
  (Metric "AdvancedMattesMutualInformation")
  (MovingImagePyramid "MovingSmoothingImagePyramid")
  (NewSamplesEveryIteration "true")
  (NumberOfResolutions 4)
  (NumberOfSamplesForExactGradient 4096)
  (NumberOfSpatialSamples 2048)
  (Optimizer "AdaptiveStochasticGradientDescent")
  (Registration "MultiResolutionRegistration")
  (ResampleInterpolator "FinalBSplineInterpolator")
  (Resampler "DefaultResampler")
  (Transform "BSplineTransform")
  (WriteIterationInfo "false")
  (WriteResultImage "true")

For more information on the available parameters and their meaning, see the Elastix Parameters documentation.

A database of parameter sets that work well for specific types of data can be found in the elastix parameter file database.

An Elastix parameter modelzoo, in sync with it's Github repo, which will eventually replace the database is currently being made. For progress check out https://github.com/ViktorvdValk/ElastixModelZoo

Registration can either be done in one line with the registration function...

In [8]:

# Call registration function
result_image, result_transform_parameters = itk.elastix_registration_method(
    fixed_image, moving_image,
    parameter_object=parameter_object,
    log_to_console=False)

.. or by initiating an elastix image filter object.

In [9]:

# Load Elastix Image Filter Object
elastix_object = itk.ElastixRegistrationMethod.New()
elastix_object.SetFixedImage(fixed_image)
elastix_object.SetMovingImage(moving_image)
elastix_object.SetParameterObject(parameter_object)

# Set additional options
elastix_object.SetLogToConsole(False)

# Update filter object (required)
elastix_object.UpdateLargestPossibleRegion()

# Results of Registration
result_image = elastix_object.GetOutput()
result_transform_parameters = elastix_object.GetTransformParameterObject()