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

Input images

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

Output image

Output registered image.

Image difference before registration.

Image difference after registration.

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

Code

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>;


  using FixedImageReaderType = itk::ImageFileReader<FixedImageType>;
  auto fixedImageReader = FixedImageReaderType::New();
  fixedImageReader->SetFileName(fixedImageFile);

  using MovingImageReaderType = itk::ImageFileReader<MovingImageType>;
  auto movingImageReader = MovingImageReaderType::New();
  movingImageReader->SetFileName(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(fixedImageReader->GetOutput());
  registration->SetMovingImage(movingImageReader->GetOutput());

  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 (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(movingImageReader->GetOutput());
  resampler->SetTransform(outputCompositeTransform);
  auto fixedImage = fixedImageReader->GetOutput();
  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());

  using WriterType = itk::ImageFileWriter<OutputImageType>;
  auto writer = WriterType::New();
  writer->SetFileName(outputImageFile);
  writer->SetInput(caster->GetOutput());
  writer->Update();

  using DifferenceFilterType = itk::SubtractImageFilter<FixedImageType, FixedImageType, FixedImageType>;
  auto difference = DifferenceFilterType::New();
  difference->SetInput1(fixedImageReader->GetOutput());
  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);

  writer->SetInput(intensityRescaler->GetOutput());
  writer->SetFileName(differenceImageAfterFile);
  writer->Update();

  resampler->SetTransform(identityTransform);
  writer->SetFileName(differenceImageBeforeFile);
  writer->Update();

  return EXIT_SUCCESS;
}

Python

#!/usr/bin/env python
import itk

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)

if len(sys.argv) != 6:
    print("Usage: " + sys.argv[0] + " <fixedImageFile> <movingImageFile> "
          "<outputImagefile> <differenceImageAfter> <differenceImageBefore>")
    sys.exit(1)

fixedImageFile = sys.argv[1]
movingImageFile = sys.argv[2]
outputImageFile = sys.argv[3]
differenceImageAfterFile = sys.argv[4]
differenceImageBeforeFile = sys.argv[5]

PixelType = itk.ctype('float')

fixedImage = itk.imread(fixedImageFile, PixelType)
movingImage = itk.imread(movingImageFile, 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=outputImageFile)
writer.SetFileName(outputImageFile)
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(differenceImageAfterFile)
writer.Update()

resampler.SetTransform(identityTransform)
writer.SetFileName(differenceImageBeforeFile)
writer.Update()

Classes demonstrated