非刚性医学图像配准是医学图像研究领域的热门专题之一,具有重要的临床应用价值。本文提出一种改进的Demons算法,将灰度守恒模型和局部结构张量守恒模型结合起来,构造一个新的能量函数处理多模态配准问题,然后采用L-BFGS算法优化能量函数,解决复杂三维数据的优化问题,并采用多尺度分层细化的思想解决大形变配准。实验表明,本文算法对大形变和多模态三维医学图像配准都有较好的效果。
引用本文： 郝培博, 陈震, 江少锋, 汪洋. 基于Demons算法的多模态医学图像非刚性配准研究. 生物医学工程学杂志, 2014, 31(1): 161-165. doi: 10.7507/1001-5515.20140032 复制
引言
医学图像的配准已经成为研究领域的热门之一。医学图像配准是指通过寻找某种空间变换，使两幅图像的对应点达到空间位置和解剖结构上的完全一致，要求配准的结果能使两幅图像上所有的解剖点，或至少是所有具有诊断意义以及手术区域的点都达到匹配[1]。医学图像配准具有很重要的临床应用价值。
随着医学图像学的不断发展，刚性配准已经不能满足临床的需要，越来越多的学者开始研究非刚性医学图像配准。研究人员提出了大量的非刚性医学图像配准算法。这些算法可分为二大类：基于特征点的算法、基于物理模型的算法[2]。
基于特征点的算法首先要在每个需要配准的数据集上建立代表不同解剖元素的几何模型，然后与目标图像中的对应解剖元素进行配准，最后利用解剖元素之间配准得到的对应关系引导图像其余部分之间的配准。Likar等[3]提出将图像划分为若干子块，取每个子块的中心作为标记对应点，但没有旋转变换，提取的标记点有偏差，配准精度有待提高。冯林等[4]提出了一种基于分层互信息和薄板样条自动确定标记对应点的选取方法，但是由于该算法依赖灰度统计相关性，无法取得足够多的标记点。基于物理模型的算法通常采用来自数学或统计学的相似性准则，以匹配图像每个区域的灰度分布模式为目标，通常会在形变的模板图像和目标图像之间定义一个灰度相似性测度，最后通过物理模型变换得到配准图像。D’Agostino等[5]用粘滞流体模型控制形变，该算法可以实现任何复杂形状的形变，但是所需计算量大。Lester等[6]在流体模型中，把粘滞常数看成可以变化的，因此对于图像不同的部分可以有不同的形变特性。
Demons算法是一种小形变无参数非刚性配准方法[7]。它具有较高的运算效率和良好的配准效果。Hellier等[8]对常用的六种配准算法的回顾性对比研究表明，Demons算法在多个评价指标中都优于其他五种算法。但是Demons算法在处理大形变配准和多模态配准时候，会产生较大的误配准和物理上的不合理变形。因此，本研究旨在对Demons算法进行一定的改进，使其能在大形变和多模态下具有较高的准确性。
1 Demons算法原理及改进
Demons配准驱动力是以光流方程式驱动图像小形变变化。对静态图像F一个给定点p，用f表示静态图像的p点的强度，m表示浮动图像的p点的强度，那么浮动图像到静态图像配准点p的估计偏移量u为
| $~u=\frac{\left( m-f \right)\nabla f}{{{\left| \nabla f \right|}^{2}}+{{\left( m-f \right)}^{2}}},$ |
其中u=(ux，uy)是浮动图像上任意点p到静态图像上对应点的偏移量，∇f是静态图像的梯度值，(m-f)2是为了梯度方程更加稳定地添加到方程中。Wang 等[9]增加了一个方程使用浮动图像边缘作为内在驱动力，改进了算法的收敛速度和稳定性，
| $u=\frac{\left( m-f \right)\nabla f}{{{\left| \nabla f \right|}^{2}}+{{\partial }^{2}}{{\left( m-f \right)}^{2}}}+\frac{\left( m-f \right)\nabla m}{{{\left| \nabla m \right|}^{2}}+{{\partial }^{2}}{{\left( m-f \right)}^{2}}}$ |
为了调整驱动力的强度,Cachier等[10]提出了归一化因子∂。
Vercauteren等[11]描绘了标准配准模型；一个基于能量的相似性测度，求得一个配准误差函数，当误差函数为0时求得最优解。相似性测度是衡量变换后图像m(x-u)好坏的重要标准。基于灰度守恒的配准中最常用的相似性测度是像素距离平方：
| $E\left( m,f,u \right)=\frac{1}{2}\int{{{\left| m\left( x-u \right)-f\left( x \right) \right|}^{2}}dx,}$ |
所以可以计算出误差函数：
| $\nabla E\left( u \right)=2{{\left( \nabla m \right)}^{T}}\left( f-m+u\nabla m \right)$ |
若直接对式(4) 求最小化会产生不惟一解，那优化问题就存在不适定性，因此需要对式(4) 增加正则项来解决不适定性，但是Thirion提出的Demons算法只是通过简单的高斯平滑来正则化，使图像结构丢失，图像出现变形。
虽然图像的梯度方向也包含了局部结构信息，但是在计算图像的局部方向时需要在每个像素的邻域做平均，如果在某个像素的局部邻域中同时包含上升边缘的梯度和下降边缘的梯度(如一条细边缘的两侧)，由于此时梯度出现正负抵消，因此无法反映该点的局部方向。局部结构张量很好地避免了梯度方向的这一不足，同时能够区分图像中的各种不同结构，近年来被广泛应用于图像处理领域。局部结构张量守恒假设是根据图像局部结构在运动方向上基本不变的特性提出的，这样就可以避免运动方向性的约束。 形式为
| $T\left( X+w,t+1 \right)-T\left( X,t \right)=0,$ |
式(5) 中T为图像的局部结构张量，包含图像的局部运动信息(任意方向运动信息)。
| $T=\left[ \begin{matrix} I_{x}^{2} & {{I}_{x}}{{I}_{y}} \\ {{I}_{x}}{{I}_{y}} & I_{y}^{2} \\ \end{matrix} \right]$ |
在设计能量函数时，可以根据图像中运动场景的不同和各守恒假设适用的范围，合理地选择其中一种或者几种守恒假设结合的方式，需注意的是数据项中选择的守恒假设越多，越有可能导致模型计算复杂度的提高和效率的降低，因此通常选择两种假设模型的结合方式构成数据项模型。
综合考虑灰度守恒假设与图像局部结构守恒假设的适用范围和优缺点，将两者结合建立一个新的能量函数[12]，得到的能量函数(相似性测度)为
| $\eqalign{ & {E_{data}} = {E_g} + {\lambda _1}{E_t} = \cr & \int_\Omega {\left( {{{\left| {I\left( {X + W} \right) - I\left( X \right)} \right|}^2} + {\lambda _1}*{{\left| {T\left( {X + W} \right) - T\left( X \right)} \right|}^2}} \right)dX} \cr & {\int_\Omega {\left( {\left| {{m_x}\left( {x - u} \right) - f\left( x \right)} \right|} \right.} ^2} + {\lambda _1}* \cr & \left[ \matrix{ {\left| {{m_x}\left( {x - u} \right) - f\left( x \right)} \right|^2} + {\left| {{m_y}\left( {x - u} \right) - {f_y}\left( x \right)} \right|^2} + \hfill \cr {\left| {{m_x}\left( {x - u} \right)*{m_y}\left( {x - u} \right) - {f_x}\left( x \right){f_y}\left( x \right)} \right|^2} \hfill \cr} \right]\left. {} \right)dX, \cr} $ |
其中λ1表示数据项中两守恒假设的权重系数。则新的误差函数为
| $\begin{align} & \nabla E\left( u \right)=2{{\left( \nabla m \right)}^{T}}\left( f-m+u\nabla m \right)+2{{\left( \nabla {{m}_{x}} \right)}^{T}} \\ & \left( {{f}_{x}}-{{m}_{x}}+u\nabla {{m}_{x}} \right)+2{{\left( \nabla {{m}_{y}} \right)}^{T}}\left( {{f}_{y}}-{{m}_{y}}+u\nabla {{m}_{y}} \right) \\ & +2{{\left( \nabla {{m}_{x}} \right)}^{T}}{{\left( \nabla {{m}_{y}} \right)}^{T}}\left( {{f}_{x}}{{f}_{y}}-{{m}_{x}}{{m}_{y}}+2u\nabla {{m}_{x}}\nabla {{m}_{y}} \right) \\ \end{align}$ |
其中mx、my为浮动图像在x和y方向的梯度，fx、fy为参考图像在x和y方向的梯度。
为了将算法应用到三维图像中，可以把二维局部结构张量扩展到三维，三维结构张量包含了三维空间中的典型结构。设I(X)为图像在点X(x,y,z)的灰度，此时的局部结构张量为
| $T=\left[ \begin{matrix} I_{x}^{2} & {{I}_{x}}{{I}_{y}} & {{I}_{x}}{{I}_{z}} \\ {{I}_{x}}{{I}_{y}} & I_{y}^{2} & {{I}_{y}}{{I}_{z}} \\ {{I}_{x}}{{I}_{z}} & {{I}_{y}}{{I}_{z}} & I_{z}^{2} \\ \end{matrix} \right]$ |
三维图像数据量大，如果不对三维数据进行优化处理，配准过程相当耗时，本文采用L-BFGS算法优化能量函数，此算法将BFGS方法与有限内存技术相结合，用来求解大规模优化问题[13]。其特点是不需要明确地形成或存储Hessian矩阵，只保存x处前m次的更新历史以及梯度，是一种有限内存拟牛顿法。
Demons算法对于大形变非刚性配准不能得到良好效果，为了改善这种情况，采用多尺度分层细化思想，在同一幅图像中，较高分辨率下可以观察到尺寸很小或者灰度级不高的物体，而在较低的分辨率下却只能观察到尺寸很大或者灰度级较高的物体，对图像采用多分辨率处理具有明显的优势。因此，这样的从粗到精的多分辨率思想在大位移非刚性运动时特别适用。
2 改进算法的实现
改进算法的具体步骤：
(1) 将浮动图像m和静态图像f进行多尺度金字塔分层。
(2) 从金字塔最高层k层开始，将mk,fk代入式(2) 中，计算出形变场uk。
(3) 由浮动图像mk和形变场uk得到新的浮动图像mk(x-u)。
(4) 将mk(x-u)和fk代入能量函数式(7) 中，判断是否满足相似性测度。若不满足需将mk(x-u)和fk重新代入误差函数式(8) 中，计算出新的形变场uk，直到满足相似性测度；若满足则将形变场uk插值得到下一层的初始形变场uk-1。
(5) 将mk-1，fk-1和uk-1重新代入步骤3中迭代，直到计算到金字塔底层。
(6) 将底层得到的形变场u0代入浮动图像m，得到最终配准后图像。
3 实验结果及讨论
实验目的是为了验证本算法在不同模态情况下和大形变情况下的配准效果优于Demons算法。
3.1 配准精度的评价方法
目前医学图像配准还没有金标准来衡量配准结果的好坏，在这里采用主观评价，最小均方误差(MSE)和相关系数(Rcc)来评价算法的优劣。
MSE的计算公式为
| $MSE = \sqrt {\sum\limits_{x = 0}^{N - 1} {\sum\limits_{y = 0}^{N - 1} {{{\left[ {M\left( {x,y} \right) - F\left( {x,y} \right)} \right]}^2}} } } ,$ |
MSE 越小表示配准效果越好。
Rcc的计算公式为
| ${{R}_{cc}}=\frac{\frac{1}{N-1}\sum\limits_{i=0}^{N-1}{\left( {{M}_{i}}-\bar{M} \right)\left( {{F}_{i}}-\bar{F} \right)}}{\sqrt{\frac{1}{N-1}\sum\limits_{i=0}^{N-1}{{{\left( {{M}_{i}}-\bar{M} \right)}^{2}}}}\sqrt{\frac{1}{N-1}\sum\limits_{i=0}^{N-1}{{{\left( {{F}_{i}}-\bar{F} \right)}^{2}}}}},$ |
其中N为像素点个数，
3.2 实验结果
实验一：验证算法在二维图像配准中也能取得良好结果，采用256×256像素的颅脑MR图像，用质子密度加权成像和T1加权成像代表不同模态的图像。实验结果如图 1所示，表 1为本文算法与Demons算法等的比较。
图1
二维图像配准结果
(a)浮动图像;(b)参考图像;(c)互信息算法结果;(d)B样条算法结果;(e)Demons算法结果;(f)本文算法结果
Figure1. 2D image registration results(a) moving image; (b) static image; (c) mutual information algorithm result; (d) B-spline algorithm result; (e) Demons algorithm result; (f) the proposed algorithm result
表1
不同算法的配准精度比较
Table1.
Registration accuracy comparison among different algorithms
从图 1中可知，(a)为质子密度加权成像图像偏移了13x17y，(b)为T1加权成像图像，(c)为互信息算法配准后图像，(d)为B样条算法配准后图像，(e)为Demons算法配准后图像，(f)为本文算法配准后
的图像。可以直观看到，B样条算法和Demons算法在处理不同模态图像配准时，出现了严重的变形和失真；而本文算法和互信息算法能很好地处理多模态情况下的配准。从表 1中可知，本文算法的MSE小于其他几种算法，而且Rcc都高于其他几种算法，说明本文算法优于其他几种算法，能够很好地处理多模态图像配准。
实验二：验证算法在三维配准中的应用，采用模拟MRT脑数据网站(http://brainweb.bic.mni.mcgill.ca/brainweb/)中的模拟MR脑数据，选取mr_T1数据为181×217×181，切片厚度为1 mm，含有3%噪声的图像序列绕垂直轴旋转10°作为浮动图像；选取质子密度加权成像mr_PD数据为181×217×181，切片厚度为1 mm，含有3%噪声的图像序列作为参考图像。配准后的结果如
图2
三维多模态配准
(a)浮动图像;(b)参考图像;(c)配准后图像;(d)配准前后第90张切片图像差值
Figure2. 3D image registration results(a) moving image; (b) static image; (c) registration result; (d) image difference of the 90th slice before and after registration
从图 2中可知，(a)为T1数据181×217×181绕垂直轴旋转10°并偏移5x5y，(b)为质子密度加权成像数据181×217×181，(c)为配准后T1数据，从直观上可以看到矢状面和冠状面与参考图像基本一致，水平面上的偏移量和旋转角度都得到了较好的配准，本文算法对多模态三维图像配准也能够得到良好的效果。目前还没有衡量图像三维配准的金标准，故采用配准前后与参考图像的差值比较来表示配准结果的好坏。如图 2(d)所示，(d)中左边是配准前图像差值，配准前有明显的偏移和误差；右边是配准后图像差值，偏移和误差有明显的改善。本文算法对不同模态的三维医学图像也能取得较好的配准结果。
实验三：选取mr_T1数据为181×217×181，切片厚度为1 mm，含有3%噪声的图像序列绕垂直轴旋转10°作为浮动图像；选取mr_T1数据为181×217×60，切片厚度为3 mm，含有3%噪声的图像序列作为参考图像。配准结果如图 3所示。
图3
不同切片厚度的三维配准
(a)浮动图像;(b)参考图像;(c)配准后图像;(d)配准前后第90张切片图像差值
Figure3. 3D registration with different slice thicknessesl(a) moving image; (b) static image; (c) registration result; (d) image difference of the 90th slice before and after registration
如图 3所示，(a)为T1数据181×217×181绕垂直轴旋转10°并偏移5x5y，(b)为T1数据181×217×60，切片厚度为3 mm，由于切片厚度不同，在配准前先进行部分体积差值，使图像尺寸相同再进行配准。(c)为配准后的数据，从直观上可以看到冠状面和矢状面的灰度值有所变化，这是由部分体积差值引起的，水平面上选择和偏移都得到了良好的改善，配准前后的差异如(d)所示，左边是配准前图像差值，右边是配准后图像差值。从(d)中可以看到，配准前两幅图像差异非常明显，有较多的伪影和偏移；配准后图像得到了较大的改善，伪影和偏移明显减少。本文算法在处理不同尺寸的体素配准时也能取得较好的结果。
4 结论
本文从多种能量守恒相结合的思想出发，把局部结构张量守恒假设和灰度守恒相结合，得到新的能量函数，使算法能适应不同模态图像的配准。将结构张量应用到能量函数中，使图像的结构得以保持，以适应大形变的图像配准。在算法实现中，采用多尺度分层细化策略，使图像运动信息更加准确。实验表明，本文算法对大形变配准、多模态配准都能得到较好的结果，并能得到比Demons算法精确度更好的结果。
引言
医学图像的配准已经成为研究领域的热门之一。医学图像配准是指通过寻找某种空间变换，使两幅图像的对应点达到空间位置和解剖结构上的完全一致，要求配准的结果能使两幅图像上所有的解剖点，或至少是所有具有诊断意义以及手术区域的点都达到匹配[1]。医学图像配准具有很重要的临床应用价值。
随着医学图像学的不断发展，刚性配准已经不能满足临床的需要，越来越多的学者开始研究非刚性医学图像配准。研究人员提出了大量的非刚性医学图像配准算法。这些算法可分为二大类：基于特征点的算法、基于物理模型的算法[2]。
基于特征点的算法首先要在每个需要配准的数据集上建立代表不同解剖元素的几何模型，然后与目标图像中的对应解剖元素进行配准，最后利用解剖元素之间配准得到的对应关系引导图像其余部分之间的配准。Likar等[3]提出将图像划分为若干子块，取每个子块的中心作为标记对应点，但没有旋转变换，提取的标记点有偏差，配准精度有待提高。冯林等[4]提出了一种基于分层互信息和薄板样条自动确定标记对应点的选取方法，但是由于该算法依赖灰度统计相关性，无法取得足够多的标记点。基于物理模型的算法通常采用来自数学或统计学的相似性准则，以匹配图像每个区域的灰度分布模式为目标，通常会在形变的模板图像和目标图像之间定义一个灰度相似性测度，最后通过物理模型变换得到配准图像。D’Agostino等[5]用粘滞流体模型控制形变，该算法可以实现任何复杂形状的形变，但是所需计算量大。Lester等[6]在流体模型中，把粘滞常数看成可以变化的，因此对于图像不同的部分可以有不同的形变特性。
Demons算法是一种小形变无参数非刚性配准方法[7]。它具有较高的运算效率和良好的配准效果。Hellier等[8]对常用的六种配准算法的回顾性对比研究表明，Demons算法在多个评价指标中都优于其他五种算法。但是Demons算法在处理大形变配准和多模态配准时候，会产生较大的误配准和物理上的不合理变形。因此，本研究旨在对Demons算法进行一定的改进，使其能在大形变和多模态下具有较高的准确性。
1 Demons算法原理及改进
Demons配准驱动力是以光流方程式驱动图像小形变变化。对静态图像F一个给定点p，用f表示静态图像的p点的强度，m表示浮动图像的p点的强度，那么浮动图像到静态图像配准点p的估计偏移量u为
| $~u=\frac{\left( m-f \right)\nabla f}{{{\left| \nabla f \right|}^{2}}+{{\left( m-f \right)}^{2}}},$ |
其中u=(ux，uy)是浮动图像上任意点p到静态图像上对应点的偏移量，∇f是静态图像的梯度值，(m-f)2是为了梯度方程更加稳定地添加到方程中。Wang 等[9]增加了一个方程使用浮动图像边缘作为内在驱动力，改进了算法的收敛速度和稳定性，
| $u=\frac{\left( m-f \right)\nabla f}{{{\left| \nabla f \right|}^{2}}+{{\partial }^{2}}{{\left( m-f \right)}^{2}}}+\frac{\left( m-f \right)\nabla m}{{{\left| \nabla m \right|}^{2}}+{{\partial }^{2}}{{\left( m-f \right)}^{2}}}$ |
为了调整驱动力的强度,Cachier等[10]提出了归一化因子∂。
Vercauteren等[11]描绘了标准配准模型；一个基于能量的相似性测度，求得一个配准误差函数，当误差函数为0时求得最优解。相似性测度是衡量变换后图像m(x-u)好坏的重要标准。基于灰度守恒的配准中最常用的相似性测度是像素距离平方：
| $E\left( m,f,u \right)=\frac{1}{2}\int{{{\left| m\left( x-u \right)-f\left( x \right) \right|}^{2}}dx,}$ |
所以可以计算出误差函数：
| $\nabla E\left( u \right)=2{{\left( \nabla m \right)}^{T}}\left( f-m+u\nabla m \right)$ |
若直接对式(4) 求最小化会产生不惟一解，那优化问题就存在不适定性，因此需要对式(4) 增加正则项来解决不适定性，但是Thirion提出的Demons算法只是通过简单的高斯平滑来正则化，使图像结构丢失，图像出现变形。
虽然图像的梯度方向也包含了局部结构信息，但是在计算图像的局部方向时需要在每个像素的邻域做平均，如果在某个像素的局部邻域中同时包含上升边缘的梯度和下降边缘的梯度(如一条细边缘的两侧)，由于此时梯度出现正负抵消，因此无法反映该点的局部方向。局部结构张量很好地避免了梯度方向的这一不足，同时能够区分图像中的各种不同结构，近年来被广泛应用于图像处理领域。局部结构张量守恒假设是根据图像局部结构在运动方向上基本不变的特性提出的，这样就可以避免运动方向性的约束。 形式为
| $T\left( X+w,t+1 \right)-T\left( X,t \right)=0,$ |
式(5) 中T为图像的局部结构张量，包含图像的局部运动信息(任意方向运动信息)。
| $T=\left[ \begin{matrix} I_{x}^{2} & {{I}_{x}}{{I}_{y}} \\ {{I}_{x}}{{I}_{y}} & I_{y}^{2} \\ \end{matrix} \right]$ |
在设计能量函数时，可以根据图像中运动场景的不同和各守恒假设适用的范围，合理地选择其中一种或者几种守恒假设结合的方式，需注意的是数据项中选择的守恒假设越多，越有可能导致模型计算复杂度的提高和效率的降低，因此通常选择两种假设模型的结合方式构成数据项模型。
综合考虑灰度守恒假设与图像局部结构守恒假设的适用范围和优缺点，将两者结合建立一个新的能量函数[12]，得到的能量函数(相似性测度)为
| $\eqalign{ & {E_{data}} = {E_g} + {\lambda _1}{E_t} = \cr & \int_\Omega {\left( {{{\left| {I\left( {X + W} \right) - I\left( X \right)} \right|}^2} + {\lambda _1}*{{\left| {T\left( {X + W} \right) - T\left( X \right)} \right|}^2}} \right)dX} \cr & {\int_\Omega {\left( {\left| {{m_x}\left( {x - u} \right) - f\left( x \right)} \right|} \right.} ^2} + {\lambda _1}* \cr & \left[ \matrix{ {\left| {{m_x}\left( {x - u} \right) - f\left( x \right)} \right|^2} + {\left| {{m_y}\left( {x - u} \right) - {f_y}\left( x \right)} \right|^2} + \hfill \cr {\left| {{m_x}\left( {x - u} \right)*{m_y}\left( {x - u} \right) - {f_x}\left( x \right){f_y}\left( x \right)} \right|^2} \hfill \cr} \right]\left. {} \right)dX, \cr} $ |
其中λ1表示数据项中两守恒假设的权重系数。则新的误差函数为
| $\begin{align} & \nabla E\left( u \right)=2{{\left( \nabla m \right)}^{T}}\left( f-m+u\nabla m \right)+2{{\left( \nabla {{m}_{x}} \right)}^{T}} \\ & \left( {{f}_{x}}-{{m}_{x}}+u\nabla {{m}_{x}} \right)+2{{\left( \nabla {{m}_{y}} \right)}^{T}}\left( {{f}_{y}}-{{m}_{y}}+u\nabla {{m}_{y}} \right) \\ & +2{{\left( \nabla {{m}_{x}} \right)}^{T}}{{\left( \nabla {{m}_{y}} \right)}^{T}}\left( {{f}_{x}}{{f}_{y}}-{{m}_{x}}{{m}_{y}}+2u\nabla {{m}_{x}}\nabla {{m}_{y}} \right) \\ \end{align}$ |
其中mx、my为浮动图像在x和y方向的梯度，fx、fy为参考图像在x和y方向的梯度。
为了将算法应用到三维图像中，可以把二维局部结构张量扩展到三维，三维结构张量包含了三维空间中的典型结构。设I(X)为图像在点X(x,y,z)的灰度，此时的局部结构张量为
| $T=\left[ \begin{matrix} I_{x}^{2} & {{I}_{x}}{{I}_{y}} & {{I}_{x}}{{I}_{z}} \\ {{I}_{x}}{{I}_{y}} & I_{y}^{2} & {{I}_{y}}{{I}_{z}} \\ {{I}_{x}}{{I}_{z}} & {{I}_{y}}{{I}_{z}} & I_{z}^{2} \\ \end{matrix} \right]$ |
三维图像数据量大，如果不对三维数据进行优化处理，配准过程相当耗时，本文采用L-BFGS算法优化能量函数，此算法将BFGS方法与有限内存技术相结合，用来求解大规模优化问题[13]。其特点是不需要明确地形成或存储Hessian矩阵，只保存x处前m次的更新历史以及梯度，是一种有限内存拟牛顿法。
Demons算法对于大形变非刚性配准不能得到良好效果，为了改善这种情况，采用多尺度分层细化思想，在同一幅图像中，较高分辨率下可以观察到尺寸很小或者灰度级不高的物体，而在较低的分辨率下却只能观察到尺寸很大或者灰度级较高的物体，对图像采用多分辨率处理具有明显的优势。因此，这样的从粗到精的多分辨率思想在大位移非刚性运动时特别适用。
2 改进算法的实现
改进算法的具体步骤：
(1) 将浮动图像m和静态图像f进行多尺度金字塔分层。
(2) 从金字塔最高层k层开始，将mk,fk代入式(2) 中，计算出形变场uk。
(3) 由浮动图像mk和形变场uk得到新的浮动图像mk(x-u)。
(4) 将mk(x-u)和fk代入能量函数式(7) 中，判断是否满足相似性测度。若不满足需将mk(x-u)和fk重新代入误差函数式(8) 中，计算出新的形变场uk，直到满足相似性测度；若满足则将形变场uk插值得到下一层的初始形变场uk-1。
(5) 将mk-1，fk-1和uk-1重新代入步骤3中迭代，直到计算到金字塔底层。
(6) 将底层得到的形变场u0代入浮动图像m，得到最终配准后图像。
3 实验结果及讨论
实验目的是为了验证本算法在不同模态情况下和大形变情况下的配准效果优于Demons算法。
3.1 配准精度的评价方法
目前医学图像配准还没有金标准来衡量配准结果的好坏，在这里采用主观评价，最小均方误差(MSE)和相关系数(Rcc)来评价算法的优劣。
MSE的计算公式为
| $MSE = \sqrt {\sum\limits_{x = 0}^{N - 1} {\sum\limits_{y = 0}^{N - 1} {{{\left[ {M\left( {x,y} \right) - F\left( {x,y} \right)} \right]}^2}} } } ,$ |
MSE 越小表示配准效果越好。
Rcc的计算公式为
| ${{R}_{cc}}=\frac{\frac{1}{N-1}\sum\limits_{i=0}^{N-1}{\left( {{M}_{i}}-\bar{M} \right)\left( {{F}_{i}}-\bar{F} \right)}}{\sqrt{\frac{1}{N-1}\sum\limits_{i=0}^{N-1}{{{\left( {{M}_{i}}-\bar{M} \right)}^{2}}}}\sqrt{\frac{1}{N-1}\sum\limits_{i=0}^{N-1}{{{\left( {{F}_{i}}-\bar{F} \right)}^{2}}}}},$ |
其中N为像素点个数，
3.2 实验结果
实验一：验证算法在二维图像配准中也能取得良好结果，采用256×256像素的颅脑MR图像，用质子密度加权成像和T1加权成像代表不同模态的图像。实验结果如图 1所示，表 1为本文算法与Demons算法等的比较。
图1
二维图像配准结果
(a)浮动图像;(b)参考图像;(c)互信息算法结果;(d)B样条算法结果;(e)Demons算法结果;(f)本文算法结果
Figure1. 2D image registration results(a) moving image; (b) static image; (c) mutual information algorithm result; (d) B-spline algorithm result; (e) Demons algorithm result; (f) the proposed algorithm result
表1
不同算法的配准精度比较
Table1.
Registration accuracy comparison among different algorithms
从图 1中可知，(a)为质子密度加权成像图像偏移了13x17y，(b)为T1加权成像图像，(c)为互信息算法配准后图像，(d)为B样条算法配准后图像，(e)为Demons算法配准后图像，(f)为本文算法配准后
的图像。可以直观看到，B样条算法和Demons算法在处理不同模态图像配准时，出现了严重的变形和失真；而本文算法和互信息算法能很好地处理多模态情况下的配准。从表 1中可知，本文算法的MSE小于其他几种算法，而且Rcc都高于其他几种算法，说明本文算法优于其他几种算法，能够很好地处理多模态图像配准。
实验二：验证算法在三维配准中的应用，采用模拟MRT脑数据网站(http://brainweb.bic.mni.mcgill.ca/brainweb/)中的模拟MR脑数据，选取mr_T1数据为181×217×181，切片厚度为1 mm，含有3%噪声的图像序列绕垂直轴旋转10°作为浮动图像；选取质子密度加权成像mr_PD数据为181×217×181，切片厚度为1 mm，含有3%噪声的图像序列作为参考图像。配准后的结果如
图2
三维多模态配准
(a)浮动图像;(b)参考图像;(c)配准后图像;(d)配准前后第90张切片图像差值
Figure2. 3D image registration results(a) moving image; (b) static image; (c) registration result; (d) image difference of the 90th slice before and after registration
从图 2中可知，(a)为T1数据181×217×181绕垂直轴旋转10°并偏移5x5y，(b)为质子密度加权成像数据181×217×181，(c)为配准后T1数据，从直观上可以看到矢状面和冠状面与参考图像基本一致，水平面上的偏移量和旋转角度都得到了较好的配准，本文算法对多模态三维图像配准也能够得到良好的效果。目前还没有衡量图像三维配准的金标准，故采用配准前后与参考图像的差值比较来表示配准结果的好坏。如图 2(d)所示，(d)中左边是配准前图像差值，配准前有明显的偏移和误差；右边是配准后图像差值，偏移和误差有明显的改善。本文算法对不同模态的三维医学图像也能取得较好的配准结果。
实验三：选取mr_T1数据为181×217×181，切片厚度为1 mm，含有3%噪声的图像序列绕垂直轴旋转10°作为浮动图像；选取mr_T1数据为181×217×60，切片厚度为3 mm，含有3%噪声的图像序列作为参考图像。配准结果如图 3所示。
图3
不同切片厚度的三维配准
(a)浮动图像;(b)参考图像;(c)配准后图像;(d)配准前后第90张切片图像差值
Figure3. 3D registration with different slice thicknessesl(a) moving image; (b) static image; (c) registration result; (d) image difference of the 90th slice before and after registration
如图 3所示，(a)为T1数据181×217×181绕垂直轴旋转10°并偏移5x5y，(b)为T1数据181×217×60，切片厚度为3 mm，由于切片厚度不同，在配准前先进行部分体积差值，使图像尺寸相同再进行配准。(c)为配准后的数据，从直观上可以看到冠状面和矢状面的灰度值有所变化，这是由部分体积差值引起的，水平面上选择和偏移都得到了良好的改善，配准前后的差异如(d)所示，左边是配准前图像差值，右边是配准后图像差值。从(d)中可以看到，配准前两幅图像差异非常明显，有较多的伪影和偏移；配准后图像得到了较大的改善，伪影和偏移明显减少。本文算法在处理不同尺寸的体素配准时也能取得较好的结果。
4 结论
本文从多种能量守恒相结合的思想出发，把局部结构张量守恒假设和灰度守恒相结合，得到新的能量函数，使算法能适应不同模态图像的配准。将结构张量应用到能量函数中，使图像的结构得以保持，以适应大形变的图像配准。在算法实现中，采用多尺度分层细化策略，使图像运动信息更加准确。实验表明，本文算法对大形变配准、多模态配准都能得到较好的结果，并能得到比Demons算法精确度更好的结果。
