CN105303547A - 一种基于网格匹配的Demons算法的多期CT图像配准方法 - Google Patents

Abstract

一种基于网格匹配的Demons算法的多期CT图像配准方法，计算机X射线断层扫描设备对人体腹部进行CT增强扫描；获取动脉期和肝静脉期的腹部ＣＴ图像，并对获取的腹部CT图像进行阈值处理，阈值处理后的静脉期图像作为参考图像，动脉期图像作为浮动图像；对多期ＣＴ图像进行基于梯度阈值掩膜的刚性预配准；将刚性预配准的浮动图像作为新的浮动图像，参考图像不变，对其进行网格匹配；将网格匹配后的参考图像作为新的参考图像，刚性预配准的浮动图像作为新的浮动图像，对其进行基于微分同胚Demons方法的形变配准。基于网格匹配的微分同胚Demons形变配准在形变配准前加入网格匹配，使得参考图像和浮动图像的灰度一致，克服灰度差这一难题。

Description

一种基于网格匹配的Demons算法的多期CT图像配准方法

技术领域

本发明属于医学图像处理技术领域，具体涉及多期CT图像配准方法。

背景技术

多时期肝部增强CT扫描目前是肝脏疾病诊断技术之中的一种实践标准。增强CT扫描的特点是扫描前静脉内注入含碘对比剂，注入后对比剂首先经下腔静脉进入腹部动脉肝脏动脉，因此腹部动脉和肝脏动脉显影，此时获取的图像称为动脉期图像；然后对比剂进入肝门静脉和肝静脉，分别获得肝门静脉和肝静脉期的图像。由于CT部分容积效应，直径较小的病灶在CT扫描上难以显示或者显示不清，而肝脏多期扫描可以发现平扫不能显示的病灶，尤其是小于1cm的小病灶。在获取到不同时期的增强CT图像后，放射科医生需要对不同时期两幅图像的信息进行综合来对肝部疾病进行识别和定位。然而，病人在两次图像获取之间的呼吸作用会使肝部的位置和形状发生微妙的变化，这大大增加了医生诊断的难度，甚至可能导致误诊情况的出现。因此，需要一种配准的方法对两个不同时期CT图像中的肝部和其他组织进行对准，从而改善医生的诊断过程。其次，在肝脏介入治疗的图像导航***中，术中的配准是极为关键的一个环节，配准的质量和速度直接影响着定位和手术效果。但是目前对脑部配准的研究比较多，且配准方法已经成熟，而对于腹部脏器的配准因其形变以及复杂性和模糊性遇到了很多困难，而对于肝脏这种软组织的配准，变形的处理是不得不考虑的因素。

医学图像配准方法分为刚性配准方法和非刚性配准方法。其中刚性配准方法由医学图像配准框架决定，包括三个部分，变换模型，测度函数，优化策略。由于刚性配准方法仅是形变配准的预配准步骤，因此简单效率高的方法更适合做刚性预配准方法。梯度下降法是最常见最早使用的优化策略算法，由于其简单、计算复杂度低且适用于任何维度的图像配准中，因此传统的刚性配准算法经常采用梯度下降法。在同模态尤其是多模态图像间的配准中，互信息是最常用的一种相似性测度函数。其主要优点是无需建立具体的依赖关系，可以建模复杂的相互映射关系。这种灵活性使得互信息成为多模态配准常用测度函数。然而互信息与两幅图像的重叠有关，因此在配准的过程中两幅图像重叠时会出现局部极值从而导致配准的失败。多分辨率方法采用由粗到精的策略，通常用来提高配准性能减少局部极值。但是对于三维图像配准来说，随着配准精确性的提高，计算复杂度也大大增加，无法同时做到精确性和高效性。

非刚性配准方法中Demons算法是一个主流算法，由于计算复杂度低，性能较好，经常用来配准由呼吸作用引起的形变。Demons方法是由Thirion在1998年受Maxwell的热动力学原理启发提出的一种可形变的图像配准方法。该Demons方法使用两幅的灰度和梯度作为外力来驱动配准的进行，并且求解过程采用迭代方式。首先选取图像所有像素或者体素作为Demons，然后迭代一次，求解demons的位移域，但由于其位移域是基于局部信息得到的，所以需要对求得的位移域使用高斯核函数进行平滑处理来达到全局配准的目的，最后使用三线性插值方法变换目标图像，这四个过程反复迭代直到达到迭代次数或者迭代条件终止。

Thirion的方法非常高效，并且在主流的配准算中各种精确度评价指标最高。其中ActiveDemons算法将浮动图像的梯度信息作为一种正内力，参考图像的梯度信息作为负内力，利用这两种力同时驱动形变，该方法能够克服原始Demons配准算法的部分缺陷，收敛速度更快，配准精确性更高。对称的Demons配准算法在原始Demons算法的基础上将形变的内力改进为对称梯度，即把参考图像和浮动图像的梯度平均化，认为它们对于图像形变的贡献等同，而不是局限于使用参考图像的梯度，该方法相对于其他Demons力的收敛速度快，匹配误差小。微分同胚的demons方法将旧的形变域和迭代域的指数映射exp(.)相加作为当前的形变域。由于计算每一次迭代后指数映射相对复杂，因此采用一种高效的计算方法-计算矩阵指数的SS(ScalingandSquaring)方法，对形变域做若干次的复合运算就可以快速地计算位移域的指数映射。位移域(迭代域)的指数映射保证了形变的微分同胚性，即指数映射后仍然是微分同胚的，整个位移域都是微分同胚的，保证了形变的可逆性。为了能将Demons算法更好的适用于计算机解剖学中微分同胚Demons方法也被扩展为对称的，微分同胚因可以保持物体的拓扑结构以及避免不合理的物理形变的特点，因此本方法选择微分同胚Demons算法作为形变配准算法。但是微分同胚Demons算法仅仅适用于无灰度差的图像配准中，不适合多期CT图像的配准。

发明内容

针对现有技术存在的问题，本发明提供一种基于网格匹配的Demons算法的多期CT图像配准方法。

本发明的技术方案是：

一种基于网格匹配的Demons算法的多期CT图像配准方法，包括以下步骤：

步骤1：计算机X射线断层扫描设备对人体腹部进行CT扫描；

步骤2：获取腹部CT图像，每组腹部CT图像包含若干张腹部断层扫描图像；

步骤3：阈值处理得到参考图像R，浮动图像M；

步骤4：基于梯度阈值掩膜的刚性预配准，得到预配准后的浮动图像；

步骤4.1：计算参考图像和浮动图像的两层图像金字塔，用简单的高斯金字塔生成图像序列，初始金字塔层数为0，变换参数默认为0；

步骤4.2：取当前层金字塔参考图像和浮动图像，平移变换参数为上一层变换参数的2倍，旋转参数和上一层变换参数一致；

步骤4.3：计算参考图像的梯度阈值掩膜图像并在掩膜上随机采样得到一个采样点集合；

步骤4.4：对采样点进行刚性变换；

步骤4.5：对刚性变换后的采样点进行线性灰度差值，并更新联合概率分布图，然后计算出图像R和M的互信息；

步骤4.6：当互信息没有达到最大时，执行步骤4.7；否则执行步骤4.8

步骤4.7：梯度优化器根据互信息的梯度值生成下一次刚想变换参数，更新变换矩阵，返回步骤4.4；

步骤4.8：如果层数达到最大值，那么循环结束，输出浮动图像Mp；如果层数没有达到最大值，那么让层数的值加1，重新执行步骤4.2，再次进入循环，直至循环到层数达到最大值，输出浮动图像Mp；

步骤5：网格匹配，得到网格匹配后的参考图像，即新的参考图像；

步骤5.1：以步骤4所得图像Mp和源图像R作为输入，计算网格联合直方图H(i,j)；

步骤5.2：对于H(i,j)中的任意i，找出概率最大的配对(i,j)，在参考图像中用j替换i；

步骤5.3：如果网格的遍历没有完成，返回步骤5.1；否则输出变换图像R_M；

步骤5.4：将源图像R和图像Mp划分网格，在对应的网格上进行直方图匹配，从而对源图像R进行变换，得到网格直方图匹配后的图像R_H；

步骤5.4：将源图像R通过网络直方图匹配所得到的变换图像和变换图像R_M进行图像融合，输出变换图像R_I；

步骤6：微分同胚Demons形变配准,得到最终的结果图像；

步骤6.1：将参考图像R_I和浮动图像Mp作为微分同胚的Demons算法的输入，进行计算；

步骤6.2：输出步骤6.1计算计算结果，得到最终结果图像M_R。

所述步骤4.1计算参考图像和浮动图像的两层图像金字塔，用简单的高斯金字塔生成图像序列，具体步骤如下：

步骤4.1.1：设金字塔层数为L，则使用方差为(L/2)²高斯平滑图像；

步骤4.1.2：相应的层级对图像进行重采样以缩小到相应大小。

所述步骤4.3计算参考图像的梯度阈值掩膜图像并得到采样点，具体步骤如下：

步骤4.3.1：首先计算参考图像的梯度图像；

步骤4.3.2：将梯度图像进行高斯平滑，即保证梯度阈值掩膜图像中包含部分梯度为0的点；

步骤4.3.3：然后设置阈值为Th，该值一般的范围是[1,10]，即阈值处理后将梯度为0的点去除；如果图像中的点的梯度值小于设定阈值Th，则将该点的梯度阈值掩膜设为0；如果图像中的点的梯度值大于等于设定阈值Th，则将该点的梯度阈值掩膜设为1，得到梯度阈值掩膜图像和采样点；

步骤4.3.4：在掩膜上随机采样若干数量的点。

所述步骤4.7梯度优化器根据互信息的梯度值生成下一次刚想变换参数具体步骤如下：

步骤4.7.1：梯度值生成下一次的变换参数的计算公式如下，对于旋转参数，下一次的变换参数X⁽ⁿ⁺¹⁾的计算公式表示为：

X⁽ⁿ⁺¹⁾＝γ_n▽F(X⁽ⁿ⁾)

对于平移参数，计算公式表示为：

X⁽ⁿ⁺¹⁾＝X⁽ⁿ⁾-γ_n▽F(X⁽ⁿ⁾)

互信息的梯度▽F(X⁽ⁿ⁾)表示为

$\▿F\left(X^{\left(n\right)}\right)=\frac{\∂S_{\mathrm{MI}}}{\∂\mathrm{\μ}}=\underset{l\∈L_M}{\mathrm{\Σ}}\underset{k\∈L_R}{\mathrm{\Σ}}\frac{\∂p(l,k;\mathrm{\μ})}{\∂\mathrm{\μ}}{\log }_2\frac{p(l,k;\mathrm{\μ})}{p(l;\mathrm{\μ})}$

其中μ表示X⁽ⁿ⁾，p(l,k；μ)表示联合概率函数，k为参考图像的灰度级，l为浮动图像的灰度级。

此处γ_n＞0为松弛因子，具体设置为0.9；梯度下降优化器搜索参数时的步长设置为第一层：(0.001,3)，第二层：(0.0001,1.5)，迭代200次。

步骤4.7.2：更新变换矩阵，返回步骤4.4；

所述步骤6.1将参考图像R_I和浮动图像Mp作为微分同胚的Demons算法的输入，进行计算，微分同胚的Demons算法的具体步骤如下：

步骤6.1.1：预先设定迭代次数N和初始形变t为零向量(t必须属于李群空间，因此初始化为0)；

步骤6.1.2：根据下面的公式可得到稠密度形变场u，u属于李代数空间，r为参考图像RI，g为浮动图像Mp；

$u=\frac{2(r-g)\·(\▿r+\▿g)}{{|\▿r+\▿g|}^2+{(r-g)}^2\·k}$

步骤6.1.3：对u做一次高斯卷积：u←K_fluid*u，产生流体映射模型的规则化效果，规则化后仍然属于李群空间；

步骤6.1.4：c←tοexp(u)，则c仍然属于李群空间；

步骤6.1.5：对c做高斯卷积并更新t，c←tοexp(u)；

步骤6.1.6：判断迭代次数是否已经为N，如小于N，则返回(2)；否则输出最优变换t^opt。

所述步骤6.1.4计算c的复合运算具体步骤：

步骤6.1.4.1：令v＝2^-nu，选择一个n值，使得max||2^-nu(p)||≤ε，ε＝0.5；

步骤6.1.4.2：对于每一个像素位置，计算v，然后令v＝Id+v；

步骤6.1.4.3：对v做n次复合运算，则exp(u)＝vοv...οv。

本发明提出了一种基于网格匹配的Demons算法的多期CT图像配准方法。在刚性预配准阶段通过改进采样方法，提高了采样点的质量，降低采样率，从而提高了算法的效率，同时改进了算法的鲁棒性精确性，加快了收敛速度。在形变配准阶段，由于Demons算法并不能直接应用于多时期图像中，仅仅适用于同模态等灰度的图像，因此受到互信息的原理启发得出网格匹配变换方法，并通过网格直方图匹配保留了原参考图像的细节纹理信息，最终将参考图像变换到和浮动图像灰度一致的图像，创新的将网格匹配与Demons方法相结合应用于多模态或者存在灰度差的腹部图像配准中。而目前的研究没有针对多时期肝脏CT图像配准，且并没有研究人员将网格匹配和Demons方法相结合应用于多时期图像的配准中，因此本发明的研究内容具有创新性，有着重要的意义。

附图说明

图1是本发明的整体多时期肝脏配准算法流程；

图2是传统刚性配准算法GMM的流程图；

图3是本发明的GTM-GMM方法流程图；

图4是本发明的多时期肝脏切片配准结果，其中(a)参考图像切片，(b)浮动图像切片，(c)梯度掩膜第一层结果切片，(d)梯度掩膜第二层结果切片，(e)配准结果图像切片；

图5是本发明的网格匹配算法流程；

图6是本发明的网格匹配示意图，其中(a)参考图像，(b)浮动图像，(c)联合直方图，(d)基于互信息的网格匹配结果图像，(e)网格直方图匹配结果图像，(f)最终网格匹配结果图像((d)与(f)的融合结果)；

图7是本发明的多时期肝脏三维维图像配准切片。

具体实施方式

下面结合附图对本发明的具体实施方式做详细说明。

本实施方式中有计算机和计算机X射线断层扫描设备；

计算机X射线断层扫描设备用于获取多起腹部CT图像

计算机用来配准图像，硬件配置为IntelCorei7CPU、3.4HZ，8GB内存

一种基于Demons算法的多期图像配准方法，如图6所示，包括以下步骤：

步骤1：计算机X射线断层扫描设备对人体腹部进行增强CT扫描；

步骤2：获取腹部CT图像，每组腹部CT图像包含若干张腹部断层扫描图像；；

所有的图像数据格式为512×512×401pixels³，间距为0.683594×0.683594×0.5mm³；

步骤3：阈值处理得到参考图像R，浮动图像M，阈值范围为(0,300)；

步骤4：基于梯度阈值掩膜的刚性预配准，得到预配准后的浮动图像；

步骤4.1：计算参考图像和浮动图像的两层图像金字塔，用简单的高斯金字塔生成图像序列。

步骤4.1.1：使用方差σ为1高斯平滑图像，高斯方程如下：

$k\left(\right||x-x_c|\left|\right)=e^{-\frac{{\left|\right|x-x_c\left|\right|}^2}{2{\mathrm{\δ}}^2}}$

步骤4.1.2：第一层图像通过对图像进行重采样缩小到原来的二分之一，第二层图像不变

步骤4.2：取第一层金字塔参考图像和浮动图像，六个变换参数初始化为0；

步骤4.3：计算参考图像的梯度阈值掩膜图像并在掩膜上进行随机采样；

步骤4.3.1：首先计算参考图像的梯度图像；

步骤4.3.2：将梯度图像进行高斯平滑，高斯核为2，即保证梯度阈值掩膜图像中包含部分梯度为0的点；

步骤4.3.3：然后设置阈值为Th，第一层金字塔配准的Th值为1，第二层为2；如果图像中的点的梯度值小于设定阈值Th，则将该点的梯度阈值掩膜设为0；如果图像中的点的梯度值大于等于设定阈值Th，则将该点的梯度阈值掩膜设为1，得到梯度阈值掩膜图像和采样点。

步骤4.3.4：在掩膜上随机采样若干数量的点，采样点个数第一层金字塔配准为20000，第二层为40000个；

步骤4.4：对采样点进行刚性变换，变换矩阵为旋转矩阵与平移矩阵相加，三个坐标轴对应的平移参数Tx、Ty、Tz和旋转参数(q₁、q₂、q₃、q₄)，计算公式如下：

$T=\left[\begin{array}{cccc}1-2({q_2}^2+{q_3}^2)& 2(q_1{q}_2+q_3{q}_4)& 2(q_1{q}_3-q_2{q}_4)& T_x\\ 2(q_1{q}_2-q_3{q}_4)& 1-2({q_1}^2+{q_3}^2)& 2(q_2{q}_3+q_1{q}_4)& T_y\\ 2(q_1{q}_3+q_2{q}_4)& 2(q_2{q}_3-q_1{q}_4)& 1-2({q_1}^2+{q_2}^2)& T_z\\ 0& 0& 0& 1\end{array}\right]$

步骤4.5：对刚性变换后的采样点进行三线性灰度差点p被包围的体素的八个点(A-H)围成的正方体内的灰度变换是线性的：

f(x,y,x)＝(1-x)(1-y)(1-y)f(A)+(1-x)y(1-z)f(B)

+(1-x)(1-y)zf(C)+(1-x)yzf(D)

+x(1-y)(1-z)f(E)+xy(1-z)f(F)

+x(1-y)zf(G)+xyzf(H)

然后更新联合概率分布图，联合概率的计算公式如下：

$p_{\mathrm{ij}}=q_{\mathrm{ij}}\÷\left({\mathrm{\Σ}}_{i-0}^{N-1}{\mathrm{\Σ}}_{j=0}^{M-1}{q}_{\mathrm{ij}}\right)$

其中，i和j表示两幅图像的联合分布图的具体坐标，i∈[0:N-1]，j∈[0:M-1]，N和M分别表示两幅图像的bin的个数，bin的个数为30。q_ij表示一个bin在联合分布图中的频数。比如任意一对灰度级(x,y),其中x是图像X中的一个像素灰度，y对应Y中的一个像素灰度，当x落在i中，y落在j中，那么频数q_ij增加1，因此p_ij是联合概率分布图中bin(i,j)的概率。最后计算出互信息，互信息的定义如下公式：

$I(X,Y)=\underset{x,y}{\mathrm{\Σ}}{p}_{\mathrm{xy}}{\log }_2\frac{p_{\mathrm{XY}}(x,y)}{p_X\left(x\right)\·p_Y\left(y\right)}$

P_x(x)表示灰度x在图像X中出现的概率，在实际计算中，x是一个灰度区间，通常称为一个bin。x,y分别表示图像X和Y的bin。图像的互信息反映了图像处理过程中的图像像素的灰度的分布情况。

步骤4.6：当互信息没有达到最大时，执行步骤4.7，否则执行步骤4.8；

步骤4.7：梯度优化器根据互信息的梯度值生成下一次刚想变换参数，更新变换矩阵，返回步骤4.4；

步骤4.7.1：梯度值生成下一次的变换参数的计算公式如下，对于旋转参数，下一次的变换参数X⁽ⁿ⁺¹⁾的计算公式表示为：

X⁽ⁿ⁺¹⁾＝γ_n▽F(X⁽ⁿ⁾)

对于平移参数，计算公式表示为：

X⁽ⁿ⁺¹⁾＝X⁽ⁿ⁾-γ_n▽F(X⁽ⁿ⁾)

互信息的梯度▽F(X⁽ⁿ⁾)表示为

$\▿F\left(X^{\left(n\right)}\right)=\frac{\∂S_{\mathrm{MI}}}{\∂\mathrm{\μ}}=\underset{l\∈L_M}{\mathrm{\Σ}}\underset{k\∈L_R}{\mathrm{\Σ}}\frac{\∂p(l,k;\mathrm{\μ})}{\∂\mathrm{\μ}}{\log }_2\frac{p(l,k;\mathrm{\μ})}{p(l;\mathrm{\μ})}$

其中μ表示X⁽ⁿ⁾，p(l,k；μ)表示联合概率函数，k为参考图像的灰度级，l为浮动图像的灰度级。

此处γ_n＞0为松弛因子，具体设置为0.9；梯度下降优化器搜索参数时的步长设置为第一层：(0.001,3)，第二层：(0.0001,1.5)，迭代200次。

步骤4.7.2更新变换矩阵，返回步骤4.4；

步骤4.8：如果是最后一层金字塔配准，那么循环结束，输出浮动图像Mp；如果层数没有达到最大值，那么让层数的值加1，重新执行步骤4.2，再次进入循环，直至循环到层数达到最大值，输出浮动图像Mp；

步骤5：网格匹配，对参考图像和浮动图像的网格内像素进行匹配。

步骤5.1：以步骤4所得图像Mp和源图像R作为输入，计算网格联合直方图H(i,j)；

设每个网格大小为16×16，在每个网格内计算两幅图像的联合概率分布，其中Bin的个数设为30，则网格的联合概率分布H(R_sub(p),M_sub(p))可表示为：

$H(R_{\mathrm{sub}}\left(p\right),M_{\mathrm{sub}}\left(p\right))=H(I\left(p\right),J\left(p\right))=\frac{q_{\mathrm{ij}}}{{16}^2}$

其中I(p)表示像素p属于Bin的具体范围，I(p)∈[0,N-1]；J(P)同样如此。

步骤5.2：对于H(i,j)中的任意i，找出概率最大的配对(i,j)，在参考图像中用j替换i；

变换的依据是在M中找出与R中像素灰度值覆盖最多的J，并用J替换R的像素，公式可表达为：

$R_T\left(p\right)=\underset{j}{\max }\left(H(I\left(p\right),j)\right)$

步骤5.3：如果网格的遍历没有完成，那么返回步骤5.1继续循环；如果网格的遍历已经完成，那么输出变换图像R_M；

步骤5.4：将源图像R通过网络直方图匹配所得到的变换图像R_H

将源图像R和浮动图像Mp划分为32×32大小的网格，在每个网格上进行直方图匹配，匹配点为6个；

步骤5.5：R_H和变换图像R_M进行图像融合，输出变换图像R_I；

步骤6：微分同胚Demons形变配准,得到最终的结果图像；

步骤6.1：将参考图像R_I和浮动图像Mp作为微分同胚的Demons算法的输入，进行计算；

步骤6.1.1：预先设定迭代次数N＝30，和初始形变t为零向量(t必须属于李群空间，因此初始化为0)；

步骤6.1.2：根据下面的公式可得到稠密度形变场u，u属于李代数空间，δ_x＝5.0；

$u=\frac{2(r-g)\·(\▿r+\▿g)}{{|\▿r+\▿g|}^2+\frac{{(r-g)}^2}{{\mathrm{\δ}}_x^2}}$

步骤6.1.3：对u做一次高斯卷积：u←K_fluid*u产生流体映射模型的规则化效果，规则化后仍然属于李群空间，高斯核为2；

步骤6.1.4：c←tοexp(u)，则c仍然属于李群空间；

其中计算c的复合运算具体步骤为：

步骤6.1.4.1：令v＝2^-nu，选择一个n值，使得max||2^-nu(p)||≤ε，ε＝0.5；

步骤6.1.4.2：对于每一个像素位置，计算v，然后令v＝Id+v；

步骤6.1.4.3：对v做n次复合运算，则exp(u)＝vοv...οv；

步骤6.1.5：对c做高斯卷积并更新t，t←K_diff*c,高斯核为2；

步骤6.1.6：判断迭代次数是否已经为N，如小于N，则返回(2)；否则输出最优变换t^opt。

步骤6.2：输出步骤6.1计算计算结果，得到最终结果图像M_R。

Claims (1)

1.一种基于网格匹配的Demons算法的多期CT图像配准方法，包括以下步骤：

步骤1：计算机X射线断层扫描设备对人体腹部进行CT扫描；

步骤2：获取腹部CT图像，每组腹部CT图像包含若干张腹部断层扫描图像；

步骤3：阈值处理得到参考图像R，浮动图像M；

步骤4：基于梯度阈值掩膜的刚性预配准，得到预配准后的浮动图像；

步骤4.1：计算参考图像和浮动图像的两层图像金字塔，用简单的高斯金字塔生成图像序列，初始金字塔层数为0，变换参数默认为0；

步骤4.2：取当前层金字塔参考图像和浮动图像，平移变换参数为上一层变换参数的2倍，旋转参数和上一层变换参数一致；

步骤4.3：计算参考图像的梯度阈值掩膜图像并在掩膜上随即采样得到一个采样点集合；

步骤4.4：对采样点进行刚性变换；

步骤4.5：对刚性变换后的采样点进行线性灰度差值，并更新联合概率分布图，然后计算出图像R和M的互信息；

步骤4.6：当互信息没有达到最大时，执行步骤4.7；否则执行步骤4.8

步骤4.7：梯度优化器根据互信息的梯度值生成下一次刚想变换参数，更新变换矩阵，返回步骤4.4；

步骤4.8：如果层数达到最大值，那么循环结束，输出浮动图像Mp；如果层数没有达到最大值，那么让层数的值加1，重新执行步骤4.2，再次进入循环，直至循环到层数达到最大值，输出浮动图像Mp；

步骤5：网格匹配，得到网格匹配后的参考图像，即新的参考图像；

步骤5.1：以步骤4所得图像Mp和源图像R作为输入，计算网格联合直方图H(i,j)；

步骤5.2：对于H(i,j)中的任意i，找出概率最大的配对(i,j)，在参考图像中用j替换i；

步骤5.3：如果网格的遍历没有完成，返回步骤5.1；否则输出变换图像R_M；

步骤5.4：将源图像R和图像Mp划分网格，在对应的网格上进行直方图匹配，从而对源图像R进行变换，得到网格直方图匹配后的图像R_H；

步骤5.4：将源图像R通过网络直方图匹配所得到的变换图像和变换图像R_M进行图像融合，输出变换图像R_I；

步骤6：微分同胚Demons形变配准,得到最终的结果图像；

步骤6.1：将参考图像R_I和浮动图像Mp作为微分同胚的Demons算法的输入，进行计算；

步骤6.2：输出步骤6.1计算计算结果，得到最终结果图像M_R。

所述步骤4.1计算参考图像和浮动图像的两层图像金字塔，用简单的高斯金字塔生成图像序列，具体步骤如下：

步骤4.1.1：设金字塔层数为L，则使用方差为(L/2)²高斯平滑图像；

步骤4.1.2：相应的层级对图像进行重采样以缩小到相应大小。

所述步骤4.3计算参考图像的梯度阈值掩膜图像并得到采样点，具体步骤如下：

步骤4.3.1：首先计算参考图像的梯度图像；

步骤4.3.2：将梯度图像进行高斯平滑，即保证梯度阈值掩膜图像中包含部分梯度为0的点；

步骤4.3.3：然后设置阈值为Th，该值一般的范围是[1,10]，即阈值处理后将梯度为0的点去除；如果图像中的点的梯度值小于设定阈值Th，则将该点的梯度阈值掩膜设为0；如果图像中的点的梯度值大于等于设定阈值Th，则将该点的梯度阈值掩膜设为1，得到梯度阈值掩膜图像和采样点；

步骤4.3.4：在掩膜上随机采样若干数量的点。

所述步骤4.7梯度优化器根据互信息的梯度值生成下一次刚想变换参数具体步骤如下：

步骤4.7.1：梯度值生成下一次的变换参数的计算公式如下，对于旋转参数，下一次的变换参数X⁽ⁿ⁺¹⁾的计算公式表示为：

X⁽ⁿ⁺¹⁾＝γ_n▽F(X⁽ⁿ⁾)

对于平移参数，计算公式表示为：

X⁽ⁿ⁺¹⁾＝X⁽ⁿ⁾-γ_n▽F(X⁽ⁿ⁾)

互信息的梯度▽F(X⁽ⁿ⁾)表示为

其中μ表示X⁽ⁿ⁾，p(l,k；μ)表示联合概率函数，k为参考图像的灰度级，l为浮动图像的灰度级。

此处γ_n＞0为松弛因子，具体设置为0.9；梯度下降优化器搜索参数时的步长设置为第一层：(0.001,3)，第二层：(0.0001,1.5)，迭代200次。

步骤4.7.2：更新变换矩阵，返回步骤4.4；

所述步骤6.1将参考图像R_I和浮动图像Mp作为微分同胚的Demons算法的输入，进行计算，微分同胚的Demons算法的具体步骤如下：

步骤6.1.1：预先设定迭代次数N和初始形变t为零向量(t必须属于李群空间，因此初始化为0)；

步骤6.1.2：根据下面的公式可得到稠密度形变场u，u属于李代数空间，r为参考图像RI，g为浮动图像Mp；

步骤6.1.3：对u做一次高斯卷积：u←K_fluid*u，产生流体映射模型的规则化效果，规则化后仍然属于李群空间；

步骤6.1.4：c←tοexp(u)，则c仍然属于李群空间；

步骤6.1.5：对c做高斯卷积并更新t，c←tοexp(u)；

步骤6.1.6：判断迭代次数是否已经为N，如小于N，则返回(2)；否则输出最优变换t^opt。

所述步骤6.1.4计算c的复合运算具体步骤：

步骤6.1.4.1：令v＝2^-nu，选择一个n值，使得max||2^-nu(p)||≤ε，ε＝0.5；

步骤6.1.4.2：对于每一个像素位置，计算v，然后令v＝Id+v；

步骤6.1.4.3：对v做n次复合运算，则exp(u)＝vοv...οv。