基于概率和引力优化模型的医学图像配准

裘意娜 李均利金林鹏

(宁波大学信息科学与工程学院，宁波 315211)

基于概率和引力优化模型的医学图像配准

裘意娜 李均利*金林鹏

(宁波大学信息科学与工程学院，宁波 315211)

基于互信息的配准方法，其目标函数经常存在许多局部极值，给配准的优化过程带来很大困难。提出一种基于概率模型的引力优化算法，在空间中随机构造参考物体与浮动物体，根据牛顿万有引力定律，搜索空间中质量最大的物体。利用该算法，实现以归一化互信息为相似性测度的医学图像配准实验。实验结果表明，这种方法能够有效地克服互信息的局部极值，在配准精度、配准时间和抗噪性方面都有较好的性能。

概率和引力优化算法;医学图像配准;互信息

Abstract:There are lots of local maximums in image registration based on mutual information，which obstruct optimization in registration process.In this paper，a new optimization algorithm，called probability and gravity optimization，was proposed.We constructed reference objects and floating objects in space，each object was located randomly，then searched the object whose quality was the heaviest according to Newton′s law of universal gravitation in the whole space.The new method was applied to medical image registration based on normalized mutual information.Experimental results showed that this registration method could efficiently restrain local maxima of mutual information function and had better performance at registration accuracy，registration rate and noise immunity.

Key words:probability and gravity optimization;medical image registration;mutual information

引言

随着医学影像学的飞速发展，CT、MR、PET等多种模态的图像已经广泛应用于临床。多模态医学图像配准在临床诊断、放疗定位、手术导航等方面有着重要的应用。当一种模态的医学图像所提供的信息不能满足需要时，可以把多种模态医学图像融合起来，提供互补的信息。最大互信息配准是目前研究较多、应用较广的一种基于像素灰度的配准方法。在 1995年被 Collignon等［1］和 Viola等［2］应用到医学图像配准后，其有效性得到广泛认可［3-6］。但是，这种方法运算量大、配准时间长，所以寻求高效的优化算法是提高配准质量、使其有效应用于医学图像配准的关键。

一些学者在优化算法方面提出了许多可行的方法。Collignon等和Fei等采用了Powell多参数优化方法，也是基于互信息的图像配准中最常用的优化算法，由于无需计算梯度，因此配准时间短，但其优化结果的好坏与初始点密切相关［7-8］;Slomka等和Radau等采用了单纯形法，是一种多维直接搜索的局部优化方法，其缺点是收敛速度过慢［9-10］;Matsopoulos等在视网膜图像的配准实验中，分别采用了模拟退火算法、遗传算法和单纯形法进行优化［11］。Wu等在配准中采用的是简单遗传算法［12］，而Silva等和Gefen等采用的是改进后的遗传算法［13-14］。Wachowiak 等提出了采用粒子群优化算法的配准方法［15］和采用禁忌搜索优化算法的配准方法［16］。这些算法各有优点，但同时也存在着容易受到局部极值干扰或“过早收敛”等不足之处，导致错误的配准结果。此外，很多学者将多种优化策略结合起来，Xu和Dony把进化策略引入到 Powell法的优化中［17］，葛培明等将 Wells提出的快速梯度计算法融入遗传算法中［18］，杨帆等人将蚁群算法和Powell法结合起来对三维图像进行了配准［19］，冯林等将粒子群算法和Powell法结合起来应用于单模和多模图像配准［20］。因而，在整个配准过程中，如何减少计算量、如何高效准确地得到配准结果，是目前研究的重点。

万有引力是自然界中最普遍的自然现象，存在于一切物体之中。最近已有一些学者把万有引力定律应用于经济、管理领域，取得了很好的效果［21-22］。笔者建立了基于概率和引力优化模型的算法，并将其应用到医学图像配准中，通过寻找最优的图像平移和旋转变换参数来配准两幅图像。

1 基于互信息的图像配准原理

对于图像配准，如何确定图像间的配准程度是一个很重要的问题。互信息是医学图像配准中常用的相似性测度，是信息论中的一个基本概念，用来描述两个随机变量间的统计相关性，是一个变量包含另一个变量的信息度量，用熵描述为

式中，H(A)和 H(B)分别为图像 A和 B的熵，H(A，B)为二者的联合熵。

如果两幅图像完全配准，那么互信息从理论上就会取得最大值。由于互信息对重叠区域的变化比较敏感，容易产生误配准，所以Studholme等人建议用归一化的互信息［23］。本研究采用如下定义的归一化互信息作为相似性测度，有

基于互信息的配准过程是一个多参数的优化过程，即搜索使两幅图像间的互信息最大的空间变换的过程，找到最优的参数。

2 配准变换模型

2.1 2D刚体变换

对待配准的两幅图像，选择一幅作为参考图像，另一幅为浮动图像。从浮动图像的空间坐标到参考图像的空间坐标，变换公式如下:式中，X=(x，y)是像素的空间位置，A是2×2的旋转矩阵，b是2×1的平移向量。在2D刚体变换中，一般包含3个变换参数，即沿 x轴，沿 y轴的平移量和绕图像中心的旋转角度。

2.2 3D刚体变换

从浮动图像的空间坐标PF到参考图像的空间坐标RF的刚体变换，用下式描述为

式中:VF、VR为3×3的对角阵，分别代表参考图像和浮动图像3个轴向上的像素大小;CF、CR分别是两图像的中心坐标;Rx、Ry、Rz是3×3的旋转矩阵，分别表示绕3个轴的旋转量，t是平移向量。

3 验证方法

采用的图像数据来源于美国Vanderbilt大学的Retrospective Registration Evaluation Project(RREP)项目，其中有一套用于研究人员进行算法初步评估的病人数据Practice组，包括一个病人的1套CT数据、6套MR数据(依次为PD、T1、T2和分别进行了几何失真校正的 PD_rectified、T1_rectified、T2_rectified)和1套PET数据。

3.1 二维图像

分别选取Practice病人3D断层图像进行2D配准实验，人为地平移旋转变换(沿x轴平移-5个像素，y轴平移4个像素，绕图像中心旋转3°)。各类医学图像成像设备原理不同，图像形成过程也不同，因此引入噪声的统计特性也不同。CT图像噪声服从高斯(Gaussian)分布，MR图像噪声服从莱斯(Rician)分布［25］。因此，在配准中给 CT图像加均值方差为(0.08，0.08)的高斯噪声，给 MR(6种)图像加了均值方差为(0.1，0.003)的莱斯噪声，如图1所示。把变化后的图像作为参考图像、源图像作为浮动图像进行2D刚体配准实验，分别显示了高斯噪声图像和莱斯噪声图像。

图2中的曲线表示在归一化互信息为测度的情况下，不同空间变换下的配准函数曲线图(Tx，Ty，Rotate分别表示沿x轴平移，沿 y轴平移，和绕图像中心旋转，归一化互信息值无单位)，根据平移旋转变换参数，测度最大值应该分别出现在-5，4，3。显然，该函数具有恒正、具有全局最大值，故引力优化中的质量可直接用互信息函数目标值表示。

图1 图像及其变换后图像。(a)原图像;(b)变换后的图像;(c)加高斯噪声的原图像;(d)变换后的高斯噪声图像;(e)加莱斯噪声的原图像;(f)变换后的莱斯噪声图像Fig.1 The original image and transformed image.(a)originalimage;(b)transformed image;(c)originalimagecontaminated with Gaussian noise;(d)transformed Gaussian noised image;(e)original image contaminated with Rician noise;(f)transformed Rician noised image.

3.2 三维图像配准

选取Practice病人3D断层CT图像和MR图像进行配准实验，RREP项目中通过基于标准点的配准方法已经得到一个近似的标准结果，通常把该结果作为评估的“金标准”，依据这个配准标准进行误差统计。以MR图像为配准参考图像、CT图像为配准浮动图像，同时采用粒子群算法、遗传算法和Powell算法对这些图像进行配准，图3显示了三维图像在6个参数变换下的归一化互信息值。

4 实验结果分析

4.1 配准精度和配准时间

在医学图像配准中，准确性是关注的重点。由于互信息测度函数计算的复杂性，使得如何提高算法的效率、尽可能快地得到配准结果也成为研究的一个重点。

本研究选用粒子群算法(PSO)［20］、遗传算法(GA)［19］、Powell算法与本研究提出的算法进行比较。文献［19］和［20］均采用的是混合算法，而本研究提出的是单一的优化算法，为公平起见，采用文献［19］和［20］中去除 Powell算法的粒子群算法和遗传算法，并且迭代次数都设为100。每个算法分别做出15组CT图、90组MR图和15组PET图，总共进行480次实验。

对沿x轴、y轴平移、绕中心旋转的配准误差和运行时间进行比较，如图4所示。为了方便，配准误差采用像素距离，距离单位为像素。配准误差小于一个像素，即达到亚像素级，被认为是配准成功;否则，成为误配准。

如果两幅图像完全配准，那么从理论上讲互信息就会取得最大值。从图4(a)～(c)中可以看出，x轴和y轴的配准误差均在1个像素之内，旋转角度在1°内，说明以上4种算法都成功配准图像，并且本研究提出的算法在配准精度方面优于遗传算法和粒子群算法，与Powell算法没有明显区别。图4(d)显示了运行时间的比较，Powell算法的执行时间最短，其次是本算法，平均在10 s左右，遗传算法对单模配准的时间性能很差。

4.2 抗噪性

在医学成像的特定条件下，考虑到放射线、示踪剂和强磁场等对人体的影响，加上成像模式本身的一些物理限制，常导致图像不清晰，并伴有噪声，所以研究算法对噪声的鲁棒性也是算法性能的重要方面。

图2(d)～(i)分别显示了高斯噪声图像和莱斯噪声图像在3种变换下的归一化互信息值，可见对加噪图像的配准是一个存在多峰值的函数优化过程，对配准算法的健壮性是一个很大的考验。分别使用本算法、PSO算法、Powell算法和GA算法进行基于归一化互信息的配准。结果如图5所示。

图2 不同类型图像在3种变换下的配准函数图(每行从左至右分别为沿x轴作平移变换、沿y轴作平移变换和绕图像中心作旋转变化)。(a)～(c)原图像;(d)～(f)高斯噪声图像;(g)～(i)莱斯噪声图像Fig.2 The registration function image by three kinds of transformation(In each line from left to right is translation along x-axis，y-axis and rotation around the center of image respectively).(a)～ (c)original image;(d)～(f)Gaussian noised image;(g)～(i)Rician noised image

图3 三维图像的配准函数图。(a)～(c)沿x轴、y轴和z轴作平移变换;(d)～(f)绕x轴、y轴和z轴旋转Fig.3 The registration function in aligning image with 3D image.(a)～ (c)translation along x-axis，y-axis and z-axis，respectively;(d)～ (f)rotation around x-axis，y-axis and z-axis，respectively

由于函数的多极性，虽然Powell算法所使用的时间均少于其他两种算法，但计算结果的准确性极差，本算法的准确性较PSO算法、遗传算法和Powell算法均有明显的提高。为了更直观地比较各算法的配准情况，验证算法的有效性，给每种图像做15个配准实验，共做了420个配准实验。图6给出了各算法误配准率的柱形图。可以发现本算法的误配准率比其他3种算法都少，其中Powell算法在对于噪声图像配准中的性能极差，其误配准率都是100%。实验证明，所提出的算法在克服局部极值方面有较大的优势。

图4 二维图像配准结果。(a)x轴平移误差的绝对值;(b)y轴平移误差的绝对值;(c)绕中心旋转角度误差的绝对值;(d)配准运行时间Fig.4 The results of 2D image registration:(a)the absolute value of error by along x-axis translation;(b)the absolute value of error by along y-axis translation;(c)the absolute value of error by around the center of image rotation;(d)the running time of registration

图5 不同算法配准噪声图像的误配准率Fig.5 The rate of mistake in noised image registration by different algorithms

4.3 3D图像的配准

把配准后浮动图像的8个顶点qi和标准结果的8个顶点pi按下面的方法对结果的准确性进行评估，有

将上述4种方法得到的配准误差对比显示在图8中。从图6中可以看出，本算法、PSO算法、遗传算法和 Powell算法的配准成功率分别为100%，80%，60%，20%。本算法的配准结果都达到了亚像素精度的要求，并且配准结果的精度普遍优于PSO算法和Powell算法。同时，比较上述3个算法的配准时间，GA算法配准时间最长，平均在10 000 s左右，其次是PSO算法，运行时间在4 000 s左右，而本算法的运行时间在2 000 s左右，比PSO算法少了近一半。虽然Powell算法在时间方面均少于其他两个算法，但结合配准精度，其成功率远远低于引力算法和PSO算法，因此，从配准时间和配准精度两方面综合比较，所提出的优化算法更胜一筹。

图6 图6 三维配准结果。(a)配准误差;(b)配准运行时间Fig.6 Results of 3D registration:(a)registration error;(b)the running time of registration

5 结论

对基于互信息的医学图像配准中的配准算法进行实验研究，分别进行了图像配准、加噪图像配准和三维图像配准的实验。结果显示，本算法获得了很好的效果，配准结果精度达到了亚像素级水平，普遍优于PSO算法、遗传算法和 Powell算法，具有精度高、鲁棒性强的优点，并且其配准时间大大缩短，同时又能使算法收敛到全局极小。由于本算法有很好的全局优化性能和时间性能，因此可用于二维、三维的医学图像配准。

附录:基于概率和引力模型的优化算法及收敛性证明

算法描述

经典物理学中的万有引力定律认为:宇宙中任何两物体之间都存在着相互吸引力，其中任一物体所受的引力大小与两物体质量的乘积成正比，而与两物间距离的平方成反比，用公式表示为

式中:G为万有引力常数;M1、M2分别为两物体的质量;d为它们之间的距离。

将万有引力定律应用到优化算法，假设优化的问题为

式中:g(X)代表目标函数，X=(x1，x2，…，xn)T∈Rn，D={X|Ii≤xi≤Ui，i=1，2，…，n}。定义 D 为搜索空间，在搜索空间中把每一个点看成是一个物体，X表示物体的坐标，m=F′(g(X))代表物体的质量，其中F′(X)是一维尺度变换函数，严格单调递增，并且对 ∀X∈ D，F′(X)＞0。

基于概率和引力优化模型中定义n维解空间中有两类物体，它们分别是参考物体群G—={Gi，1≤i≤n1}和浮动物体群 N—={Ni，1≤i≤n2}，并且 n1＜n2，定义 Gi(m)和 Gki(p)分别为参考物体Gi的质量和在第k维的位置，Ni(m)和Nki(p)分别为浮动物体Ni的质量和在第k维的位置，这里1≤k≤n。基于概率和引力优化算法的计算步骤如下:

步骤1:随机产生参考物体群G—和浮动物体群N—，令iter=0;比较所有物体的质量(归一化互信息值)，将质量较大的n1个物体保存到参考物体群中，并保持如下关系:

步骤2:参考物体Gi与Gj之间在搜索空间内按照下式构造新物体new，即

式中，C为常数，rand表示在(0，1)区间内产生的一个随机数。求出新物体的质量，然后淘汰参考物体群和新物体中质量最小的物体。

步骤3:计算每个参考物体Gi(1≤i≤n1)和浮动物体Nj(1≤j≤n2)之间的引力测度Fij和距离测度rij，即

对于每个浮动物体Nj(1≤j≤n2)选择和其引力最大的参考物体Gsj，有

对于每个有浮动物体吸附的参考物体Gi(1≤i≤n1)，选择对它吸引力最小的浮动物体Nwi，用变异算子在解空间内进行随机赋值，有

对于未被淘汰掉的浮动物体Nj(1≤j≤ n2)，即j≠ wi，则 Nj和 Gsj在搜索空间内构造新物体new，构造方法同步骤 2。用新物体替换 Nj，即。对于被淘汰的浮动物体Nwi，则用变异算子在解空间内进行随机赋值。求出更新过后浮动物体的质量，iter=iter+1。

步骤4:比较更新过后的参考物体群和浮动物体群的质量，将质量较大的n1个物体保存到参考物体中，此时亦满足式(7)。

步骤5:若满足终止条件则停机，否则转步骤2。

收敛性证明

基于概率和引力的优化算法是一种随机优化算法，关于随机优化算法以概率1收敛于全局最优解的条件，Solis和 Wets对其进行了证明［24］，其主要结论如下:

假设 1［24］若 f(D(z，ζ))≤ f(z)，ζ∈ S，则f(D(z，ζ))≤ f(ζ)

式中，D为产生问题解的函数，ζ为从概率空间(Rn，B，μk)产生的随机向量，f为目标函数，S为Rn的子集，表示问题的约束空间，μk为B上的概率度量，B为Rn子集的σ域。

假设 2［24］若对 S的任意 Borel子集 A，有ν(A)＞ 0，则

式中，ν(A)为子集A的n维闭包，μk(A)为由μk产生A的概率。

定理［24］设f为一可测函数，S为Rn的一可测子集，{zk}∞0为随机算法产生的解序列，则当满足假设1和假设2时，有，其中 Rε为全局最优点的集合。

收敛性证明

由上述结论可知:对于基于概率和引力优化算法，只要能够满足假设1和假设2，根据定理就可保证其以概率1收敛于全局最优解。下面是对该算法的寻优分析。

在基于概率和引力优化算法中，其返回的结果是第t次迭代前的参考物体XG1的位置，X(t)为第t次迭代的某个物体位置，f(X)是物体的适应度函数。因此，令假设1中的函数D为

可以看出，假设1是满足的。

对于假设2，只需证明规模为S的群体的样本空间的并集包含S，即

式中，Mi，t是第i次迭代时物体样本空间的支撑集。设迭代N次，其中第i次物体的范围为Si，即其支撑集。所以，群体的空间并集为。由于物体的位置范围是可以调节的，当其覆盖范围是从其所在位置到问题域的边界时，尽管边界所在物体很少，但却可以使得，即假设2成立。

综上所述，由定理可以基于概率和引力优化算法、以概率1收敛于全局最优解。

［1］Collignon A，Maces F，Delaere D，etal.Automated multimodality medical image registration using information theory［A］.In:Bizais Y，Barillot C，and Palola RD，Eds.Information Processing in Medical Image［C］.Dordrecht:Kluwer Academic Publishers，1995.263-274.

［2］Viola P，Wells WM.Alignment by maximization of mutual information［A］.In:Crimson E，Shafer S，Blake A，Eds.Proceedings of the 5thinternational conference on computer vision［C］.Washington D.C.:IEEE Computer Society，1995.16-23.

［3］Maes F，Collignon A，Vandermeulen D，et al.Mutimodality image registration by maximization of mutual information［J］.IEEE Transactions on Medical Imaging，1997，16(2):187-198.

［4］Maes F.Segmentation and registration of multimodal medical images:From theory，implementation and validation to a useful tool in clinical practice［D］.Leuven:Catholic University，1998.

［5］Studholme C，Hill DLG，Hawker DJ.An overlap invariant entropy measure of 3D medical image alignment［J］.Pattern Recognition，1999，32(1):71-86.

［6］Pluim JPW，Maintz JBA，Viergever MA.Mutual information based registration of medical images:a survey［J］.IEEE Transactions on Medical Imaging，2003，22(8):986-1004.

［7］Collignon A.Muti-modality medicalimage registration by maximization of mutual information［D］.Leuven，Belgium:Catholic University of Leuven，1998.

［8］Fei B，Wheaton A，Lee Z，et al.Automatic MR Volume registration and its evaluation for the pelvis and prostate［J］.Physics in Medicine and Biology，2002，47(6):823-838.

［9］Slomka PJ，Mandel J.Downey D，et al.Evaluation of voxelbased registration of 3-D power Doppler ultrasound and 3-D magnetic resonance angiographic images of carotid arteries［J］.Ultrasound in Medicine and Biology，2001，27(7):945-955.

［10］Radau PE，Slomka PJ，Julin P，et al.Wahlund.Evaluation of linear registration algorithm for brain SPECT and the errors due to hypoperfusion lesions［J］.Medical Physics，2001，28(8):1660-1668.

［11］Matsopoulos GK，Mouravliansky NA，Delibasis KK，et al.Automatic retinal image registration scheme using global optimization techniques［J］.IEEE Transactions on Information Technology in Biomedicine，1993，3(1):47-60.

［12］Wu Z，Hartov A，Paulsen K D，et al.Multimodal image reregistration via mutual information to account for initial issue motion during image-guided neurosurgery［A］.In:Proceedings of the 26thAnnual International Conference of the IEEE EMBS，San Francisco，USA，IEEE Press，2004，1:1675-1678.

［13］Silva L，Bellon O R P，Boyer K L.Precision range image registration using a robust surface interpenetration measure and enhanced genetic algorithms［J］.IEEE Transcations on Pattern Analysis and Machine Intelligence，2005，27(5):762-776.

［14］Gefen S，Bertrand L，Kriyati N，et al.Localization of sections within the brain via 2D to 3D image registration［A］.In:Proceedings of 2005 IEEE International Conference on Acoustics，Speech and Signal Processing，Philadelphia，USA IEEE Press，2005，2:733-736.

［15］Wachowiak MP，Smolikova R，Zheng Y.An approach to multimodal biomedical image registration utilizing particle swarm optimization ［J］.IEEE Transactions on Evolutionary Computation，2004，8(3):289-301.

［16］Wachowiak MP，Smolikova R，Elmaphraby AS.Hybrid optimization for unltrasound and multimodal image registration［A］.In:Proceedings of IEEE Engineering in Medicine and Biology Annual Conference.Istanbul，Turkey，2001，3:2418-2421.

［17］Xu X，Dony RD.Differential evolution with powell′s direction set method in medical image registration［A］.In:Proceedings of the 2004 IEEE International Symposium on Biomedical Imaging:From Nano to Macro，Arlington，USA，IEEE Press，2004，1:732-735.

［18］葛培明，陈虬，李尧臣.混合遗传算法在医学图像配准中的应用［J］.中国生物医学工程学报，2007，26(3):326-331.

［19］杨帆，张汗灵.蚁群算法和 Powell法结合的多分辨率三维图像配准［J］.电子与信息学报，2007，29(3):622-625.

［20］冯林，严亮.PSO和Powell混合算法在医学图像配准中的应用研究［J］.北京生物医学工程，2005，24(1):8-12.

［21］顾丽敏，吴小俊.基于万有引力定律的人脸识别方法［J］.计算机工程，2008，34(11):181-183.

［22］李熙.基于万有引力的遥感影像分类［J］.遥感信息，2004(3):36-38.

［23］Studholme C，Hill DLG，Hawkes DJ.An overlap invariant entropy measure of 3D medical image alignment［J］.Pattern Recongnition，1999，32(1):71-86.

［24］SolisFJ，Wets RTB.Minimization by Random Search Techniques［J］.Mathematics of Operations Research，1981，06(1):19-30.

［25］Pierre Gravel，Gilles Beaudoin，Jacques A.De Guise.A Method for Modeling Noise in Medical Images［J］.IEEE Transactions on Medical Imaging，2004，23(1):1221-1223.

A Medical Image Registration Method Based on Probability and Gravity Optimization Model

QIU Yi-Na LI Jun-Li*JIN Lin-Peng
(School of Information Science and Engineering，Ningbo University，Ningbo 315211，China)

TP391

A

0258-8021(2010)03-0345-08

10.3969/j.issn.0258-8021.2010.03.005

2009-07-14，

2010-03-03

国家自然科学基金资助项目(60672072);浙江省自然科学基金(Y106505);宁波市自然科学基金(2009A610089);宁波大学王宽诚基金

*通讯作者。 E-mail:li.junli@vip.163.com