itk registration 7

 参考：https://itk.org/ITKSoftwareGuide/html/Book2/ITKSoftwareGuide-Book2ch3.html#x26-1740003.18

3.7 Multi-Resolution Registration

使用多分辨率方法进行图像配准被广泛用于提高速度、精度和鲁棒性。其基本思想是，首先在图像具有较少像素的粗尺度上进行配准。然后，在粗级别上确定的空间映射用于在下一个更细级别上初始化注册。这个过程不断重复，直到达到可能的最佳规模。这种由粗到细的策略极大地提高了配准成功率，并且通过消除粗尺度上的局部最优值而增强了鲁棒性。通过在粗尺度上平滑图像，鲁棒性可以得到进一步提高。

在前面的所有示例中，我们以单一分辨率运行配准过程。然而，ITKv4配准框架的结构提供了一种多分辨率的配准方法。为此，我们只需要定义层数以及每层输入图像的分辨率和平滑度。配准滤波器根据用户定义的ShrinkFactor和SmoothingSigma向量对图像进行平滑和子采样。

3.7.1 Fundamentals

在ITK中，可以使用ITK::MultiResolutionPyramidImageFilter从输入图像创建一个降低分辨率的图像序列。下采样是根据用户定义的多分辨率schedule进行的。schedule被指定为整数itk::Array2D，包含每个维度(列)的每个多分辨率级别(行)的收缩因子。例如

8 4 4（不同列对应不同的维度[x,y,z]，不同行表示不同采样等级）  
4 4 2

是一个调度的三维图像为两个多分辨率水平。在第一个(最粗糙的)级别，图像在column列维上减少了8倍，在row行维上减少了4倍，在slice切片维上减少了4倍。在第二层中，图像在column列维上减少了4倍，在row行维上减少了4倍，在slice切片维上减少了2倍。

方法SetNumberOfLevels()用于设置金字塔中的分辨率级别的数量。该方法将为调度分配内存，并生成一个默认表，其中所有维度的起始(最粗)收缩因子设置为2(M-1)，其中M是级别的数量。在随后的所有级别中，所有因素都被减半。例如，如果我们将级别数设置为4，那么默认的计划是：

8 8 8  （不同列对应不同的维度[x,y,z]，不同行表示不同采样等级）  
4 4 4  
2 2 2  
1 1 1

用户可以使用方法GetSchedule()获取调度的副本，进行修改，并使用方法SetSchedule()重置它。或者，用户可以使用方法SetStartingShrinkFactors()指定起始(最粗)收缩因子来创建默认表。后续级别的因子是通过将因子减半或设置为1来生成的，具体取决于哪个更大。例如，对于4层金字塔和8、8和4的起始因子，生成的schedule将是：

8 8 4  （不同列对应不同的维度[x,y,z]，不同行表示不同采样等级）  
4 4 2  
2 2 1  
1 1 1

当这个过滤器被Update()触发时，会产生M个输出，其中第M个输出对应于金字塔的第M层。为了生成这些图像，首先使用itk::DiscreteGaussianImageFilter进行高斯平滑，其方差设置为(s/2)^2，其中s为收缩因子。然后使用itk::ShrinkImageFilter对平滑后的图像进行次采样。

3.7.2 Fundamentals

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

这与以前的ITK版本相比有很大的优势，因为在ITKv3中，我们必须使用不同的过滤器(ITK::MultiResolutionImageRegistrationMethod)来运行多分辨率进程。此外，我们必须使用图像金字塔过滤器(itk::MultiResolutionPyramidImageFilter)来创建下降采样图像的序列。

图3.33 多分辨率注册过程的概念表示。

ImageRegistrationMethodv4以一种由粗到细的方式解决了注册问题，如图3.33所示。首先使用固定和移动图像金字塔的第一层图像在最粗的层次上进行配准。由配准决定的变换参数，然后使用金字塔的第二层的图像在下一个更精细的层次上初始化配准。这个过程会重复进行，直到我们达到图像分辨率的最佳水平。

在典型的注册场景中，用户将调整组件设置，甚至在多分辨率级别之间交换组件。例如，当优化在一个粗分辨率，它可能采取更激进的步长和有一个更宽松的收敛标准。

可以使用ITK的Command/Observer设计模式实现在不同分辨率级别之间调整组件。在开始在每个分辨率级别配准之前，ImageRegistrationMethodv4调用MultiResolutionIterationEvent()。可以通过实现响应事件的itk::Command来更改配准组件。在前面的3.4节中对events和commands之间的交互进行了简要描述。

#include "itkImageRegistrationMethodv4.h"//头文件
#include "itkTranslationTransform.h"
#include "itkMattesMutualInformationImageToImageMetricv4.h"
#include "itkRegularStepGradientDescentOptimizerv4.h"

#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"

#include "itkResampleImageFilter.h"
#include "itkCastImageFilter.h"
#include "itkCheckerBoardImageFilter.h"

#include "itkCommand.h"#包含Command

template <typename TRegistration>
class RegistrationInterfaceCommand : public itk::Command
{ //我们新的接口command类称为RegistrationInterfaceCommand。它源自Command，并在多分辨率配准类型上进行模板化。
public:
//类似于CommandIterationUpdate类的方式定义Self，Superclass，Pointer，New（）
  using Self = RegistrationInterfaceCommand;
  using Superclass = itk::Command;
  using Pointer = itk::SmartPointer<Self>;
  itkNewMacro(Self);

protected:
  RegistrationInterfaceCommand() = default;//构造函数
  
  using RegistrationType = TRegistration;//为了方便起见，我们在Execute（）方法中声明对转换指针有用的类型。
  using RegistrationPointer = RegistrationType *;
  using OptimizerType = itk::RegularStepGradientDescentOptimizerv4<double>;
  using OptimizerPointer = OptimizerType *;
  
  void
  Execute(itk::Object * object, const itk::EventObject & event) override
  {//有两个参数传递给Execute（）方法：第一个是指向调用事件的对象的指针，第二个是调用的事件。
   
    if (!(itk::MultiResolutionIterationEvent().CheckEvent(&event)))//我们验证调用的事件是正确的类型，即itk :: MultiResolutionIterationEvent（）。如果没有，我们将不采取任何进一步行动而返回。
    {
      return;
    }
    
    auto registration = static_cast<RegistrationPointer>(object);//然后我们将输入object指针转换为RegistrationPointer。注意，这里没有进行错误检查来验证dynamic_cast是否成功，因为我们知道实际的对象是一个配准方法。然后我们从配准方法中请求优化器对象。
    auto optimizer =
      static_cast<OptimizerPointer>(registration->GetModifiableOptimizer());

    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 level = " << currentLevel << std::endl;
    std::cout << "    shrink factor = " << shrinkFactors << std::endl;
    std::cout << "    smoothing sigma = ";
    std::cout << smoothingSigmas[currentLevel] << std::endl;
    std::cout << std::endl;
   
    if (registration->GetCurrentLevel() == 0)//如果这是第一个分辨率级别，我们将学习率(表示第一步大小)和最小步长(表示收敛准则)设置为较大的值。在随后的每个分辨率级别上，我们将把最小步长减少到原来的5倍，以便优化器能够逐步关注更小的区域。学习率设置为当前步长。这样，当优化器在下一级别的注册过程开始时被重新初始化时，步骤长度将从上一级别使用的最后一个值开始。这将保证优化器通过参数空间所采用的路径的连续性。
    {
      optimizer->SetLearningRate(16.00);
      optimizer->SetMinimumStepLength(2.5);
    }
    else
    {
      optimizer->SetLearningRate(optimizer->GetCurrentStepLength());
      optimizer->SetMinimumStepLength(optimizer->GetMinimumStepLength() *
                                      0.2);
    }
  }
 
  void
  Execute(const itk::Object *, const itk::EventObject &) override
  {//还需要另一个接受const输入对象​​的Execute（）方法版本，因为此方法在基类中定义为纯虚函数。此版本只是返回而无需采取任何措施。
    return;
  }
};

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::RegularStepGradientDescentOptimizerv4<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() << "   ";
    std::cout << m_CumulativeIterationIndex++ << std::endl;
  }

private:
  unsigned int m_CumulativeIterationIndex{ 0 };
};


int
main(int argc, const 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;

  const std::string fixedImageFile = argv[1];
  const std::string movingImageFile = argv[2];
  const std::string outImagefile = argv[3];
  const PixelType backgroundGrayLevel = (argc > 4) ? std::stoi(argv[4]) : 100;
  const std::string checkerBoardBefore = (argc > 5) ? argv[5] : "";
  const std::string checkerBoardAfter = (argc > 6) ? argv[6] : "";
  const int         numberOfBins = (argc > 7) ? std::stoi(argv[7]) : 0;

  using FixedImageType = itk::Image<PixelType, Dimension>;
  using MovingImageType = itk::Image<PixelType, Dimension>;
  using TransformType = itk::TranslationTransform<double, Dimension>;

  using OptimizerType = itk::RegularStepGradientDescentOptimizerv4<double>;

  using MetricType =
    itk::MattesMutualInformationImageToImageMetricv4<FixedImageType,
                                                     MovingImageType>;
  using RegistrationType = itk::
    ImageRegistrationMethodv4<FixedImageType, MovingImageType, TransformType>;

  TransformType::Pointer    transform = TransformType::New();
  OptimizerType::Pointer    optimizer = OptimizerType::New();
  MetricType::Pointer       metric = MetricType::New();
  RegistrationType::Pointer registration = RegistrationType::New();

  registration->SetOptimizer(optimizer);
  registration->SetMetric(metric);

  using FixedImageReaderType = itk::ImageFileReader<FixedImageType>;
  using MovingImageReaderType = itk::ImageFileReader<MovingImageType>;

  FixedImageReaderType::Pointer fixedImageReader =
    FixedImageReaderType::New();
  MovingImageReaderType::Pointer movingImageReader =
    MovingImageReaderType::New();

  fixedImageReader->SetFileName(fixedImageFile);
  movingImageReader->SetFileName(movingImageFile);

  registration->SetFixedImage(fixedImageReader->GetOutput());
  registration->SetMovingImage(movingImageReader->GetOutput());


  using ParametersType = OptimizerType::ParametersType;
  ParametersType initialParameters(transform->GetNumberOfParameters());

  initialParameters[0] = 0.0; // Initial offset in mm along X
  initialParameters[1] = 0.0; // Initial offset in mm along Y

  transform->SetParameters(initialParameters);

  registration->SetInitialTransform(transform);
  registration->InPlaceOn();

  metric->SetNumberOfHistogramBins(24);

  if (argc > 7)
  {

    metric->SetNumberOfHistogramBins(numberOfBins);
  }

  optimizer->SetNumberOfIterations(200);//要设置优化器参数，请注意，LearningRate和MinimumStepLength是在每个分辨率级别开始时在obsever中设置的。
  optimizer->SetRelaxationFactor(0.5);

  CommandIterationUpdate::Pointer observer = CommandIterationUpdate::New();
  optimizer->AddObserver(itk::IterationEvent(), observer);
//我们将多分辨率级别设置为三个，并为每个分辨率级别设置相应的shrink factor和smoothing sigma value。在低分辨率的下采样图像中使用平滑可以避免metric函数的大波动，从而防止优化器陷入局部极小值。在这个简单的例子中，我们没有平滑，我们在前两个分辨率级别中使用了小收缩。
  constexpr unsigned int numberOfLevels = 3;

  RegistrationType::ShrinkFactorsArrayType shrinkFactorsPerLevel;
  shrinkFactorsPerLevel.SetSize(3);
  shrinkFactorsPerLevel[0] = 3;
  shrinkFactorsPerLevel[1] = 2;
  shrinkFactorsPerLevel[2] = 1;

  RegistrationType::SmoothingSigmasArrayType smoothingSigmasPerLevel;
  smoothingSigmasPerLevel.SetSize(3);
  smoothingSigmasPerLevel[0] = 0;
  smoothingSigmasPerLevel[1] = 0;
  smoothingSigmasPerLevel[2] = 0;
//fixed和moving图像类型如前面的示例中定义。根据ShrinkFactor和SmoothingSigma向量提供的值，通过配准方法在内部创建了不同分辨率级别的下采样图像。
  registration->SetNumberOfLevels(numberOfLevels);
  registration->SetShrinkFactorsPerLevel(shrinkFactorsPerLevel);
  registration->SetSmoothingSigmasPerLevel(smoothingSigmasPerLevel);

  using CommandType = RegistrationInterfaceCommand<RegistrationType>;//一旦所有配准组件都就位，我们就可以创建接口命令的实例，并使用AddObserver（）方法将其连接到配准对象。
  CommandType::Pointer command = CommandType::New();

  registration->AddObserver(itk::MultiResolutionIterationEvent(), command);

  try
  {
    registration->Update();//通过调用Update（）触发配准过程。
    std::cout << "Optimizer stop condition: "
              << registration->GetOptimizer()->GetStopConditionDescription()
              << std::endl;
  }
  catch (const itk::ExceptionObject & err)
  {
    std::cout << "ExceptionObject caught !" << std::endl;
    std::cout << err << std::endl;
    return EXIT_FAILURE;
  }

  ParametersType finalParameters = transform->GetParameters();

  double TranslationAlongX = finalParameters[0];
  double TranslationAlongY = finalParameters[1];

  unsigned int numberOfIterations = optimizer->GetCurrentIteration();

  double bestValue = optimizer->GetValue();

  std::cout << "Result = " << std::endl;
  std::cout << " Translation X = " << TranslationAlongX << std::endl;
  std::cout << " Translation Y = " << TranslationAlongY << std::endl;
  std::cout << " Iterations    = " << numberOfIterations << std::endl;
  std::cout << " Metric value  = " << bestValue << std::endl;

  using ResampleFilterType =
    itk::ResampleImageFilter<MovingImageType, FixedImageType>;

  ResampleFilterType::Pointer resample = ResampleFilterType::New();

  resample->SetTransform(transform);
  resample->SetInput(movingImageReader->GetOutput());

  FixedImageType::Pointer fixedImage = fixedImageReader->GetOutput();


  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(outImagefile);


  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);

  TransformType::Pointer identityTransform = TransformType::New();
  identityTransform->SetIdentity();
  resample->SetTransform(identityTransform);

  for (int q = 0; q < argc; ++q)
  {
    std::cout << q << " " << argv[q] << std::endl;
  }
  if (checkerBoardBefore != std::string(""))
  {
    writer->SetFileName(checkerBoardBefore);
    writer->Update();
  }


  resample->SetTransform(transform);
  if (checkerBoardAfter != std::string(""))
  {
    writer->SetFileName(checkerBoardAfter);
    writer->Update();
  }

  return EXIT_SUCCESS;
}