Perform Registration on Vector Images

Synopsis

Register images that have multi-component pixels like an itk::VectorImage or an itk::Image< itk::Vector, Dimension >.

Results

Fixed image.	Moving image.	Resampled, registered moving image.

Code

C++

#include "itkMeanSquaresImageToImageMetricv4.h"
#include "itkGradientDescentOptimizerv4.h"
#include "itkRegistrationParameterScalesFromPhysicalShift.h"
#include "itkVectorImageToImageMetricTraitsv4.h"

#include "itkGaussianSmoothingOnUpdateDisplacementFieldTransform.h"

#include "itkCastImageFilter.h"

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

#include <iomanip>

int
main(int argc, char * argv[])
{

  if (argc < 4)
  {
    std::cerr << "Missing Parameters " << std::endl;
    std::cerr << "Usage: " << argv[0];
    std::cerr << " fixedImageFile movingImageFile ";
    std::cerr << " resampledMovingImageFile ";
    std::cerr << " [numberOfAffineIterations numberOfDisplacementIterations] ";
    std::cerr << std::endl;
    return EXIT_FAILURE;
  }
  const char * fixedImageFile = argv[1];
  const char * movingImageFile = argv[2];
  const char * resampledMovingImageFile = argv[3];

  unsigned int numberOfAffineIterations = 2;
  unsigned int numberOfDisplacementIterations = 2;
  if (argc >= 5)
  {
    numberOfAffineIterations = std::stoi(argv[4]);
  }
  if (argc >= 6)
  {
    numberOfDisplacementIterations = std::stoi(argv[5]);
  }

  constexpr unsigned int Dimension = 2;

  // The input images have red, blue, and green pixel components.
  constexpr unsigned int NumberOfPixelComponents = 3;
  using PixelComponentType = float;
  using InputPixelType = itk::Vector<PixelComponentType, NumberOfPixelComponents>;
  using InputImageType = itk::Image<InputPixelType, Dimension>;

  using ParametersValueType = double;

  using ReaderType = itk::ImageFileReader<InputImageType>;
  ReaderType::Pointer fixedImageReader = ReaderType::New();
  fixedImageReader->SetFileName(fixedImageFile);
  ReaderType::Pointer movingImageReader = ReaderType::New();
  movingImageReader->SetFileName(movingImageFile);

  // Get the input images
  try
  {
    fixedImageReader->Update();
    movingImageReader->Update();
  }
  catch (itk::ExceptionObject & error)
  {
    std::cerr << "Error while reading images: " << error << std::endl;
    return EXIT_FAILURE;
  }
  InputImageType::Pointer fixedImage = fixedImageReader->GetOutput();
  InputImageType::Pointer movingImage = movingImageReader->GetOutput();

  using AffineTransformType = itk::AffineTransform<ParametersValueType, Dimension>;
  AffineTransformType::Pointer affineTransform = AffineTransformType::New();

  using DisplacementTransformType =
    itk::GaussianSmoothingOnUpdateDisplacementFieldTransform<ParametersValueType, Dimension>;
  DisplacementTransformType::Pointer displacementTransform = DisplacementTransformType::New();

  using DisplacementFieldType = DisplacementTransformType::DisplacementFieldType;
  DisplacementFieldType::Pointer displacementField = DisplacementFieldType::New();

  // Set the displacement field to be the same as the fixed image region, which will
  // act by default as the virtual domain in this example.
  displacementField->SetRegions(fixedImage->GetLargestPossibleRegion());
  // make sure the displacement field has the same spatial information as the image
  displacementField->CopyInformation(fixedImage);
  std::cout << "fixedImage->GetLargestPossibleRegion(): " << fixedImage->GetLargestPossibleRegion() << std::endl;
  displacementField->Allocate();
  // Fill it with 0's
  DisplacementTransformType::OutputVectorType zeroVector;
  zeroVector.Fill(0.0);
  displacementField->FillBuffer(zeroVector);
  // Assign to transform
  displacementTransform->SetDisplacementField(displacementField);
  displacementTransform->SetGaussianSmoothingVarianceForTheUpdateField(5.0);
  displacementTransform->SetGaussianSmoothingVarianceForTheTotalField(6.0);

  // Identity transform for fixed image
  using IdentityTransformType = itk::IdentityTransform<ParametersValueType, Dimension>;
  IdentityTransformType::Pointer identityTransform = IdentityTransformType::New();

  // The metric
  using VirtualImageType = itk::Image<double, Dimension>;
  using MetricTraitsType =
    itk::VectorImageToImageMetricTraitsv4<InputImageType, InputImageType, VirtualImageType, InputPixelType::Length>;
  using MetricType =
    itk::MeanSquaresImageToImageMetricv4<InputImageType, InputImageType, VirtualImageType, double, MetricTraitsType>;
  using PointSetType = MetricType::FixedSampledPointSetType;
  MetricType::Pointer metric = MetricType::New();

  using PointType = PointSetType::PointType;
  PointSetType::Pointer                             pointSet = PointSetType::New();
  itk::ImageRegionIteratorWithIndex<InputImageType> fixedImageIt(fixedImage, fixedImage->GetLargestPossibleRegion());
  itk::SizeValueType                                index = 0;
  itk::SizeValueType                                count = 0;
  for (fixedImageIt.GoToBegin(); !fixedImageIt.IsAtEnd(); ++fixedImageIt)
  {
    // take every N^th point
    if (count % 2 == 0)
    {
      PointType point;
      fixedImage->TransformIndexToPhysicalPoint(fixedImageIt.GetIndex(), point);
      pointSet->SetPoint(index, point);
      ++index;
    }
    ++count;
  }
  metric->SetFixedSampledPointSet(pointSet);
  metric->SetUseSampledPointSet(true);


  // Assign images and transforms.
  // By not setting a virtual domain image or virtual domain settings,
  // the metric will use the fixed image for the virtual domain.
  metric->SetFixedImage(fixedImage);
  metric->SetMovingImage(movingImage);
  metric->SetFixedTransform(identityTransform);
  metric->SetMovingTransform(affineTransform);
  const bool gaussian = false;
  metric->SetUseMovingImageGradientFilter(gaussian);
  metric->SetUseFixedImageGradientFilter(gaussian);
  metric->Initialize();

  using RegistrationParameterScalesFromShiftType = itk::RegistrationParameterScalesFromPhysicalShift<MetricType>;
  RegistrationParameterScalesFromShiftType::Pointer shiftScaleEstimator =
    RegistrationParameterScalesFromShiftType::New();
  shiftScaleEstimator->SetMetric(metric);


  //
  // Affine registration
  //
  std::cout << "First do an affine registration..." << std::endl;
  using OptimizerType = itk::GradientDescentOptimizerv4;
  OptimizerType::Pointer optimizer = OptimizerType::New();
  optimizer->SetMetric(metric);
  optimizer->SetNumberOfIterations(numberOfAffineIterations);
  optimizer->SetScalesEstimator(shiftScaleEstimator);
  optimizer->StartOptimization();
  std::cout << "After optimization affine parameters are: " << affineTransform->GetParameters() << std::endl;

  //
  // Deformable registration
  //
  std::cout << "Now, add the displacement field to the composite transform..." << std::endl;
  using CompositeType = itk::CompositeTransform<ParametersValueType, Dimension>;
  CompositeType::Pointer compositeTransform = CompositeType::New();
  compositeTransform->AddTransform(affineTransform);
  compositeTransform->AddTransform(displacementTransform);
  // Optimize only the displacement field, but still using the
  // previously computed affine transformation.
  compositeTransform->SetOnlyMostRecentTransformToOptimizeOn();
  metric->SetMovingTransform(compositeTransform);
  metric->SetUseSampledPointSet(false);
  metric->Initialize();

  // Optimizer
  optimizer->SetScalesEstimator(shiftScaleEstimator);
  optimizer->SetMetric(metric);
  optimizer->SetNumberOfIterations(numberOfDisplacementIterations);
  try
  {
    if (numberOfDisplacementIterations > 0)
    {
      optimizer->StartOptimization();
    }
    else
    {
      std::cout << "*** Skipping displacement field optimization.\n";
    }
  }
  catch (itk::ExceptionObject & error)
  {
    std::cerr << "Error: " << error << std::endl;
    return EXIT_FAILURE;
  }

  // Warp the image with the displacement field
  using ResampleFilterType = itk::ResampleImageFilter<InputImageType, InputImageType>;
  ResampleFilterType::Pointer resampler = ResampleFilterType::New();
  resampler->SetTransform(compositeTransform);
  resampler->SetInput(movingImageReader->GetOutput());
  resampler->SetReferenceImage(fixedImage);
  resampler->SetUseReferenceImage(true);

  // Write the warped image into a file
  using OutputPixelType = itk::Vector<unsigned char, NumberOfPixelComponents>;
  using OutputImageType = itk::Image<OutputPixelType, Dimension>;
  using CastFilterType = itk::CastImageFilter<InputImageType, OutputImageType>;
  CastFilterType::Pointer caster = CastFilterType::New();
  caster->SetInput(resampler->GetOutput());

  using WriterType = itk::ImageFileWriter<OutputImageType>;
  WriterType::Pointer writer = WriterType::New();
  writer->SetFileName(resampledMovingImageFile);
  writer->SetInput(caster->GetOutput());
  try
  {
    writer->UpdateLargestPossibleRegion();
  }
  catch (itk::ExceptionObject & error)
  {
    std::cerr << "Error: " << error << std::endl;
    return EXIT_FAILURE;
  }

  return EXIT_SUCCESS;
}

Classes demonstrated

template<typename TFixedImageType, typename TMovingImageType, typename TVirtualImageType, unsigned int NumberOfComponents, typename TCoordRep = typename ObjectToObjectMetricBase::CoordinateRepresentationType>
class VectorImageToImageMetricTraitsv4

A simple structure holding type information for ImageToImageMetricv4 classes.

This class provides type information for class members and methods used in gradient calculation. This class is used for images with vector pixel types, including VectorImage. For images with scalar pixel types, see itkDefaultImageToImageMetricTraitsv4.

See

itkDefaultImageToImageMetricTraitsv4

See itk::VectorImageToImageMetricTraitsv4 for additional documentation.