Loading [MathJax]/extensions/TeX/AMSmath.js
ITK 6.0.0
Insight Toolkit
 
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages
Loading...
Searching...
No Matches
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
*
* https://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);
const unsigned int currentLevel = registration->GetCurrentLevel();
const 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 };
};
namespace
{
int
ExampleMain(int argc, const char * const 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
auto transOptimizer = TOptimizerType::New();
auto transMetric = MetricType::New();
auto 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
auto 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>;
auto compositeTransform = CompositeTransformType::New();
compositeTransform->AddTransform(movingInitTx);
// Software Guide : EndCodeSnippet
const auto fixedImage = itk::ReadImage<FixedImageType>(argv[1]);
const auto movingImage = itk::ReadImage<MovingImageType>(argv[2]);
transRegistration->SetFixedImage(fixedImage);
transRegistration->SetMovingImage(movingImage);
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 aggressive 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.
//
auto observer1 = CommandIterationUpdate::New();
transOptimizer->AddObserver(itk::IterationEvent(), observer1);
// Create the Registration interface observer and register it with the
// registration method.
//
using TranslationCommandType =
RegistrationInterfaceCommand<TRegistrationType>;
auto 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
transRegistration->Update();
std::cout
<< "Optimizer stop condition: "
<< transRegistration->GetOptimizer()->GetStopConditionDescription()
<< std::endl;
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
auto affineOptimizer = AOptimizerType::New();
auto affineMetric = MetricType::New();
auto 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(fixedImage);
affineRegistration->SetMovingImage(movingImage);
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 registration method, so first
// we instantiate an object of the output transform type.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
auto 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;
using RegionType = FixedImageType::RegionType;
using SizeType = FixedImageType::SizeType;
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 =
auto 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 beginning 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.
//
auto 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>;
auto 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
affineRegistration->Update();
std::cout
<< "Optimizer stop condition: "
<< affineRegistration->GetOptimizer()->GetStopConditionDescription()
<< std::endl;
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 =
auto resample = ResampleFilterType::New();
resample->SetTransform(compositeTransform);
resample->SetInput(movingImage);
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 =
auto caster = CastFilterType::New();
caster->SetInput(resample->GetOutput());
itk::WriteImage(caster->GetOutput(), argv[3]);
// 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>;
auto checker = CheckerBoardFilterType::New();
checker->SetInput1(fixedImage);
checker->SetInput2(resample->GetOutput());
caster->SetInput(checker->GetOutput());
resample->SetDefaultPixelValue(0);
// Write out checkerboard outputs
// Before registration
TransformType::Pointer identityTransform;
identityTransform = TransformType::New();
identityTransform->SetIdentity();
resample->SetTransform(identityTransform);
if (argc > 5)
{
itk::WriteImage(caster->GetOutput(), argv[5]);
}
// After registration
resample->SetTransform(compositeTransform);
if (argc > 6)
{
itk::WriteImage(caster->GetOutput(), argv[6]);
}
return EXIT_SUCCESS;
}
} // namespace
int
main(int argc, char * argv[])
{
try
{
return ExampleMain(argc, argv);
}
catch (const itk::ExceptionObject & exceptionObject)
{
std::cerr << "ITK exception caught:\n" << exceptionObject << '\n';
}
catch (const std::exception & stdException)
{
std::cerr << "std exception caught:\n" << stdException.what() << '\n';
}
catch (...)
{
std::cerr << "Unhandled exception!\n";
}
return EXIT_FAILURE;
}
Casts input pixels to output pixel type.
Combines two images in a checkerboard pattern.
Superclass for callback/observer methods.
Definition itkCommand.h:46
virtual void Execute(Object *caller, const EventObject &event)=0
This class contains a list of transforms and concatenates them by composition.
Conjugate gradient descent optimizer with a golden section line search for nonlinear optimization.
Abstraction of the Events used to communicating among filters and with GUIs.
Standard exception handling object.
Implementation of an Identity Transform.
Interface method for the current registration framework.
Templated n-dimensional image class.
Definition itkImage.h:89
Computes the mutual information between two images to be registered using the method of Mattes et al.
Base class for most ITK classes.
Definition itkObject.h:62
virtual const std::string & GetObjectName() const
Registration helper class for estimating scales of transform parameters a step sizes,...
Resample an image via a coordinate transform.
Implements transparent reference counting.
Translation transformation of a vector space (e.g. space coordinates)
ITK_TEMPLATE_EXPORT void WriteImage(TImagePointer &&image, const std::string &filename, bool compress=false)
TOutputImage::Pointer ReadImage(const std::string &filename)
const SizeValueType * GetSize() const
Definition itkSize.h:167