【ITK库学习】使用itk库进行图像配准：“Hello World”配准（一）

基本组成：参考图像、待配准图像、传递函数、路劲选择、校对机、优化器
传递函数：从参考图像上的点到待配准图像上点的空间映射关系
校对机：评估待配准图象在非网格位置的程度
路劲选择：提供一种参考图像呗待配准图像配准的程度
优化器：路径选择形成数量上的标准，优化器通过寻找被传递参数定义的空间去最大优化此标准

1、itkImageRegistrationMethod / itkImageRegistrationMethodv4

该类为图像配准方法的基类。

该类定义配准方法的通用接口。它根据要配准的两个图像的类型进行模板化，使用通用 Transform，允许在运行时选择用于配准图像的特定转换类型。该方法使用通用度量来比较两个图像，配准方法的最终目标是找到优化度量的转换参数集。

此类中使用术语：固定图像和移动图像来指示变换正在映射什么图像。

此类使用固定图像的坐标系作为参考，并搜索将点从固定图像的空间映射到移动图像的空间变换。

为此，将连续应用度量来将固定图像与变换后的移动图像进行比较，此过程还需要从移动图像中插入值。

itkImageRegistrationMethodv是基于传统的优化方法，如梯度下降法等，算法复杂度相对简单一些，更容易上手。支持可以在运行时选择通用优化器，优化器的唯一限制是它应能在单值成本函数中运行，因为用于比较图像的指标提供单个值作为输出。

itkImageRegistrationMethodv4是基于图像类的优化方法，使用图像类提供的优化接口进行配准，算法更复杂，如多分辨率策略、局部变形场模型等，可以处理更复杂的配准问题，参数设置和使用稍微复杂一些，需要一定的图像处理和算法知识。允许多阶段配准，其中每个阶段的特征在于可能不同的变换和不同的图像度量。例如，许多执行线性配准，然后执行变形配准，其中两个阶段都在多个级别中执行。

对算法和参数设置要求较高，可以选择itkImageRegistrationMethodv4，如果只需要简单的图像配准，可以选择itkImageRegistrationMethodv。
共同的成员函数：

SetFixedImage()：设置/获取固定图像（参考图像）
Set/GetMovingImage()：设置/获取移动图像（待配准图像）
SetOptimizer()/GetModifiableOptimizer()：设置/获取优化器
SetMetric()/GetModifiableMetric()：设置/获取度量(校对机)

itkImageRegistrationMethodv常用的成员函数：

SetTransform()/GetModifiableTransform：设置/获取变换（传递函数）
Set/GetFixedImageRegion()：将固定图像的区域设置/获取为在配准期间被视为感兴趣区域，该区域将被传递给ImageMetric，以限制度量计算仅考虑该区域
SetFixedImageRegionDefined()：打开/关闭固定图像区域的使用，ImageMetric将限制其计算
GetFixedImageRegionDefined()：如已为固定图像定义了ImageMetric将限制其计算的区域，则为 True
Set/GetInitialTransformParameters()：设置/获取初始转换参数
SetInterpolator()GetModifiableInterpolator()：设置/获取插值器
Set/GetLastTransformParameters()：为派生类提供设置此私有变量的能力
StartOptimization()：初始化
GetMTime()：返回该对象或其任何缓存的 ivar 的最新修改时间的方法

itkImageRegistrationMethodv4常用的成员函数：

Set/GetOptimizerWeights()：设置/获取优化器权重，允许设置每个局部参数的加权数组，如果未设置，权重将被视为同一性，权重用于在优化期间屏蔽特定参数以保持其恒定，或者它们可以用来应用另一种先验知识，权重的大小必须等于局部变换参数的数量
Set/GetFixedPointSet()：设置/获取固定点集
Set/GetFixedInitialTransformInput()：设置/获取初始固定变换
SetMovingPointSet()：设置移动点集
Set/GetMovingInitialTransformInput()：设置/获取初始移动变换
Set/GetNumberOfLevels()：设置/获取多分辨率级别的数量，在设置级别数时，需要为每个级别设置以下内容：a）虚拟域的收缩因子，b）Sigma平滑参数，c）使用指定级别的特定参数转换适配器
Set/GetMetricSamplingStrategy()：设置/获取指标采样策略
SetMetricSamplingPercentagePerLevel()：设置指标采样百分比，有效值范围[0.0,1.0]
SetMetricSamplePoints()：设置度量采样点
SetInPlace()：请求将 InitialTransform 移植到输出上，而不创建副本
Set/GetInitialTransformInput()：设置/获取要优化的初始变换，该变换与MovingInitialTransform 一起指定从运动图像到虚拟图像的初始变换，它用于默认参数，可用于指定变换类型，如果过滤器设置了 “InPlace” ，则此变换将是输出变换对象或“嫁接”到输出，否则，此 InitialTransform 将被深度复制或“克隆”到输出，如果未设置此参数，则使用默认构造的输出变换
SetInitializeCenterOfLinearOutputTransform()：使用队列中前一个变换的中心初始化要优化的当前线性变换，这提供了比默认原点更好的初始化
Set/GetTransformParametersAdaptorsPerLevel()：设置/获取转换适配器
Set/GetSmoothingSigmasPerLevel()：设置/获取每个级别的平滑Sigma值，在每个分辨率级别，都会应用高斯平滑滤波器（ itkDiscreteGaussianImageFilter ），Sigma值根据选项 m_SmoothingSigmasAreSpecifiedInPhysicalUnits 指定
Set/GetSmoothingSigmasAreSpecifiedInPhysicalUnits()：设置/获取是否以物理单位（默认）或体素指定每个级别的平滑Sigma
SetShrinkFactorsPerLevel()：设置每个级别的收缩因子，其中每个级别的每个维度都有一个恒定的收缩因子，例如，输入到 Factors = [4,2,1] 函数将在第一级将每个维度的图像缩小 4，然后在第二级缩小 2，最后一级的1原始分辨率（使用 itkShrinkImageFilter ）
SetShrinkFactorsPerDimension()：为每个维度设置特定级别的收缩因子

2、itkTranslationTransform

向量空间的平移变换（例如空间坐标）

使用映射变换可以获得相同的功能，但性能差异很大。

常用的成员函数：

Set/GetParameters()：设置/获取指定的变换值的参数
SetIdentity()：将参数设置为 IdentityTransform
SetFixedParameters()：设置固定参数并更新内部变换，Translation Transform 不需要固定参数，因此该方法的实现是空操作
Set/GetOffset()：设置平移变换的偏移量，此方法将 TranslationTransform 的偏移量设置为用户指定的值
IsLinear()：表明该变换是线性的，也就是说，给定两个点 P 和 Q，以及标量系数 a 和 b，则T(aP+bQ)=aT§+bT(Q)
GetInverseTransform()：返回该变换的逆变换

3、itkMeanSquaresImageToImageMetric / itkMeanSquaresImageToImageMetric4

这两个类用于计算图像间均方差的度量类。

指标越小，表示两幅图像的相似度越高。

使用该指标需要调用SetFixedImage()和SetMovingImage()函数设置待比较的两幅图像，并可以使用SetFixedImageRegion()函数进一步指定比较区域，然后调用Initialize()函数进行初始化，然后可以通过调用GetValue()函数计算出指标值。

itkMeanSquaresImageToImageMetric是ITK中较早的版本的度量类，它使用的是传统的像素级别的计算方法。它使用的计算公式为：sum((pixel1 - pixel2)^2) / N，其中pixel1和pixel2分别表示两个图像中的像素值，N表示图像中的像素总数。这种计算方法简单直观，但是对于图像配准等应用可能不够灵活，无法适应复杂的数据类型、图像的变换等。

itkMeanSquaresImageToImageMetricv4在ITK在新版本中被引入，它使用的是基于像素间距离场的计算方法。它使用的计算公式为：sum(weight * (distance1 - distance2)^2) / N，其中weight是像素间的权重值，distance1和distance2分别表示两个图像中的像素距离值，N表示图像中的像素总数。这种计算方法可以更灵活地处理图像间的变换和重采样，并且可以用于处理多通道和向量图像等复杂数据。

常用的成员函数：

GetValueThreadProcessSample()：获取匹配测量值
GetValueAndDerivative()：获取单值优化器的值和导数
GetValue()：获取值
GetDerivative()：获取匹配度量的导数

4、itkRegularStepGradientDescentOptimizerv / itkRegularStepGradientDescentOptimizerv4

二者均为梯度下降优化器，ITK新旧版本的不同。

itkRegularStepGradientDescentOptimizerv基于固定步长的梯度下降算法，它有固定的步长参数，需要手动设置。

常用的成员函数：

StepAlongGradient()：考虑由因子表示的步长，沿着校正的梯度前进一步，该方法由 AdvanceOneStep 调用，它预计会被非向量空间中的优化方法覆盖

itkRegularStepGradientDescentOptimizerv4基于线搜索的梯度下降算法，它的步长在每次迭代时根据目标函数的变化情况进行动态调整，在许多情况下也具有更好的优化性能和更快的收敛速度。为防止采取太大的步骤，在每次迭代中，该优化器将沿着度量导数的方向采取一步，每次导数的方向突然改变时，优化器都会假设已通过局部极值，并通过将步长减小默认设置为 0.5 的松弛因子来做出反应，初始步长的默认值为 1，并且该值只能通过 SetLearningRate() 手动更改，因为此优化器不使用 ScaleEstimator来自动估计学习率，另请注意，与itkRegularStepGradientDescentOptimizerv不同，它没有“最大化/最小化”选项来修改度量导数的效果，假定指定的度量返回“改进”优化的参数导数结果。

常用的成员函数：

Set/GetRelaxationFactor()：设置/获取松弛因子值
Set/GetMinimumStepLength()：设置/获取用于收敛检查的最小步长（学习率）值，当通过大步越过局部极小值时，将通过松弛因子调整（减小）步长，以便朝着最小值点（收敛）采取更小的步，当步长值达到较小值时，视为收敛，默认值设置为 1e-4 以通过所有测试
Set/GetGradientMagnitudeTolerance()：设置/获取梯度幅度容差值以进行收敛检查
Set/GetCurrentLearningRateRelaxation()：设置/获取学习率的当前比例
EstimateLearningRate()：根据当前梯度估计学习率
AdvanceOneStep()：沿着梯度方向前进一步，包括变换更新
StartOptimization()：启动并运行优化

示例代码：

#include "itkImageRegistrationMethodv4.h"
#include "itkTranslationTransform.h"
#include "itkMeanSquaresImageToImageMetricv4.h"
#include "itkRegularStepGradientDescentOptimizerv4.h"
#include "itkImage.h"using namespace itk;const unsigned int  Dim2D = 2;       //数据的Dimension
typedef itk::Image<float, Dim2D> FixedImageType;
typedef itk::Image<float, Dim2D> MovingImageType;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*;voidExecute(itk::Object* caller, const itk::EventObject& event) override{Execute((const itk::Object*)caller, event);}voidExecute(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::endl;}
};bool imageRegistration1(FixedImageType* fixImage, MovingImageType* moveImage, bool useEstimator)
{typedef itk::LinearInterpolateImageFunction<FixedImageType, double> FixedLinearInterpolatorType;typedef itk::LinearInterpolateImageFunction<MovingImageType, double>  MovingLinearInterpolatorType;typedef itk::TranslationTransform<double, Dim2D> TransformType;typedef itk::MeanSquaresImageToImageMetricv4<FixedImageType, MovingImageType> MetricType;typedef itk::RegularStepGradientDescentOptimizerv4<double> OptimizerType;typedef itk::ImageRegistrationMethodv4<FixedImageType, MovingImageType, TransformType> RegistrationType;typename TransformType::Pointer movingInitialTransform = TransformType::New();typename TransformType::ParametersType initialParameters(movingInitialTransform->GetNumberOfParameters());initialParameters[0] = 0.0;   // Initial offset in mm along XinitialParameters[1] = 0.0;   // Initial offset in mm along Y//initialParameters[2] = 0.0;   // Initial offset in mm along ZmovingInitialTransform->SetParameters(initialParameters);typename TransformType::Pointer identityTransform = TransformType::New();identityTransform->SetIdentity();typename FixedLinearInterpolatorType::Pointer fixedInterpolator = FixedLinearInterpolatorType::New();typename MovingLinearInterpolatorType::Pointer movingInterpolator = MovingLinearInterpolatorType::New();typename MetricType::Pointer metric = MetricType::New();metric->SetFixedInterpolator(fixedInterpolator);metric->SetMovingInterpolator(movingInterpolator);typename OptimizerType::Pointer optimizer = OptimizerType::New();optimizer->SetLearningRate(4);optimizer->SetMinimumStepLength(0.001);optimizer->SetRelaxationFactor(0.5);if (useEstimator){using ScalesEstimatorType =itk::RegistrationParameterScalesFromPhysicalShift<MetricType>;auto scalesEstimator = ScalesEstimatorType::New();scalesEstimator->SetMetric(metric);scalesEstimator->SetTransformForward(true);optimizer->SetScalesEstimator(scalesEstimator);optimizer->SetDoEstimateLearningRateOnce(true);}optimizer->SetNumberOfIterations(200);auto observer = CommandIterationUpdate::New();optimizer->AddObserver(itk::IterationEvent(), observer);constexpr unsigned int numberOfLevels = 1;typename RegistrationType::Pointer registration = RegistrationType::New();typename RegistrationType::ShrinkFactorsArrayType shrinkFactorsPerLevel;shrinkFactorsPerLevel.SetSize(1);shrinkFactorsPerLevel[0] = 1;typename RegistrationType::SmoothingSigmasArrayType smoothingSigmasPerLevel;smoothingSigmasPerLevel.SetSize(1);smoothingSigmasPerLevel[0] = 0;registration->SetFixedImage(fixImage);              //参考图像registration->SetMovingImage(moveImage);            //待配准图像registration->SetMetric(metric);                    //校对机registration->SetOptimizer(optimizer);              //优化器registration->SetMovingInitialTransform(movingInitialTransform);    //初始移动变换registration->SetFixedInitialTransform(identityTransform);     //初始固定变换registration->SetNumberOfLevels(numberOfLevels);    //多分辨率级别的数量 registration->SetSmoothingSigmasPerLevel(smoothingSigmasPerLevel);   //每个级别的平滑Sigma值registration->SetShrinkFactorsPerLevel(shrinkFactorsPerLevel);       //每个级别的收缩因子try{registration->Update();}catch (itk::ExceptionObject& ex){//读取过程发生错误std::cerr << "Error: " << ex << std::endl;return false;}TransformType::ConstPointer transform = registration->GetTransform();TransformType::ParametersType finalParameters = transform->GetParameters();const double  TranslationAlongX = finalParameters[0];const double  TranslationAlongY = finalParameters[1];const unsigned int numberOfIterations = optimizer->GetCurrentIteration();    //抵达收敛的迭代的实际次数const double bestValue = optimizer->GetValue();       //参数集合的图像量规值return true;
}