CN107871325B - 基于Log-Euclidean协方差矩阵描述符的图像非刚性配准方法 - Google Patents

Abstract

本发明公开了一种基于Log‑Euclidean协方差矩阵描述符的图像非刚性配准方法，将能够有效表示图像局部结构特性且对大的旋转、缩放、照度等变化具有不变性的特征描述符结合进经典Log‑Demons配准模型的目标函数中，构建了一种新配准算法模型。新模型将先验结构信息结合进配准过程，避免了形变驱动力仅依赖灰度差和梯度信息的弊端，有效减少了形变误差，获得了优于基于B样条的FFD模型、Log‑Demons模型等当前主流非刚性配准算法的配准精度和鲁棒性，增强了配准鲁棒性，对在临床中的实际应用有着重要的参考价值。

Description

基于Log-Euclidean协方差矩阵描述符的图像非刚性配准 方法

技术领域

本发明涉及医学图像处理及应用的技术领域，尤其是指一种基于Log-Euclidean协方差矩阵描述符的图像非刚性配准方法。

背景技术

由于不同时间采集到的肺部CT会受到呼吸运动、脏器蠕动、重力作用、体位变动等因素使软组织产生复杂形变，而为观察解剖结构的变化情况，满足临床病情诊断、康复治疗效果评估、治疗计划的确定或调整等，对图像进行非刚性配准具有重要的理论意义和实用价值。

在配准过程中，形变模型的选择是非常重要的，因为它反映了几何变换的类型，很大程度上限定了适用范围。当前，最常用的形变模型有Demons模型和自由形变(Free-FormDeformation,FFD)模型.FFD模型具有平滑的几何变换函数，通过对控制网格位置的移动来获得任何形变，但却不能够确保拓扑保持性。Log-Demons模型是Demons模型的一个变体，通过迭代地计算力矢量场以在适当的方向上驱动形变，然后通过卷积来平滑力矢量场以便更新形变。另外，它具有微分同胚的空间变换，能够保持图像的拓扑结构在配准前后不发生改变。与FFD模型相比，Demons模型具有较高的灵活性、配准精度和计算效率，更适合于医学图像的配准。

经典Demons配准算法的形变驱动力受限于图像梯度的局部范围,使其局限于小形变、弱鲁棒性、低精度等。针对这些局限，许多学者积极开展基于Demons模型的图像配准算法研究，提出很多改进的Demons算法。改进策略主要分为两种：一是在驱动力中添加新的信息项，克服了仅依赖灰度的梯度信息出现的驱动力不足问题，如Wang等提出了ActiveDemons算法；二是在配准能量表达式中添加新的信息项，旨在获得具有必要属性的形变场，如Mansi等提出的Log-Demons配准算法具有拓扑保持性的形变场，符合真实材料的形变特性，可以保持物体拓扑结构以避免不合理的物理形变的优势。但是，Log-Demons算法对具有灰度偏差的图像间配准呈现出弱鲁棒性，对大形变的配准表现出低配准精度的问题。

本发明提供一种基于Log-Euclidean协方差矩阵描述符的图像非刚性配准方法，将能够有效表示图像局部结构特性且对大的旋转、缩放、照度等变化具有不变性的特征描述符结合进经典Log-Demons配准模型的目标函数中，构建了一种新配准算法模型。新模型将先验结构信息结合进配准过程，避免了形变驱动力仅依赖灰度差和梯度信息的弊端，有效减少了形变误差，增强了配准鲁棒性。

发明内容

本发明的目的在于克服现有技术的缺点与不足，提出了一种基于Log-Euclidean协方差矩阵描述符的图像非刚性配准方法，获得了优于基于B样条的FFD模型、Log-Demons模型等当前主流非刚性配准算法的配准精度和鲁棒性，对在临床中的实际应用有着重要的参考价值。

为实现上述目的，本发明所提供的技术方案为：基于Log-Euclidean协方差矩阵描述符的图像非刚性配准方法，包括以下步骤：

1)对参考图像和浮动图像采用基于B样条的FFD配准算法进行刚性预配准，获得预配准后的参考图像和浮动图像；

2)采用向导图像滤波作用于预配准后的参考图像和浮动图像，获得滤波后的参考图像和浮动图像；

3)在向导图像滤波的基础上，采用多分辨率策略获得3层金字塔图像，初始金字塔层数为0；

4)对当前层金字塔参考图像和浮动图像的每一像素点提取如下7个特征：空间坐标x和y、归一化后灰度值、x和y方向一阶和二阶梯度的范数，将提取的特征分别映射为行数×列数×7维的特征空间F₁和F₂，再对F₁和F₂通过使用积分图像方法来快速计算参考图像和浮动图像的Log-Euclidean协方差矩阵描述符C_R(x,y)和

5)将Log-Euclidean协方差矩阵描述符C_R(x,y)和

取对数后的差作为新的匹配项添加到Log-Demons算法的目标能量函数中，构建新目标能量函数；

6)根据所提出新目标能量函数计算图像的变形驱动力，得到瞬时位移向量；

7)利用得到瞬时位移向量更新形变矩阵；

8)判断是否达到最大迭代次数，若没有达到，则返回步骤3)，继续迭代；否则判断是否达到最大分解级数，若没有达到，则对低分辨率下得到的形变矩阵进行超采样，将其作为上一级高分辨率的初始变换，返回步骤3)，否则配准结束，进入步骤9)；

9)将最终获得的形变矩阵通过双三次插值方法作用于刚性配准后的浮动图像，完成浮动图像的最终配准。

在步骤2)中，对预配准后的参考图像和浮动图像进行边缘和细节保存的平滑运算处理，采用的是向导图像滤波：

q_i＝a_kI_i+b_k,i∈ω_k

式中，q_i是输出图，I_i是向导图，i是像素索引，a_k和b_k是常数，ω_k是以向导图I_i中一点k为中心半径为r的局部窗。

在步骤3)中，在配准算法实现过程中采用了多分辨率策略，包括如下步骤：

3.1)通过重采样的方法对参考图像和浮动图像计算3层金字塔模型；

3.2)在低分辨率下利用较少的时间进行粗配准，得到本尺度的形变矩阵；

3.3)对在低分辨率下得到的形变矩阵进行超采样，将其作为上一级高分辨率的初始变换；

3.4)逐层的迭代，最终实现整幅图像的配准。

在步骤4)中，所述图像的特征描述符计算，包括以下步骤：

4.1)对每一像素提取的特征映射为行数×列数×7维的特征空间F：

F(x,y)＝[x,y,I(x,y),|I_x(x,y)|,|I_y(x,y)|,|I_xx(x,y)|,|I_yy(x,y)|]

式中，(x,y)是像素点位置，I(x,y)是(x,y)像素点的灰度值，|I_x(x,y)|和|I_y(x,y)|是像素点(x,y)的一阶梯度范数，|I_xx(x,y)|和|I_yy(x,y)|是像素点(x,y)的二阶梯度范数；

4.2)对特征空间F采用积分图像方法计算得到每一像素的Log-Euclidean协方差矩阵描述符C(x,y)：

式中，μ是C上中心点在(x,y)的r邻域内的均值，N_r是r邻域内点的总个数，x'和y'满足不等式|x-x'|≤r，|y-y'|≤r。

在步骤5)中，在Log-Demons算法基础上，采用了新目标能量函数：

式中，R和M分别表示参考图像和浮动图像的像素值，符号

表示复合运算符，exp和log分别表示指数函数和对数函数，ν是第n次迭代时未正则化的速度场，ν_c是第n-1次迭代后正则化的速度场，C_R和

分别表示参考图像和浮动图像Log-Euclidean协方差矩阵描述符，φ(p)表示当前估计的变换，符号▽表示梯度，α、β和γ分别是影响着图像噪声、光滑程度和正则化在目标能量函数中的权重，ζ为新添加匹配项的所占权重，p是像素位置。

在步骤6)中，对所提出新目标能量函数进行形变驱动力的计算，得到瞬时位移向量：

式中，R(p)和M(p)分别表示参考图像和浮动图像在位置p处的像素值，符号

表示M与函数φ(p)的复合运算，φ(p)表示当前估计的变换，α、β和γ分别是影响着图像噪声、光滑程度和正则化在目标能量函数中的权重，ζ为新添加匹配项的所占权重，log表示对数函数，C_R和

分别表示参考图像和浮动图像Log-Euclidean协方差矩阵描述符，J_p和K_p分别表示新目标能量函数第一项中函数

和第三项中函数

的Jacobian矩阵，其表达式分别为：

在步骤7)中，利用得到瞬时位移向量更新形变矩阵，包括以下步骤：

7.1)高斯滤波平滑瞬时位移向量：

式中，

为高斯滤波，*表示卷积；

7.2)更新形变矩阵：

式中，ν_n表示第n次迭代后正则化的速度场，符号

表示复合函数运算。

本发明与现有技术相比，具有如下优点与有益效果：

1、本发明的配准精度在经典Log-Demons算法基础上平均改善75％以上。

2、本发明的配准鲁棒性优于当前主流的非刚性配准算法如FFD算法和Log-Demons算法。

3、本发明的配准收敛速度胜于经典的Log-Demons算法。

附图说明

图1为本发明逻辑流程示意图。

图2为本发明由第1组人工合成进行试验得到的配准精度比较图。

图3为本发明由第2组人工合成进行试验得到的配准精度比较图。

图4为本发明由第1组肺CT图像(有肿瘤)进行试验得到的配准精度比较图。

图5为本发明由第2组肺CT图像(无肿瘤)进行试验得到的配准精度比较图。

图6a为两组人工合成图像配准的目标能量函数随迭代次数变化的曲线图。

图6b为三组肺CT图像配准的目标能量函数随迭代次数变化的曲线图。

具体实施方式

下面结合具体实施例对本发明作进一步说明。

如图1至图6b所示，本实施例所提供的基于Log-Euclidean协方差矩阵描述符的图像非刚性配准方法，使用的实验平台是Matlab7.11，计算机主机配置：Intel(R)处理器，CPU主频3.0GHz，内存4.0GB，其包括以下步骤：

1)对参考图像和浮动图像采用基于B样条的FFD配准算法进行刚性预配准，获得预配准后的参考图像和浮动图像，确保在非刚性配准过程中的参考图像和浮动图像是空间对齐的。

2)采用向导图像滤波分别作用于预配准后的参考图像和浮动图像，获得了边缘和细节结构保存的平滑后图像。向导图像滤波能够保存图像边缘和细节结构,同时对图像平坦区域进行近似Gaussian平滑处理,有效避免重采样过程细节丢失问题,增加特征描述符的可区分性。

3)在配准算法实现过程中采用了多分辨率策略，包括如下步骤：

3.1)通过重采样的方法对参考图像和浮动图像计算3层金字塔模型；

3.2)在低分辨率下利用较少的时间进行粗配准，得到本尺度的形变矩阵；

3.3)对在低分辨率下得到的形变矩阵进行超采样，将其作为上一级高分辨率的初始变换；

3.4)逐层的迭代，最终实现整幅图像的配准。

多分辨率策略中，从高分辨率到低分辨率的最大迭代次数分别为100次、50次和20次，目标能量函数值不小于1×10^-4为迭代停止准则，最大分解级数设为3，初始层从0开始。多分率策略限制了局部极值的产生,有助于提高配准鲁棒性,并降低计算代价。

4)对当前层金字塔参考图像和浮动图像的特征描述符计算，包括以下步骤：

4.1)对每一像素提取的特征映射为行数×列数×7维的特征空间F：

F(x,y)＝[x,y,I(x,y),|I_x(x,y)|,|I_y(x,y)|,|I_xx(x,y)|,|I_yy(x,y)|]

式中，(x,y)是像素点位置，I(x,y)是(x,y)像素点的灰度值，|I_x(x,y)|和|I_y(x,y)|是像素点(x,y)的一阶梯度范数，|I_xx(x,y)|和|I_yy(x,y)|是像素点(x,y)的二阶梯度范数；

4.2)对特征空间F采用积分图像方法计算得到每一像素的Log-Euclidean协方差矩阵描述符C(x,y)：

式中，μ是C上中心点在(x,y)的r邻域内的均值，N_r是r邻域内点的总个数，x'和y'满足不等式|x-x'|≤r，|y-y'|≤r。

所获得的特征描述符不仅能够精确显示解剖变化，而且对大的旋转、缩放、照度等变化具有不变性，有效地描述了图像局部特性。

5)在Log-Demons算法基础上，配准实现过程采用新目标能量函数：

式中，R和M分别表示参考图像和浮动图像的像素值，符号

表示复合运算符，exp和log分别表示指数函数和对数函数，ν是第n次迭代时未正则化的速度场，ν_c是第n-1次迭代后正则化的速度场，C_R和

分别表示参考图像和浮动图像的Log-Euclidean协方差矩阵描述符，ν是第n次迭代时未正则化的速度场，ν_c是第n-1次迭代后正则化的速度场，φ(p)表示当前估计的变换，符号▽表示梯度，α、β和γ分别是影响着图像噪声、光滑程度和正则化在目标能量函数中的权重，ζ为新添加匹配项的所占权重，p是像素位置。

新目标能量函数结合了一个微分同胚的平滑约束项

确保所提出的目标能量函数的可微性，进一步强化了目标能量函数的平滑性，有助于提高配准精度。

6)对所提出新目标能量函数进行形变驱动力的计算，得到瞬时位移向量：

式中，R(p)和M(p)分别表示参考图像和浮动图像在位置p处的像素值，符号

表示M与函数φ(p)的复合运算，φ(p)表示当前估计的变换，α、β和γ分别是影响着图像噪声、光滑程度和正则化在目标能量函数中的权重，ζ为新添加匹配项的所占权重，log表示对数函数，C_R和

分别表示参考图像和浮动图像Log-Euclidean协方差矩阵描述符，J_p和K_p分别表示新目标能量函数第一项中函数

和第三项中函数

的Jacobian矩阵，其表达式分别为：

所提出的形变驱动力避免了仅依靠灰度差和梯度信息的弊端，有效的减少了形变误差，有助于增强配准的鲁棒性。

为了分析配准鲁棒性和配准精度，本发明以32例术后随访检查的胸部CT为配准实验对象，配准结果定量评估分析如表1所示：

表1配准结果评估分析

评估方法	NMI	NCC	DSC	RMSE	RSSD
						配准前	1.2811±0.12	0.7696±0.09	0.8070±0.13	707.75±426	-
FFD算法	1.3495±0.15	0.9098±0.03	0.6085±0.19	495.46±211	0.6602±0.18
						Log-Demons	1.3749±0.14	0.9273±0.03	0.9523±0.02	403.90±280	0.3533±0.13
所提算法	1.5628±0.12	0.9945±0.006	0.9894±0.004	99.37±104	0.0405±0.03

从表中可以看出，本发明的配准后图像与参考图像的归一化互信息NMI均值最大且标准差最小，互相关系数NCC均值最大高达99％以上且标准差最小，Dice相似度度量DSC均值最大且高达98％以上且标准差最小，均方根误差RMSE均值最小且标准差最小，相对差方和RSSD均值最小且标准差最小，因而所提算法对大形变和灰度变化的配准鲁棒性最好。

7)利用得到瞬时位移向量更新形变矩阵，包括以下步骤：

7.1)高斯滤波平滑瞬时位移向量：

式中，

为高斯滤波，*表示卷积；

7.2)更新形变矩阵：

式中，ν_n表示第n次迭代后正则化的速度场，符号

表示复合函数运算。

8)判断是否达到最大迭代次数，若没有达到，则返回步骤3)，继续迭代；否则判断是否达到最大分解级数，若没有达到，则对低分辨率下得到的形变矩阵进行超采样，将其作为上一级高分辨率的初始变换，返回步骤3)，否则配准结束，进入步骤9)；

9)将最终获得的形变矩阵通过双三次插值方法作用于刚性配准后的浮动图像，完成浮动图像的最终配准。

综上所述，本发明将先验结构信息结合进配准过程，避免了形变驱动力仅依赖灰度差和梯度信息的弊端，有效减少了形变误差，增强了配准鲁棒性，对在临床中的实际应用有着重要的参考价值，值得推广。

以上所述实施例只为本发明之较佳实施例，并非以此限制本发明的实施范围，故凡依本发明之形状、原理所作的变化，均应涵盖在本发明的保护范围内。

Claims (4)

1.基于Log-Euclidean协方差矩阵描述符的图像非刚性配准方法，其特征在于，包括以下步骤：

1)对参考图像和浮动图像采用基于B样条的FFD配准算法进行刚性预配准，获得预配准后的参考图像和浮动图像；

2)采用向导图像滤波作用于预配准后的参考图像和浮动图像，获得滤波后的参考图像和浮动图像；

3)在向导图像滤波的基础上，采用多分辨率策略获得3层金字塔图像，初始金字塔层数为0；

其中，在配准算法实现过程中采用多分辨率策略，包括如下步骤：

3.1)通过重采样的方法对参考图像和浮动图像计算3层金字塔模型；

3.2)在低分辨率下进行粗配准，得到本尺度的形变矩阵；

3.3)对在低分辨率下得到的形变矩阵进行超采样，将其作为上一级高分辨率的初始变换；

3.4)逐层的迭代，最终实现整幅图像的配准；

4)对当前层金字塔参考图像和浮动图像的每一像素点提取如下7个特征：空间坐标x和y、归一化后灰度值、x和y方向一阶和二阶梯度的范数，将提取的特征分别映射为行数×列数×7维的特征空间F₁和F₂，再对F₁和F₂通过使用积分图像方法来快速计算参考图像和浮动图像的Log-Euclidean协方差矩阵描述符C_R(x,y)和

5)将Log-Euclidean协方差矩阵描述符C_R(x,y)和

取对数后的差作为新的匹配项添加到Log-Demons算法的目标能量函数中，构建新目标能量函数；

6)根据所提出新目标能量函数计算图像的变形驱动力，得到瞬时位移向量；

7)利用得到瞬时位移向量更新形变矩阵；

8)判断是否达到最大迭代次数，若没有达到，则返回步骤3)，继续迭代；否则判断是否达到最大分解级数，若没有达到，则对低分辨率下得到的形变矩阵进行超采样，将其作为上一级高分辨率的初始变换，返回步骤3)，否则配准结束，进入步骤9)；

9)将最终获得的形变矩阵通过双三次插值方法作用于刚性配准后的浮动图像，完成浮动图像的最终配准。

2.根据权利要求1所述的基于Log-Euclidean协方差矩阵描述符的图像非刚性配准方法，其特征在于：在步骤2)中，对预配准后的参考图像和浮动图像进行边缘和细节保存的平滑运算处理，采用的是向导图像滤波：

q_i＝a_kI_i+b_k,i∈ω_k

式中，q_i是输出图，I_i是向导图，i是像素索引，a_k和b_k是常数，ω_k是以向导图I_i中一点k为中心半径为r的局部窗。

3.根据权利要求1所述的基于Log-Euclidean协方差矩阵描述符的图像非刚性配准方法，其特征在于，在步骤4)中，所述图像的Log-Euclidean协方差矩阵描述符计算，包括以下步骤：

4.1)对每一像素提取的特征映射为行数×列数×7维的特征空间F：

F(x,y)＝[x,y,I(x,y),|I_x(x,y)|,|I_y(x,y)|,|I_xx(x,y)|,|I_yy(x,y)|]

式中，(x,y)是像素点位置，I(x,y)是(x,y)像素点的灰度值，|I_x(x,y)|和|I_y(x,y)|是像素点(x,y)的一阶梯度范数，|I_xx(x,y)|和|I_yy(x,y)|是像素点(x,y)的二阶梯度范数；

4.2)对特征空间F采用积分图像方法计算得到每一像素的Log-Euclidean协方差矩阵描述符C(x,y)：

式中，μ是C上中心点在(x,y)的r邻域内的均值，N_r是r邻域内点的总个数，x'和y'满足不等式|x-x'|≤r，|y-y'|≤r。

4.根据权利要求1所述的基于Log-Euclidean协方差矩阵描述符的图像非刚性配准方法，其特征在于，在步骤5)中，在Log-Demons算法基础上，采用新目标能量函数：

式中，R和M分别表示参考图像和浮动图像的像素值，符号

表示复合运算符，exp和log分别表示指数函数和对数函数，ν是第n次迭代时未正则化的速度场，ν_c是第n-1次迭代后正则化的速度场，C_R和

分别表示参考图像和浮动图像Log-Euclidean协方差矩阵描述符，φ(p)表示当前估计的变换，符号▽表示梯度，α、β和γ分别是影响着图像噪声、光滑程度和正则化在目标能量函数中的权重，ζ为新添加匹配项的所占权重，p是像素位置；

在步骤6)中，对所提出新目标能量函数进行形变驱动力的计算，得到瞬时位移向量δν_pro：

式中，R(p)和M(p)分别表示参考图像和浮动图像在位置p处的像素值，符号

表示M与函数φ(p)的复合运算，φ(p)表示当前估计的变换，α、β和γ分别是影响着图像噪声、光滑程度和正则化在目标能量函数中的权重，ζ为新添加匹配项的所占权重，log表示对数函数，C_R和

分别表示参考图像和浮动图像Log-Euclidean协方差矩阵描述符，J_p和K_p分别表示新目标能量函数第一项中函数

和第三项中函数

的Jacobian矩阵，

其表达式分别为：

在步骤7)中，利用得到瞬时位移向量更新形变矩阵，包括以下步骤：

7.1)高斯滤波平滑瞬时位移向量：

式中，

为高斯滤波，*表示卷积；

7.2)更新形变矩阵：

式中，ν_n表示第n次迭代后正则化的速度场，符号

表示复合函数运算。

Images