Global Registration of Two Images#
Synopsis#
A basic global registration of two images.
Results#
fixing.png#
moving.png#
output.png#
Output:
Result =
Translation X = 15.0103
Translation Y = -1.12679
Iterations = 15
Metric value = 2374.3
Code#
C++#
#include "itkCastImageFilter.h"
#include "itkEllipseSpatialObject.h"
#include "itkImage.h"
#include "itkImageRegistrationMethod.h"
#include "itkLinearInterpolateImageFunction.h"
#include "itkImageFileWriter.h"
#include "itkMeanSquaresImageToImageMetric.h"
#include "itkRegularStepGradientDescentOptimizer.h"
#include "itkResampleImageFilter.h"
#include "itkRescaleIntensityImageFilter.h"
#include "itkSpatialObjectToImageFilter.h"
#include "itkTranslationTransform.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::TranslationTransform<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());
initialParameters[0] = 0.0; // Initial offset along X
initialParameters[1] = 0.0; // Initial offset along Y
registration->SetInitialTransformParameters(initialParameters);
optimizer->SetMaximumStepLength(4.00);
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();
// In the case of the \doxygen{TranslationTransform}, there is a
// straightforward interpretation of the parameters. Each element of the
// array corresponds to a translation along one spatial dimension.
const double TranslationAlongX = finalParameters[0];
const double TranslationAlongY = finalParameters[1];
// The optimizer can be queried for the actual number of iterations
// performed to reach convergence. The \code{GetCurrentIteration()}
// method returns this value. A large number of iterations may be an
// indication that the maximum step length has been set too small, which
// is undesirable since it results in long computational times.
const unsigned int numberOfIterations = optimizer->GetCurrentIteration();
// 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 << " Translation X = " << TranslationAlongX << std::endl;
std::cout << " Translation Y = " << TranslationAlongY << std::endl;
std::cout << " Iterations = " << numberOfIterations << 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}. Please
// refer to Section~\ref{sec:ResampleImageFilter} for details on the use
// of this filter. First, a ResampleImageFilter type is instantiated
// using the image types. It is convenient to use the fixed image type as
// the output type since it is likely that the transformed moving image
// will be compared with the fixed image.
using ResampleFilterType = itk::ResampleImageFilter<ImageType, ImageType>;
// A resampling filter is created and the moving image is connected as its input.
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");
/*
// The fixed image and the transformed moving image can easily be compared
// using the \doxygen{SubtractImageFilter}. This pixel-wise filter computes
// the difference between homologous pixels of its two input images.
using DifferenceFilterType = itk::SubtractImageFilter<
FixedImageType,
FixedImageType,
FixedImageType >;
auto difference = DifferenceFilterType::New();
difference->SetInput1( fixedImageReader->GetOutput() );
difference->SetInput2( resampler->GetOutput() );
*/
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 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.
-
template<typename TParametersValueType = double, unsigned int NDimensions = 3>
class TranslationTransform : public itk::Transform<TParametersValueType, NDimensions, NDimensions> Translation transformation of a vector space (e.g. space coordinates)
The same functionality could be obtained by using the Affine transform, but with a large difference in performance.
- ITK Sphinx Examples:
