ITK multi-stage registration

Blog address ：https://blog.csdn.net/fanre/article/details/109118273

1、 First translation registration rough registration , Then affine transform is used for multi-resolution registration

#include "itkImageRegistrationMethodv4.h"// The header file 
 
#include "itkMattesMutualInformationImageToImageMetricv4.h"
 
#include "itkRegularStepGradientDescentOptimizerv4.h"
#include "itkConjugateGradientLineSearchOptimizerv4.h"
 
#include "itkTranslationTransform.h"
#include "itkAffineTransform.h"
#include "itkCompositeTransform.h"
 
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
 
#include "itkResampleImageFilter.h"
#include "itkCastImageFilter.h"
#include "itkCheckerBoardImageFilter.h"
 
#include "itkCommand.h"
 
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;
 
  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;
  }
};
 
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>;
 
  using TTransformType = itk::TranslationTransform<double, Dimension>;// First , We need to configure the registration component in the initial phase .transform The instantiation of a type requires only the dimension of the space and the type used to represent the spatial coordinates .
  
  using TOptimizerType = itk::RegularStepGradientDescentOptimizerv4<double>;// The types of other registration components are defined here . Use itk::RegularStepGradientDescentOptimizerv4 As the first stage optimizer .
  using MetricType =
    itk::MattesMutualInformationImageToImageMetricv4<FixedImageType,
                                                     MovingImageType>;// Besides , We use itk::MattesMutualInformationImageToImageMetricv4 As a measure , Because it is suitable for multi-mode registration .
  using TRegistrationType = itk::ImageRegistrationMethodv4<FixedImageType,
                                                           MovingImageType,
                                                           TTransformType>;
// All components use their New（） Method to instantiate and connect to the registration object ,
  TOptimizerType::Pointer    transOptimizer = TOptimizerType::New();
  MetricType::Pointer        transMetric = MetricType::New();
  TRegistrationType::Pointer transRegistration = TRegistrationType::New();
 
  transRegistration->SetOptimizer(transOptimizer);
  transRegistration->SetMetric(transMetric);
 
  TTransformType::Pointer movingInitTx = TTransformType::New();// The output of the registration process transform Will be constructed inside the registration filter , Because it's related TransformType It has been passed to the registration method as a template parameter . however , if necessary , We should provide an initial shift transformation for the registration method .
 
  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);
 
  transRegistration->SetMovingInitialTransform(movingInitTx);// After setting the initial parameters , Can pass SetMovingInitialTransform（） Method passes the initial transformation to the registration filter .
 
  using CompositeTransformType = itk::CompositeTransform<double, Dimension>;// We can use itk::CompositeTransform To stack all outputs from multiple stages transforms. This composite transformation should also include the initial transformation of the move ( If it exists ), Because as 3.6.1 As explained in section , The output of each registration stage does not include the initial conversion of input to that stage .
  CompositeTransformType::Pointer compositeTransform =
    CompositeTransformType::New();
  compositeTransform->AddTransform(movingInitTx);
 
  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");
 
  constexpr unsigned int numberOfLevels1 = 1;// In the case of this simple example , The first stage runs in one registration level only at coarse resolution .
 
  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);
 
  transMetric->SetNumberOfHistogramBins(24);
 
  if (argc > 7)
  {
    transMetric->SetNumberOfHistogramBins(std::stoi(argv[7]));
  }
 
  transOptimizer->SetNumberOfIterations(200);// Again , For this initial stage , We can use a more positive set of parameters for the optimizer by taking a larger step size and relaxing the stop condition .
  transOptimizer->SetRelaxationFactor(0.5);
 
  transOptimizer->SetLearningRate(16);
  transOptimizer->SetMinimumStepLength(1.5);
 
  CommandIterationUpdate::Pointer observer1 = CommandIterationUpdate::New();
  transOptimizer->AddObserver(itk::IterationEvent(), observer1);
 
  using TranslationCommandType =
    RegistrationInterfaceCommand<TRegistrationType>;
  TranslationCommandType::Pointer command1 = TranslationCommandType::New();
  transRegistration->AddObserver(itk::MultiResolutionIterationEvent(),
                                 command1);
  try
  {
    transRegistration->Update();// Once all registration components are in place , We call Update() To trigger the registration process , And output the result transform Add to the final composite transform , Therefore, the composite transformation can be used to initialize the next registration stage .
    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());
 
  using ATransformType = itk::AffineTransform<double, Dimension>;// Now we can upgrade to an affine transformation as the second stage of the registration process . Affine transformation is a linear transformation that maps a line to a line . It can be used to represent translation 、 rotate 、 Anisotropic scaling 、 Cut or any combination of them . Details about affine transformation can be found in 3.9.16 See in section . The instantiation of the transformation type only needs the dimension of the space and the type used to represent the spatial coordinates .
  
  using AOptimizerType =
    itk::ConjugateGradientLineSearchOptimizerv4Template<double>;// Different optimizers are used in the configuration of the second stage , At the same time, the measure remains the same as before .
  using ARegistrationType = itk::ImageRegistrationMethodv4<FixedImageType,
                                                           MovingImageType,
                                                           ATransformType>;
// Again , All components use their New（） Method to instantiate and connect to the registration object , Just like the previous stage .
  AOptimizerType::Pointer    affineOptimizer = AOptimizerType::New();
  MetricType::Pointer        affineMetric = MetricType::New();
  ARegistrationType::Pointer affineRegistration = ARegistrationType::New();
 
  affineRegistration->SetOptimizer(affineOptimizer);
  affineRegistration->SetMetric(affineMetric);
// The current phase can be initialized using the initial transformation of the registration and the result transformation of the previous phase , To combine both into a composite transform .
  affineRegistration->SetMovingInitialTransform(compositeTransform);
 
  affineRegistration->SetFixedImage(fixedImageReader->GetOutput());
  affineRegistration->SetMovingImage(movingImageReader->GetOutput());
  affineRegistration->SetObjectName("AffineRegistration");
 
  affineMetric->SetNumberOfHistogramBins(24);
 
  if (argc > 7)
  {
    affineMetric->SetNumberOfHistogramBins(std::stoi(argv[7]));
  }
  ATransformType::Pointer affineTx = ATransformType::New();// The user must set the registration transformation in addition to the registration method fixed Parameters , So first we instantiate the output transform Object of type .
 
  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;// In affine transformation , The center of rotation consists of fixed Parameter set definitions , The default setting is [0,0]. however , Consider a situation like this , That is, run the registered virtual The origin of space is far away from the zero origin . under these circumstances , Using the center of rotation as the default makes the optimization process unstable . therefore , We are always keen to set the rotation center at virtual Center of space ,virtual Space is usually fixed Image space .
// Be careful , Any center of gravity or geometric center can be used as the center of rotation . In this case , The rotation center is set to fixed The geometric center of the image . We can also use itk::ImageMomentsCalculator  Filter to calculate the center of mass .
// then , We calculated fixed Physical center of image , And set it as the center of the output affine transformation .
 
  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);
 
  affineRegistration->SetInitialTransform(affineTx);// then , You should use SetInitialTransform（） Method will initialize the output transform Connect to the registration object . It is important to distinguish SetInitialTransform() and SetMovingInitialTransform(), The latter is used to initialize the registration stage according to the results of the previous stages . You can assume that the first one is used to directly manipulate the optimizable transformation in the current registration process .
 
  using ScalesEstimatorType =
    itk::RegistrationParameterScalesFromPhysicalShift<MetricType>;
  ScalesEstimatorType::Pointer scalesEstimator = ScalesEstimatorType::New();// choice physical Filter to estimate parameters scales.scales estimator Then pass it to optimizer.
  scalesEstimator->SetMetric(affineMetric);
  scalesEstimator->SetTransformForward(true);
// However , We need not worry about the above considerations . because ScalesEstimator, The initial step size can also be estimated automatically , It can be done at each iteration or only at the first iteration . In this case , We choose to estimate the learning rate at the beginning of the registration process .
  affineOptimizer->SetScalesEstimator(scalesEstimator);
  
  affineOptimizer->SetDoEstimateLearningRateOnce(true);
  affineOptimizer->SetDoEstimateLearningRateAtEachIteration(false);
 
  affineOptimizer->SetLowerLimit(0);
  affineOptimizer->SetUpperLimit(2);
  affineOptimizer->SetEpsilon(0.2);
  affineOptimizer->SetNumberOfIterations(200);
  affineOptimizer->SetMinimumConvergenceValue(1e-6);
  affineOptimizer->SetConvergenceWindowSize(5);
 
  CommandIterationUpdate::Pointer observer2 = CommandIterationUpdate::New();
  affineOptimizer->AddObserver(itk::IterationEvent(), observer2);
 
  constexpr unsigned int numberOfLevels2 = 2;// We run two registration levels , The second level runs at full resolution , Here we make the final adjustment to the output parameters .
 
  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);
 
  using AffineCommandType = RegistrationInterfaceCommand<ARegistrationType>;
  AffineCommandType::Pointer command2 = AffineCommandType::New();
  affineRegistration->AddObserver(itk::MultiResolutionIterationEvent(),
                                  command2);
  try
  {
// Last , We call Update() Trigger the registration process , And add the output transformation of the last stage to the composite transformation . This composite transformation will be regarded as the final transformation of the multi-level registration process , And will be used by the resampler to moving Image resampled to virtual Domain space ( without fixed Initial transformation of , Then for fixed Image space ).
    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());
 
  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;
 
  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();
 
  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);
 
  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();
  }
 
  resample->SetTransform(compositeTransform);
  if (argc > 6)
  {
    writer->SetFileName(argv[6]);
    writer->Update();
  }
 
  return EXIT_SUCCESS;
}

2、itk::CompositeTransform Function is used to fuse translation transform and affine transform ;

The registration function uses v4, Can be multi-resolution settings ;

 

Multistage registration , The output of the first stage is directly linked to the second stage , The whole registration process only calls after the final stage Update() Trigger once .

https://github.com/InsightSoftwareConsortium/ITK/blob/master/Examples/RegistrationITKv4/MultiStageImageRegistration2.cxx

#include "itkImageRegistrationMethodv4.h"
 
#include "itkMattesMutualInformationImageToImageMetricv4.h"
 
#include "itkRegularStepGradientDescentOptimizerv4.h"
#include "itkConjugateGradientLineSearchOptimizerv4.h"
 
#include "itkTranslationTransform.h"
#include "itkAffineTransform.h"
#include "itkCompositeTransform.h"
 
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
 
#include "itkImageMomentsCalculator.h"
#include "itkResampleImageFilter.h"
#include "itkCastImageFilter.h"
#include "itkCheckerBoardImageFilter.h"
 
#include "itkCommand.h"
 
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;
 
  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);
    if (registration == nullptr)
    {
      itkExceptionMacro(<< "Dynamic cast failed, object of type "
                        << object->GetNameOfClass());
    }
 
    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;
  }
};
 
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 (optimizer == nullptr)
    {
      return;
    }
    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>;
 
  using TTransformType = itk::TranslationTransform<double, Dimension>;// Define the different types of the first phase .
  using TOptimizerType = itk::RegularStepGradientDescentOptimizerv4<double>;
  using MetricType =
    itk::MattesMutualInformationImageToImageMetricv4<FixedImageType,
                                                     MovingImageType>;
  using TRegistrationType =
    itk::ImageRegistrationMethodv4<FixedImageType, MovingImageType>;// In multi-stage registration , The type definition is the same as the previous example , But there is an important subtle change : The conversion type is no longer passed to the registration method as a template parameter . under these circumstances , The registration filter will convert the itk:: transform As the type of output conversion .
 
  TOptimizerType::Pointer    transOptimizer = TOptimizerType::New();
  MetricType::Pointer        transMetric = MetricType::New();
  TRegistrationType::Pointer transRegistration = TRegistrationType::New();
 
  transRegistration->SetOptimizer(transOptimizer);
  transRegistration->SetMetric(transMetric);
 
  TTransformType::Pointer translationTx = TTransformType::New();// We didn't deliver transform type , Instead, it is created outside the registration filter transform Explicit instantiation of objects , And use SetInitialTransform() Method to connect it to the registration object . in addition , By calling InPlaceOn() Method , This transformation object will become the output transformation of the registration filter , Or it will be grafted to the output .
 
  transRegistration->SetInitialTransform(translationTx);
  transRegistration->InPlaceOn();
 
  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");
 
  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);
 
  transMetric->SetNumberOfHistogramBins(24);
 
  if (argc > 7)
  {
    transMetric->SetNumberOfHistogramBins(std::stoi(argv[7]));
  }
 
  transOptimizer->SetNumberOfIterations(200);
  transOptimizer->SetRelaxationFactor(0.5);
  transOptimizer->SetLearningRate(16);
  transOptimizer->SetMinimumStepLength(1.5);
 
  CommandIterationUpdate::Pointer observer1 = CommandIterationUpdate::New();
  transOptimizer->AddObserver(itk::IterationEvent(), observer1);
 
  using TranslationCommandType =
    RegistrationInterfaceCommand<TRegistrationType>;
  TranslationCommandType::Pointer command1 = TranslationCommandType::New();
  transRegistration->AddObserver(itk::MultiResolutionIterationEvent(),
                                 command1);
// At the coarse resolution level , The first phase runs with only one registration level . however , Please note that , In this step, we do not need to update the translation registration filter , Because the output of this phase will be directly connected to the initial input of the next phase . because ITK Pipe structure , When we call at the last stage Update() when , The first phase will also be updated . Now we upgrade to an affine transformation as the second stage of the registration process , Same as before , We initially define and instantiate the different components of the current registration phase . We used a new optimizer , But the same metrics are used in the new configuration .
// Again, pay attention to ,TransformType There is no type definition passed to the registration filter . This is important , Because when the registration filter will convert the base class itk::transform As the type of its output transform , It prevents type mismatches when two phases are cascaded .
  using ATransformType = itk::AffineTransform<double, Dimension>;
  using AOptimizerType =
    itk::ConjugateGradientLineSearchOptimizerv4Template<double>;
  using ARegistrationType =
    itk::ImageRegistrationMethodv4<FixedImageType, MovingImageType>;
 
  AOptimizerType::Pointer    affineOptimizer = AOptimizerType::New();
  MetricType::Pointer        affineMetric = MetricType::New();
  ARegistrationType::Pointer affineRegistration = ARegistrationType::New();
 
  affineRegistration->SetOptimizer(affineOptimizer);
  affineRegistration->SetMetric(affineMetric);
 
  affineMetric->SetNumberOfHistogramBins(24);
  if (argc > 7)
  {
    affineMetric->SetNumberOfHistogramBins(std::stoi(argv[7]));
  }
 
  fixedImageReader->Update();
  FixedImageType::Pointer fixedImage = fixedImageReader->GetOutput();
 
  using FixedImageCalculatorType =
    itk::ImageMomentsCalculator<FixedImageType>;
 
  FixedImageCalculatorType::Pointer fixedCalculator =
    FixedImageCalculatorType::New();// then , Use their New() Method to instantiate all components , And connect to the registration object in the conversion type . In addition to the previous examples , Here we use fixed The centroid of the image to initialize the affine transformation fixed Parameters .ImageMomentsCalculator Filters are used for this purpose .
  fixedCalculator->SetImage(fixedImage);
  fixedCalculator->Compute();
 
  FixedImageCalculatorType::VectorType fixedCenter =
    fixedCalculator->GetCenterOfGravity();
 
  ATransformType::Pointer affineTx = ATransformType::New();// then , We are Affine Initialization in transformation fixed Parameters （ Center of rotation ）, And connect it to the registration object .
 
  const unsigned int numberOfFixedParameters =
    affineTx->GetFixedParameters().Size();
  ATransformType::ParametersType fixedParameters(numberOfFixedParameters);
  for (unsigned int i = 0; i < numberOfFixedParameters; ++i)
  {
    fixedParameters[i] = fixedCenter[i];
  }
  affineTx->SetFixedParameters(fixedParameters);
 
  affineRegistration->SetInitialTransform(affineTx);
  affineRegistration->InPlaceOn();
 
  affineRegistration->SetFixedImage(fixedImageReader->GetOutput());
  affineRegistration->SetMovingImage(movingImageReader->GetOutput());
  affineRegistration->SetObjectName("AffineRegistration");
 
  affineRegistration->SetMovingInitialTransformInput(
    transRegistration->GetTransformOutput());// Now? , The output of the first stage passes itk::DataObjectDecorator For packaging , And pass SetMovingInitialTransformInput() Method as a moving The input passed from the initial transformation to the second stage . Be careful , This API In another initialization method SetMovingInitialTransform() The name of is appended with a “Input” word , This method has been used in the previous example . This extension means the following API Data object decorator type required .
 
  using ScalesEstimatorType =
    itk::RegistrationParameterScalesFromPhysicalShift<MetricType>;
  ScalesEstimatorType::Pointer scalesEstimator = ScalesEstimatorType::New();
  scalesEstimator->SetMetric(affineMetric);
  scalesEstimator->SetTransformForward(true);
 
  affineOptimizer->SetScalesEstimator(scalesEstimator);
  affineOptimizer->SetDoEstimateLearningRateOnce(true);
  affineOptimizer->SetDoEstimateLearningRateAtEachIteration(false);
  affineOptimizer->SetLowerLimit(0);
  affineOptimizer->SetUpperLimit(2);
  affineOptimizer->SetEpsilon(0.2);
  affineOptimizer->SetNumberOfIterations(200);
  affineOptimizer->SetMinimumConvergenceValue(1e-6);
  affineOptimizer->SetConvergenceWindowSize(10);
 
  CommandIterationUpdate::Pointer observer2 = CommandIterationUpdate::New();
  affineOptimizer->AddObserver(itk::IterationEvent(), observer2);
 
  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);
 
  using AffineCommandType = RegistrationInterfaceCommand<ARegistrationType>;
  AffineCommandType::Pointer command2 = AffineCommandType::New();
  affineRegistration->AddObserver(itk::MultiResolutionIterationEvent(),
                                  command2);
  try
  {
// Once all registration components are ready , Last , We call Update（） To trigger the entire registration process , It consists of two cascaded registration stages .
    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;
  }
 
  using CompositeTransformType = itk::CompositeTransform<double, Dimension>;// Last , Composite transformations are used to link the results of all phases together , This will be seen as the final output of the multistage process , Will be passed to resample resample image moving To virtual Domain space ( without fxied Initial fixed Image space transformation ).
  CompositeTransformType::Pointer compositeTransform =
    CompositeTransformType::New();
  compositeTransform->AddTransform(translationTx);
  compositeTransform->AddTransform(affineTx);
 
  std::cout << " Translation transform parameters after registration: "
            << std::endl
            << transOptimizer->GetCurrentPosition() << std::endl
            << " Last LearningRate: "
            << transOptimizer->GetCurrentStepLength() << std::endl;
 
  std::cout << " Affine transform parameters after registration: "
            << std::endl
            << affineOptimizer->GetCurrentPosition() << std::endl
            << " Last LearningRate: " << affineOptimizer->GetLearningRate()
            << std::endl;
 
  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();
 
  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);
 
  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();
  }
 
  resample->SetTransform(compositeTransform);
  if (argc > 6)
  {
    writer->SetFileName(argv[6]);
    writer->Update();
  }
 
  return EXIT_SUCCESS;
}