Download this example

Perform 2D Translation Registration With Mean Squares

Synopsis

This example illustrates the use of the image registration framework in Insight. It should be read as a Hello World for ITK registration. Instead of means to an end, this example should be read as a basic introduction to the elements typically involved when solving a problem of image registration.

A registration method requires the following set of components: two input images, a transform, a metric and an optimizer. Some of these components are parameterized by the image type for which the registration is intended. The following header files provide declarations of common types used for these components.

This example corresponds to the ImageRegistration1.cxx example from the ITK software guide.

Results

Input images (fixed image left/ moving image right).

Output registered image.

Image difference between moving and fixed image (left before registration, right after registration).

Code

Python

#!/usr/bin/env python

import sys
import itk
import argparse

from distutils.version import StrictVersion as VS

if VS(itk.Version.GetITKVersion()) < VS("4.9.0"):
    print("ITK 4.9.0 is required.")
    sys.exit(1)

parser = argparse.ArgumentParser(
    description="Perform 2D Translation Registration With Mean Squares."
)
parser.add_argument("fixed_input_image")
parser.add_argument("moving_input_image")
parser.add_argument("output_image")
parser.add_argument("difference_image_after")
parser.add_argument("difference_image_before")
args = parser.parse_args()

PixelType = itk.ctype("float")

fixedImage = itk.imread(args.fixed_input_image, PixelType)
movingImage = itk.imread(args.moving_input_image, PixelType)

Dimension = fixedImage.GetImageDimension()
FixedImageType = itk.Image[PixelType, Dimension]
MovingImageType = itk.Image[PixelType, Dimension]

TransformType = itk.TranslationTransform[itk.D, Dimension]
initialTransform = TransformType.New()

optimizer = itk.RegularStepGradientDescentOptimizerv4.New(
    LearningRate=4,
    MinimumStepLength=0.001,
    RelaxationFactor=0.5,
    NumberOfIterations=200,
)

metric = itk.MeanSquaresImageToImageMetricv4[FixedImageType, MovingImageType].New()

registration = itk.ImageRegistrationMethodv4.New(
    FixedImage=fixedImage,
    MovingImage=movingImage,
    Metric=metric,
    Optimizer=optimizer,
    InitialTransform=initialTransform,
)

movingInitialTransform = TransformType.New()
initialParameters = movingInitialTransform.GetParameters()
initialParameters[0] = 0
initialParameters[1] = 0
movingInitialTransform.SetParameters(initialParameters)
registration.SetMovingInitialTransform(movingInitialTransform)

identityTransform = TransformType.New()
identityTransform.SetIdentity()
registration.SetFixedInitialTransform(identityTransform)

registration.SetNumberOfLevels(1)
registration.SetSmoothingSigmasPerLevel([0])
registration.SetShrinkFactorsPerLevel([1])

registration.Update()

transform = registration.GetTransform()
finalParameters = transform.GetParameters()
translationAlongX = finalParameters.GetElement(0)
translationAlongY = finalParameters.GetElement(1)

numberOfIterations = optimizer.GetCurrentIteration()

bestValue = optimizer.GetValue()

print("Result = ")
print(" Translation X = " + str(translationAlongX))
print(" Translation Y = " + str(translationAlongY))
print(" Iterations    = " + str(numberOfIterations))
print(" Metric value  = " + str(bestValue))

CompositeTransformType = itk.CompositeTransform[itk.D, Dimension]
outputCompositeTransform = CompositeTransformType.New()
outputCompositeTransform.AddTransform(movingInitialTransform)
outputCompositeTransform.AddTransform(registration.GetModifiableTransform())

resampler = itk.ResampleImageFilter.New(
    Input=movingImage,
    Transform=outputCompositeTransform,
    UseReferenceImage=True,
    ReferenceImage=fixedImage,
)
resampler.SetDefaultPixelValue(100)

OutputPixelType = itk.ctype("unsigned char")
OutputImageType = itk.Image[OutputPixelType, Dimension]

caster = itk.CastImageFilter[FixedImageType, OutputImageType].New(Input=resampler)

writer = itk.ImageFileWriter.New(Input=caster, FileName=args.output_image)
writer.Update()

difference = itk.SubtractImageFilter.New(Input1=fixedImage, Input2=resampler)

intensityRescaler = itk.RescaleIntensityImageFilter[
    FixedImageType, OutputImageType
].New(
    Input=difference,
    OutputMinimum=itk.NumericTraits[OutputPixelType].min(),
    OutputMaximum=itk.NumericTraits[OutputPixelType].max(),
)

resampler.SetDefaultPixelValue(1)
writer.SetInput(intensityRescaler.GetOutput())
writer.SetFileName(args.difference_image_after)
writer.Update()

resampler.SetTransform(identityTransform)
writer.SetFileName(args.difference_image_before)
writer.Update()

C++

#include "itkImageRegistrationMethodv4.h"
#include "itkTranslationTransform.h"
#include "itkMeanSquaresImageToImageMetricv4.h"
#include "itkRegularStepGradientDescentOptimizerv4.h"

#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"

#include "itkResampleImageFilter.h"
#include "itkCastImageFilter.h"
#include "itkRescaleIntensityImageFilter.h"
#include "itkSubtractImageFilter.h"


int
main(int argc, char * argv[])
{
  if (argc != 6)
  {
    std::cerr << "Missing Parameters " << std::endl;
    std::cerr << "Usage: " << argv[0];
    std::cerr << " fixedImageFile  movingImageFile ";
    std::cerr << "outputImagefile [differenceImageAfter]";
    std::cerr << "[differenceImageBefore]" << std::endl;
    return EXIT_FAILURE;
  }

  const char * fixedImageFile = argv[1];
  const char * movingImageFile = argv[2];
  const char * outputImageFile = argv[3];
  const char * differenceImageAfterFile = argv[4];
  const char * differenceImageBeforeFile = argv[5];

  constexpr unsigned int Dimension = 2;
  using PixelType = float;

  using FixedImageType = itk::Image<PixelType, Dimension>;
  using MovingImageType = itk::Image<PixelType, Dimension>;

  const auto fixedImage = itk::ReadImage<FixedImageType>(fixedImageFile);
  const auto movingImage = itk::ReadImage<MovingImageType>(movingImageFile);


  using TransformType = itk::TranslationTransform<double, Dimension>;
  auto initialTransform = TransformType::New();

  using OptimizerType = itk::RegularStepGradientDescentOptimizerv4<double>;
  auto optimizer = OptimizerType::New();
  optimizer->SetLearningRate(4);
  optimizer->SetMinimumStepLength(0.001);
  optimizer->SetRelaxationFactor(0.5);
  optimizer->SetNumberOfIterations(200);

  using MetricType = itk::MeanSquaresImageToImageMetricv4<FixedImageType, MovingImageType>;
  auto metric = MetricType::New();

  using RegistrationType = itk::ImageRegistrationMethodv4<FixedImageType, MovingImageType>;
  auto registration = RegistrationType::New();
  registration->SetInitialTransform(initialTransform);
  registration->SetMetric(metric);
  registration->SetOptimizer(optimizer);
  registration->SetFixedImage(fixedImage);
  registration->SetMovingImage(movingImage);

  auto                          movingInitialTransform = TransformType::New();
  TransformType::ParametersType initialParameters(movingInitialTransform->GetNumberOfParameters());
  initialParameters[0] = 0.0;
  initialParameters[1] = 0.0;
  movingInitialTransform->SetParameters(initialParameters);
  registration->SetMovingInitialTransform(movingInitialTransform);

  auto identityTransform = TransformType::New();
  identityTransform->SetIdentity();
  registration->SetFixedInitialTransform(identityTransform);

  constexpr unsigned int numberOfLevels = 1;
  registration->SetNumberOfLevels(numberOfLevels);

  RegistrationType::ShrinkFactorsArrayType shrinkFactorsPerLevel;
  shrinkFactorsPerLevel.SetSize(1);
  shrinkFactorsPerLevel[0] = 1;
  registration->SetShrinkFactorsPerLevel(shrinkFactorsPerLevel);

  RegistrationType::SmoothingSigmasArrayType smoothingSigmasPerLevel;
  smoothingSigmasPerLevel.SetSize(1);
  smoothingSigmasPerLevel[0] = 0;
  registration->SetSmoothingSigmasPerLevel(smoothingSigmasPerLevel);

  try
  {
    registration->Update();
    std::cout << "Optimizer stop condition: " << registration->GetOptimizer()->GetStopConditionDescription()
              << std::endl;
  }
  catch (const itk::ExceptionObject & err)
  {
    std::cerr << "ExceptionObject caught !" << std::endl;
    std::cerr << err << std::endl;
    return EXIT_FAILURE;
  }

  auto transform = registration->GetTransform();
  auto finalParameters = transform->GetParameters();
  auto translationAlongX = finalParameters[0];
  auto translationAlongY = finalParameters[1];

  auto numberOfIterations = optimizer->GetCurrentIteration();

  auto bestValue = optimizer->GetValue();

  std::cout << "Result = " << std::endl;
  std::cout << " Translation X = " << translationAlongX << std::endl;
  std::cout << " Translation Y = " << translationAlongY << std::endl;
  std::cout << " Iterations    = " << numberOfIterations << std::endl;
  std::cout << " Metric value  = " << bestValue << std::endl;

  using CompositeTransformType = itk::CompositeTransform<double, Dimension>;
  auto outputCompositeTransform = CompositeTransformType::New();
  outputCompositeTransform->AddTransform(movingInitialTransform);
  outputCompositeTransform->AddTransform(registration->GetModifiableTransform());

  using ResampleFilterType = itk::ResampleImageFilter<MovingImageType, FixedImageType>;
  auto resampler = ResampleFilterType::New();
  resampler->SetInput(movingImage);
  resampler->SetTransform(outputCompositeTransform);
  resampler->SetUseReferenceImage(true);
  resampler->SetReferenceImage(fixedImage);
  resampler->SetDefaultPixelValue(100);

  using OutputPixelType = unsigned char;
  using OutputImageType = itk::Image<OutputPixelType, Dimension>;

  using CastFilterType = itk::CastImageFilter<FixedImageType, OutputImageType>;
  auto caster = CastFilterType::New();
  caster->SetInput(resampler->GetOutput());

  itk::WriteImage(caster->GetOutput(), outputImageFile);

  using DifferenceFilterType = itk::SubtractImageFilter<FixedImageType, FixedImageType, FixedImageType>;
  auto difference = DifferenceFilterType::New();
  difference->SetInput1(fixedImage);
  difference->SetInput2(resampler->GetOutput());

  using RescalerType = itk::RescaleIntensityImageFilter<FixedImageType, OutputImageType>;
  auto intensityRescaler = RescalerType::New();
  intensityRescaler->SetInput(difference->GetOutput());
  intensityRescaler->SetOutputMinimum(itk::NumericTraits<OutputPixelType>::min());
  intensityRescaler->SetOutputMaximum(itk::NumericTraits<OutputPixelType>::max());

  resampler->SetDefaultPixelValue(1);
  itk::WriteImage(intensityRescaler->GetOutput(), differenceImageAfterFile);

  resampler->SetTransform(identityTransform);
  itk::WriteImage(intensityRescaler->GetOutput(), differenceImageBeforeFile);

  return EXIT_SUCCESS;
}

Classes demonstrated

template<typename TFixedImage, typename TMovingImage, typename TOutputTransform = Transform<double, TFixedImage::ImageDimension, TFixedImage::ImageDimension>, typename TVirtualImage = TFixedImage, typename TPointSet = PointSet<unsigned int, TFixedImage::ImageDimension>>
class ImageRegistrationMethodv4 : public itk::ProcessObject

Interface method for the current registration framework.

This interface method class encapsulates typical registration usage by incorporating all the necessary elements for performing a simple image registration between two images. This method also allows for multistage registration whereby each stage is characterize by possibly different transforms of and different image metrics. For example, many users will want to perform a linear registration followed by deformable registration where both stages are performed in multiple levels. Each level can be characterized by:

• the resolution of the virtual domain image (see below)

• smoothing of the fixed and moving images

• the coarseness of the current transform via transform adaptors (see below)

Multiple stages are handled by linking multiple instantiations of this class where the output transform is added to the optional composite transform input.

Transform adaptors: To accommodate new changes to the current ITK registration framework, we introduced the concept of transform adaptors. Whereas each stage is associated with a moving and, possibly, fixed transform, each level of each stage is defined by a transform adaptor which describes how to adapt the transform to the current level. For example, if one were to use the B-spline transform during a deformable registration stage, common practice is to increase the resolution of the B-spline mesh (or, analogously, the control point grid size) at each level. At each level, one would define the parameters of the B-spline transform adaptor at that level which increases the resolution from the previous level. For many transforms, such as affine, this concept of an adaptor may be nonsensical. For this reason, the base transform adaptor class does not do anything to the transform but merely passes it through. Each level of each stage must define a transform adaptor but, by default, the base adaptor class is assigned which, again, does not do anything to the transform. A special mention should be made of the transform adaptor at level 0 of any stage. Most likely, the user will not want to do anything to the transform as it enters into the given stage so typical use will be to assign the base adaptor class to level 0 of all stages but we leave that open to the user.

Output: The output is the updated transform.

Author

Nick Tustison

Author

Brian Avants

Subclassed by itk::SyNImageRegistrationMethod< TFixedImage, TMovingImage, TOutputTransform, TVirtualImage, TPointSet >, itk::TimeVaryingBSplineVelocityFieldImageRegistrationMethod< TFixedImage, TMovingImage, TOutputTransform, TVirtualImage, TPointSet >, itk::TimeVaryingVelocityFieldImageRegistrationMethodv4< TFixedImage, TMovingImage, TOutputTransform, TVirtualImage, TPointSet >

See itk::ImageRegistrationMethodv4 for additional documentation.