itk registration 16

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

3.16 Model Based Registration

本节介绍了几何模型与图像配准的概念。我们将这一概念称为基于模型的配准，但这可能不是最广泛的术语。该方法首先建立几何模型，并在模型中识别出若干参数。这些参数的变化使模型能够适应特定病人的形态。配准的任务是找到模型参数的最佳组合，使该模型能很好地表示图像中所包含的解剖结构。

例如，我们假设在大脑图像的axial轴向视图中，我们可以用一个椭圆来近似颅骨。椭圆成为我们的简化几何模型，配准的任务是找到椭圆的最佳中心，测量它的轴线长度和它在平面上的方向。如图3.57所示。如果我们将此方法与图像到图像的配准问题进行比较，我们可以看出这里的主要区别在于，除了映射模型的空间位置外，我们还可以定制内部参数来改变其形状。

图3.57:Model-to-Image配准的基本概念。通过应用空间变换和修改模型的内部参数，将简化的几何模型(椭圆)配准解剖结构(颅骨)。该图像并不是实际配准的结果，这里展示的目的只是为了说明模型对图像配准的概念。

图3.56说明了在配置基于模型的配准问题时ITK中配准框架的主要组件。配准的基本输入数据由itk::Image中的像素数据和存储在itk::SpatialObject中的几何数据提供。为了评估model和image之间的适合度，必须定义一个度量。这个适应度值可以通过在SpatialObject的空间定位中引入变化和/或通过改变其内部参数来提高。优化器的搜索空间现在是transform参数和shape内部参数的组合。

图3.56:基于模型配准的基本组成部分是image、spatial object、transform、metric、interpolator和optimizer。

同样的方法可以被认为是一种分割技术，因为一旦模型被最优地叠加在图像上，我们就可以根据它们与模型特定部分的关联来标记像素。模型在图像配准/分割中的应用是无穷无尽的。这种方法的主要优点可能是，与图像到图像的配准不同，它实际上提供了对图像中包含的解剖结构的Insight。经过修改的模型变成了对解剖结构的基本元素的本质表现。

ITK提供了一个类的层次结构，旨在支持形状模型的构造。这个层次结构以SpatialObject作为它的基类。在这个级别上定义了许多基本功能，包括评估给定点是在模型inside还是outside的能力，通过创建基本形状的层次组合形成复杂形状，以及支持基本空间参数化，如scale、orientation和position。

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

这个例子演示了如何使用itk::SpatialObject作为配准框架的一个组件，以便执行基于模型的配准。当前示例创建了一个由几个椭圆组成的几何模型。然后，利用该模型生成椭圆的合成二值图像。然后，它引入了对模型位置和形状的扰动，最后它使用扰动版本作为配准问题的输入。定义了一个度量来评价几何模型与图像之间的适合度。 

 

#include "itkEllipseSpatialObject.h"//让我们先看看支持SpatialObject所需的类。在本例中，我们使用itk::EllipseSpatialObject作为基本形状组件，并使用itk::GroupSpatialObject将它们组合在一起，作为更复杂形状的表示。它们各自的标题包括在下面。
#include "itkGroupSpatialObject.h"

#include "itkSpatialObjectToImageFilter.h"//为了生成椭圆的初始合成图像，我们使用itk::SpatialObjectToImageFilter对图像中的每个像素进行测试（无论像素（以及空间对象）在几何模型的inside还是outside）。

#include "itkImageToSpatialObjectRegistrationMethod.h"

#include "itkImageToSpatialObjectMetric.h"//定义了一个度量来评估SpatialObject和Image之间的适合度。这种度量标准的基类是itk::ImageToSpatialObjectMetric，其标头包含在下面。

#include "itkLinearInterpolateImageFunction.h"//与先前的配准问题一样，我们必须评估非网格位置的图像强度。为此，使用了itk::LinearInterpolateImageFunction。

#include "itkEuler2DTransform.h"//通过使用itk::Transform将SpatialObject从其自己的空间映射到图像空间。在此示例中，我们使用itk::Euler2DTransform。

#include "itkOnePlusOneEvolutionaryOptimizer.h"//配准从根本上说是一个优化问题。这里我们包含了用于搜索参数空间和确定将形状模型映射到图像顶部的最佳变换的优化器。本例中使用的优化器是itk::OnePlusOneEvolutionaryOptimizer，它实现了一个进化算法。

#include "itkDiscreteGaussianImageFilter.h"
#include "itkNormalVariateGenerator.h"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkCastImageFilter.h"
#include "itkRescaleIntensityImageFilter.h"

#include "itkCommand.h"//在前面的配准示例中，跟踪优化器在参数空间中的发展非常重要。这是通过使用Command/Observer范式完成的。下面的代码实现了监视配准过程的itk::Command observer。该代码与我们在之前的注册示例中使用的代码非常相似。
//该命令将在优化程序的每次迭代时调用，并将打印出当前转换参数的组合。

template <class TOptimizer>
class IterationCallback : public itk::Command
{
public:
  using Self = IterationCallback;
  using Superclass = itk::Command;
  using Pointer = itk::SmartPointer<Self>;
  using ConstPointer = itk::SmartPointer<const Self>;

  itkTypeMacro(IterationCallback, Superclass);
  itkNewMacro(Self);

  /** Type defining the optimizer. */
  using OptimizerType = TOptimizer;

  /** Method to specify the optimizer. */
  void
  SetOptimizer(OptimizerType * optimizer)
  {
    m_Optimizer = optimizer;
    m_Optimizer->AddObserver(itk::IterationEvent(), this);
  }

  /** Execute method will print data at each iteration */
  void
  Execute(itk::Object * caller, const itk::EventObject & event) override
  {
    Execute((const itk::Object *)caller, event);
  }

  void
  Execute(const itk::Object *, const itk::EventObject & event) override
  {
    if (typeid(event) == typeid(itk::StartEvent))
    {
      std::cout << std::endl << "Position              Value";
      std::cout << std::endl << std::endl;
    }
    else if (typeid(event) == typeid(itk::IterationEvent))
    {
      std::cout << m_Optimizer->GetCurrentIteration() << "   ";
      std::cout << m_Optimizer->GetValue() << "   ";
      std::cout << m_Optimizer->GetCurrentPosition() << std::endl;
    }
    else if (typeid(event) == typeid(itk::EndEvent))
    {
      std::cout << std::endl << std::endl;
      std::cout << "After " << m_Optimizer->GetCurrentIteration();
      std::cout << "  iterations " << std::endl;
      std::cout << "Solution is    = " << m_Optimizer->GetCurrentPosition();
      std::cout << std::endl;
    }
  }

protected:
  IterationCallback() = default;
  itk::WeakPointer<OptimizerType> m_Optimizer;
};

template <typename TFixedImage, typename TMovingSpatialObject>
class SimpleImageToSpatialObjectMetric
  : public itk::ImageToSpatialObjectMetric<TFixedImage, TMovingSpatialObject>
{//现在考虑一下这种新配准方法最关键的组成部分:metric。这个组件评估SpatialObject和图像之间的匹配。该指标的smoothness和regularity决定了分配给优化器的任务的难度。在这种情况下，我们使用一个非常robust的优化器，应该能够找到它的方式，即使在最不连续的成本函数。要实现的度量应该派生自ImageToSpatialObjectMetric类。
//下面的代码实现了一个简单的度量，它计算空间对象内部像素的总和。实际上，当模型与图像对齐时，度量值最大。度量标准模板化在空间对象的类型和图像的类型上。
public:
  /** Standard class type aliases. */
  using Self = SimpleImageToSpatialObjectMetric;
  using Superclass =
    itk::ImageToSpatialObjectMetric<TFixedImage, TMovingSpatialObject>;
  using Pointer = itk::SmartPointer<Self>;
  using ConstPointer = itk::SmartPointer<const Self>;

  using PointType = itk::Point<double, 2>;
  using PointListType = std::list<PointType>;
  using MovingSpatialObjectType = TMovingSpatialObject;
  using ParametersType = typename Superclass::ParametersType;
  using DerivativeType = typename Superclass::DerivativeType;
  using MeasureType = typename Superclass::MeasureType;

  /** Method for creation through the object factory. */
  itkNewMacro(Self);

  /** Run-time type information (and related methods). */
  itkTypeMacro(SimpleImageToSpatialObjectMetric, ImageToSpatialObjectMetric);

  static constexpr unsigned int ParametricSpaceDimension = 3;

  /** Specify the moving spatial object. */
  void
  SetMovingSpatialObject(const MovingSpatialObjectType * object) override
  {
    if (!this->m_FixedImage)
    {
      std::cout << "Please set the image before the moving spatial object"
                << std::endl;
      return;
    }
    this->m_MovingSpatialObject = object;
    m_PointList.clear();
    using myIteratorType =
      itk::ImageRegionConstIteratorWithIndex<TFixedImage>;

    myIteratorType it(this->m_FixedImage,
                      this->m_FixedImage->GetBufferedRegion());

    itk::Point<double, 2> point;

    while (!it.IsAtEnd())
    {
      this->m_FixedImage->TransformIndexToPhysicalPoint(it.GetIndex(), point);

      if (this->m_MovingSpatialObject->IsInsideInWorldSpace(point, 99999))
      {
        m_PointList.push_back(point);
      }
      ++it;
    }

    std::cout << "Number of points in the metric = "
              << static_cast<unsigned long>(m_PointList.size()) << std::endl;
  }

  /** Get the Derivatives of the Match Measure */
  void
  GetDerivative(const ParametersType &, DerivativeType &) const override
  {
    return;
  }
  
  /** Get the value for SingleValue optimizers. */
  MeasureType  GetValue(const ParametersType & parameters) const override
  {
//度量的基本操作是它的GetValue()方法。适应度值就是在这个方法中计算出来的。在我们当前的例子中，适应度是在SpatialObject的点上计算的。对于每个点，通过变换将其坐标映射到图像空间中。生成的点用于对图像进行评估，生成的值以总和的形式累积。由于我们不允许缩放变化，当所有空间对象点被映射到图像的白色区域时，总和的最优值就会产生。注意，GetValue()方法的参数是变换的参数数组。
    double value;
    this->m_Transform->SetParameters(parameters);

    value = 0;
    for (auto it : m_PointList)
    {
      PointType transformedPoint = this->m_Transform->TransformPoint(it);
      if (this->m_Interpolator->IsInsideBuffer(transformedPoint))
      {
        value += this->m_Interpolator->Evaluate(transformedPoint);
      }
    }
    return value;
  }

  /** Get Value and Derivatives for MultipleValuedOptimizers */
  void
  GetValueAndDerivative(const ParametersType & parameters,
                        MeasureType &          Value,
                        DerivativeType &       Derivative) const override
  {
    Value = this->GetValue(parameters);
    this->GetDerivative(parameters, Derivative);
  }

private:
  PointListType m_PointList;
};

int
main(int argc, char * argv[])
{
  if (argc > 1)
  {
    std::cerr << "Too many parameters " << std::endl;
    std::cerr << "Usage: " << argv[0] << std::endl;
  }
//定义完所有配准组件后，我们准备将各个部分放在一起并执行配准过程。
//首先，我们实例化GroupSpatialObject和EllipseSpatialObject。这两个对象由空间的尺寸参数化。在我们当前的示例中，创建了2D实例化。
  using GroupType = itk::GroupSpatialObject<2>;
  using EllipseType = itk::EllipseSpatialObject<2>;
  
  using ImageType = itk::Image<float, 2>;//下面几行使用像素类型和空间维度实例化图像。此图像使用浮动像素类型，因为我们计划平滑它，以增加优化器的捕获半径。真实像素类型的图像比整数像素类型的图像在平滑条件下表现更好。

  EllipseType::Pointer ellipse1 = EllipseType::New();//使用它们的New（）方法创建EllipseSpatialObjects，并将结果分配给SmartPointers。这些线将创建三个椭圆。
  EllipseType::Pointer ellipse2 = EllipseType::New();
  EllipseType::Pointer ellipse3 = EllipseType::New();

  ellipse1->SetRadiusInObjectSpace(10.0);//每个从SpatialObject派生的类都有特定的参数，使用户可以定制其形状。对于EllipseSpatialObject, SetRadius()用于定义椭圆的大小。另外一个SetRadius(Array)方法允许用户独立地定义椭圆轴。
  ellipse2->SetRadiusInObjectSpace(10.0);
  ellipse3->SetRadiusInObjectSpace(10.0);
 
  // Place each ellipse at the right position to form a triangle

  EllipseType::TransformType::OffsetType offset;//默认情况下，椭圆在空间中居中。我们使用以下代码行将椭圆排列成三角形。与对象内在关联的空间转换由GetTransform()方法访问。该转换可以使用SetOffset()方法在空间中定义转换。我们利用这个特性将椭圆放置在空间的特定点上。

  offset[0] = 100.0;
  offset[1] = 40.0;

  ellipse1->GetModifiableObjectToParentTransform()->SetOffset(offset);
  ellipse1->Update();

  offset[0] = 40.0;
  offset[1] = 150.0;
  ellipse2->GetModifiableObjectToParentTransform()->SetOffset(offset);
  ellipse2->Update();

  offset[0] = 150.0;
  offset[1] = 150.0;
  ellipse3->GetModifiableObjectToParentTransform()->SetOffset(offset);
  ellipse3->Update();
  //注意，在转换中进行更改之后，SpatialObject将调用ComputeGlobalTransform()方法来更新其全局变换。这样做的原因是SpatialObjects可以按照层次结构进行排列。然后就可以通过moving组的父对象来改变一组空间对象的位置。

//现在我们将三个EllipseSpatialObjects添加到一个GroupSpatialObject中，该对象随后将被传递给配准方法。GroupSpatialObject可以方便地将三个椭圆作为表示复杂形状的高级结构进行管理。组可以嵌套任意数量的级别，以便以更高的细节表示形状。
  GroupType::Pointer group = GroupType::New();
  group->AddChild(ellipse1);
  group->AddChild(ellipse2);
  group->AddChild(ellipse3);
  
  using SpatialObjectToImageFilterType =
    itk::SpatialObjectToImageFilter<GroupType, ImageType>;//准备好几何模型后，我们继续生成代表椭圆占据的空间印记的二值图像。SpatialObjectToImageFilter用于此目的。请注意，此过滤器是在所使用的空间对象和要生成的图像类型上实例化的。
  
  SpatialObjectToImageFilterType::Pointer imageFilter =
    SpatialObjectToImageFilterType::New();//对于定义的类型，我们使用New（）方法构造一个过滤器。新创建的过滤器已分配给SmartPointer。
 
  imageFilter->SetInput(group);//GroupSpatialObject作为输入传递到过滤器。
 
  ImageType::SizeType size;//itk::SpatialObjectToImageFilter充当重采样过滤器。因此，它要求用户定义所需输出图像的尺寸。这是通过SetSize（）方法指定的。
  size[0] = 200;
  size[1] = 200;
  imageFilter->SetSize(size);
 
  imageFilter->Update();//最后，我们通过调用Update（）方法来触发过滤器的执行。
 
  using GaussianFilterType =
    itk::DiscreteGaussianImageFilter<ImageType, ImageType>;//为了获得更平滑的度量，我们使用itk::DiscreteGaussianImageFilter对图像进行平滑处理。这扩展了度量的捕获半径，并产生了一个更连续的成本函数来优化。下面几行实例化高斯过滤器，并使用New()方法创建此类型的一个对象。
  GaussianFilterType::Pointer gaussianFilter = GaussianFilterType::New();

  gaussianFilter->SetInput(imageFilter->GetOutput());//SpatialObjectToImageFilter的输出作为输入连接到DiscreteGaussianImageFilter。
  
  constexpr double variance = 20;//为了增加捕获半径，将滤波器的方差定义为较大的值。最后，使用Update（）方法触发过滤器的执行。
  gaussianFilter->SetVariance(variance);
  gaussianFilter->Update();
 
  using RegistrationType =
    itk::ImageToSpatialObjectRegistrationMethod<ImageType, GroupType>;//下面，我们实例化itk::ImageToSpatialObjectRegistrationMethod方法的类型，并用New()方法实例化一个配准对象。注意，配准类型模板化在图像和SpatialObject类型上。在本例中，空间对象是一组空间对象。
  RegistrationType::Pointer registration = RegistrationType::New();
  
  using MetricType = SimpleImageToSpatialObjectMetric<ImageType, GroupType>;//我们实例化在图像类型和空间对象类型上模板化的度量。与往常一样，使用New（）方法创建一个对象。
  MetricType::Pointer metric = MetricType::New();
  
  using InterpolatorType =
    itk::LinearInterpolateImageFunction<ImageType, double>;//需要一个插值器来评估非网格位置的图像。在这里，我们实例化了线性内插器类型。
  InterpolatorType::Pointer interpolator = InterpolatorType::New();
 
  using OptimizerType = itk::OnePlusOneEvolutionaryOptimizer;//实例化了evolutionary optimizer。
  OptimizerType::Pointer optimizer = OptimizerType::New();

  using TransformType = itk::Euler2DTransform<>;//接下来，我们实例化转换类。在这种情况下，我们使用在2D空间中实现刚性变换的Euler2DTransform。
  TransformType::Pointer transform = TransformType::New();
 
  itk::Statistics::NormalVariateGenerator::Pointer generator =
    itk::Statistics::NormalVariateGenerator::New();//Evolutionary algorithms是基于测试参数的随机变化。为了支持随机值的计算，ITK提供了一组随机数生成器。在本例中，我们使用itk::NormalVariateGenerator生成具有正态分布的值。
  
  generator->Initialize(12345);//随机数生成器必须使用种子初始化。
 
  optimizer->SetNormalVariateGenerator(generator);//通过指定随机数生成器，初始总体的样本数和最大迭代数，可以初始化OnePlusOneEvolutionaryOptimizer。
  optimizer->Initialize(10);
  optimizer->SetMaximumIteration(400);
  
  TransformType::ParametersType parametersScale;//在前面的配准示例中，我们注意了对不同变换参数的动态范围进行规格化。特别地，我们必须对Euler2DTransform的角度和平移的范围进行补偿。为了实现这一目标，我们向优化器提供了一组scales。
  parametersScale.set_size(3);
  parametersScale[0] = 1000; // angle scale

  for (unsigned int i = 1; i < 3; i++)
  {
    parametersScale[i] = 2; // offset scale
  }
  optimizer->SetScales(parametersScale);
 
  using IterationCallbackType = IterationCallback<OptimizerType>;//这里我们实例化了Command对象，它将充当配准方法的观察者，并在每次迭代时打印出参数。在前面，我们将此命令定义为在优化器类型上模板化的类。一旦用New()方法创建了它，我们就将优化器连接到命令。
  IterationCallbackType::Pointer callback = IterationCallbackType::New();
  callback->SetOptimizer(optimizer);
  
  registration->SetFixedImage(gaussianFilter->GetOutput());//所有组件都插入到ImageToSpatialObjectRegistrationMethod对象中。这里使用的是典型的Set()方法。注意使用SetMovingSpatialObject()方法连接空间对象。我们提供原始合成二值图像的模糊版本作为输入图像。
  registration->SetMovingSpatialObject(group);
  registration->SetTransform(transform);
  registration->SetInterpolator(interpolator);
  registration->SetOptimizer(optimizer);
  registration->SetMetric(metric);
  
  TransformType::ParametersType initialParameters(
    transform->GetNumberOfParameters());//使用SetInitialTransformParameters()方法将转换参数的初始集传递给配准方法。注意，由于我们的原始模型已经与合成图像进行了配准，我们引入了一个人为的错误配准，以便使优化初始化在远离最优值的某个点上。

  initialParameters[0] = 0.2; // Angle
  initialParameters[1] = 7.0; // Offset X
  initialParameters[2] = 6.0; // Offset Y
  registration->SetInitialTransformParameters(initialParameters);

  std::cout << "Initial Parameters  : " << initialParameters << std::endl;
  optimizer->MaximizeOn();//由于用于评估空间对象和图像之间适合度的度量的特性，我们必须告诉优化器，我们感兴趣的是找到度量的最大值。有些指标将低数值与良好匹配关联在一起，而另一些指标将高数值与良好匹配关联在一起。MaximizeOn()和MaximizeOff()方法允许用户处理这两种类型的指标。
  
  try//最后，我们使用Update（）方法触发注册过程的执行。我们将这个调用放在try / catch块中，以防在此过程中引发任何异常。
  {
    registration->Update();
    std::cout << "Optimizer stop condition: "
              << registration->GetOptimizer()->GetStopConditionDescription()
              << std::endl;
  }
  catch (const itk::ExceptionObject & exp)
  {
    std::cerr << "Exception caught ! " << std::endl;
    std::cerr << exp << std::endl;
  }
 
  RegistrationType::ParametersType finalParameters =
    registration->GetLastTransformParameters();//配准产生的转换参数集可以用GetLastTransformParameters()方法恢复。此方法返回转换参数数组，应该根据每个转换的实现对其进行解释。在我们当前的例子中，Euler2DTransform有三个参数:旋转角度、x上的平移量和y上的平移量。

  std::cout << "Final Solution is : " << finalParameters << std::endl;
  return EXIT_SUCCESS;
}