ITK  5.2.0
Insight Toolkit
Examples/RegistrationITKv4/MultiStageImageRegistration1.cxx
/*=========================================================================
*
* Copyright NumFOCUS
*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0.txt
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*
*=========================================================================*/
// Software Guide : BeginCommandLineArgs
// INPUTS: {BrainT1SliceBorder20.png}
// INPUTS: {BrainProtonDensitySliceR10X13Y17.png}
// OUTPUTS: {MultiStageImageRegistration1Output.png}
// ARGUMENTS: 100
// OUTPUTS: {MultiStageImageRegistration1CheckerboardBefore.png}
// OUTPUTS: {MultiStageImageRegistration1CheckerboardAfter.png}
// Software Guide : EndCommandLineArgs
// Software Guide : BeginLatex
//
// This example illustrates the use of more complex components of the
// registration framework. In particular, it introduces a multistage,
// multi-resolution approach to run a multi-modal registration process
// using two linear \doxygen{TranslationTransform} and
// \doxygen{AffineTransform}. Also, it shows the use of \emph{Scale
// Estimators} for fine-tuning the scale parameters of the optimizer when an
// Affine transform is used. The
// \doxygen{RegistrationParameterScalesFromPhysicalShift} filter is used for
// automatic estimation of the parameters scales.
//
// \index{itk::ImageRegistrationMethodv4!AffineTransform}
// \index{itk::ImageRegistrationMethodv4!Scaling parameter space}
// \index{itk::AffineTransform!Image Registration}
// \index{itk::ImageRegistrationMethodv4!Multi-Stage}
// \index{itk::ImageRegistrationMethodv4!Multi-Resolution}
// \index{itk::ImageRegistrationMethodv4!Multi-Modality}
//
// To begin the example, we include the headers of the registration
// components we will use.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
// Software Guide : EndCodeSnippet
#include "itkCommand.h"
// The following section of code implements a Command observer
// that will monitor the configurations of the registration process
// at every change of stage and resolution level.
template <typename TRegistration>
class RegistrationInterfaceCommand : public itk::Command
{
public:
using Self = RegistrationInterfaceCommand;
using Superclass = itk::Command;
using Pointer = itk::SmartPointer<Self>;
itkNewMacro(Self);
protected:
RegistrationInterfaceCommand() = default;
public:
using RegistrationType = TRegistration;
// The Execute function simply calls another version of the \code{Execute()}
// method accepting a \code{const} input object
void
Execute(itk::Object * object, const itk::EventObject & event) override
{
Execute((const itk::Object *)object, event);
}
void
Execute(const itk::Object * object, const itk::EventObject & event) override
{
if (!(itk::MultiResolutionIterationEvent().CheckEvent(&event)))
{
return;
}
std::cout << "\nObserving from class " << object->GetNameOfClass();
if (!object->GetObjectName().empty())
{
std::cout << " \"" << object->GetObjectName() << "\"" << std::endl;
}
const auto * registration = static_cast<const RegistrationType *>(object);
unsigned int currentLevel = registration->GetCurrentLevel();
typename RegistrationType::ShrinkFactorsPerDimensionContainerType
shrinkFactors =
registration->GetShrinkFactorsPerDimension(currentLevel);
typename RegistrationType::SmoothingSigmasArrayType smoothingSigmas =
registration->GetSmoothingSigmasPerLevel();
std::cout << "-------------------------------------" << std::endl;
std::cout << " Current multi-resolution level = " << currentLevel
<< std::endl;
std::cout << " shrink factor = " << shrinkFactors << std::endl;
std::cout << " smoothing sigma = " << smoothingSigmas[currentLevel]
<< std::endl;
std::cout << std::endl;
}
};
// The following section of code implements an observer
// that will monitor the evolution of the registration process.
class CommandIterationUpdate : public itk::Command
{
public:
using Self = CommandIterationUpdate;
using Superclass = itk::Command;
using Pointer = itk::SmartPointer<Self>;
itkNewMacro(Self);
protected:
CommandIterationUpdate() = default;
public:
using OptimizerPointer = const OptimizerType *;
void
Execute(itk::Object * caller, const itk::EventObject & event) override
{
Execute((const itk::Object *)caller, event);
}
void
Execute(const itk::Object * object, const itk::EventObject & event) override
{
auto optimizer = static_cast<OptimizerPointer>(object);
if (!(itk::IterationEvent().CheckEvent(&event)))
{
return;
}
std::cout << optimizer->GetCurrentIteration() << " ";
std::cout << optimizer->GetValue() << " ";
std::cout << optimizer->GetCurrentPosition() << " "
<< m_CumulativeIterationIndex++ << std::endl;
}
private:
unsigned int m_CumulativeIterationIndex{ 0 };
};
int
main(int argc, char * argv[])
{
if (argc < 4)
{
std::cerr << "Missing Parameters " << std::endl;
std::cerr << "Usage: " << argv[0];
std::cerr << " fixedImageFile movingImageFile ";
std::cerr << " outputImagefile [backgroundGrayLevel]";
std::cerr << " [checkerboardbefore] [CheckerBoardAfter]";
std::cerr << " [numberOfBins] " << std::endl;
return EXIT_FAILURE;
}
constexpr unsigned int Dimension = 2;
using PixelType = float;
using FixedImageType = itk::Image<PixelType, Dimension>;
using MovingImageType = itk::Image<PixelType, Dimension>;
// Software Guide : BeginLatex
//
// In a multistage scenario, each stage needs an individual instantiation
// of the \doxygen{ImageRegistrationMethodv4}, so each stage can possibly
// have a different transform, a different optimizer, and a different image
// metric and can be performed in multiple levels.
// The configuration of the registration method at each stage closely
// follows the procedure in the previous section.
//
// In early stages we can use simpler transforms and more aggressive
// optimization parameters to take big steps toward the optimal value.
// Then, at the final stage we can have a more complex transform to do fine
// adjustments of the final parameters.
//
// A possible scheme is to use a simple translation transform for initial
// coarse registration levels and upgrade to an affine transform at the
// finer level.
// Since we have two different types of transforms, we can use a multistage
// registration approach as shown in the current example.
//
// First we need to configure the registration components of the initial
// stage. The instantiation of the transform type requires only the
// dimension of the space and the type used for representing space
// coordinates.
//
// \index{itk::TranslationTransform!Instantiation}
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// The types of other registration components are defined here.\newline
// \doxygen{RegularStepGradientDescentOptimizerv4} is used as the
// optimizer of the first stage. Also, we use
// \doxygen{MattesMutualInformationImageToImageMetricv4} as the metric
// since it is fitted for a multi-modal registration.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
using MetricType =
MovingImageType>;
using TRegistrationType = itk::ImageRegistrationMethodv4<FixedImageType,
MovingImageType,
TTransformType>;
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// Then, all the components are instantiated using their \code{New()}
// method and connected to the registration object as in previous examples.
//
// Software Guide : EndLatex
TOptimizerType::Pointer transOptimizer = TOptimizerType::New();
MetricType::Pointer transMetric = MetricType::New();
TRegistrationType::Pointer transRegistration = TRegistrationType::New();
transRegistration->SetOptimizer(transOptimizer);
transRegistration->SetMetric(transMetric);
// Software Guide : BeginLatex
//
// The output transform of the registration process will be constructed
// internally in the registration filter since the related
// \emph{TransformType} is already passed to the registration method as a
// template parameter. However, we should provide an initial moving
// transform for the registration method if needed.
//
// \index{itk::TranslationTransform!New()}
// \index{itk::TranslationTransform!Pointer}
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
TTransformType::Pointer movingInitTx = TTransformType::New();
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// After setting the initial parameters, the initial transform can be
// passed to the registration filter by \code{SetMovingInitialTransform()}
// method.
//
// \index{itk::Image\-Registration\-Methodv4!SetMovingInitialTransform()}
//
// Software Guide : EndLatex
using ParametersType = TOptimizerType::ParametersType;
ParametersType initialParameters(movingInitTx->GetNumberOfParameters());
initialParameters[0] = 3.0; // Initial offset in mm along X
initialParameters[1] = 5.0; // Initial offset in mm along Y
movingInitTx->SetParameters(initialParameters);
// Software Guide : BeginCodeSnippet
transRegistration->SetMovingInitialTransform(movingInitTx);
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// We can use a \doxygen{CompositeTransform} to stack all the output
// transforms resulted from multiple stages. This composite
// transform should also hold the moving initial transform (if it exists)
// because as explained in section \ref{sec:RigidRegistrationIn2D},
// the output of each registration stage does not include the input
// initial transform to that stage.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
using CompositeTransformType = itk::CompositeTransform<double, Dimension>;
CompositeTransformType::Pointer compositeTransform =
CompositeTransformType::New();
compositeTransform->AddTransform(movingInitTx);
// Software Guide : EndCodeSnippet
using FixedImageReaderType = itk::ImageFileReader<FixedImageType>;
using MovingImageReaderType = itk::ImageFileReader<MovingImageType>;
FixedImageReaderType::Pointer fixedImageReader =
FixedImageReaderType::New();
MovingImageReaderType::Pointer movingImageReader =
MovingImageReaderType::New();
fixedImageReader->SetFileName(argv[1]);
movingImageReader->SetFileName(argv[2]);
transRegistration->SetFixedImage(fixedImageReader->GetOutput());
transRegistration->SetMovingImage(movingImageReader->GetOutput());
transRegistration->SetObjectName("TranslationRegistration");
// Software Guide : BeginLatex
//
// In the case of this simple example, the first stage is run only
// in one level of registration at a coarse resolution.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
constexpr unsigned int numberOfLevels1 = 1;
TRegistrationType::ShrinkFactorsArrayType shrinkFactorsPerLevel1;
shrinkFactorsPerLevel1.SetSize(numberOfLevels1);
shrinkFactorsPerLevel1[0] = 3;
TRegistrationType::SmoothingSigmasArrayType smoothingSigmasPerLevel1;
smoothingSigmasPerLevel1.SetSize(numberOfLevels1);
smoothingSigmasPerLevel1[0] = 2;
transRegistration->SetNumberOfLevels(numberOfLevels1);
transRegistration->SetShrinkFactorsPerLevel(shrinkFactorsPerLevel1);
transRegistration->SetSmoothingSigmasPerLevel(smoothingSigmasPerLevel1);
// Software Guide : EndCodeSnippet
transMetric->SetNumberOfHistogramBins(24);
if (argc > 7)
{
// optionally, override the values with numbers taken from the command
// line arguments.
transMetric->SetNumberOfHistogramBins(std::stoi(argv[7]));
}
transOptimizer->SetNumberOfIterations(200);
transOptimizer->SetRelaxationFactor(0.5);
// Software Guide : BeginLatex
//
// Also, for this initial stage we can use a more agressive parameter
// set for the optimizer by taking a big step size and relaxing
// stop criteria.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
transOptimizer->SetLearningRate(16);
transOptimizer->SetMinimumStepLength(1.5);
// Software Guide : EndCodeSnippet
// Create the Command observer and register it with the optimizer.
//
CommandIterationUpdate::Pointer observer1 = CommandIterationUpdate::New();
transOptimizer->AddObserver(itk::IterationEvent(), observer1);
// Create the Registration interface observer and register it with the
// registration method.
//
using TranslationCommandType =
RegistrationInterfaceCommand<TRegistrationType>;
TranslationCommandType::Pointer command1 = TranslationCommandType::New();
transRegistration->AddObserver(itk::MultiResolutionIterationEvent(),
command1);
// Software Guide : BeginLatex
//
// Once all the registration components are in place, we trigger the
// registration process by calling \code{Update()} and add the result
// output transform to the final composite transform, so this composite
// transform can be used to initialize the next registration stage.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
try
{
transRegistration->Update();
std::cout
<< "Optimizer stop condition: "
<< transRegistration->GetOptimizer()->GetStopConditionDescription()
<< std::endl;
}
catch (const itk::ExceptionObject & err)
{
std::cout << "ExceptionObject caught !" << std::endl;
std::cout << err << std::endl;
return EXIT_FAILURE;
}
compositeTransform->AddTransform(
transRegistration->GetModifiableTransform());
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// Now we can upgrade to an Affine transform as the second stage of
// registration process. The AffineTransform is a linear transformation
// that maps lines into lines. It can be used to represent translations,
// rotations, anisotropic scaling, shearing or any combination of them.
// Details about the affine transform can be seen in
// Section~\ref{sec:AffineTransform}. The instantiation of the transform
// type requires only the dimension of the space and the type used for
// representing space coordinates.
//
// \index{itk::AffineTransform!Instantiation}
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// We also use a different optimizer in configuration of the second stage
// while the metric is kept the same as before.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
using AOptimizerType =
using ARegistrationType = itk::ImageRegistrationMethodv4<FixedImageType,
MovingImageType,
ATransformType>;
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// Again all the components are instantiated using their \code{New()}
// method and connected to the registration object like in previous stages.
//
// Software Guide : EndLatex
AOptimizerType::Pointer affineOptimizer = AOptimizerType::New();
MetricType::Pointer affineMetric = MetricType::New();
ARegistrationType::Pointer affineRegistration = ARegistrationType::New();
affineRegistration->SetOptimizer(affineOptimizer);
affineRegistration->SetMetric(affineMetric);
// Software Guide : BeginLatex
//
// The current stage can be initialized using the initial transform of the
// registration and the result transform of the previous stage, so that
// both are concatenated into the composite transform.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
affineRegistration->SetMovingInitialTransform(compositeTransform);
// Software Guide : EndCodeSnippet
affineRegistration->SetFixedImage(fixedImageReader->GetOutput());
affineRegistration->SetMovingImage(movingImageReader->GetOutput());
affineRegistration->SetObjectName("AffineRegistration");
affineMetric->SetNumberOfHistogramBins(24);
if (argc > 7)
{
// optionally, override the values with numbers taken from the command
// line arguments.
affineMetric->SetNumberOfHistogramBins(std::stoi(argv[7]));
}
// Software Guide : BeginLatex
//
// In Section \ref{sec:InitializingRegistrationWithMoments} we showed
// the importance of center of rotation in the registration process.
// In Affine transforms, the center of rotation is defined by the fixed
// parameters set, which are set by default to [0, 0].
// However, consider a situation where the
// origin of the virtual space, in which the registration is run, is far
// away from the zero origin. In such cases, leaving the center of rotation
// as the default value can make the optimization process unstable.
// Therefore, we are always interested to set the center of rotation to the
// center of virtual space which is usually the fixed image space.
//
// Note that either center of gravity or geometrical center can be used
// as the center of rotation. In this example center of rotation is set
// to the geometrical center of the fixed image. We could also use
// \doxygen{ImageMomentsCalculator} filter to compute the center of mass.
//
// Based on the above discussion, the user must set the fixed parameters of
// the registration transform outside of the registraton method, so first
// we instantiate an object of the output transform type.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
ATransformType::Pointer affineTx = ATransformType::New();
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// Then, we compute the physical center of the fixed image and set
// that as the center of the output Affine transform.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
using SpacingType = FixedImageType::SpacingType;
using OriginType = FixedImageType::PointType;
FixedImageType::Pointer fixedImage = fixedImageReader->GetOutput();
const SpacingType fixedSpacing = fixedImage->GetSpacing();
const OriginType fixedOrigin = fixedImage->GetOrigin();
const RegionType fixedRegion = fixedImage->GetLargestPossibleRegion();
const SizeType fixedSize = fixedRegion.GetSize();
ATransformType::InputPointType centerFixed;
centerFixed[0] = fixedOrigin[0] + fixedSpacing[0] * fixedSize[0] / 2.0;
centerFixed[1] = fixedOrigin[1] + fixedSpacing[1] * fixedSize[1] / 2.0;
const unsigned int numberOfFixedParameters =
affineTx->GetFixedParameters().Size();
ATransformType::ParametersType fixedParameters(numberOfFixedParameters);
for (unsigned int i = 0; i < numberOfFixedParameters; ++i)
{
fixedParameters[i] = centerFixed[i];
}
affineTx->SetFixedParameters(fixedParameters);
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// Then, the initialized output transform should be connected to
// the registration object by using \code{SetInitialTransform()} method.
//
// It is important to distinguish between the \code{SetInitialTransform()}
// and \code{SetMovingInitialTransform()} that was used to initialize the
// registration stage based on the results of the previous stages.
// You can assume that the first one is used for direct manipulation of the
// optimizable transform in current registration process.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
affineRegistration->SetInitialTransform(affineTx);
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// The set of optimizable parameters in the Affine transform have different
// dynamic ranges. Typically the parameters associated with the matrix
// have values around $[-1:1]$, although they are not restricted to this
// interval. Parameters associated with translations, on the other hand,
// tend to have much higher values, typically on the order of $10.0$ to
// $100.0$. This difference in dynamic range negatively affects the
// performance of gradient descent optimizers. ITK provides some mechanisms
// to compensate for such differences in values among the parameters when
// they are passed to the optimizer.
//
// The first mechanism consists of providing an
// array of scale factors to the optimizer. These factors re-normalize the
// gradient components before they are used to compute the step of the
// optimizer at the current iteration.
// These scales are estimated by the user intuitively as shown in previous
// examples of this chapter. In our particular case, a common choice
// for the scale parameters is to set all those associated
// with the matrix coefficients to $1.0$, that is, the first $N \times N$
// factors. Then, we set the remaining scale factors to a small value.
//
// Software Guide : EndLatex
// Software Guide : BeginLatex
//
// Here the affine transform is represented by the matrix $\bf{M}$ and the
// vector $\bf{T}$. The transformation of a point $\bf{P}$ into $\bf{P'}$
// is expressed as
//
// \begin{equation}
// \left[
// \begin{array}{c}
// {P'}_x \\ {P'}_y \\ \end{array}
// \right]
// =
// \left[
// \begin{array}{cc}
// M_{11} & M_{12} \\ M_{21} & M_{22} \\ \end{array}
// \right]
// \cdot
// \left[
// \begin{array}{c}
// P_x \\ P_y \\ \end{array}
// \right]
// +
// \left[
// \begin{array}{c}
// T_x \\ T_y \\ \end{array}
// \right]
// \end{equation}
//
//
// Software Guide : EndLatex
// Software Guide : BeginLatex
//
// Based on the above discussion, we need much smaller scales for
// translation parameters of vector $\bf{T}$ ($T_x$, $T_y$) compared to the
// parameters of matrix $\bf{M}$ ($M_{11}$, $M_{12}$, $M_{21}$, $M_{22}$).
// However, it is not easy to have an intuitive estimation of all parameter
// scales when we have to deal with a large parameter space.
//
// Fortunately, ITKv4 provides a framework for automated parameter scaling.
// \doxygen{RegistrationParameterScalesEstimator} vastly reduces the
// difficulty of tuning parameters for different transform/metric
// combinations. Parameter scales are estimated by analyzing the result of
// a small parameter update on the change in the magnitude of physical
// space deformation induced by the transformation.
//
// The impact from a unit change of a parameter may be defined in multiple
// ways, such as the maximum shift of voxels in index or physical space, or
// the average norm of transform Jacobian. Filters
// \doxygen{RegistrationParameterScalesFromPhysicalShift} and
// \doxygen{RegistrationParameterScalesFromIndexShift} use the first
// definition to estimate the scales, while the
// \doxygen{RegistrationParameterScalesFromJacobian} filter estimates
// scales based on the later definition. In all methods, the goal is to
// rescale the transform parameters such that a unit change of each
// \emph{scaled parameter} will have the same impact on deformation.
//
// In this example the first filter is chosen to estimate the parameter
// scales. The scales estimator will then be passed to optimizer.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
using ScalesEstimatorType =
ScalesEstimatorType::Pointer scalesEstimator = ScalesEstimatorType::New();
scalesEstimator->SetMetric(affineMetric);
scalesEstimator->SetTransformForward(true);
affineOptimizer->SetScalesEstimator(scalesEstimator);
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// The step length has to be proportional to the expected values of the
// parameters in the search space. Since the expected values of the matrix
// coefficients are around $1.0$, the initial step of the optimization
// should be a small number compared to $1.0$. As a guideline, it is
// useful to think of the matrix coefficients as combinations of
// $cos(\theta)$ and $sin(\theta)$. This leads to use values close to the
// expected rotation measured in radians. For example, a rotation of $1.0$
// degree is about $0.017$ radians.
//
// However, we need not worry about the above considerations.
// Thanks to the \emph{ScalesEstimator}, the initial step size can also be
// estimated automatically, either at each iteration or only at the first
// iteration. In this example we choose to estimate learning rate
// once at the begining of the registration process.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
affineOptimizer->SetDoEstimateLearningRateOnce(true);
affineOptimizer->SetDoEstimateLearningRateAtEachIteration(false);
// Software Guide : EndCodeSnippet
// Set the other parameters of optimizer
affineOptimizer->SetLowerLimit(0);
affineOptimizer->SetUpperLimit(2);
affineOptimizer->SetEpsilon(0.2);
affineOptimizer->SetNumberOfIterations(200);
affineOptimizer->SetMinimumConvergenceValue(1e-6);
affineOptimizer->SetConvergenceWindowSize(5);
// Create the Command observer and register it with the optimizer.
//
CommandIterationUpdate::Pointer observer2 = CommandIterationUpdate::New();
affineOptimizer->AddObserver(itk::IterationEvent(), observer2);
// Software Guide : BeginLatex
//
// At the second stage, we run two levels of registration, where the second
// level is run in full resolution in which we do the final adjustments
// of the output parameters.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
constexpr unsigned int numberOfLevels2 = 2;
ARegistrationType::ShrinkFactorsArrayType shrinkFactorsPerLevel2;
shrinkFactorsPerLevel2.SetSize(numberOfLevels2);
shrinkFactorsPerLevel2[0] = 2;
shrinkFactorsPerLevel2[1] = 1;
ARegistrationType::SmoothingSigmasArrayType smoothingSigmasPerLevel2;
smoothingSigmasPerLevel2.SetSize(numberOfLevels2);
smoothingSigmasPerLevel2[0] = 1;
smoothingSigmasPerLevel2[1] = 0;
affineRegistration->SetNumberOfLevels(numberOfLevels2);
affineRegistration->SetShrinkFactorsPerLevel(shrinkFactorsPerLevel2);
affineRegistration->SetSmoothingSigmasPerLevel(smoothingSigmasPerLevel2);
// Software Guide : EndCodeSnippet
// Create the Registration interface observer and register it with the
// registration object.
using AffineCommandType = RegistrationInterfaceCommand<ARegistrationType>;
AffineCommandType::Pointer command2 = AffineCommandType::New();
affineRegistration->AddObserver(itk::MultiResolutionIterationEvent(),
command2);
// Software Guide : BeginLatex
//
// Finally we trigger the registration process by calling \code{Update()}
// and add the output transform of the last stage to the composite
// transform. This composite transform will be considered as the final
// transform of this multistage registration process and will be used by
// the resampler to resample the moving image in to the virtual domain
// space (fixed image space if there is no fixed initial transform).
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
try
{
affineRegistration->Update();
std::cout
<< "Optimizer stop condition: "
<< affineRegistration->GetOptimizer()->GetStopConditionDescription()
<< std::endl;
}
catch (const itk::ExceptionObject & err)
{
std::cout << "ExceptionObject caught !" << std::endl;
std::cout << err << std::endl;
return EXIT_FAILURE;
}
compositeTransform->AddTransform(
affineRegistration->GetModifiableTransform());
// Software Guide : EndCodeSnippet
std::cout << "\nInitial parameters of the registration process: "
<< std::endl
<< movingInitTx->GetParameters() << std::endl;
std::cout << "\nTranslation parameters after registration: " << std::endl
<< transOptimizer->GetCurrentPosition() << std::endl
<< " Last LearningRate: "
<< transOptimizer->GetCurrentStepLength() << std::endl;
std::cout << "\nAffine parameters after registration: " << std::endl
<< affineOptimizer->GetCurrentPosition() << std::endl
<< " Last LearningRate: " << affineOptimizer->GetLearningRate()
<< std::endl;
// Software Guide : BeginLatex
//
// Let's execute this example using the following multi-modality images:
//
// \begin{itemize}
// \item BrainT1SliceBorder20.png
// \item BrainProtonDensitySliceR10X13Y17.png
// \end{itemize}
//
// The second image is the result of intentionally rotating the first
// image by $10$ degrees and then translating by $(-13,-17)$. Both images
// have unit-spacing and are shown in Figure
// \ref{fig:FixedMovingMultiStageImageRegistration1}.
//
// The registration converges after $5$ iterations in the translation
// stage. Also, in the second stage, the registration converges after $46$
// iterations in the first level, and $6$ iterations in the second level.
// The final results when printed as an array of parameters are:
//
// \begin{verbatim}
// Initial parameters of the registration process:
// [3, 5]
//
// Translation parameters after first registration stage:
// [9.0346, 10.8303]
//
// Affine parameters after second registration stage:
// [0.9864, -0.1733, 0.1738, 0.9863, 0.9693, 0.1482]
// \end{verbatim}
//
// As it can be seen, the translation parameters after the first stage
// compensate most of the offset between the fixed and moving images.
// When the images are close to each other, the affine registration is
// run for the rotation and the final match.
// By reordering the Affine array of parameters as coefficients of matrix
// $\bf{M}$ and vector $\bf{T}$ they can now be seen as
//
// \begin{equation}
// M =
// \left[
// \begin{array}{cc}
// 0.9864 & -0.1733 \\ 0.1738 & 0.9863 \\ \end{array}
// \right]
// \mbox{ and }
// T =
// \left[
// \begin{array}{c}
// 0.9693 \\ 0.1482 \\ \end{array}
// \right]
// \end{equation}
//
// In this form, it is easier to interpret the effect of the
// transform. The matrix $\bf{M}$ is responsible for scaling, rotation and
// shearing while $\bf{T}$ is responsible for translations.
//
// The second component of the matrix values is usually associated with
// $\sin{\theta}$. We obtain the rotation through SVD of the affine
// matrix. The value is $9.975$ degrees, which is approximately the
// intentional misalignment of $10.0$ degrees.
//
// Also, let's compute the total translation values resulting from initial
// transform, translation transform, and the Affine transform together.
//
// In $X$ direction:
// \begin{equation}
// 3 + 9.0346 + 0.9693 = 13.0036
// \end{equation}
// In $Y$ direction:
// \begin{equation}
// 5 + 10.8303 + 0.1482 = 15.9785
// \end{equation}
//
// It can be seen that the translation values closely match the true
// misalignment introduced in the moving image.
//
// It is important to note that once the images are registered at a
// sub-pixel level, any further improvement of the registration relies
// heavily on the quality of the interpolator. It may then be reasonable to
// use a coarse and fast interpolator in the lower resolution levels and
// switch to a high-quality but slow interpolator in the final resolution
// level. However, in this example we used a linear interpolator for all
// stages and different registration levels since it is so fast.
//
// \begin{figure}
// \center
// \includegraphics[width=0.44\textwidth]{BrainProtonDensitySliceBorder20}
// \includegraphics[width=0.44\textwidth]{BrainProtonDensitySliceR10X13Y17}
// \itkcaption[AffineTransform registration]{Fixed and moving images
// provided as input to the registration method using the AffineTransform.}
// \label{fig:FixedMovingMultiStageImageRegistration1}
// \end{figure}
//
// Software Guide : EndLatex
using ResampleFilterType =
ResampleFilterType::Pointer resample = ResampleFilterType::New();
resample->SetTransform(compositeTransform);
resample->SetInput(movingImageReader->GetOutput());
PixelType backgroundGrayLevel = 100;
if (argc > 4)
{
backgroundGrayLevel = std::stoi(argv[4]);
}
resample->SetSize(fixedImage->GetLargestPossibleRegion().GetSize());
resample->SetOutputOrigin(fixedImage->GetOrigin());
resample->SetOutputSpacing(fixedImage->GetSpacing());
resample->SetOutputDirection(fixedImage->GetDirection());
resample->SetDefaultPixelValue(backgroundGrayLevel);
using OutputPixelType = unsigned char;
using OutputImageType = itk::Image<OutputPixelType, Dimension>;
using CastFilterType =
WriterType::Pointer writer = WriterType::New();
CastFilterType::Pointer caster = CastFilterType::New();
writer->SetFileName(argv[3]);
caster->SetInput(resample->GetOutput());
writer->SetInput(caster->GetOutput());
writer->Update();
// Software Guide : BeginLatex
//
// \begin{figure}
// \center
// \includegraphics[width=0.32\textwidth]{MultiStageImageRegistration1Output}
// \includegraphics[width=0.32\textwidth]{MultiStageImageRegistration1CheckerboardBefore}
// \includegraphics[width=0.32\textwidth]{MultiStageImageRegistration1CheckerboardAfter}
// \itkcaption[Multistage registration input images]{Mapped moving image
// (left) and composition of fixed and moving images before (center) and
// after (right) registration.}
// \label{fig:MultiStageImageRegistration1Outputs}
// \end{figure}
//
// The result of resampling the moving image is presented in the left image
// of Figure \ref{fig:MultiStageImageRegistration1Outputs}. The center and
// right images of the figure depict a checkerboard composite of the fixed
// and moving images before and after registration.
//
// Software Guide : EndLatex
// Generate checkerboards before and after registration
using CheckerBoardFilterType = itk::CheckerBoardImageFilter<FixedImageType>;
CheckerBoardFilterType::Pointer checker = CheckerBoardFilterType::New();
checker->SetInput1(fixedImage);
checker->SetInput2(resample->GetOutput());
caster->SetInput(checker->GetOutput());
writer->SetInput(caster->GetOutput());
resample->SetDefaultPixelValue(0);
// Write out checkerboard outputs
// Before registration
TransformType::Pointer identityTransform;
try
{
identityTransform = TransformType::New();
}
catch (const itk::ExceptionObject & err)
{
err.Print(std::cerr);
return EXIT_FAILURE;
}
identityTransform->SetIdentity();
resample->SetTransform(identityTransform);
if (argc > 5)
{
writer->SetFileName(argv[5]);
writer->Update();
}
// After registration
resample->SetTransform(compositeTransform);
if (argc > 6)
{
writer->SetFileName(argv[6]);
writer->Update();
}
return EXIT_SUCCESS;
}
itk::Object::GetObjectName
virtual const std::string & GetObjectName() const
itk::CastImageFilter
Casts input pixels to output pixel type.
Definition: itkCastImageFilter.h:104
itk::CompositeTransform
This class contains a list of transforms and concatenates them by composition.
Definition: itkCompositeTransform.h:87
itkRegularStepGradientDescentOptimizerv4.h
itk::IdentityTransform
Implementation of an Identity Transform.
Definition: itkIdentityTransform.h:50
itk::GradientDescentOptimizerv4Template
Gradient descent optimizer.
Definition: itkGradientDescentOptimizerv4.h:78
itk::GTest::TypedefsAndConstructors::Dimension2::PointType
ImageBaseType::PointType PointType
Definition: itkGTestTypedefsAndConstructors.h:51
itkImageFileReader.h
itk::CheckerBoardImageFilter
Combines two images in a checkerboard pattern.
Definition: itkCheckerBoardImageFilter.h:46
itk::GTest::TypedefsAndConstructors::Dimension2::SizeType
ImageBaseType::SizeType SizeType
Definition: itkGTestTypedefsAndConstructors.h:49
itk::SmartPointer< Self >
itkCastImageFilter.h
itkAffineTransform.h
itk::AffineTransform
Definition: itkAffineTransform.h:101
itkImageRegistrationMethodv4.h
itkTranslationTransform.h
itk::ImageFileReader
Data source that reads image data from a single file.
Definition: itkImageFileReader.h:75
itk::RegularStepGradientDescentOptimizerv4
Regular Step Gradient descent optimizer.
Definition: itkRegularStepGradientDescentOptimizerv4.h:47
itk::Command
Superclass for callback/observer methods.
Definition: itkCommand.h:45
itkCheckerBoardImageFilter.h
itk::ImageFileWriter
Writes image data to a single file.
Definition: itkImageFileWriter.h:87
itk::Command
class ITK_FORWARD_EXPORT Command
Definition: itkObject.h:43
itk::GTest::TypedefsAndConstructors::Dimension2::RegionType
ImageBaseType::RegionType RegionType
Definition: itkGTestTypedefsAndConstructors.h:54
itk::TranslationTransform
Translation transformation of a vector space (e.g. space coordinates)
Definition: itkTranslationTransform.h:43
itk::Command::Execute
virtual void Execute(Object *caller, const EventObject &event)=0
itkImageFileWriter.h
itk::ImageRegistrationMethodv4
Interface method for the current registration framework.
Definition: itkImageRegistrationMethodv4.h:117
itk::ResampleImageFilter
Resample an image via a coordinate transform.
Definition: itkResampleImageFilter.h:90
itk::Object
Base class for most ITK classes.
Definition: itkObject.h:62
itk::Math::e
static constexpr double e
Definition: itkMath.h:54
itk::RegistrationParameterScalesFromPhysicalShift
Registration helper class for estimating scales of transform parameters a step sizes,...
Definition: itkRegistrationParameterScalesFromPhysicalShift.h:35
itk::Image
Templated n-dimensional image class.
Definition: itkImage.h:86
itk::EventObject
Abstraction of the Events used to communicating among filters and with GUIs.
Definition: itkEventObject.h:57
itkCompositeTransform.h
itkResampleImageFilter.h
itk::ConjugateGradientLineSearchOptimizerv4Template
Conjugate gradient descent optimizer with a golden section line search for nonlinear optimization.
Definition: itkConjugateGradientLineSearchOptimizerv4.h:48
itk::GTest::TypedefsAndConstructors::Dimension2::Dimension
constexpr unsigned int Dimension
Definition: itkGTestTypedefsAndConstructors.h:44
itkCommand.h
itkMattesMutualInformationImageToImageMetricv4.h
itkConjugateGradientLineSearchOptimizerv4.h
itk::MattesMutualInformationImageToImageMetricv4
Computes the mutual information between two images to be registered using the method of Mattes et al.
Definition: itkMattesMutualInformationImageToImageMetricv4.h:103
itk::Size::GetSize
const SizeValueType * GetSize() const
Definition: itkSize.h:169