Download this example

Register Image to Another Using Landmarks

Synopsis

Rigidly register one image to another using manually specified landmarks.

Results

fixed.png

moving.png

output.png

Code

Python

#!/usr/bin/env python3

import itk

Dimension = 2
PixelType = itk.ctype("unsigned char")
ImageType = itk.Image[PixelType, Dimension]


def create_fixed_image():
    start = [
        0,
    ] * Dimension
    size = [
        100,
    ] * Dimension

    region = itk.ImageRegion[Dimension](start, size)

    ImageType = itk.Image[PixelType, Dimension]
    image = ImageType.New()
    image.SetRegions(region)
    image.Allocate()
    image.FillBuffer(0)

    image[11:20, 11:20] = 255

    itk.imwrite(image, "fixedPython.png")

    return image


def create_moving_image():
    start = [
        0,
    ] * Dimension
    size = [
        100,
    ] * Dimension

    region = itk.ImageRegion[Dimension](start, size)

    ImageType = itk.Image[PixelType, Dimension]
    image = ImageType.New()
    image.SetRegions(region)
    image.Allocate()
    image.FillBuffer(0)

    image[51:60, 51:60] = 100

    itk.imwrite(image, "movingPython.png")

    return image


fixed_image = create_fixed_image()
moving_image = create_moving_image()

LandmarkPointType = itk.Point[itk.D, Dimension]
LandmarkContainerType = itk.vector[LandmarkPointType]

fixed_landmarks = LandmarkContainerType()
moving_landmarks = LandmarkContainerType()

fixed_point = LandmarkPointType()
moving_point = LandmarkPointType()

fixed_point[0] = 10
fixed_point[1] = 10
moving_point[0] = 50
moving_point[1] = 50
fixed_landmarks.push_back(fixed_point)
moving_landmarks.push_back(moving_point)

fixed_point[0] = 20
fixed_point[1] = 10
moving_point[0] = 60
moving_point[1] = 50
fixed_landmarks.push_back(fixed_point)
moving_landmarks.push_back(moving_point)

fixed_point[0] = 20
fixed_point[1] = 20
moving_point[0] = 60
moving_point[1] = 60
fixed_landmarks.push_back(fixed_point)
moving_landmarks.push_back(moving_point)

TransformInitializerType = itk.LandmarkBasedTransformInitializer[
    itk.Transform[itk.D, Dimension, Dimension]
]
transform_initializer = TransformInitializerType.New()

transform_initializer.SetFixedLandmarks(fixed_landmarks)
transform_initializer.SetMovingLandmarks(moving_landmarks)

transform = itk.Rigid2DTransform[itk.D].New()
transform_initializer.SetTransform(transform)
transform_initializer.InitializeTransform()

output = itk.resample_image_filter(
    moving_image,
    transform=transform,
    use_reference_image=True,
    reference_image=fixed_image,
    default_pixel_value=200,
)

itk.imwrite(output, "outputPython.png")

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)[])
{
  auto fixedImage = ImageType::New();
  CreateFixedImage(fixedImage);

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

  using TransformType = itk::Rigid2DTransform<double>;
  using LandmarkBasedTransformInitializerType =
    itk::LandmarkBasedTransformInitializer<TransformType, 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);

  auto transform = TransformType::New();

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

  using ResampleFilterType = itk::ResampleImageFilter<ImageType, ImageType, double>;
  auto resampleFilter = ResampleFilterType::New();
  resampleFilter->SetInput(movingImage);
  resampleFilter->SetTransform(transform);
  resampleFilter->SetUseReferenceImage(true);
  resampleFilter->SetReferenceImage(fixedImage);
  resampleFilter->SetDefaultPixelValue(200);

  itk::WriteImage(resampleFilter->GetOutput(), "output.png");

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

  itk::WriteImage(image, "fixed.png");
}


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

  itk::WriteImage(image, "moving.png");
}

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

• Register Image To Another Using Landmarks

See itk::LandmarkBasedTransformInitializer for additional documentation.