基于空间加权互信息的非刚性医学图像配准

基于空间加权互信息的非刚性医学图像配准
王腾飞;杨艳
【摘要】As the similarity metric criterion,mutual information is widely used in medical image registration field,in particular,in rigidregistration research,its higher accuracy and good robustness in registration effect have been demonstrated.However,due to the complexity ofnonlinear deformation,the registration methods based on mutual information must be thoroughly studied when used in non-rigid registrationarea.The method we proposed in this article uses B-spline based free-form deformation model to simulate nonlinear deformation of anatomicalstructure in medical image,and meanwhile according to this deformation model,and taking the influence of spatial information on registrationeffect into consideration,we improve the mutual information method in the way of spatial weight,and apply LBFGS method,which has goodeffect in large scale parameters optimisation,to optimising the registration parameters,we implement the program through programming,andvalidate and analyse the effect by
experiments.Experimental results show that the proposed scheme obviously outperforms the traditionalmutual information method in registration accuracy.%互信息作为相似性测度标准在医学图像配准领域应用广泛，尤其是在刚性配准研究中，能实现较高配准精度和非常稳健的配准效果。

然而由于非线性形变的复杂性，在非刚性配准领域，对基于互信息的配准方法进行深入研究势在必行。

提出的方案采用基于B 样条的自由变换模型模拟医学图像中解剖
结构的非线性形变，同时根据变换模型，考虑空间信息对配准效果的影响，以空间加权的方式对互信息方法进行改进。

使用对大规模参数优化效果较好的LBFGS方法对配准参数进行优化，编程实现程序，并通过实验对效果进行验证和分析。

实验结果显示，该方案配准精度明显优于传统的互信息方法。

【期刊名称】《计算机应用与软件》
【年(卷),期】2014(000)009
【总页数】5页(P213-216,282)
【关键词】非刚性;互信息;B样条;自由变换;空间加权
【作者】王腾飞;杨艳
【作者单位】武汉大学物理科学与技术学院湖北武汉430072;武汉大学物理科学与技术学院湖北武汉430072
【正文语种】中文
【中图分类】TP391.41
0 引言
医学影像技术为各种疾病的诊断和治疗提供必要的信息依据，因而在临床诊断中应用广泛。

在医学领域，图像配准是指对在同一生理解剖区域(或者近邻的区域)不同时间点或者采用不同模式获取的两幅图像通过线性和非线性的几何变形实现解剖区域完全对准的过程。

图像配准的基本方法分为基于特征的方法和基于像素强度的方法两类，其中基于特征的方法需要首先对待配准图像进行特征提取，配准的精确性很大程度上依赖于特征提取的精确性，而基于像素强度的方法(以基于互信息的方
法为例)则不需要对图像的结构特征和像素强度做任何先验性的假设，具有更高的
通用性。

然而在非刚性配准领域，完全基于互信息的方法很难满足要求，对基于互信息的方法进行改进依然是一个值得深入研究的课题。

引入空间信息是现阶段基于像素强度医学图像配准领域的研究方向之一，尤其是在对非刚性配准的研究中。

对待配准图像进行分块，计算每个图像块对的局部互信息，对每个图像块对分别进行配准是最基本的空间信息引入方式之一，采用这种方式的配准方法效果要明显优于基于全局互信息的配准方法[1]。

然而，引入空间互信息
的方式并不仅限于此，Studholme等在文献[2]中将空间位置信息作为与参考图像和浮动图像地位等同的第三维信息引入到图像配准中，提出基于区域互信息
RMI(Region Mutual Information)的配准方法，Loeckx等在文献[3]中对Studholme的方案进行了改进，提出了条件互信息cMI(conditional Mutual Information)的概念。

Zhuang等[4]采用了与Loeckx类似的思路，但是不同的实现方法,Lu等[5]在统计直方图时进行空间延拓，在像素层面引入空间信息。

上述方法都取得了一定的效果。

本研究针对非刚性配准，考虑医学图像中存在非线性形变的情形，采用基于B样
条的自由变换模型表征待配准图像间的几何形变，使用通过空间加权方式改进的互信息作为配准测度函数，分块实现非刚性医学图像配准，并通过实验对本文方案的效果进行验证和分析。

1 配准方案研究
1.1 基于B样条的自由变换模型
医学图像上生理解剖结构的空间几何形变分为线性形变和非线性形变两类，其中线性形变包括平移、旋转、缩放和剪切这些规则的形变方式，在图像配准领域可以使用仿射变换模型完全表征。

而非线性形变则是一些不规则的形变，比线性形变更加复杂更加难以量化，一般是由人的一些正常生理活动，如呼吸、心跳等导致的。

本
文提出的方案采用基于B样条的自由变换模型[6，7]FFD(Free Form Deformation)来表征非线性形变。

基于B样条的自由变换模型曾被用于研究心脏运动产生的形变，是一种非常有效的非线性形变工具，其基本原理是通过操纵均匀分布于图像上的控制点坐标来实现对图像上解剖结构几何形变的控制，最终的变形结果是每个控制点坐标变化对图像形变影响的总和，是一种二阶导连续的光滑形变。

考虑二维图像，将待配准图像的空间坐标定义在Ω={(x,y)|0≤x<X,0≤y<Y}上,定义一个大小为nx×ny、间隔为Δ的控制点集Φ={φi,j}，则基于B样条的自由变换模型可以表示为：
(1)
其中，i=⎣⎣⎣⎣表示第l个3阶B样条曲线的基函数：
(2)
控制点集Φ是自由变换模型的自由变量集，共有2×nx×ny个自由度，每一种控制点位置变化组合对应唯一的一种图像形变。

与其他一些非线性形变模型，如薄板样条模型[8]和弹性体模型[9]等相比，基于B样条的自由变换模型每个控制点对图像形变的影响都是局部性的，都限制在其近邻的4×4个控制点网格内，如图1所示，因而即使在控制点规模非常大的情形下，也会有相对较高的计算效率。

图1 每个控制点位置的变换只会对其近邻的16个网格内的点的坐标产生影响1.2 基于互信息的图像配准原理
互信息(Mutual Information，MI)是信息论中的一个概念，原本用以表征两个随机变量的统计相关性，在1995年被引入到图像配准领域，是基于像素强度的图像配准领域最常用的相似性测度标准之一[10]。

在图像配准领域，两幅图像的互信息表达的是一幅图像中包含有另一幅图像的信息的大小。

将两幅图像分别定义为参考图像R和浮动图像F，对浮动图像F进行空间变形来实现与参考图像的配准，其互信息I(R,F)可表示为：
I(R,F)=H(R)+H(F)-H(R,F)
(3)
其中，H(R)、H(F)代表参考图像和浮动图像的边缘信息熵，H(R,F)代表两幅图像的联合信息熵，可以根据其强度概率分布函数分别计算得到：
(4)
(5)
(6)
互信息可以作为图像配准相似测度标准的依据是，当两幅图像完全对准时，其互信息值达到最大。

在图像配准过程中，由于几何变换的存在，待配准图像的重叠面积可能会发生改变，而互信息的大小与图像重叠面积相关联，为了规避互信息值对图像重叠面积的依赖性，Studholme[11]提出了使用归一化的互信息作为图像配准的相似性度量标准，可用公式表述为：
(7)
1.3 空间加权互信息方案
将空间信息引入到图像配准中是该领域目前的研究趋势之一。

在前面介绍基于B
样条的自由变换模型时我们注意到，每个控制点φi,j对图像形变的影响都是局部
性的，其范围仅限于以其为中心的4×4个网格以内，且影响程度由近及远逐渐减小。

因此，本文方案考虑将浮动图像进行空间分块处理，针对每一个控制点φi,j，将受其影响的最大范围，也就是一个以φi,j为中心面积为4Δ×4Δ的图像块单独处理，计算其与参考图像中对应区域的归一化互信息。

对大小为nx×ny的控制点网格，共要将待配准图像分成nx×ny个相互重叠的图像块，分别计算对应的局部区域归一化互信息，总的互信息根据每个图像块面积占图像总面积的比重pi,j依概率加权得到。

(8)
同时，考虑到控制点坐标变化对其周围图像形变的影响程度随距离衰减，且衰减程度受B样条基函数控制，因此在计算浮动图像上每一个图像块与对应参考图像块的互信息时，与当前控制点距离不同的点对待配准图像的直方图和联合直方图的贡献不等同考虑。

本文方案采用与形变模型相对应的二维3阶B样条曲线对不同位置像素点的贡献加权。

考虑参考点φi,j对应的图像块，设φi,j的坐标为(xi,j,yi,j)，对浮动图像上的强度为f的点(xf,yf)，经过函数T(x,y)变形后，坐标为(T(xf),T(yf))，对应的参考图像上的点的强度值为r，则当前待配准浮动图像块与对应参考图像块的联合直方图可表示为：
(9)
其中w(x,y)为二维空间权函数，即本方案所选用的二维3阶B样条曲线，其一维形式在[-2,2]上定义为：
(10)
假设待配准图像总的像素点数为N，则当前局部互信息对整幅图像互信息的贡献率：
(11)
每个控制点对不同距离像素点的影响及本文方案使用的加权方案如图2所示。

图2 加权方案
1.4 参数优化方案选择和梯度信息计算
定义惩罚函数C=-I(R,F)，图像配准的过程实际就是函数C对控制点参数Φ寻优的过程，对于该过程，本文提出的方案在实现过程中选用的是对大规模参数优化问题效果比较好的LBFGS优化算法[12，13]。

LBFGS算法是一种基于有限存储单元的类牛顿算法，其名称中的BFGS是其四位发明者姓氏首字母的组合。

与传统的牛顿法相比，这种方法不直接计算当前变量的黑塞矩阵，而是根据其前有限次迭代得到的黑塞矩阵和当前的梯度信息进行构造，对一些黑塞矩阵计算比较困难或者变形参数密度较大的优化情形效果很好。

本文方案实现过程中使用的梯度信息是采用差分的方法计算得到，其计算公式为：▽
(12)
得益于基于B样条自由变换模型的局部形变特性，互不影响的参考点的梯度信息可以同时计算，从而提高梯度计算效率。

1.5 算法流程图
图3 本文配准方案的算法流程图
2 配准实验设计
本文方案在Windows环境下基于Matlab7.8.0平台编程实现，分别对标准图像数据和医学图像数据加入特定形变之后进行配准实验，形变采用文献[14，15]中
介绍的方法。

为了提高配准精度和效率，方便实验的进行，方案实现过程中对方案做了以下处理：1)网格划分时，采用x和y方向间隔相同的均匀网格；2)计算控制点梯度信息时，相互之间无直接影响的控制点同时处理；3)基于本文方案的在进行直方图统计前，对待配准图像的强度值进行二次采样，将8位的灰度图像采样为7位；4)在配准过程开始前，根据图像控制点间隔和图像像素点数，生成一个二维3阶B样条权重表，在直方图加权统计时重复使用，以提高配准效率。

本文首先使用标准图像数据对方案效果进行验证，同时增加了基于全局互信息的方案作为对比项，实验结果如图4所示。

图4 标准图像实验结果
基于本文提出的方案，在配准后的图像的差值图像上几乎看不到原图像的任何结构，相比于未配准的图像和基于全局互信息配准方法的结果，效果有明显的提升。

为了更加科学地评价配准效果，本文还选用均方根误差：
(13)
和图像相干系数：
(14)
量化评价图像配准的效果[6]。

对标准图像配准效果的SSD和CC计算结果见表1
所示。

表1 标准图像配准效果评价未配准全局MI配准本文方案配准
SSD0.03190.02510.0150CC0.82710.92560.9878
在进行医学图像数据实验时，采用的是人体左肾的磁共振图像，为了更加清晰地看到本文方案的效果，增加了只进行强度二次采样的全局互信息配准方案和二次采样后分块计算互信息的方案作为对比项。

配准结果如图5所示(为更加清晰地看到配
准效果，对四组配准后的差值图像进行了增强处理)。

图5 左肾磁共振图像配准结果
与图5对应的均方根误差和相干系数计算结果见表2所示。

表2 左肾图像配准效果评价未配准全局MI全局MI+强度采样分块MI+强度采样
本文方案
SSD0.02660.01220.01090.00940.0083CC0.83870.96990.97660.98580.9911
两组数据配准后的差值图像和数据评价结果显示，本文提出的方案对普通图像和磁共振图像的配准效果要优于传统的基于互信息的方法和基于分块互信息的方法。

3 实验结果分析和讨论
本文针对基于B样条的非刚性图像形变模型，考虑控制点对与其距离不同的像素
点的影响程度不同，提出了相适应的分块空间加权互信息方案，并在设计的实验中取得了很好的配准效果，本文方案是一个值得肯定的研究方案。

根据空间加权互信息方法的原理和一些实验结果判断，控制点间隔的设置对不同形变程度的图像配准效果是不同的，实验中发现，最大形变程度小于控制点间隔的一半时，图像配准都能取得很好的效果，因此，在具体应用中应根据实际情况设定控制点间隔的大小，对一些形变程度比较大的图像进行配准时，还可以考虑分层次配准的方法，即先设置较大的控制点间隔，实现图像的粗配准，然后细分控制点网格，实现精细配准。

在配准效率方面，由于采用生成空间权重表的形式，空间加权计算方面并不会对配准效率产生较大的影响，但在进行图像配准实验时发现，配准过程需要的时间依然要超过传统的全局互信息配准过程，这主要是因为要重复计算大量相互重叠的图像块的互信息和每个控制点的差分梯度。

然而我们还注意到，每个图像块的互信息和每个控制点的梯度信息在计算时是互不影响的，因此在进一步的研究过程中，尤其是在基于硬件平台实现时，可以优先采用并行计算的方式，同时计算各个图像块的
互信息和每个控制点的梯度信息，以提高计算效率。

4 结语
针对医学图像的非刚性配准，经过对相关理论的深入研究，本文以自己的理解在传统的基于互信息的方法中加入空间信息，提出了一个基于空间加权互信息的分块图像配准方案。

实验结果表明，本文提出方案能很好地完成一般医学图像非刚性配准的任务，且具有好高的配准精度。

同时注意到，本文只是一个初步的研究和验证，在后续的研究中，使用不同的空间插值方式进行直方图统计、算法的并行处理、分层配准方法的使用、梯度计算方法的研究、本文方案在多模态配准和超声图像配准中的应用等很多课题都值得深入地探讨。

参考文献
[1] Zhao Yi,Soatto S.Non-rigid registration combining global and local statistics[C]. IEEE-Computer-Society Conference on Computer Vision and Pattern Recognition Workshops,[s.l.]:IEEE,2009:2200-2207.
[2] Studholme C,Drapaca C,Iordanova B,et al.Deformation-based mapping of volume change from serial brain MRI in the presence of local tissue contrast change [J].IEEE Transactions on Medical Imaging,2006,25(5):626-639.
[3] D Loe Loeckx D,Slagmolen P,Vandermeulen D,et al.Non-rigid image registration using conditional mutual information [J].IEEE Transactions on Medical Imaging,2010,29(1):19-29.
[4] Zhuang Xiahai,Arridge S,Hawkes D J,et al.A non-rigid registration framework using spatially encoded mutual information and free-form deformations [J].IEEE Transactions on Medical Imaging,2011,30(10):1819-
1828.
[5] Lu Xuesong,Zhang Su,Su He,et al.Mutual information based multimodal image registration using a novel joint histogram estimation
[J].Computerized Medical Imaging and Graphics,2008,32(3):202-209. [6] Rueckert D,Sonoda L I,Hayes C,et al.Non-rigid registration using free-form deformations:application to breast MR images [J].IEEE Transactions on Medical Imaging,1999,18(8):712-721.
[7] Tustison N J,Avants B B,Gee J C.Directly manipulated free-form deformation image registration [J].IEEE Transactions on Image Processing,2009,18(3):624-635.
[8] Bookstein F L.Principal warps:Thin-plate splines and the decomposition of deformations [J].IEEE Transactions on Pattern Analysis and Machine Intelligence,1989,11(6):567-589.
[9] Davis M H,Khotanzad A,Flamig D P,et al.A physics-based coordinate transformation for 3-D image matching [J].IEEE Transactions on Medical Imaging,1997,16(3):317-328.
[10] Pluim J,Maintz J,Viergever M.Mutual Information Based Registration of Medical Images:A Survey [J].IEEE Transactions on Medical
Imaging,2003,22(8):986-1004.
[11] Studholme C,Hill D L,Hawkes D J.An overlap invariant entropy measure of 3D medical image alignment [J].Pattern Recognition,1998,32(1):71-86.
[12] Byrd R H,Lu Peihua,Nocedal J,et al.A limited memory algorithm for bound constrained optimization [J].SIAM Journal on Scientific Computing,1995,16(5):1190-1208.
[13] Klein S,Staring M,Pluim J.Evaluation of optimization methods for non-rigid medical image registration using mutual information and B-splines [J].IEEE Transactions on Image Processing,2007,16(12):2879-2890. [14] Lu Xuesong,Suo Shiteng,Liu Haihua,et al.Three-dimensional multimodal image non-rigid registration and fusion in a high intensity focused ultrasound system [J].Computer Aided Surgery,2012,17(1):1-12.
[15] Lu Xuesong,Zhang Su,Zha Yunfei,et al.Three-dimensional non-rigid registration and fusion for image-guided surgery navigation
system[C].IEEE:International Conference of Medical Image Analysis and Clinical Application (MIACA),[S.l.]:IEEE,2010:135-138.。