Global Registration Two Images (Affine)#
Synopsis#
A basic global registration of two images.
Results#
fixed.png#
moving.png#
output.png#
Output:
Final parameters: [1.1921533496320087, 0.10372902569023787, -0.18132002016416124, 1.1464158834351523, 0.0021451859042650244, -0.003039195975788232]
Result =
Metric value = 1836.41
Code#
C++#
#include "itkCastImageFilter.h"
#include "itkEllipseSpatialObject.h"
#include "itkImage.h"
#include "itkImageRegistrationMethod.h"
#include "itkLinearInterpolateImageFunction.h"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkMeanSquaresImageToImageMetric.h"
#include "itkRegularStepGradientDescentOptimizer.h"
#include "itkResampleImageFilter.h"
#include "itkRescaleIntensityImageFilter.h"
#include "itkSpatialObjectToImageFilter.h"
#include "itkAffineTransform.h"
constexpr unsigned int Dimension = 2;
using PixelType = unsigned char;
using ImageType = itk::Image<PixelType, Dimension>;
static void
CreateEllipseImage(ImageType::Pointer image);
static void
CreateSphereImage(ImageType::Pointer image);
int
main()
{
// The transform that will map the fixed image into the moving image.
using TransformType = itk::AffineTransform<double, Dimension>;
// An optimizer is required to explore the parameter space of the transform
// in search of optimal values of the metric.
using OptimizerType = itk::RegularStepGradientDescentOptimizer;
// The metric will compare how well the two images match each other. Metric
// types are usually parameterized by the image types as it can be seen in
// the following type declaration.
using MetricType = itk::MeanSquaresImageToImageMetric<ImageType, ImageType>;
// Finally, the type of the interpolator is declared. The interpolator will
// evaluate the intensities of the moving image at non-grid positions.
using InterpolatorType = itk::LinearInterpolateImageFunction<ImageType, double>;
// The registration method type is instantiated using the types of the
// fixed and moving images. This class is responsible for interconnecting
// all the components that we have described so far.
using RegistrationType = itk::ImageRegistrationMethod<ImageType, ImageType>;
// Create components
auto metric = MetricType::New();
auto transform = TransformType::New();
auto optimizer = OptimizerType::New();
auto interpolator = InterpolatorType::New();
auto registration = RegistrationType::New();
// Each component is now connected to the instance of the registration method.
registration->SetMetric(metric);
registration->SetOptimizer(optimizer);
registration->SetTransform(transform);
registration->SetInterpolator(interpolator);
// Get the two images
auto fixedImage = ImageType::New();
auto movingImage = ImageType::New();
CreateSphereImage(fixedImage);
CreateEllipseImage(movingImage);
// Write the two synthetic inputs
itk::WriteImage(fixedImage, "fixed.png");
itk::WriteImage(movingImage, "moving.png");
// Set the registration inputs
registration->SetFixedImage(fixedImage);
registration->SetMovingImage(movingImage);
registration->SetFixedImageRegion(fixedImage->GetLargestPossibleRegion());
// Initialize the transform
using ParametersType = RegistrationType::ParametersType;
ParametersType initialParameters(transform->GetNumberOfParameters());
// rotation matrix
initialParameters[0] = 1.0; // R(0,0)
initialParameters[1] = 0.0; // R(0,1)
initialParameters[2] = 0.0; // R(1,0)
initialParameters[3] = 1.0; // R(1,1)
// translation vector
initialParameters[4] = 0.0;
initialParameters[5] = 0.0;
registration->SetInitialTransformParameters(initialParameters);
optimizer->SetMaximumStepLength(.1); // If this is set too high, you will get a
//"itk::ERROR: MeanSquaresImageToImageMetric(0xa27ce70): Too many samples map outside moving image buffer: 1818 /
// 10000" error
optimizer->SetMinimumStepLength(0.01);
// Set a stopping criterion
optimizer->SetNumberOfIterations(200);
// Connect an observer
// auto observer = CommandIterationUpdate::New();
// optimizer->AddObserver( itk::IterationEvent(), observer );
try
{
registration->Update();
}
catch (const itk::ExceptionObject & err)
{
std::cerr << "ExceptionObject caught !" << std::endl;
std::cerr << err << std::endl;
return EXIT_FAILURE;
}
// The result of the registration process is an array of parameters that
// defines the spatial transformation in an unique way. This final result is
// obtained using the \code{GetLastTransformParameters()} method.
ParametersType finalParameters = registration->GetLastTransformParameters();
std::cout << "Final parameters: " << finalParameters << std::endl;
// The value of the image metric corresponding to the last set of parameters
// can be obtained with the \code{GetValue()} method of the optimizer.
const double bestValue = optimizer->GetValue();
// Print out results
//
std::cout << "Result = " << std::endl;
std::cout << " Metric value = " << bestValue << std::endl;
// It is common, as the last step of a registration task, to use the
// resulting transform to map the moving image into the fixed image space.
// This is easily done with the \doxygen{ResampleImageFilter}.
using ResampleFilterType = itk::ResampleImageFilter<ImageType, ImageType>;
auto resampler = ResampleFilterType::New();
resampler->SetInput(movingImage);
// The Transform that is produced as output of the Registration method is
// also passed as input to the resampling filter. Note the use of the
// methods \code{GetOutput()} and \code{Get()}. This combination is needed
// here because the registration method acts as a filter whose output is a
// transform decorated in the form of a \doxygen{DataObject}. For details in
// this construction you may want to read the documentation of the
// \doxygen{DataObjectDecorator}.
resampler->SetTransform(registration->GetOutput()->Get());
// As described in Section \ref{sec:ResampleImageFilter}, the
// ResampleImageFilter requires additional parameters to be specified, in
// particular, the spacing, origin and size of the output image. The default
// pixel value is also set to a distinct gray level in order to highlight
// the regions that are mapped outside of the moving image.
resampler->SetSize(fixedImage->GetLargestPossibleRegion().GetSize());
resampler->SetOutputOrigin(fixedImage->GetOrigin());
resampler->SetOutputSpacing(fixedImage->GetSpacing());
resampler->SetOutputDirection(fixedImage->GetDirection());
resampler->SetDefaultPixelValue(100);
// The output of the filter is passed to a writer that will store the
// image in a file. An \doxygen{CastImageFilter} is used to convert the
// pixel type of the resampled image to the final type used by the
// writer. The cast and writer filters are instantiated below.
using CastFilterType = itk::CastImageFilter<ImageType, ImageType>;
auto caster = CastFilterType::New();
caster->SetInput(resampler->GetOutput());
itk::WriteImage(caster->GetOutput(), "output.png");
return EXIT_SUCCESS;
}
void
CreateEllipseImage(ImageType::Pointer image)
{
using EllipseType = itk::EllipseSpatialObject<Dimension>;
using SpatialObjectToImageFilterType = itk::SpatialObjectToImageFilter<EllipseType, ImageType>;
auto imageFilter = SpatialObjectToImageFilterType::New();
ImageType::SizeType size;
size[0] = 100;
size[1] = 100;
imageFilter->SetSize(size);
ImageType::SpacingType spacing;
spacing.Fill(1);
imageFilter->SetSpacing(spacing);
auto ellipse = EllipseType::New();
EllipseType::ArrayType radiusArray;
radiusArray[0] = 10;
radiusArray[1] = 20;
ellipse->SetRadiusInObjectSpace(radiusArray);
using TransformType = EllipseType::TransformType;
auto 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
CreateSphereImage(ImageType::Pointer image)
{
using EllipseType = itk::EllipseSpatialObject<Dimension>;
using SpatialObjectToImageFilterType = itk::SpatialObjectToImageFilter<EllipseType, ImageType>;
auto imageFilter = SpatialObjectToImageFilterType::New();
ImageType::SizeType size;
size[0] = 100;
size[1] = 100;
imageFilter->SetSize(size);
ImageType::SpacingType spacing;
spacing.Fill(1);
imageFilter->SetSpacing(spacing);
auto ellipse = EllipseType::New();
EllipseType::ArrayType radiusArray;
radiusArray[0] = 10;
radiusArray[1] = 10;
ellipse->SetRadiusInObjectSpace(radiusArray);
using TransformType = EllipseType::TransformType;
auto 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>
class AffineTransform : public itk::MatrixOffsetTransformBase<TParametersValueType, NDimensions, NDimensions> Affine transformation of a vector space (e.g. space coordinates)
This class allows the definition and manipulation of affine transformations of an n-dimensional affine space (and its associated vector space) onto itself. One common use is to define and manipulate Euclidean coordinate transformations in two and three dimensions, but other uses are possible as well.
An affine transformation is defined mathematically as a linear transformation plus a constant offset. If A is a constant n x n matrix and b is a constant n-vector, then y = Ax+b defines an affine transformation from the n-vector x to the n-vector y.
The difference between two points is a vector and transforms linearly, using the matrix only. That is, (y1-y2) = A*(x1-x2).
The AffineTransform class determines whether to transform an object as a point or a vector by examining its type. An object of type Point transforms as a point; an object of type Vector transforms as a vector.
One common use of affine transformations is to define coordinate conversions in two- and three-dimensional space. In this application, x is a two- or three-dimensional vector containing the “source” coordinates of a point, y is a vector containing the “target” coordinates, the matrix A defines the scaling and rotation of the coordinate systems from the source to the target, and b defines the translation of the origin from the source to the target. More generally, A can also define anisotropic scaling and shearing transformations. Any good textbook on computer graphics will discuss coordinate transformations in more detail. Several of the methods in this class are designed for this purpose and use the language appropriate to coordinate conversions.
Any two affine transformations may be composed and the result is another affine transformation. However, the order is important. Given two affine transformations T1 and T2, we will say that “precomposing T1 with T2” yields the transformation which applies T1 to the source, and then applies T2 to that result to obtain the target. Conversely, we will say that “postcomposing T1 with T2” yields the transformation which applies T2 to the source, and then applies T1 to that result to obtain the target. (Whether T1 or T2 comes first lexicographically depends on whether you choose to write mappings from right-to-left or vice versa; we avoid the whole problem by referring to the order of application rather than the textual order.)
There are two template parameters for this class:
TParametersValueType The type to be used for scalar numeric values. Either float or double.
NDimensions The number of dimensions of the vector space.
This class provides several methods for setting the matrix and vector defining the transform. To support the registration framework, the transform parameters can also be set as an Array<double> of size (NDimension + 1) * NDimension using method SetParameters(). The first (NDimension x NDimension) parameters defines the matrix in row-major order (where the column index varies the fastest). The last NDimension parameters defines the translation in each dimensions.
This class also supports the specification of a center of rotation (center) and a translation that is applied with respect to that centered rotation. By default the center of rotation is set to the origin.
Subclassed by itk::AzimuthElevationToCartesianTransform< TParametersValueType, NDimensions >, itk::CenteredAffineTransform< TParametersValueType, NDimensions >, itk::ScalableAffineTransform< TParametersValueType, NDimensions >
-
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.
