Global Registration Of Two Images (BSpline)

Synopsis

A global registration of two images.

Results

Output:

Intial Parameters =
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]

Starting Registration
*************************************************
N=50   NUMBER OF CORRECTIONS=5       INITIAL VALUES F= 5410.08   GNORM= 98.7811
*************************************************
I   NFN    FUNC        GNORM       STEPLENGTH
1    2      5173.917      47.318       0.010
2    8      1312.526      35.022      67.448
3   10      1241.036      21.952       0.095
4   11      1214.603      20.777       0.500
5   12      1169.045      47.965       0.500
6   13      1066.160      73.507       0.500
7   14       750.612      84.979       0.500
8   16       427.232      44.597       0.235
9   17       371.768      44.133       0.500
10   19       320.897      30.003       0.054
11   20       319.018      23.934       0.500
12   21       276.782      16.253       0.500
13   22       190.496      16.864       0.500
14   23       121.434      17.719       0.500
15   24        80.920      23.984       0.500
16   25        50.736      11.243       0.500
17   26        40.862       5.633       0.500
THE MINIMIZATION TERMINATED WITHOUT DETECTING ERRORS.
Optimizer stop condition = LBFGSOptimizer: Failure
Last Transform Parameters
[0.000004242652759130258, -0.09801257315348906, -0.3773136074775823, 0.04057356343246512, 0.011460974084144235, 0.046750769129838214, 5.442980649652609, -1.3233672220673571, 10.138325385009857, 1.2593858799352904, 0.1744432989808139, 27.255697794222773, 18.42787109104701, 27.82714519423285, 3.443535460389487, 0.04083243642873029, 5.114318164332667, 0.5697948754499333, 5.4282475197978695, 0.7721289215254623, -0.00001989578361620887, -0.034703497612920686, -0.08550784338556144, 0.006239176991931418, 0.001367593141535802, -4.808966670313369e-7, -0.2187505995413243, -1.8087191436710883, -0.9922884310554664, -0.014742349905146031, 0.003053794283557083, -11.874136737768657, -75.22296974556679, -33.96218955913907, -0.42061613449053764, 0.011561263304423593, -0.7922226379398286, -7.185598670736005, 0.7735213860631797, 0.1316131069607225, 0.0033284552078211094, 10.085406006964824, 56.37568206497135, 25.246189532623205, 0.3387536225030382, 0.000015865649566500348, 0.11028526291334259, 0.6117750781918082, 0.26063979972963747, 0.0029171257653189493]
../../../../_images/GlobalRegistrationTwoImagesBSpline.png

Output Image

Code

C++

#include "itkImageRegistrationMethod.h"
#include "itkMeanSquaresImageToImageMetric.h"
#include "itkTimeProbesCollectorBase.h"
#include "itkSpatialObjectToImageFilter.h"
#include "itkEllipseSpatialObject.h"

#if ITK_VERSION_MAJOR < 4
#  include "itkBSplineDeformableTransform.h"
#else
#  include "itkBSplineTransform.h"
#endif
#include "itkLBFGSOptimizer.h"
#include "itkImageFileWriter.h"
#include "itkResampleImageFilter.h"
#include "itkCastImageFilter.h"
#include "itkSquaredDifferenceImageFilter.h"

#ifdef ENABLE_QUICKVIEW
#  include "QuickView.h"
#endif

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

using ImageType = itk::Image<PixelType, ImageDimension>;

static void
CreateEllipseImage(ImageType::Pointer image);
static void
CreateCircleImage(ImageType::Pointer image);

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

  const unsigned int     SpaceDimension = ImageDimension;
  constexpr unsigned int SplineOrder = 3;
  using CoordinateRepType = double;

#if ITK_VERSION_MAJOR < 4
  using TransformType = itk::BSplineDeformableTransform<CoordinateRepType, SpaceDimension, SplineOrder>;
#else
  using TransformType = itk::BSplineTransform<CoordinateRepType, SpaceDimension, SplineOrder>;
#endif
  using OptimizerType = itk::LBFGSOptimizer;


  using MetricType = itk::MeanSquaresImageToImageMetric<ImageType, ImageType>;

  using InterpolatorType = itk::LinearInterpolateImageFunction<ImageType, double>;

  using RegistrationType = itk::ImageRegistrationMethod<ImageType, ImageType>;

  MetricType::Pointer       metric = MetricType::New();
  OptimizerType::Pointer    optimizer = OptimizerType::New();
  InterpolatorType::Pointer interpolator = InterpolatorType::New();
  RegistrationType::Pointer registration = RegistrationType::New();


  // The old registration framework has problems with multi-threading
  // For now, we set the number of work units to 1
  registration->SetNumberOfWorkUnits(1);

  registration->SetMetric(metric);
  registration->SetOptimizer(optimizer);
  registration->SetInterpolator(interpolator);

  TransformType::Pointer transform = TransformType::New();
  registration->SetTransform(transform);

  // Create the synthetic images
  ImageType::Pointer fixedImage = ImageType::New();
  CreateCircleImage(fixedImage);

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

  // Setup the registration
  registration->SetFixedImage(fixedImage);
  registration->SetMovingImage(movingImage);

  ImageType::RegionType fixedRegion = fixedImage->GetBufferedRegion();
  registration->SetFixedImageRegion(fixedRegion);

  //  Here we define the parameters of the BSplineDeformableTransform grid.  We
  //  arbitrarily decide to use a grid with $5 \times 5$ nodes within the image.
  //  The reader should note that the BSpline computation requires a
  //  finite support region ( 1 grid node at the lower borders and 2
  //  grid nodes at upper borders). Therefore in this example, we set
  //  the grid size to be $8 \times 8$ and place the grid origin such that
  //  grid node (1,1) coincides with the first pixel in the fixed image.

#if ITK_VERSION_MAJOR < 4
  using RegionType = TransformType::RegionType;
  RegionType           bsplineRegion;
  RegionType::SizeType gridSizeOnImage;
  RegionType::SizeType gridBorderSize;
  RegionType::SizeType totalGridSize;

  gridSizeOnImage.Fill(5);
  gridBorderSize.Fill(3); // Border for spline order = 3 ( 1 lower, 2 upper )
  totalGridSize = gridSizeOnImage + gridBorderSize;

  bsplineRegion.SetSize(totalGridSize);

  using SpacingType = TransformType::SpacingType;
  SpacingType spacing = fixedImage->GetSpacing();

  using OriginType = TransformType::OriginType;
  OriginType origin = fixedImage->GetOrigin();

  ImageType::SizeType fixedImageSize = fixedRegion.GetSize();

  for (unsigned int r = 0; r < ImageDimension; r++)
  {
    spacing[r] *= static_cast<double>(fixedImageSize[r] - 1) / static_cast<double>(gridSizeOnImage[r] - 1);
  }

  ImageType::DirectionType gridDirection = fixedImage->GetDirection();
  SpacingType              gridOriginOffset = gridDirection * spacing;

  OriginType gridOrigin = origin - gridOriginOffset;

  transform->SetGridSpacing(spacing);
  transform->SetGridOrigin(gridOrigin);
  transform->SetGridRegion(bsplineRegion);
  transform->SetGridDirection(gridDirection);
#else
  TransformType::PhysicalDimensionsType fixedPhysicalDimensions;
  TransformType::MeshSizeType           meshSize;
  for (unsigned int i = 0; i < ImageDimension; i++)
  {
    fixedPhysicalDimensions[i] =
      fixedImage->GetSpacing()[i] * static_cast<double>(fixedImage->GetLargestPossibleRegion().GetSize()[i] - 1);
  }
  unsigned int numberOfGridNodesInOneDimension = 5;
  meshSize.Fill(numberOfGridNodesInOneDimension - SplineOrder);
  transform->SetTransformDomainOrigin(fixedImage->GetOrigin());
  transform->SetTransformDomainPhysicalDimensions(fixedPhysicalDimensions);
  transform->SetTransformDomainMeshSize(meshSize);
  transform->SetTransformDomainDirection(fixedImage->GetDirection());
#endif

  using ParametersType = TransformType::ParametersType;

  const unsigned int numberOfParameters = transform->GetNumberOfParameters();

  ParametersType parameters(numberOfParameters);

  parameters.Fill(0.0);

  transform->SetParameters(parameters);

  //  We now pass the parameters of the current transform as the initial
  //  parameters to be used when the registration process starts.

  registration->SetInitialTransformParameters(transform->GetParameters());

  std::cout << "Intial Parameters = " << std::endl;
  std::cout << transform->GetParameters() << std::endl;

  //  Next we set the parameters of the LBFGS Optimizer.

  optimizer->SetGradientConvergenceTolerance(0.05);
  optimizer->SetLineSearchAccuracy(0.9);
  optimizer->SetDefaultStepLength(.5);
  optimizer->TraceOn();
  optimizer->SetMaximumNumberOfFunctionEvaluations(1000);

  std::cout << std::endl << "Starting Registration" << std::endl;

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

  OptimizerType::ParametersType finalParameters = registration->GetLastTransformParameters();

  std::cout << "Last Transform Parameters" << std::endl;
  std::cout << finalParameters << std::endl;

  transform->SetParameters(finalParameters);

  using ResampleFilterType = itk::ResampleImageFilter<ImageType, ImageType>;

  ResampleFilterType::Pointer resample = ResampleFilterType::New();

  resample->SetTransform(transform);
  resample->SetInput(movingImage);

  resample->SetSize(fixedImage->GetLargestPossibleRegion().GetSize());
  resample->SetOutputOrigin(fixedImage->GetOrigin());
  resample->SetOutputSpacing(fixedImage->GetSpacing());
  resample->SetOutputDirection(fixedImage->GetDirection());
  resample->SetDefaultPixelValue(100);

  using OutputPixelType = unsigned char;

  using OutputImageType = itk::Image<OutputPixelType, ImageDimension>;

  using CastFilterType = itk::CastImageFilter<ImageType, OutputImageType>;

  using OutputWriterType = itk::ImageFileWriter<OutputImageType>;

  OutputWriterType::Pointer writer = OutputWriterType::New();
  CastFilterType::Pointer   caster = CastFilterType::New();


  writer->SetFileName("output.png");

  caster->SetInput(resample->GetOutput());
  writer->SetInput(caster->GetOutput());

  try
  {
    writer->Update();
  }
  catch (itk::ExceptionObject & err)
  {
    std::cerr << "ExceptionObject caught !" << std::endl;
    std::cerr << err << std::endl;
    return EXIT_FAILURE;
  }

#ifdef ENABLE_QUICKVIEW
  QuickView viewer;
  viewer.AddImage(fixedImage.GetPointer(), true, "Fixed Image");
  viewer.AddImage(movingImage.GetPointer(), true, "Moving Image");
  viewer.AddImage(resample->GetOutput(), true, "Resampled Moving Image");

  viewer.Visualize();
#endif

  return EXIT_SUCCESS;
}

void
CreateEllipseImage(ImageType::Pointer image)
{
  using EllipseType = itk::EllipseSpatialObject<ImageDimension>;

  using SpatialObjectToImageFilterType = itk::SpatialObjectToImageFilter<EllipseType, ImageType>;

  SpatialObjectToImageFilterType::Pointer imageFilter = SpatialObjectToImageFilterType::New();

  ImageType::SizeType size;
  size[0] = 100;
  size[1] = 100;

  imageFilter->SetSize(size);

  ImageType::SpacingType spacing;
  spacing.Fill(1);
  imageFilter->SetSpacing(spacing);

  EllipseType::Pointer   ellipse = EllipseType::New();
  EllipseType::ArrayType radiusArray;
  radiusArray[0] = 10;
  radiusArray[1] = 20;
  ellipse->SetRadiusInObjectSpace(radiusArray);

  using TransformType = EllipseType::TransformType;
  TransformType::Pointer transform = TransformType::New();
  transform->SetIdentity();

  TransformType::OutputVectorType translation;
  translation[0] = 65;
  translation[1] = 45;
  transform->Translate(translation, false);

  ellipse->SetObjectToParentTransform(transform);

  imageFilter->SetInput(ellipse);

  ellipse->SetDefaultInsideValue(255);
  ellipse->SetDefaultOutsideValue(0);
  imageFilter->SetUseObjectValue(true);
  imageFilter->SetOutsideValue(0);

  imageFilter->Update();

  image->Graft(imageFilter->GetOutput());
}

void
CreateCircleImage(ImageType::Pointer image)
{
  using EllipseType = itk::EllipseSpatialObject<ImageDimension>;

  using SpatialObjectToImageFilterType = itk::SpatialObjectToImageFilter<EllipseType, ImageType>;

  SpatialObjectToImageFilterType::Pointer imageFilter = SpatialObjectToImageFilterType::New();

  ImageType::SizeType size;
  size[0] = 100;
  size[1] = 100;

  imageFilter->SetSize(size);

  ImageType::SpacingType spacing;
  spacing.Fill(1);
  imageFilter->SetSpacing(spacing);

  EllipseType::Pointer   ellipse = EllipseType::New();
  EllipseType::ArrayType radiusArray;
  radiusArray[0] = 10;
  radiusArray[1] = 10;
  ellipse->SetRadiusInObjectSpace(radiusArray);

  using TransformType = EllipseType::TransformType;
  TransformType::Pointer transform = TransformType::New();
  transform->SetIdentity();

  TransformType::OutputVectorType translation;
  translation[0] = 50;
  translation[1] = 50;
  transform->Translate(translation, false);

  ellipse->SetObjectToParentTransform(transform);

  imageFilter->SetInput(ellipse);

  ellipse->SetDefaultInsideValue(255);
  ellipse->SetDefaultOutsideValue(0);
  imageFilter->SetUseObjectValue(true);
  imageFilter->SetOutsideValue(0);

  imageFilter->Update();

  image->Graft(imageFilter->GetOutput());
}

Classes demonstrated

template<typename TParametersValueType = double, unsigned int NDimensions = 3, unsigned int VSplineOrder = 3>
class BSplineDeformableTransform : public itk::BSplineBaseTransform<TParametersValueType, NDimensions, VSplineOrder>

Deformable transform using a BSpline representation.

This class encapsulates a deformable transform of points from one N-dimensional space to another N-dimensional space. The deformation field is modelled using B-splines. A deformation is defined on a sparse regular grid of control points

\vec{\lambda}_j and is varied by defining a deformation \vec{g}(\vec{\lambda}_j) of each control point. The deformation D(\vec{x}) at any point \vec{x} is obtained by using a B-spline interpolation kernel.
Note

BSplineTransform is a newer version of this class, and it is preferred.

The deformation field grid is defined by a user specified GridRegion, GridSpacing and GridOrigin. Each grid/control point has associated with it N deformation coefficients \vec{\delta}_j, representing the N directional components of the deformation. Deformation outside the grid plus support region for the BSpline interpolation is assumed to be zero.

Additionally, the user can specified an addition bulk transform B such that the transformed point is given by:

\vec{y} = B(\vec{x}) + D(\vec{x})

The parameters for this transform is an N x N-D grid of spline coefficients. The user specifies the parameters as one flat array: each N-D grid is represented by an array in the same way an N-D image is represented in the buffer; the N arrays are then concatenated together on form a single array.

For efficiency, this transform does not make a copy of the parameters. It only keeps a pointer to the input parameters and assumes that the memory is managed by the caller.

The following illustrates the typical usage of this class:

using TransformType = BSplineDeformableTransform<double,2,3>;
TransformType::Pointer transform = TransformType::New();

transform->SetGridRegion( region );
transform->SetGridSpacing( spacing );
transform->SetGridOrigin( origin );

// NB: the region must be set first before setting the parameters

TransformType::ParametersType parameters( transform->GetNumberOfParameters() );

// Fill the parameters with values

transform->SetParameters( parameters )

outputPoint = transform->TransformPoint( inputPoint );

An alternative way to set the B-spline coefficients is via array of images. The grid region, spacing and origin information is taken directly from the first image. It is assumed that the subsequent images are the same buffered region. The following illustrates the API:

TransformType::ImageConstPointer images[2];

// Fill the images up with values

transform->SetCoefficientImages( images );
outputPoint = transform->TransformPoint( inputPoint );

Warning: use either the SetParameters() or SetCoefficientImages() API. Mixing the two modes may results in unexpected results.

The class is templated coordinate representation type (float or double), the space dimension and the spline order.

See

BSplineTransform

ITK Sphinx Examples:

See itk::BSplineDeformableTransform for additional documentation.
template<typename TFixedImage, typename TMovingImage>
class ImageRegistrationMethod : public itk::ProcessObject

Base class for Image Registration Methods.

This Class define the generic interface for a registration method.

This class is templated over the type of the two image to be registered. A generic Transform is used by this class. That allows to select at run time the particular type of transformation that is to be applied for registering the images.

This method use a generic Metric in order to compare the two images. the final goal of the registration method is to find the set of parameters of the Transformation that optimizes the metric.

The registration method also support a generic optimizer that can be selected at run-time. The only restriction for the optimizer is that it should be able to operate in single-valued cost functions given that the metrics used to compare images provide a single value as output.

The terms : Fixed image and Moving image are used in this class to indicate what image is being mapped by the transform.

This class uses the coordinate system of the Fixed image as a reference and searches for a Transform that will map points from the space of the Fixed image to the space of the Moving image.

For doing so, a Metric will be continuously applied to compare the Fixed image with the Transformed Moving image. This process also requires to interpolate values from the Moving image.

ITK Sphinx Examples:

See itk::ImageRegistrationMethod for additional documentation.