ITK 多阶段配准

博客地址：https://blog.csdn.net/fanre/article/details/109118273

1、先平移配准粗配准，再仿射变换进行多分辨率配准

#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 "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>;//首先，我们需要配置初始阶段的配准组件。transform类型的实例化仅需要空间的维数和用于表示空间坐标的类型。
  
  using TOptimizerType = itk::RegularStepGradientDescentOptimizerv4<double>;//这里定义了其他配准组件的类型。使用itk::RegularStepGradientDescentOptimizerv4作为第一阶段的优化器。
  using MetricType =
    itk::MattesMutualInformationImageToImageMetricv4<FixedImageType,
                                                     MovingImageType>;//此外，我们使用itk::MattesMutualInformationImageToImageMetricv4作为度量，因为它适合多模式配准。
  using TRegistrationType = itk::ImageRegistrationMethodv4<FixedImageType,
                                                           MovingImageType,
                                                           TTransformType>;
//所有组件都使用它们的New（）方法实例化并连接到配准对象，
  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();//配准过程的输出transform将在配准过滤器内部构造，因为相关的TransformType已经作为模板参数传递给注册方法。但是，如果需要，我们应该为配准方法提供一个初始移动变换。
 
  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);//设置初始参数后，可以通过SetMovingInitialTransform（）方法将初始变换传递到配准过滤器。
 
  using CompositeTransformType = itk::CompositeTransform<double, Dimension>;//我们可以使用itk::CompositeTransform来堆叠来自多个阶段的所有输出transforms。这个复合转换还应该包含移动的初始转换(如果存在的话)，因为正如3.6.1节中解释的那样，每个配准阶段的输出不包括到该阶段的输入初始转换。
  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;//在此简单示例的情况下，第一阶段仅以粗略分辨率在一个配准级别中运行。
 
  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);
  try
  {
    transRegistration->Update();//一旦所有配准组件就位，我们通过调用Update()来触发配准过程，并将结果输出transform添加到最终的复合变换，因此该复合变换可用于初始化下一个配准阶段。
    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>;//现在我们可以升级到一个仿射变换作为配准过程的第二阶段。仿射变换是将直线映射为直线的线性变换。它可以用来表示平移、旋转、各向异性缩放、剪切或它们的任何组合。关于仿射变换的细节可以在3.9.16节中看到。转换类型的实例化只需要空间的维度和用于表示空间坐标的类型。
  
  using AOptimizerType =
    itk::ConjugateGradientLineSearchOptimizerv4Template<double>;//第二阶段的配置中使用了不同的优化程序，同时该度量标准保持与以前相同。
  using ARegistrationType = itk::ImageRegistrationMethodv4<FixedImageType,
                                                           MovingImageType,
                                                           ATransformType>;
//同样，所有组件都使用它们的New（）方法实例化并连接到配准对象，就像之前的阶段一样。
  AOptimizerType::Pointer    affineOptimizer = AOptimizerType::New();
  MetricType::Pointer        affineMetric = MetricType::New();
  ARegistrationType::Pointer affineRegistration = ARegistrationType::New();
 
  affineRegistration->SetOptimizer(affineOptimizer);
  affineRegistration->SetMetric(affineMetric);
//可以使用配准的初始变换和前一阶段的结果变换来初始化当前阶段，以便将两者都合并为复合变换。
  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();//用户必须在配准方法之外设置配准变换的fixed参数，因此首先我们实例化输出transform类型的对象。
 
  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;//在仿射变换中，旋转中心由fixed参数集定义，默认设置为[0,0]。但是，考虑这样一种情况，即运行配准的virtual空间的原点远离零原点。在这种情况下，将旋转中心作为默认值会使优化过程不稳定。因此，我们一直热衷于将旋转中心设置在virtual空间的中心，virtual空间通常是fixed的图像空间。
//注意，任何一个重心或几何中心都可以用作旋转中心。在这个例子中，旋转中心被设置为fixed图像的几何中心。我们也可以使用itk::ImageMomentsCalculator 过滤器来计算质心。
//然后，我们计算fixed图像的物理中心，并将其设置为输出仿射变换的中心。
 
  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);//然后，应使用SetInitialTransform（）方法将初始化后的输出transform连接到配准对象。重要的是要区分SetInitialTransform()和SetMovingInitialTransform()，后者用于根据前几个阶段的结果初始化配准阶段。您可以假设第一个用于直接操作当前配准过程中的可优化变换。
 
  using ScalesEstimatorType =
    itk::RegistrationParameterScalesFromPhysicalShift<MetricType>;
  ScalesEstimatorType::Pointer scalesEstimator = ScalesEstimatorType::New();//选择physical滤波器以估计参数scales。scales estimator然后将传递给optimizer。
  scalesEstimator->SetMetric(affineMetric);
  scalesEstimator->SetTransformForward(true);
//然而，我们不需要担心上述考虑。由于ScalesEstimator，初始步长也可以自动估计，可以在每次迭代或仅在第一次迭代。在这个例子中，我们选择在注册过程开始时估计一次学习率。
  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;//我们运行两个配准级别，其中第二个级别以全分辨率运行，在此我们对输出参数进行最终调整。
 
  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
  {
//最后，我们通过调用Update()触发配准过程，并将最后阶段的输出转换添加到复合转换中。此复合变换将被视为此多级配准过程的最终变换，并将被重采样器用于将moving图像重采样到virtual域空间(如果没有fixed的初始变换，则为fixed的图像空间)。
    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函数用来将平移变换和仿射变换两个变换进行融合；

其中注册函数用的是v4，都可以进行多分辨率设置；

 

多阶段配准，第一阶段的输出直接链接到第二阶段，整个配准过程只在最后阶段之后通过调用Update()触发一次。

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>;//定义第一阶段的不同类型。
  using TOptimizerType = itk::RegularStepGradientDescentOptimizerv4<double>;
  using MetricType =
    itk::MattesMutualInformationImageToImageMetricv4<FixedImageType,
                                                     MovingImageType>;
  using TRegistrationType =
    itk::ImageRegistrationMethodv4<FixedImageType, MovingImageType>;//多阶配准中，类型定义与前面的示例相同，但有一个重要的细微变化:转换类型不再作为模板参数传递给配准方法。在这种情况下，配准过滤器将把转换基类itk:: transform作为输出转换的类型。
 
  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();//我们没有传递transform类型，而是在配准过滤器外部创建transform对象的显式实例化，并使用SetInitialTransform()方法将其连接到配准对象。另外，通过调用InPlaceOn()方法，这个转换对象将成为配准过滤器的输出转换，或者将被嫁接到输出。
 
  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);
//在粗分辨率级别上，第一个阶段仅使用一个配准级别运行。但是，请注意，在这一步我们不需要更新平移配准过滤器，因为这一阶段的输出将直接连接到下一阶段的初始输入。由于ITK的管道结构，当我们在最后阶段调用Update()时，第一个阶段也将被更新。现在我们升级到一个仿射变换作为配准过程的第二阶段，和之前一样，我们最初定义并实例化当前配准阶段的不同组件。我们使用了新的优化器，但在新的配置中使用了相同的度量。
//再次注意，TransformType没有传递给配准过滤器的类型定义。这一点很重要，因为当配准过滤器将转换基类itk::transform作为其输出变换的类型时，它可以防止在两个阶段级联时出现类型不匹配。
  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();//然后，使用它们的New()方法实例化所有组件，并在转换类型中连接到配准对象。除了前面的例子，这里我们使用fixed图像的质心来初始化仿射变换的fixed参数。ImageMomentsCalculator过滤器用于此目的。
  fixedCalculator->SetImage(fixedImage);
  fixedCalculator->Compute();
 
  FixedImageCalculatorType::VectorType fixedCenter =
    fixedCalculator->GetCenterOfGravity();
 
  ATransformType::Pointer affineTx = ATransformType::New();//然后，我们在Affine变换中初始化fixed参数（旋转中心），并将其连接到配准对象。
 
  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());//现在，第一阶段的输出通过itk::DataObjectDecorator进行包装，并通过SetMovingInitialTransformInput()方法作为moving初始变换传递给第二阶段的输入。注意，这个API在另一个初始化方法SetMovingInitialTransform()的名称后面附加了一个“Input”字，该方法在前面的示例中已经使用过了。此扩展意味着以下API需要数据对象装饰器类型。
 
  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
  {
//一旦所有配准组件都准备就绪，最后，我们通过在最后一个阶段的配准过滤器上调用Update（）来触发整个配准过程，包括两个级联的配准阶段。
    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>;//最后,复合变换用于所有阶段的结果连接在一起,这将被视为最终的输出的多级过程,将传递给重新取样重新取样图像moving到virtual域空间(如果没有fxied的初始fixed图像空间变换)。
  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;
}