Register Image to Another Using Landmarks¶

Synopsis¶

Rigidly register one image to another using manually specified landmarks.

Results¶

Code¶

C++¶

#include "itkImageFileWriter.h"
#include "itkImage.h"
#include "itkVector.h"
#include "itkResampleImageFilter.h"
#include "itkLandmarkBasedTransformInitializer.h"
#include "itkRigid2DTransform.h"

constexpr unsigned int Dimension = 2;
using PixelType = unsigned char;
using ImageType = itk::Image<PixelType, Dimension>;

static void
CreateFixedImage(ImageType::Pointer image);
static void
CreateMovingImage(ImageType::Pointer image);

int
main(int itkNotUsed(argc), char * itkNotUsed(argv)[])
{
  ImageType::Pointer fixedImage = ImageType::New();
  CreateFixedImage(fixedImage);

  ImageType::Pointer movingImage = ImageType::New();
  CreateMovingImage(movingImage);

  using Rigid2DTransformType = itk::Rigid2DTransform<double>;
  using LandmarkBasedTransformInitializerType =
    itk::LandmarkBasedTransformInitializer<Rigid2DTransformType, ImageType, ImageType>;

  LandmarkBasedTransformInitializerType::Pointer landmarkBasedTransformInitializer =
    LandmarkBasedTransformInitializerType::New();
  //  Create source and target landmarks.
  using LandmarkContainerType = LandmarkBasedTransformInitializerType::LandmarkPointContainer;
  using LandmarkPointType = LandmarkBasedTransformInitializerType::LandmarkPointType;

  LandmarkContainerType fixedLandmarks;
  LandmarkContainerType movingLandmarks;

  LandmarkPointType fixedPoint;
  LandmarkPointType movingPoint;

  fixedPoint[0] = 10;
  fixedPoint[1] = 10;
  movingPoint[0] = 50;
  movingPoint[1] = 50;
  fixedLandmarks.push_back(fixedPoint);
  movingLandmarks.push_back(movingPoint);

  fixedPoint[0] = 10;
  fixedPoint[1] = 20;
  movingPoint[0] = 50;
  movingPoint[1] = 60;
  fixedLandmarks.push_back(fixedPoint);
  movingLandmarks.push_back(movingPoint);

  fixedPoint[0] = 20;
  fixedPoint[1] = 10;
  movingPoint[0] = 60;
  movingPoint[1] = 50;
  fixedLandmarks.push_back(fixedPoint);
  movingLandmarks.push_back(movingPoint);

  fixedPoint[0] = 20;
  fixedPoint[1] = 20;
  movingPoint[0] = 60;
  movingPoint[1] = 60;
  fixedLandmarks.push_back(fixedPoint);
  movingLandmarks.push_back(movingPoint);

  landmarkBasedTransformInitializer->SetFixedLandmarks(fixedLandmarks);
  landmarkBasedTransformInitializer->SetMovingLandmarks(movingLandmarks);

  Rigid2DTransformType::Pointer transform = Rigid2DTransformType::New();

  transform->SetIdentity();
  landmarkBasedTransformInitializer->SetTransform(transform);
  landmarkBasedTransformInitializer->InitializeTransform();

  using ResampleFilterType = itk::ResampleImageFilter<ImageType, ImageType, double>;
  ResampleFilterType::Pointer resampleFilter = ResampleFilterType::New();
  resampleFilter->SetInput(movingImage);
  resampleFilter->SetTransform(transform);
  resampleFilter->SetSize(fixedImage->GetLargestPossibleRegion().GetSize());
  resampleFilter->SetOutputOrigin(fixedImage->GetOrigin());
  resampleFilter->SetOutputSpacing(fixedImage->GetSpacing());
  resampleFilter->SetOutputDirection(fixedImage->GetDirection());
  resampleFilter->SetDefaultPixelValue(200);
  resampleFilter->GetOutput();

  // Write the output
  using WriterType = itk::ImageFileWriter<ImageType>;
  WriterType::Pointer writer = WriterType::New();
  writer->SetInput(resampleFilter->GetOutput());
  writer->SetFileName("output.png");
  writer->Update();

  return EXIT_SUCCESS;
}

void
CreateFixedImage(ImageType::Pointer image)
{
  // Create a black image with a white square
  ImageType::IndexType start;
  start.Fill(0);

  ImageType::SizeType size;
  size.Fill(100);

  ImageType::RegionType region;
  region.SetSize(size);
  region.SetIndex(start);

  image->SetRegions(region);
  image->Allocate();
  image->FillBuffer(0);

  itk::ImageRegionIterator<ImageType> imageIterator(image, region);

  while (!imageIterator.IsAtEnd())
  {
    if (imageIterator.GetIndex()[0] > 10 && imageIterator.GetIndex()[0] < 20 && imageIterator.GetIndex()[1] > 10 &&
        imageIterator.GetIndex()[1] < 20)
    {
      imageIterator.Set(255);
    }
    ++imageIterator;
  }

  // Write the deformation field
  using WriterType = itk::ImageFileWriter<ImageType>;
  WriterType::Pointer writer = WriterType::New();
  writer->SetInput(image);
  writer->SetFileName("fixed.png");
  writer->Update();
}


void
CreateMovingImage(ImageType::Pointer image)
{
  // Create a black image with a white square
  ImageType::IndexType start;
  start.Fill(0);

  ImageType::SizeType size;
  size.Fill(100);

  ImageType::RegionType region;
  region.SetSize(size);
  region.SetIndex(start);

  image->SetRegions(region);
  image->Allocate();
  image->FillBuffer(0);

  itk::ImageRegionIterator<ImageType> imageIterator(image, region);

  while (!imageIterator.IsAtEnd())
  {
    if (imageIterator.GetIndex()[0] > 50 && imageIterator.GetIndex()[0] < 60 && imageIterator.GetIndex()[1] > 50 &&
        imageIterator.GetIndex()[1] < 60)
    {
      imageIterator.Set(100);
    }
    ++imageIterator;
  }

  // Write the deformation field
  using WriterType = itk::ImageFileWriter<ImageType>;
  WriterType::Pointer writer = WriterType::New();
  writer->SetInput(image);
  writer->SetFileName("moving.png");
  writer->Update();
}

Classes demonstrated¶

template<typename TTransform, typename TFixedImage = itk::ImageBase<TTransform::InputSpaceDimension>, typename TMovingImage = itk::ImageBase<TTransform::OutputSpaceDimension>> class LandmarkBasedTransformInitializer : public itk::Object

This class computes the transform that aligns the fixed and moving images given a set of pair landmarks. The class is templated over the Transform type as well as fixed image and moving image types. The transform computed gives the best fit transform that maps the fixed and moving images in a least squares sense. The indices are taken to correspond, so point 1 in the first set will get mapped close to point 1 in the second set, etc.

Currently, the following transforms are supported by the class: VersorRigid3DTransform Rigid2DTransform AffineTransform BSplineTransform

An equal number of fixed and moving landmarks need to be specified using SetFixedLandmarks() and SetMovingLandmarks(). Any number of landmarks may be specified. In the case of the Affine transformation the number of landmarks must be greater than the landmark dimensionality. If this is not the case an exception is thrown. In the case of the VersorRigid3DTransform and Rigid2DTransform the number of landmarks must be equal or greater than the landmark dimensionality. If this is not the case, only the translational component of the transformation is computed and the rotation is the identity. In the case of using Affine or BSpline transforms, each landmark pair can contribute in the final transform based on its defined weight. Number of weights should be equal to the number of landmarks and can be specified using SetLandmarkWeight(). By defaults are weights are set to one. Call InitializeTransform() to initialize the transform.

The class is based in part on Hybrid/vtkLandmarkTransform originally implemented in python by David G. Gobbi.

The solution is based on Berthold K. P. Horn (1987), “Closed-form solution of absolute orientation

using unit quaternions,”

http://people.csail.mit.edu/bkph/papers/Absolute_Orientation.pdf

The Affine Transform

initializer is based on an algorithm by H Spaeth, and is described in the Insight Journal Article “Affine Transformation for Landmark Based Registration Initializer

in ITK” by Kim E.Y., Johnson H., Williams N. available at

http://midasjournal.com/browse/publication/825

ITK Sphinx Examples:

All ITK Sphinx Examples

See itk::LandmarkBasedTransformInitializer for additional documentation.