基于多目标交叉变异粒子群算法的多模医学图像非刚性配准

许鸿奎 江铭炎 杨明强（山东大学信息科学与工程学院，济南 5000）（山东建筑大学，济南 500）

引言

不同模态的医学图像提供不同的信息，如CT（computer-assisted tomography）和MRI（magnetic resonance image）图像提供清晰的解剖结构信息，而PET（positron emission tomography）和 fMRI（functional MRI）则提供各种功能信息。图像融合可以充分利用多模态医学图像的各种相互补充的信息，在医疗诊断、治疗方案的制定等方面发挥重要的作用。

图像融合的前提是图像配准，而以互信息为相似性测度的图像配准不需要对图像进行分割或其他预处理，特别适合于多模态医学图像的配准。基于互信息的配准最早由Maes提出［1］，基于互信息的非刚性配准最早由Meyer［2］采用薄板样条变换实现［2］；文献［3-5］则基于互信息采用B样条变换完成了非刚性配准，其中文献［5］相对［3］提高了配准速度，但是配准精度有所下降。虽然基于互信息的配准已经得到了很大发展和应用，但在提高配准精度上还面临许多挑战，主要问题在于传统的优化方法如梯度下降法、拟牛顿法（limited-memory Broyden-Fletcher-Goldfarb-Shanno，LBFGS）等很容易陷入互信息的局部极值，不能很好地搜索到最佳配准参数。对图像进行预滤波处理，改善互信息曲面，使其更为平滑，一定程度上可以提高配准精度。然而实验证明，使用不同的滤波器得到的配准结果差别很大，而不同的图像对滤波器的要求也不一样，这就限制了基于互信息的配准方法的应用。为此，文献［6-7］分别使用遗传算法、粒子群算法去全局搜索最大互信息以获得最佳刚性变换参数，取得了良好的结果；文献［8］使用海明窗滤波对图像进行预处理，再利用粒子群算法进行全局搜索。然而，这些方法都是针对刚性配准的。

基于互信息的非刚性配准（如基于B样条的配准），特别适合于具有局部差异的多模态医学图像之间的配准。相对刚性配准，更多的配准参数及其相互作用，使得互信息曲面更加复杂；各种因素产生的互信息局部极值对搜索策略提出了更高的要求，若采用传统的以及上述优化方法，难以得到理想的配准结果。针对这一问题，本研究提出一种新的配准参数的全局搜索方法，用于基于互信息的多模态医学图像的B样条配准。它以LBFGS的优化结果为参考来构造初始种群，初始种群更优良；采用粒子群算法进行多目标优化，结合交叉变异策略来增加种群的多样性，克服粒子群算法受初始值的选取等因素的影响易陷入局部最优的缺点。实验证明，不必对图像进行预平滑处理也能得到理想的配准精度，提高了互信息配准的鲁棒性。

1 基于互信息的图像配准

尽管多模医学图像的灰度并不相似，甚至差别很大，但对同一人体组织而言，对应像素点之间的灰度在统计学上并非独立而是相关的；而互信息用于描述两个系统间的统计相关性，或者一个系统中所包含的另一个系统中信息的多少，因而以互信息为测度进行配准是可行的［2］。

给定图像A和B，它们的互信息量定义为：

其中，H（A）和H（B）分别为图像A和B的熵，H（A，B）为二者的联合熵。互信息的计算可以通过两幅图像的联合直方图统计来获得：

式中，pAB（a，b）为两图像的灰度对（a，b）在两图像相同坐标位置的出现概率，pA（a）是图像A中像素取灰度值a的概率，pB（b）是图像B中像素取灰度值b的概率。

基于互信息的图像配准可以描述为：

式中，A为参考图像，B为浮动图像，T表示对浮动图像施加的某种变换，采用三次B样条变换。配准过程就是对互信息函数的最优求解过程。

1.1 基于B样条的图像变换

基于B样条自由变形的图像配准是通过操控覆盖在图像上的控制网格来实现的。B样条函数的局部特性使得一个控制点的移动只会影响它附近的像素点，而不是整个图像所有的像素点都发生移动，因而特别适合于具有局部变形的图像之间非刚性配准。

定义域Ω为包含图像所有像素的最小平面矩形域，用Φ表示由nx×ny个控制点Φij（0≤i≤nx，0≤j≤ny）所组成的网格，网格步长为δ。将网格Φ覆盖在平面矩形域Ω上，Ω上的一个点P（x，y）经过三次B样条变换后的位置用下式表示：

某点的位移由其周围的16个网格控制点及其位移决定。0≤u＜1，表示网格控制点对T贡献的权重，该权重基于该点到各网格顶点的距离。

基于互信息的B样条配准，就是找到一组最合适的网格控制点Φij对浮动图像进行变换，使得变换后的浮动图像与参考图像之间的互信息最大，而最优变换参数由搜索策略完成。

1.2 传统优化算法对于互信息配准尚存在的问题

传统的优化算法如牛顿法、梯度下降法具有理论完备、算法成熟、效率较高和结果稳定性较好等诸多优点，因而在实际应用过程中获得了广泛应用。但是它们对函数的连续性、导数的存在性都有相对较高的要求，采用梯度下降法、LBFGS法寻找上述最优配准参数，需要预先对图像进行高斯平滑处理，以减少互信息的局部极值点，否则无法收敛到最优，配准效果不理想。

互信息局部极值的产生原因是多方面的［9］，比如图像本身固有的特性，图像变换时插值的影响，图像重叠区域大小的变化，噪声的影响等等。抑制局部极值有多种方法，比如可以采用更高阶的插值方法，或对图像先进行图像分割等预处理，再进行配准。高阶的插值意味着计算量的增大，采用分割预处理与图像本身的特性有较大关系，不具有普遍适用性。对于基于B样条的非刚性配准，互信息局部极值的产生原因更加复杂，更多的配准参数及其相互作用使得互信息函数更加复杂。互信息函数的非凸性、不规则性使得传统优化算法在寻优过程中很容易陷入局部极值。实验证明，若不进行高斯平滑预处理，梯度下降法和LBFGS找不到最优配准参数；而不同的图像，不同阶次的B样条，甚至不同的B样条网格尺寸对高斯平滑滤波器的参数要求也都各不相同。这种对图像预处理的要求，影响了互信息配准的鲁棒性。为此，拟采用全局搜索的智能优化方法寻找最优配准参数。

2 粒子群搜索策略

从20世纪60年代开始，以遗传算法（genetic algorithm，GA）和粒子群优化方法（particle swarm optimization，PSO）为代表的智能算法得到了很大发展，其主要特点是全局随机搜索策略，不用考虑目标函数本身是否连续或者可微。而粒子群优化算法相对遗传算法又具有本质的并行性，具有更高的搜索效率。PSO算法起源于对一个简化社会模型的仿真，它和鸟类或鱼类的群集现象有明显的联系。作为种群进化方法，每个优化问题的解被认为是搜索空间的一个飞行的粒子，所有的粒子都有一个由目标函数决定的适应度值，每个粒子由一个速度决定其飞行方向和距离。PSO算法初始化为一群随机粒子，粒子们根据对个体和群体的飞行经验的综合分析来动态调整各自的飞行速度，在解空间中通过迭代寻找最优解。每一次迭代中，粒子通过跟踪两个“极值”来更新自己，一个是个体极值pbest，即粒子目前自身所找到的最优解；另一个是全局极值gbest，即整个种群目前找到的最优解。在基本粒子群算法中，速度和位置的更新按下面公式进行［10］：

式中，下标i代表第i个粒子，下标j代表速度或位置的第j维，上标k代表迭代次数；c1和c2是学习因子，通常c1、c2∈［0，4］；r1和r2是介于［0，1］之间的随机数；是粒子Pi在第j维的个体极值的坐标；是群体在第j维的全局极值的坐标。

粒子群算法这种全局搜索策略，不必考虑目标函数本身是否连续或者可微，因而用于基于互信息的多模医学图像非刚性配准时，就没有对图像进行平滑预处理的特别要求。实验表明，高斯平滑滤波对其收敛性影响很小，因而采用粒子群优化算法来搜索B样条配准参数。但是实验中也发现，即使采用带惯性权重因子的粒子群算法［11］，受初始值随机选取的影响，结果也不稳定，出现了早熟现象，配准结果不甚理想，为此本研究提出多目标交叉变异粒子群优化算法加以改进。

3 多目标交叉变异粒子群优化算法

随着基本粒子群算法的不断应用，出现了很多改进方法。除了上述引入惯性权重和压缩因子外，典型的还有自适应模糊粒子群算法、杂交微粒群优化算法、免疫粒子群算法、小生境粒子群算法［12］等等，在不同的应用领域取得了良好的效果。

针对基于互信息的多模态医学图像的非刚性配准，采用粒子群算法进行多目标优化［13］，结合交叉变异策略［14］来跳出局部极值，克服粒子群算法早熟的缺点；同时对种群的初始化加以改进，以LBFGS的优化结果作为参考值构造初始种群，这样构造的初始种群离目标函数的全局最优更近，陷入局部最优的机会大大减少。虽然这样会消耗一定的时间，但可以进一步提高配准精度。

3.1 多目标粒子群优化算法

多目标优化最早由意大利经济学家V.F.Pareto提出，通常一个多目标优化问题可以表示如下［13］：

式中，S⊂Rn称为可行解区域，x为决策变量，E=｛F（x）|x∈Rn｝为目标向量，fi（x）i=1，2，…n是第i个目标函数。在大多数情况下，各个目标函数可能是有冲突的，这就使得多目标优化问题不存在惟一的全局最优解，使所有目标函数同时最优。

由于误差的平方和（sum of square differences，SSD）也是配准图像的常用的相似性测度，采用互信息为第一目标函数f1（x），而SSD为第二目标函数f2（x）。采用下述交叉变异的方法，同时优化两个目标函数，以期增加种群的多样性。

3.2 子群间的交叉变异

在两个子群中按一定的杂交概率选取指定数量的粒子放入繁殖池中随机地两两进行杂交，产生相同数目的子代粒子。子代粒子的位置按式（9）和式（10）根据父代粒子的位置计算得到，子代粒子的速度按式（11）和（12）根据父代粒子的速度计算得到［15］。杂交后的子代粒子回归到各自的子群并代替父代粒子，从而保持种群的规模不变。

式中，x是粒子的位置矢量，v是粒子的速度矢量；Ci（x）和Pi（x）i=1，2是子代粒子和父代粒子的位置；Ci（v）和Pi（v）i=1，2是子代粒子和父代粒子的速度；p是杂交概率。在迭代过程中，后代不断地继承来自两个子群的父代的特征，然后又寻找使自己的适应度函数最优的位置。这种交叉变异机制，保证了种群的多样化，使得算法在解空间搜索的遍历性得到改善，可以有效地避免求解陷入局部最优。

3.3 算法流程

步骤1：子群划分：根据目标函数f1（x）和f2（x）把微粒群平均分成2个子群ChiledS1和ChiledS2。

步骤2：设置子群参数：设置粒子的杂交概率p，设置其它参数学习因子c1、c2，惯性参量ω、速度限制vmax等。

步骤3：以LBFGS的优化结果为参考，初始化两个子群ChiledS1和ChiledS2中各粒子的初始位置。

步骤4：按照惯性权重粒子群算法更新子群中粒子的速度和位置。

步骤5：子群间的交叉变异：从两个子群中按参数选择一定的粒子进行杂交，根据式（9）和式（10）以及式（11）和式（12）分别产生子代的速度和位置。

步骤6：满足迭代中止的条件，则停止迭代，否则返回步骤4继续迭代。

由两个目标函数的搜索结果取均值作为最终配准参数，然后对浮动图像进行B样条变换，完成和参考图像的配准。

4 实验结果与分析

脑部实验图像来自BrainWeb［16］。图1（a）是217像素×181像素的MR-T2图像，作为配准的参考图像；图1（b）是相对应的MR-PD图像，在左上方和右下方经过了局部变形，作为待配准的浮动图像。此变形是对两个网格控制点进行一定的移位产生的。

首先以互信息为相似性测度，采用LBFGS优化方法进行三次B样条变换实现配准，B样条变换的网格大小为32像素×32像素。正如前面1.2节分析，由于互信息曲面存在局部极值，配准前需要预滤波平滑处理。图1（c）是使用模板大小为20像素×20像素和尺度为5的高斯滤波器得到配准图像，表1列出了使用不同参数的高斯平滑滤波器得到的配准结果。可见，使用不同的滤波器配准的结果差别很大。LBFGS优化方法对图像预处理的依赖，影响了基于互信息配准方法的鲁棒性。

粒子群算法是一种全局优化方法，不需要图像预处理，也能得到较高的配准精度。图2是采用惯性权重粒子群优化方法的配准结果，配准前没有进行高斯滤波；其中，图2（a）是图1（a）和图1（b）配准前的差，图2（b）是图1（a）和图1（b）配准后的差。对20次配准的结果进行平均计算，互信息和误差的平均值分别是2.049和0.045，这是不错的配准结果。然而从图2（c）所示控制网格的变化可以看出仍有早熟的迹象，所得配准参数距离真正的网格点变化值仍有一定的差距。

图1 待配准图像及采用LBFGS优化的配准结果。（a）MR-T2参考图像；（b）MR-PD浮动图像；（c）配准后的浮动图像Fig.1 Images to be made registration and the result using LBFGS optimization.（a）the MR-T2 reference image；（b）the MR-PD float image；（c）the float image after registration

表1 不同参数的高斯平滑滤波器的配准结果Tab.1 Registration results using different filters

图2 采用惯性权重粒子群优化算法的配准结果。（a）配准前的差；（b）配准后的差；（c）控制网格的变化Fig.2 Registration result using inertia weight PSO.（a）the difference before registration；（b）the difference after registration；（c）the changes of control grid

图3 采用多目标交叉变异粒子群算法的配准结果。（a）配准后的浮动图像；（b）配准后的差；（c）控制网格的变化Fig.3 Registration result using proposed optimization in the paper.（a）the float image after registration；（b）the difference after registration；（c）the changes of control grid

图3是采用所提出的多目标交叉变异粒子群优化算法得到的配准结果；其中，图3（a）是配准后的输出图像，图3（b）是两图像配准后的差。对比图3（b）和图2（b），配准效果比惯性权重粒子群算法得到了明显改善。对20次配准后的结果进行平均计算，互信息和误差的平均值分别是2.075和0.019，相比惯性权重粒子群算法的2.049和0.045，提高了配准精度。图3（c）是控制网格的变化，对比图2（c）可以看出，所得配准参数更接近真正的网格点变化值。

将分别采用梯度下降法［3，5］、惯性权重粒子群优化方法以及本研究优化方法得到的配准结果数据列于表2加以比较，本方法的配准精度优于另外两种方法。

除了上述脑部MR-T2图像和MR-PD图像，从BrainWeb中又取了50个MR-T2图像及序列号对应的50个MR-PD图像进行配准实验。按前述相同的方法，人为地对MR-PD图像进行局部变形。假设将配准误差大于和等于0.02定义为误配准，依此统计各种方法的配准率，如表3所示。可见，本方法取得了较高的配准率。

为了进一步验证本算法的有效性，还采用其它模态的医学图像进行配准实验。图4是CT和CBCT（Cone beam computer tomography）图像的配准实验结果，其中图4（a）参考图像，图4（b）是浮动图像，图4（c）是采用改进的多目标交叉变异粒子群优化算法的配准输出，图4（g）是配准输出与参考图像的差；作为比较，给出其它方法配准后两图像的差，图4（d）是刚性配准后浮动图像与参考图像的差，图4（e）是采用LBFGS优化方法进行非刚性配准后的浮动图像与参考图像的差，图4（f）是采用惯性权重粒子群算法进行非刚性配准后的浮动图像与参考图像的差。很明显，采用本方法配准后的浮动图像与参考图像的差最小，配准效果最好。配准前，两图像的互信息及距离（像素值差的平法和，以下相同）分别是0.527和11 583，采用本方法配准后互信息及距离分别是1.082和7 155，配准效果显著。

表2 不同优化方法的配准结果比较Tab.2 Comparison of registration results using different optimizations

表3 不同优化方法的配准率比较Tab.3 Comparison of registration rate using different optimizations

图4 对CT和CBCT图像采用不同配准方法进行配准的结果比较。（a）CT参考图像；（b）CBCT浮动图像；（c）使用本文方法配准后的浮动图像；（d）刚性配准后的差；（e）LBFGS B样条配准后的差；（f）惯性权重PSO B样条配准后的差；（g）本方法配准后的差Fig.4 Comparison of CT and CBCT image registration results using different registration methods.（a）the CT reference image；（b）the CBCT float image；（c）the float image after registration using the method proposed in this paper；（d）the difference after rigid registration；（e）the difference after LBFGS B-spline registration；（f）the difference after B-spline registration with inertia weight PSO；（g）the difference after registration proposed in this paper

图5是CT和PET图像的配准实验结果。其中，图5（a）参考图像，图5（b）是浮动图像，图5（c）是采用改进的多目标交叉变异粒子群优化算法的配准输出。配准前，两图像的互信息及距离分别是0.884和48 254，采用本文方法配准后互信息及距离分别是1.004和47 527，也取得了不错的配准效果。

图5 CT和PET图像的配准。（a）CT参考图像；（b）PET浮动图像；（c）使用本方法配准后的浮动图像Fig.5 Registration results of CT and PET image using the method proposed in this paper.（a）the CT reference image；（b）the PET float image；（c）the float image after registration

然而，配准精度提高的代价是配准速度的下降。采用图1中的图像进行配准实验时，本方法相比惯性权重粒子群算法，配准后两图像的互信息改善了0.026，而配准时间增加了近10％，大约19s。这是因为多目标操作及交叉变异都需要额外占用一定的时间。

5 结论

基于多目标交叉变异粒子群算法，以互信息为相似性测度采用B样条变换对多模态医学图像进行非刚性配准，克服了传统优化方法如梯度下降法、LBFGS依赖于图像预平滑处理的缺点，提高了基于互信息的非刚性配准的鲁棒性；该方法以LBFGS的优化结果为参考构造的初始种群更优良，充分继承了传统优化方法和智能优化方法各自的优点，所构造的初始种群距离全局最优更近；针对基本粒子群算法易陷入局部最优的缺点，采用多目标优化方法引入交叉变异机制对粒子群算法加以改进，使得算法在解空间搜索的遍历性得到改善，优化结果更接近全局最优，从而减少了配准误差，提高了配准的精度。如何在减少配准误差的同时进一步提高配准的速度将是以后研究的重点。

［1］Maes F.Multimodality image registration by maximization of mutual information［J］.IEEE Trans on Medical Imaging，1997，16（2）：187-198.

［2］Meyer CR，Leichtman GS，Brunberg JA，et al.Demonstration of accuracy and clinical versatility of mutual information for automatic multimodality image fusion using affine and thin plane spline warped geometric deformations［J］.Medical Image Analysis，1997，1（3）：195-206.

［3］Rueckert D，Sonoda L，Hayes C，et al.Nonrigid registration using free-form deformations-application to breast MR images［J］.IEEE Trans on Medical Imaging，1999，18（8）：712-721.

［4］Kybic J.Fast parametric elastic image registration［J］.IEEE transaction on image processing，2003，12（11）：1427-1441.

［5］彭晓明，陈武凡，马茜.基于B样条的快速弹性图像配准方法［J］.计算机工程与应用，2006，11：186-189.

［6］Ingole VT，Deshmukh CN，Anjali Joshi，et al.Medical image registration using genetic algorithm［C］∥Ajith A，Norway.2009 2nd International Conference on Emerging Trends in Engineering and Technology.USA：IEEE computer Society，2009：63-66.

［7］Chen Yenwei，Aya Mimori.Hybrid particle swarm optimization for medical image registration［C］∥Hua Zhang.2009 Fifth International Conference on Natural Computation.Tianjin：Institute of Electrical and Electronics Engineers（IEEE），2009：26-30.

［8］裴继红，田剑豪，杨烜.基于海明窗滤波及粒子群优化搜索的医学图像配准［J］.生物医学工程学杂志，2007，24（02）：262-267.

［9］Ardizzone E，Gallea R.Effective and efficient interpolation for mutual information based multimodality elastic image registration［C］∥Takashi M.2009 IEEE 12th International Conference on Computer Vision Workshops.Washington：IEEE Computer Society，2009：376-381.

［10］Kennedy J，Eberhart R.Particle swarm optimization［C］∥Proceedings of IEEE International Conference on Neural Networks.New Jersey：IEEE Service Center，1995：1942-1948.

［11］Ememipour J，Nejad MS，Ebadzadeh MM，et.al.Introduce a new inertia weight for particle swarm optimization［C］∥Sungwon S.2009 Fourth International Conference on Computer Sciences and Convergence Information Technology.Seoul：IEEE Computer Society，2009：1650-1653.

［12］王万良，唐宇.微粒群算法的研究现状与展望［J］.浙江工业大学学报，2007，35（2）：136-141.

［13］Xiang Feng，Francic CM.A parallel evolutionary approach to multi-objective optimization［C］∥Kay CT.2007 IEEE Congress on Evolutionary Computation.Singapore：IEEE Computational Intelligence Society，2007：1199-1206.

［14］Coelho Ls，Krohling Ra.Predictive controller tuning using modified particle swarm optimization based on cauchy and gaussian distributions［C］∥Mario K.8th On-Line World Conference on Soft Computing in Industrial Applications.Berlin：Springer，2005：287-298.

［15］张敏慧.改进的粒子群计算智能算法及其多目标优化的应用研究［D］.杭州：浙江大学，2005.

［16］EVANS AC.Brain Web：online simulated brain database［DB/OL］.http://www.bic.mni.mcgill.ca/brainweb，2006-06-08/2007-03-01.