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
#include "itkImageRegistrationMethodv4.h"
 
#include "itkMattesMutualInformationImageToImageMetricv4.h"
 
#include "itkRegularStepGradientDescentOptimizerv4.h"
#include "itkConjugateGradientLineSearchOptimizerv4.h"
 
#include "itkTranslationTransform.h"
#include "itkAffineTransform.h"
#include "itkCompositeTransform.h"
// Software Guide : EndCodeSnippet
 
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
 
#include "itkResampleImageFilter.h"
#include "itkCastImageFilter.h"
#include "itkCheckerBoardImageFilter.h"
 
#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 OptimizerType = itk::GradientDescentOptimizerv4Template<double>;
  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
  using TTransformType = itk::TranslationTransform<double, Dimension>;
  // 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 TOptimizerType = itk::RegularStepGradientDescentOptimizerv4<double>;
  using MetricType =
    itk::MattesMutualInformationImageToImageMetricv4<FixedImageType,
                                                     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
  using ATransformType = itk::AffineTransform<double, Dimension>;
  // 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 =
    itk::ConjugateGradientLineSearchOptimizerv4Template<double>;
  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;
  using RegionType = FixedImageType::RegionType;
  using SizeType = FixedImageType::SizeType;
 
  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 =
    itk::RegistrationParameterScalesFromPhysicalShift<MetricType>;
  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 =
    itk::ResampleImageFilter<MovingImageType, FixedImageType>;
  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 =
    itk::CastImageFilter<FixedImageType, OutputImageType>;
  using WriterType = itk::ImageFileWriter<OutputImageType>;
 
  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
  using TransformType = itk::IdentityTransform<double, Dimension>;
  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;
}