CN103345741A - 一种非刚性多模医学图像精确配准方法 - Google Patents

Abstract

本发明公开了一种非刚性多模医学图像精确配准方法，包括以下步骤：分别对参考图像和浮动图像进行归一化处理并使用脉冲发送皮层模型进行处理得出两幅点火映射图；然后采用索贝尔算子分别提取两幅点火映射图的边缘特征得到对应的梯度图，并计算两梯度图的差值平方和；最后基于自由网格变形模型及差值平方和确定目标函数，采用拟牛顿法对目标函数进行迭代优化，获得最佳形变参数。本发明将脉冲发送皮层模型和差值平方和用于非刚性多模医学图像配准，大大提高了配准精度。

Description

一种非刚性多模医学图像精确配准方法

技术领域

本发明属于图像分析与处理中的图像配准领域，更具体地，涉及一种非刚性多模医学图像精确配准方法。

背景技术

随着计算机技术和图形图像学的发展，医学成像技术也不断得到提高。由于不同医学成像设备存在各自的使用范围和局限，单独使用某一类图像难以获得正确的诊断结论，将多种模式医学图像融合，充分利用图像自身的特点做到信息互补，对提高医学影像诊断的准确性具有重要意义。因待融合的多模图像其成像方位、角度和分辨率等可能存在较大差异，要实现有效的医学图像融合，需先对多模图像进行配准。图像配准是图像融合的关键技术，它通过寻找一种空间变换把一幅图像映射到另一幅图像，使两者在几何上能够匹配对应，已在外科手术计划的制定、放射治疗计划的制定、病理变化的跟踪和治疗效果的评价等方面得到了广泛应用。

多模医学图像配准主要包括刚性配准和非刚性配准两类，前者目前已基本趋于成熟，后者作为近年来图像配准问题的研究热点则有待进一步发展。针对非刚性多模图像配准问题，目前提出的方法主要包括基于灰度和基于特征的配准方法。在基于灰度的配准方法中，典型例子是基于互信息测度的配准方法，如传统互信息法、基于归一化互信息（Normalized MutualInformation,简称NMI）的方法及基于条件互信息（Conditional MutualInformation,简称CMI）的方法，这类方法其计算复杂度高，易陷入局部最优解，且可能引入图像变换误差。基于特征的配准方法中，典型例子包括基于局部方向梯度、基于局部频率、基于尺度不变特征变换(Scale InvariantFeature Transform，SIFT)及基于熵图（Entropy image）的配准方法。其中，基于熵图的配准方法先计算待匹配图像的熵图，然后计算熵图的差值平方和（Sum of Squared Differences,简称SSD），这种将Entropy image和SSD相结合的ESSD方法其配准精度优于传统互信息法，但该方法中熵图无法有效表征复杂医学图像的结构特征，因此难以实现多模医学图像的精确配准。

发明内容

针对现有技术的上述缺陷，本发明提供了一种非刚性多模医学图像精确配准方法，其利用自由形变模型和待配准图像的梯度图的差值平方和确定目标函数，利用拟牛顿法对目标函数进行迭代优化，获得最佳形变参数。本发明提出的方法能有效提取待配准图像的特征信息，克服现有方法其配准精度欠佳的不足。

为实现上述目的，按照本发明的一个方面，提供了一种非刚性多模医学图像精确配准方法，包括以下步骤：

步骤1对参考图像

进行归一化处理，以得到图像I_R；

步骤2基于脉冲发送皮层模型获得图像I_R的点火映射图I_R′；

步骤3基于索贝尔算子获得点火映射图I_R′的梯度图本步骤具体包括以下子步骤：

步骤3-1分别用水平方向的索贝尔算子

和竖直方向的索贝尔算子在点火映射图I_R′上滑动，使其中心都位于I_R′的像素点(i,j)，以进行卷积运算，并选取其中的较大值作为(i,j)处的梯度值

$G_R^{w_1}(i,j)=\max \{S_x^{w_1}\⊗{I_R}^{\′},S_y^{w_1}\⊗{I_R}^{\′}\}$

其中，表示卷积运算；w₁×w₁表示索贝尔算子的窗口大小；max{}表示取极大值运算。

步骤3-2重复索贝尔算子的滑动操作，直到求出图像所有的像素点的取值，从而得到参考图像对应点火映射图的梯度图

步骤4对浮动图像重复步骤1至3，得到其对应点火映射图的梯度图

步骤5采用以下公式计算两梯度图和

的差值平方和

$\mathrm{SSD}(G_R^{w_1},G_F^{w_1})=\underset{c=1}{\overset{N}{\mathrm{\Σ}}}{\left|\right|G_R^{w_1}\left(x_c\right)-G_F^{w_1}\left(x_c\right)\left|\right|}_2^2$

其中，x_c表示任意单个像素点；‖·‖₂表示欧几里得距离；

步骤6改变索贝尔算子窗口大小为w₂×w₂，重复步骤3至5，以得到两梯度图

的差值平方和

步骤7确定目标函数并利用拟牛顿法进行优化迭代，得到最佳形变参数，本步骤具体包括以下子步骤：

步骤7-1由步骤5和步骤6的输出结果计算两点火映射图I_R′和I_F′的相似性：

$\mathrm{sobSSD}({I_R}^{\′},{I_F}^{\′})=\frac{\mathrm{SSD}(G_R^{w_1},G_F^{w_1})+\mathrm{SSD}(G_R^{w_2},G_F^{w_2})}{2}$

步骤7-2基于步骤7-1的相似性度量公式，结合自由网格变形模型T确定目标函数f：

f(T;I_R′,I_F′)＝sobSSD(T;I_R′,I_F′)+γC_smooth(T)

其中γ为权重参数，C_smooth(T)为正则化项。

步骤7-3使用布洛伊登-弗莱彻-戈德法布-香农（Broyden-Fletcher-Goldfarb-Shannon，简称BFGS）拟牛顿法对目标函数f进行迭代，得到最佳形变模型

$\hat{T}=\arg \min f(T;{I_R}^{\′},{I_F}^{\′})$

2、根据权利要求1所述的一种非刚性多模医学图像精确配准方法，其特征在于，步骤（2）具体包括以下子步骤：

步骤2-1初始化迭代n=1，设置最大迭代次数为N，内部活动项U(n)、活动项阈值E(n)、脉冲输出Y(n)和点火映射图I_R′(n)的初始矩阵U(0)、E(0)、Y(0)和I_R′(0)为零矩阵；

步骤2-2计算图像I_R在像素点(i,j)时的脉冲输出Y_ij(n)：

其中，Y_ij(n)为脉冲发送皮层模型的发放脉冲输出；S_ij为神经元所受的外部激励大小，表示待匹配图像在(i,j)处的灰度值的归一化结果；U_ij(n)为神经元(i,j)在网络运行中第n次迭代时的内部活动项；f为神经元内部活动的衰减系数，g为阈值衰减系数，h为阈值放大系数；W_ijkl为神经元(i,j)与神经元(k,l)间的连接突触权重。

步骤2-3对于图像I_R的所有像素点，重复步骤2-2，以得到整幅图像所有像素点的脉冲输出，从而得出第n次迭代时的点火图Y(n)；

步骤2-4判断n是否小于预定的最大迭代次数N，若是则设置n＝n+1，并返回步骤2-2，否则将N次迭代得到的所有点火图相加，得到图像I_R对应的点火映射图I_R′并进入步骤3；

3、根据权利要求1所述的一种非刚性多模医学图像精确配准方法，其特征在于，w₁×w₁代表3×3的窗口，其索贝尔算子为：

$S_x^3=\left[\begin{array}{ccc}-1& 0& +1\\ -2& 0& +2\\ -1& 0& +1\end{array}\right]$ 和 $S_y^3=\left[\begin{array}{ccc}+1& +2& +1\\ 0& 0& 0\\ -1& -2& -1\end{array}\right].$

4、根据权利要求1所述的非刚性多模医学图像精确配准方法，其特征在于，w₂×w₂代表5×5的窗口，其索贝尔算子为：

$S_x^5=\left[\begin{array}{ccccc}+2& +3& 0& -3& -2\\ +3& +4& 0& -4& -3\\ +6& +6& 0& -6& -6\\ +3& +4& 0& -4& -3\\ +2& +3& 0& -3& -2\end{array}\right]$ 和 $S_y^5=\left[\begin{array}{ccccc}+2& +3& +6& +3& +2\\ +3& +4& +6& +4& +3\\ 0& 0& 0& 0& 0\\ -3& -4& -6& -4& -3\\ -2& -3& -6& -3& -2\end{array}\right].$

5、根据权利要求1所述的非刚性多模医学图像精确配准方法，其特征在于，正则化项C_smooth(T)是通过以下公式计算：

$C_{\mathrm{smooth}}\left(T\right)=\frac{1}{V}{\∫}_0^x{\∫}_0^y[{\left(\frac{{\∂}^{2}T}{\∂x^2}\right)}^2+2\left(\frac{{\∂}^{2}T}{\∂\mathrm{xy}}\right)+{\left(\frac{{\∂}^{2}T}{\∂y^2}\right)}^2]\mathrm{dxdy}$

其中，V为图像域的大小。

总体而言，本发明提出的上述技术方案与现有技术相比，能够取得下列有益效果：

（1）本发明将脉冲发送皮层模型应用于图像配准，充分利用了该模型对待配准图像其特征信息的表征能力，为实现基于特征的精确图像配准提供了基础。

（2）本发明在利用脉冲发送皮层模型获得的点火映射图的基础上，利用索贝尔算子获得的梯度图的差值平方和作为相似性度量，结合BFGS拟牛法获得最佳形变参数，可有效降低目标配准误差，为实现多模医学图像精确配准提供了有效手段。

附图说明

图1是本发明非刚性多模医学图像配准方法的流程图。

图2是本发明的参考图像和浮动图像所采用的T1-MRI、PD-MRI、T2-MRI三种模式的解剖标记图，图上数字标注的点的位置为解剖学标记点。其中，图2a为T1-MRI的图像（浮动图像），图2b为PD-MRI的图像（浮动图像），图2c为T2-MRI的图像（参考图像）。

图3是本发明的方法与现有NMI和ESSD方法的配准结果对比，其中，图3a为T1-MRI图（浮动图像），图3b为PD-MRI图（浮动图像）,图3c为T2-MRI图（参考图像），图3d为本发明的算法的配准结果图（T1-T2），图3e为ESSD算法的配准结果图（T1-T2），图3f为CMI算法的配准结果图（T1-T2），图3g为NMI的算法的配准结果图（T1-T2），图3h为本发明的算法的配准结果图（PD-T2），图3i为ESSD算法的配准结果图（PD-T2），图3j为CMI算法的配准结果图（PD-T2）,图3k为NMI算法的配准结果图（PD-T2）

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅用以解释本发明，并不用于限定本发明。此外，下面所描述的本发明各个实施方式中所涉及到的技术特征只要彼此之间未构成冲突就可以相互组合。

如图1所示，本发明非刚性多模医学图像精确配准方法包括以下步骤：

步骤1对参考图像

进行归一化处理得到图像I_R；

步骤2基于脉冲发送皮层模型获得图像I_R的点火映射图I_R′，本步骤具体包括以下子步骤：

步骤2-1初始化迭代次数n=1，在实例中设置最大迭代次数N为40次，内部活动项U(n)、活动项阈值E(n)、脉冲输出Y(n)和点火映射图I_R′(n)的初始矩阵U(0)、E(0)、Y(0)和I_R′(0)为零矩阵；

步骤2-2计算图像I_R在(i,j)时的脉冲输出Y_ij(n)：

其中，Y_ij(n)为脉冲发送皮层模型的发放脉冲输出；S_ij为神经元所受的外部激励大小，在数字图像处理中表示像素点(i,j)位置灰度图；U_ij(n)为神经元(i,j)在网络运行中第n次迭代时的内部活动项；f为神经元内部活动的衰减系数，实例中具体取值0.7；g为阈值衰减系数，实例中具体取值为0.8；h为阈值放大系数，实例中具体取值为20；W_ijkl为神经元(i,j)与神经元(k,l)间的连接突触权重，且(k,l)≠(i,j)，实例中取值如下：

$W_{\mathrm{ijkl}}=\left[\begin{array}{ccc}0.1091& 0.1409& 0.1091\\ 0.1409& 0& 0.1409\\ 0.1091& 0.1409& 0.1091\end{array}\right]$

步骤2-3对于图像I_R的所有像素点，重复步骤2-2，以得到整幅图像所有像素点的脉冲输出，从而得出第n次迭代时的点火图Y(n)；

步骤2-4判断n是否小于预定的最大迭代次数N，若是则设置n＝n+1，并返回步骤2-2，否则将N次迭代得到的所有点火图相加，得到图像I_R对应的点火映射图I_R′并进入步骤3；

步骤3基于索贝尔算子获得点火映射图I_R′的梯度图

本步骤具体包括以下子步骤：

步骤3-1分别用水平方向的索贝尔算子

和竖直方向的索贝尔算子

在点火映射图I_R′上滑动，使其中心都位于I_R′的像素点(i,j)，以进行卷积运算，并选取其中的较大值作为梯度图在(i,j)处的取值

$G_R^{w_1}(i,j)=\max \{S_x^{w_1}\⊗{I_R}^{\′},S_y^{w_1}\⊗{I_R}^{\′}\}$

其中，

表示卷积运算；max{}表示取极大值运算；w₁×w₁表示索贝尔算子的窗口大小，实例中w₁×w₁代表3×3的窗口，此时3×3的索贝尔算子为：

$S_x^3=\left[\begin{array}{ccc}-1& 0& +1\\ -2& 0& +2\\ -1& 0& +1\end{array}\right]$ 和 $S_y^3=\left[\begin{array}{ccc}+1& +2& +1\\ 0& 0& 0\\ -1& -2& -1\end{array}\right]$

步骤3-2重复索贝尔算子的滑动操作，直到求出图像所有的像素点的取值，从而得到参考图像对应点火映射图的梯度图

步骤4对浮动图像

重复步骤1至3，得到其对应的梯度图

步骤5采用以下公式计算两梯度图和

的差值平方和

$\mathrm{SSD}(G_R^{w_1},G_F^{w_1})=\underset{c=1}{\overset{N}{\mathrm{\Σ}}}{\left|\right|G_R^{w_1}\left(x_c\right)-G_F^{w_1}\left(x_c\right)\left|\right|}_2^2$

其中，x_c表示任意单个像素点；‖·‖₂表示欧几里得距离；

步骤6改变索贝尔算子窗口大小为w₂×w₂，重复步骤3至5，以得出两梯度图

的差值平方和为

其中，w₂×w₂在实例中代表5×5的窗口，5×5的索贝尔算子为：

$S_x^5=\left[\begin{array}{ccccc}+2& +3& 0& -3& -2\\ +3& +4& 0& -4& -3\\ +6& +6& 0& -6& -6\\ +3& +4& 0& -4& -3\\ +2& +3& 0& -3& -2\end{array}\right]$ 和 $S_y^5=\left[\begin{array}{ccccc}+2& +3& +6& +3& +2\\ +3& +4& +6& +4& +3\\ 0& 0& 0& 0& 0\\ -3& -4& -6& -4& -3\\ -2& -3& -6& -3& -2\end{array}\right].$

步骤7确定目标函数并进行迭代，以得出最佳形变参数，本步骤具体包括以下子步骤：

步骤7-1由步骤5和步骤6的输出结果计算两点火映射图I_R′和I_F′的相似性：

$\mathrm{sobSSD}({I_R}^{\′},{I_F}^{\′})=\frac{\mathrm{SSD}(G_R^{w_1},G_F^{w_1})+\mathrm{SSD}(G_R^{w_2},G_F^{w_2})}{2}$

步骤7-2基于步骤7-1的相似性度量公式，结合自由网格变形模型T确定目标函数f：

f(T;I_R′,I_F′)＝sobSSD(T;I_R′,I_F′)+γC_smooth(T)

其中γ为权重参数，实例中取值为0.01；C_smooth(T)为正则化项，其计算公式为：

$C_{\mathrm{smooth}}\left(T\right)=\frac{1}{V}{\∫}_0^x{\∫}_0^y[{\left(\frac{{\∂}^{2}T}{\∂x^2}\right)}^2+2\left(\frac{{\∂}^{2}T}{\∂\mathrm{xy}}\right)+{\left(\frac{{\∂}^{2}T}{\∂y^2}\right)}^2]\mathrm{dxdy}$

其中，V为图像域的大小。

步骤7-3使用布洛伊登-弗莱彻-戈德法布-香农（Broyden-Fletcher-Goldfarb-Shannon，简称BFGS）牛顿优化算法对目标函数f进行迭代得出最佳形变模型

$\hat{T}=\arg \min f(T;{I_R}^{\′},{I_F}^{\′})$

实例

本实例采用图像来自于BrainWeb数据库，大小为256×212的T1-MRI,T2-MRI,和PD-MRI图像，其中T1-MRI和PD-MRI为形变图像，即浮动图像，T2-MRI为原始图像，即参考图像，如图2所示。其配准精度由以模拟变形作为金标准的目标配准误差TRE_d和以解剖标志点作为金标准的目标配准误差TRE_l进行评价。

以模拟变形作为金标准的目标配准误差TRE_d可以表示为：

${\mathrm{TRE}}_d=\frac{\left|\right|T_s\left(M\right)-T_c\left(M\right){\left|\right|}_2}{\left|M\right|}$

其中T_s表示模拟形变，T_c表示程序运行生成的形变，M表示整个图像域的像素个数。

对于一个由程序运行自动生成的形变T_c和一组由专家标记出来的解剖标志点{I_LF，I_LR}={(x_i，y_i)，(x_i′，y_i′)}(i=1,2,···,k,其中k是解剖标志点的数目)。以解剖标志点作为金标准的配准误差目标TRE_l由下式计算：

${\mathrm{TRE}}_l=\frac{{\left|\right|T_c\left(I_{\mathrm{LF}}\right)-I_{\mathrm{LR}}\left|\right|}_2}{k}$

基于TRE_l的评估需要参考图像和浮动图像的解剖学标记点，由相关领域专家为经过形变的T1-MRI,PD-MRI和原始的T2-MRI图像进行人工选择标记点，为每一图像选择15个标记点，如图2所示。

我们将本发明的方法与ESSD方法、CMI方法和NMI方法在配准精度方面进行比较。如表1和表2所示，可明显的看出本发明的方法能够较其它方法实现更低的目标配准误差TRE_d和TRE_l，也就是说本发明的方法能后较其它方法获得更高的配准精度。例如表1所示，本发明的方法以模拟变形作为金标准的目标配准误差TRE_d的均值可达到1.1mm，明显低于配准前的3.2mm，并明显低于其它配准方法。

表1本发明方法SCM-sobSSD与NMI，ESSD方法的TRE_d(mm)对比

例如表2所示，本发明的方法以解剖标志点作为金标准的配准误差目标TREl的均值可达7.2mm，明显低于配准前的18.6mm，也低于其它配准方法，从而明显的证明了本发明在配准精度上的有效性。

表2本发明方法SCM-sobSSD与NMI，ESSD方法的TRE_l(mm)对比

另外，从图3中我们可以更加直观的看出本发明提出的方法比ESSD方法和NMI方法相比可以取得更好的配准结果。从图3中我们可以直观的看出T1-MRI或PD-MRI浮动图像经配准后较其它两种方法更接近于参考图像T2-MRI，充分证明了本发明的方法的实践性较强，精确度较高。

本领域的技术人员容易理解，以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

Claims (5)

1.一种非刚性多模医学图像精确配准方法，其特征在于，包括以下步骤： 

步骤1对参考图像进行归一化处理，以得到图像I_R； 

步骤2基于脉冲发送皮层模型获得图像I_R的点火映射图I_R′； 

步骤3基于索贝尔算子获得点火映射图I_R′的梯度图

本步骤具体包括以下子步骤： 

步骤3-1分别用水平方向的索贝尔算子

和竖直方向的索贝尔算子 

在点火映射图I_R′上滑动，使其中心都位于I_R′的像素点(i,j)，以进行卷积运算，并选取其中的较大值作为(i,j)处的梯度值

其中，

表示卷积运算；w₁×w₁表示索贝尔算子的窗口大小；max{}表示取极大值运算。 

步骤3-2重复索贝尔算子的滑动操作，直到求出图像所有的像素点的取值，从而得到参考图像对应点火映射图的梯度图

步骤4对浮动图像

重复步骤1至3，得到其对应点火映射图的梯度图

步骤5采用以下公式计算两梯度图

和

的差值平方和 

其中，x_c表示任意单个像素点；‖·‖₂表示欧几里得距离； 

步骤6改变索贝尔算子窗口大小为w₂×w₂，重复步骤3至5，以得到 两梯度图的差值平方和

步骤7确定目标函数并利用拟牛顿法进行优化迭代，得到最佳形变参数，本步骤具体包括以下子步骤： 

步骤7-1由步骤5和步骤6的输出结果计算两点火映射图I_R′和I_F′的相似性： 

步骤7-2基于步骤7-1的相似性度量公式，结合自由网格变形模型T确定目标函数f： 

f(T;I_R′,I_F′)＝sobSSD(T;I_R′,I_F′)+γC_smooth(T) 

其中γ为权重参数，C_smooth(T)为正则化项。 

步骤7-3使用布洛伊登-弗莱彻-戈德法布-香农（Broyden-Fletcher-Goldfarb-Shannon，简称BFGS）拟牛顿法对目标函数f进行迭代，得到最佳形变模型

。

2.根据权利要求1所述的一种非刚性多模医学图像精确配准方法，其特征在于，步骤（2）具体包括以下子步骤： 

步骤2-1初始化迭代n=1，设置最大迭代次数为N，内部活动项U(n)、活动项阈值E(n)、脉冲输出Y(n)和点火映射图I_R′(n)的初始矩阵U(0)、E(0)、Y(0)和I_R′(0)为零矩阵； 

步骤2-2计算图像I_R在像素点(i,j)时的脉冲输出Y_ij(n)： 

其中，Y_ij(n)为脉冲发送皮层模型的发放脉冲输出；S_ij为神经元所受的外部激励大小，表示待匹配图像在(i,j)处的灰度值的归一化结果；U_ij(n)为神经元(i,j)在网络运行中第n次迭代时的内部活动项；f为神经元内部活动的衰减系数，g为阈值衰减系数，h为阈值放大系数；W_ijkl为神经元(i,j)与神经元(k,l)间的连接突触权重。 

步骤2-3对于图像I_R的所有像素点，重复步骤2-2，以得到整幅图像所有像素点的脉冲输出，从而得出第n次迭代时的点火图Y(n)； 

步骤2-4判断n是否小于预定的最大迭代次数N，若是则设置n＝n+1，并返回步骤2-2，否则将N次迭代得到的所有点火图相加，得到图像I_R对应的点火映射图I_R′并进入步骤3。 

3.根据权利要求1所述的一种非刚性多模医学图像精确配准方法，其特征在于，w₁×w₁代表3×3的窗口，其索贝尔算子为： 

和

4.根据权利要求1所述的非刚性多模医学图像精确配准方法，其特征在于，w₂×w₂代表5×5的窗口，其索贝尔算子为： 

和

5.根据权利要求1所述的非刚性多模医学图像精确配准方法，其特征在于，正则化项C_smooth(T)是通过以下公式计算： 

其中，V为图像域的大小。 

Images