【ITK库学习】使用itk库进行图像配准:“Hello World”配准(一)

2023-12-20 18:50:35

图像配准

将一幅图上的点映射到另一幅图像上同源点的空间转换过程;是同意吗目标的两幅图想在空间位置上的对准。

左图点P经过传递函数T,得到在右图的配准点位Q:

配准框架

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

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,其中pixel1pixel2分别表示两个图像中的像素值,N表示图像中的像素总数。这种计算方法简单直观,但是对于图像配准等应用可能不够灵活,无法适应复杂的数据类型、图像的变换等。

itkMeanSquaresImageToImageMetricv4在ITK在新版本中被引入,它使用的是基于像素间距离场的计算方法。它使用的计算公式为:sum(weight * (distance1 - distance2)^2) / N,其中weight是像素间的权重值,distance1distance2分别表示两个图像中的像素距离值,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*;

	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::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 X
	initialParameters[1] = 0.0;   // Initial offset in mm along Y
	//initialParameters[2] = 0.0;   // Initial offset in mm along Z
	movingInitialTransform->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;
}

注:调试过程中报错:itkStatisticsAlgorithm.hxx的533行报: ”InsertSort”: 找不到标识符,是4.13版本的InsertSort引用在函数定义前,修改文件内容,将InsertSort函数定义放在引用函数前bug解决

bug解决办法参考自博客:https://blog.csdn.net/assjaa/article/details/111317338

文章来源:https://blog.csdn.net/weixin_41804847/article/details/135090037
本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。