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]
Output Image#
Code#
C++#
#include "itkImageRegistrationMethod.h"
#include "itkMeanSquaresImageToImageMetric.h"
#include "itkTimeProbesCollectorBase.h"
#include "itkSpatialObjectToImageFilter.h"
#include "itkEllipseSpatialObject.h"
#include "itkBSplineTransform.h"
#include "itkLBFGSOptimizer.h"
#include "itkImageFileWriter.h"
#include "itkResampleImageFilter.h"
#include "itkCastImageFilter.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(const ImageType::Pointer & image);
static void
CreateCircleImage(const ImageType::Pointer & image);
int
main(int itkNotUsed(argc), char * itkNotUsed(argv)[])
{
const unsigned int SpaceDimension = ImageDimension;
constexpr unsigned int SplineOrder = 3;
using CoordinateRepType = double;
using TransformType = itk::BSplineTransform<CoordinateRepType, SpaceDimension, SplineOrder>;
using OptimizerType = itk::LBFGSOptimizer;
using MetricType = itk::MeanSquaresImageToImageMetric<ImageType, ImageType>;
using InterpolatorType = itk::LinearInterpolateImageFunction<ImageType, double>;
using RegistrationType = itk::ImageRegistrationMethod<ImageType, ImageType>;
auto metric = MetricType::New();
auto optimizer = OptimizerType::New();
auto interpolator = InterpolatorType::New();
auto 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);
auto transform = TransformType::New();
registration->SetTransform(transform);
// Create the synthetic images
auto fixedImage = ImageType::New();
CreateCircleImage(fixedImage);
auto 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.
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());
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 (const 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>;
auto 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>;
auto caster = CastFilterType::New();
caster->SetInput(resample->GetOutput());
try
{
itk::WriteImage(caster->GetOutput(), "output.png");
}
catch (const 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(const ImageType::Pointer & image)
{
using EllipseType = itk::EllipseSpatialObject<ImageDimension>;
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
CreateCircleImage(const ImageType::Pointer & image)
{
using EllipseType = itk::EllipseSpatialObject<ImageDimension>;
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, 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
and is varied by defining a deformation 
of each control point. The deformation 
at any point 
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
, 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
such that the transformed point is given by: 
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:
-
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.
