超声和MR/CT图像的非刚性配准技术

刘雪莉，宋志坚

（复旦大学数字医学研究中心，上海市医学图像处理与计算机辅助手术重点实验室，上海200032）

图像配准就是寻找待配准图像之间的变换关系，使图像达到空间上一致的过程。图像配准的目标函数可以表示为公式（1）[1]。其中，S代表原图像，T代表待配准目标图像，W代表图像的空间变换参数，M代表两幅图像的相似度，R为正则项。

超声成像实时性高、成本低、无辐射，临床上常用超声引导穿刺或肿瘤切除。但是，超声图像质量不高，很难正确的显示病变和定位病灶。MR有多种成像模式，图像质量高，利于准确的定位肿瘤，但MR引导手术成本高、实时性差。针对目前超声引导和MR引导手术存在的问题，研究人员提出了一种术前MR/CT图像和术中或术后超声图像融合引导手术的新方法，该方法不仅能够保证超声引导的实时性，而且因术前MR/CT图像的引入保证了引导准确性。术前MR/CT图像和术中实时超声图像融合引导手术需要解决的首要问题是多模态非刚性图像配准问题。

1 应用于医学图像非刚性配准技术归纳

非刚性图像配准可以分为如下三个基本过程：（1）构建形变模型；（2）确定目标函数；（3）选择优化算法，得到最优变换参数。

1.1 形变模型

用于非刚性配准的形变模型可以分别为四类：（1）基于物理模型的形变模型；（2）基于插值理论的形变模型；（3）基于已有知识的形变模型；（4）满足特殊约束的形变模型。

基于物理模型的形变模型主要有弹性体模型、粘性流体模型、扩散模型、曲率配准和微分同胚映射流动模型等；基于插值理论的形变模型主要有径向基函数模型、弹性体样条模型、自由形变模型和局部仿射变换模型等；基于已有知识的形变模型包括统计学约束模型和生物材料或生物物理模型等；满足特殊约束的形变模型包括拓扑特性约束模型、体积约束模型和刚性约束模型等。

1.2 目标函数

根据衡量图像匹配程度的标准的不同，可以将目标函数分为两类：（1）几何法，构建两幅图像标记点之间的相关性。标记点可以人工选取，也可以通过SIFT[2]、SURF[3]和GLOH[4]等检测算法得到。（2）图像法，通过灰度、梯度和像素信息计算两幅图像的相关性。

近年来，基于几何法的配准技术中比较常用的有GMMs（高斯混合模型）[5]和ICP（迭代最近点）算法及相关的改进，如EM-ICP、T分布模型等。

图像法是超声和MR/CT图像配准的一般的方法，可以分为两类：第一类是单模态图像配准，将MR/CT图像转化为模拟超声图像后进行单模态图像配准，或将两种模态图像转化为中间模态进行配准，如目前常用的特征描述子、熵图和拉普拉斯图法；另一类是多模态图像配准，目前常用的方法可以分为互信息及其改进方法和相关率及其改进方法两类。

1.3 优化方法

根据目标函数的不同，需要采用不同的优化方法。优化方法可以分为两类：第一类为连续优化，例如梯度下降算法、共轭梯度算法、powell算法、QN算法、LM算法和随机梯度下降算法等；第二类为离散优化，例如马尔可夫场、基于图的方法、BP法和LP法等。

2 超声和MR/CT图像配准技术

本章是本文的重点，总结近年来应用于超声和MR/CT图像的非刚性配准方法，归纳其原理、创新性和优缺点。本章主要分为三部分：第一部分是多模态图像转化为单模态图像的配准方法，基本思路是将多模态图像转化为单模态图像应用于超声和MR/CT图像配准；第二部分是多模态图像配准方法，主要分为互信息及其改进和相关率及其改进两类展开；第三部分是基于特征点和图像分割的多模态图像配准方法，目前腹部的超声和MR/CT图像配准主要采用这类方法。

2.1 多模态图像转化为单模态图像

在超声和MR/CT图像的配准中，将多模态图像转化为单模态图像有两种方法，一种是将MR/CT图像转化为模拟超声图像，另一种是将两种模态图像转化为中间模态图像。本章2.1.1节属于第一种方法，2.1.2和2.1.3节属于第二种方法。

2.1.1MR/CT图像模拟超声图像

以CT图像模拟超声图像为例。根据超声成像原理，超声信号通过不同密度介质时会产生折射、反射和衍射现象。模拟超声图像主要根据反射和衰减两个特性构建[8～11]。

超声波在两个声阻值不同的组织交界面通过时将会发生反射如图1。将组织的声阻定义为Z=ρc，其中ρ为组织密度，c为超声在组织中传播速度，不考虑骨骼和空气。根据超声的物理性质，超声在组织界面发生反射时的衰减比率如公式（2）：

为了进一步提高图像质量，研究人员通过对数压缩和组织图谱等方法对图像进行进一步修正。

由于图像处理部分在术前进行，该类方法可以降低术中或术后图像配准所需要的时间。但是，该方法是在假设的理想条件下进行，模拟出的图像不能体现组织内部的非均质性和散斑。

图 1超声入射组织界面示意图

2.1.2基于熵图和拉普拉斯图的图像配准方法

熵值计算来源于信息论，基于熵图的图像配准方法的基本思想如图2。在原始图像中取每个像素的邻域块，计算邻域块的熵（即该点的邻域结构信息）赋值该像素[13]，通过这种方法不同模态的图像转化为相同模态的熵图。计算熵值的图像块大小选择取决于配准的需求，图像块越小局部信息越多，图像块越大统计信息越充分。

图 2基于熵图的图像配准方法原理图

另一种结构描述的方法为拉普拉斯结构映射，即拉普拉斯图。原理如图3[14]，将RN空间中的k个点（a1,a2,a3, …ak）转化为高阶流形M，流形通过构建邻域图进行降维（low dimensional embedding），得到Rn（n≪N）空间中的点（b1,b2,b3, …bk），然后对降维后的图像进行配准。低维空间中点的距离可以表示为公式（5），其中，二阶可微函数m：M→Rn，ai→bi，表示降维前后的空间变换关系。

和熵图相比，拉普拉斯图在体现局部信息方面更有优势，并且计算过程中涉及流形变换、邻域图和空间降维，去除了图像中的离群值和错误点，使该方法的抗噪能力更强。但是，拉普拉斯图的变换计算复杂度比熵图高，且特征降维也损失了图像信息。

2.1.3基于特征描述的图像配准方法

这类方法基于自相似原理。由于同一成像体的不同模态图像的结构相同，因此，内部图像块之间的结构相似性能够反映图像的相似性。对比单一图像特征或者多个图像特征权重加和的方法，这类描述子的优势在于配准过程能够用到图像的点，边和纹理等多个信息，且对噪声不敏感。

图 3基于拉普拉斯图的图像配准方法原理图

文献[16]使用模态独立邻域描述子（MIND）描述结构信息，计算方法如公式（7）。其中，r代表邻域块，D代表邻域块之间的相似性。

相似性一般用SSD距离（Sum of Squared Differences）进行计算，也可以直接用滤波简化计算。计算得到图像中各个像素块的MIND后，仍然用其的SSD距离构建配准目标函数，如公式（8）：

文献[15]使用自相似上下文（SSC）对结构信息进行描述。比MIND交互性更好。SSC和MIND计算示意图如图4。另外，文献[15]改进了文献[16]的优化方法，采用汉明距离和离散优化。

图 4MIND和SSD图像块之间关系计算示意图

文献[17]采用梯度方向（GOA）描述图像信息，图像配准的目标是最小化两幅待配准图像的梯度夹角，目标函数如公式（9）：

该方法的优势在于梯度的方向能够反映几何边界的方向信息，因此，配准过程能够用到图像的结构信息。

基于特征描述的配准方法能够反映图像的结构和纹理信息，计算速度快，抗噪能力强。但是，基于图像块的运算丢失了部分局部信息，影响了配准精度。并且，该类方法适用于较小变形，随着变形程度的增加，配准误差增大。

2.2 基于改进互信息和相关率的配准方法

互信息（MI）和相关率（CR）是较早应用于图像配准尤其是多模态非刚性配准的方法之一。但是传统的MI和CR方法仅基于灰度信息的图像配准方法鲁棒性不高，不能反映图像的局部特征和结构特征，计算量大。因此，研究人员提出各种改进互信息和相关率的方法。

2.2.1互信息及其改进方法

传统的互信息（MI）方法应用于医学图像配准存在很多问题。因此，研究人员不断对互信息配准方法进行改进，由于传统互信息不能保证图像的重叠不变性，人们提出了标准互信息（NMI）[20]；为了使互信息能够反映局部信息，人们提出了局部互信息（LMI）[21]；为了在互信息计算中使用更多的图像信息，人们提出了多参数互信息（α-MI）等等。但是，这些方法都不能解决互信息计算过程中的非旋转不变性和结构相似度描述等问题。因此，基于上述方法，研究人员又提出了进一步的改进方法。

文献[19]对此进行改进，提出了一种自相似权重互信息（SeSaMI），该方法相对于α-MI做出了两点改进：第一，该方法使用了基于kNN图的α-MI，在目标函数中加入了局部信息；第二，该方法在特征距离的计算中加入了自相似权重，反映图像的结构特征。该方法配准的基本过程是：第一步，利用Pearson相关系数得到结构特征明显的区域；第二步，计算图像的自相似性，计算方法如图5。首先，将图像块表示为二维灰度直方图，这种直方图的表示形式保证了图像的旋转不变性。然后，用EMD距离（Earth Mover Distance）[23]计算灰度直方图的相似性（即直方图所表示的图像块之间的相似性）；第三步，计算基于k临近距离（kNN图）的α-MI，如公式（10）:

第五步，采用改进的随机梯度下降法[24]对目标函数进行优化。

文献[18]中的方法和文献[19]的方法类似，该方法是基于局部互信息（LMI）的改进，提出条件上下文互信息（CoCoMI）的配准方法，该方法相对于LMI作了两点改进。第一，保证配准过程的旋转不变性；第二，在目标函数中加入了结构信息。LMI的计算方法如公式（13），CoCoMI的计算方法如公式（14），CoCoMI相对LMI的改进如图3-7。

从公式（13）、（14）和图3-7可以看出两种方法的不同在于图像块的选取方式。LMI的计算过程是随机选取M个种子点，在每个种子点周围选取邻域，在邻域内随机选取N个点，计算两幅图像中以这N个点为中心的图像块互信息的平均值。而CoCoMI方法邻域内的N个点并不像LMI那样随机选取，选取的是和种子点结构相似的N个点，然后计算以它们为中心的图像块互信息（图6中，N=5）。结构相似点的衡量方法和SeSaMI的自相似权重的计算方法相同。

该节提到的两种改进互信息方法在保证图像信息的完整的同时，加入了局部信息和结构信息，进一步提高配准精度。文献[18]对各个配准方法的配准误差进行了对比，结果显示SeSaMI和CoCoMI比MIND等误差小，并且随着变形程度的增加，配准精度更稳定。但是，这类方法比其他图像配准方法计算复杂度高，用时长。

图5图像块之间的相似性计算

图 6LMI和CoCoMI比较

2.2.2相关率及其改进方法

互相关法是最基本的基于灰度统计的图像配准的方法，通常被用于进行模板匹配和模式识别。它通过计算模板图像和搜索窗口之间的互相关值确定匹配的程度，互相关最大时的搜索窗口位置决定了模板图像在待配准图像中的位置。图像X和图像Y的相关率（CR）计算如公式（15）[25]，式中参数计算如公式（16）、（17），推导过程参考文献[29]：

由公式可知，CR计算不需要像MI那样统计联合分布直方图，因此基于CR的图像配准方法的计算量小于MI。但是传统CR和MI一样，不能够反映图像的结构信息和局部信息。

文献[26]对传统CR方法进行了改进，提出了一种基于鲁棒性块相关率（RaPTOR）的配准方法。该方法存在三个方面的优势：首先，这种方法是对浮动图像X灰度分布直方图中的分布块进行处理而不是直接处理E(Y|X)，加入了灰度局部信息；其次，对块进行处理而非对点进行处理的方法降低了计算复杂度；再次，这种方法允许使用窗函数对直方图进行处理，减少了量化误差对于算法精度的影响。块相关率的计算方法如图3-7，如公式（18），公式中μj计算如公式（19）:

公式中，I表示超声图像，J表示MR/CT图像。该方法认为超声图像中某点的灰度不仅和MR图像中该点的灰度值相关，而且和该点的梯度值相关。因此，可以表示为MR图像在x点的灰度pi和梯度gi的线性组合，如公式（21）：

系数c=｛α,β,γ｝可以用最小二乘形式表示，如公式（22）：

图7RAPTOR方法中CR计算示意图

综合以上结果，基于的图像配准方法的相似度计算如公式（23）：

同时，文献基于的图像配准又做出了进一步的改进，和局部互信息相对于互信息的改进相似，文献计算局部图像块的相关率后求平均，在目标函数中加入图像局部信息。但是，该方法在配准过程中仅用到梯度信息，灰度信息和局部信息，相对于SeSaMI和RaPTOR等方法，该方法没有体现图像的更高层次的结构信息，配准精度没有SeSaMI和RaPTOR等方法高。但是该方法的优势在于比SeSaMI和RaPTOR等方法配准速度快，且比传统的MI和CR方法配准精度高。

2.3 基于特征点和图像分割的配准方法

本章2.1和2.2节所总结的超声和MR/CT图像的配准方法除LC2外，主要应用于脑部图像。应用于腹部图像的配准方法大多基于特征点和图像分割。

文献[6]将TPS-RPM算法和ICP算法结合应用于前列腺超声和MR图像配准，这是一种基于特征点的配准方法，特征点由医生手动选取，这种配准方法保证了ICP对于整体线性变换和TPS-RPM算法对于局部非线性变换的优势，取得了鲁棒性较高的配准结果。

文献[12]提出了一种利用肝脏的下腔静脉和表面在超声和MR图像中都较为清晰的特点，在配准过程中医生只需在下腔静脉处点击一次鼠标，通过点扩张分割出下腔静脉，通过马尔可夫随机场模型提取肝脏表面点云，然后利用PCA算法结合ICP算法进行配准，得到了快速、鲁棒性高的配准结果。

文献[30]是基于MR和US图像分割的配准方法，该方法的配准步骤如下：第一步，使用迭代优化算法通过已有前列腺图谱的变形和仿射变换分割前列腺，然后将分割得到的前列腺图像网格化；第二步，将前列腺MR图像网格化后的顶点作为前列腺US图像网格的初始顶点，构建能量函数，迭代分割出前列腺；第三步，构建线弹性能量函数作为配准的目标函数。这种配准方法在术前分割MR图像，节省了术中和术后图像配准时间。以分割后的MR图像作为对照，通过迭代优化，保证了超声图像分割的质量。

文献[31]是一种基于前列腺表面的配准方法， MR图像通过多分辨率的基于图分割的方法分割出前列腺，US图像是医生分割，然后用EM-ICP算法进行配准。该文献的主要优点体现在其MR图像的分割方法是采用关系图优化法。

但是，基于特征点和图像分割的配准技术首先要确定标记点或分割图像。由于超声图像质量不高，自动选取标记点或者图像分割的误差都比较较大，人工选取标记点或人工图像分割的实时性不高。

3 总结与思考

目前，超声和MR/CT图像配准并没有一种普适的方法。要实现图像的快速的配准，需要舍弃一部分的配准精度，如特征描述和法等。要实现较高的精度，就需要较长的配准时间，如SeSaMI和RaPTOR等。另外，特征描述子比熵图和拉普拉斯图的运算速度快是由于特征描述的改进；SSC相对于MIND配准速度的提高来自于优化方法的改进；SeSaMI、RaPTOR、和CoCoMI相对于传统MI和CR精度的提高，来自于更多的图像信息。因此，如何在配准方法中用到更多的图像信息，如何在提高配准精度的同时保证配准的实时性，是超声和MR/CT图像的非刚性配准领域仍需探索的问题。

[1] Sotiras A , Davatzikos C, Paragios N. Deformable Medical Image Registration: A Survey. Medical Imaging, IEEE Transactions on.2013, 32: 1153-1190.

[2] Ke Y , Sukthankar R. PCA-sift: A more distinctive representation for local image descriptors. Proc. Int. Conf. Comput. Vis. Pattern Recognit, 2003, 2: 506-513.

[3] Bay H , Ess A , Tuytelaars T ,et al.Speeded-up robust features(SURF). Comput. Vis. Image Understand., 2008, 110(3): 346-359.

[4] Mikolajczyk K, Schmid C . A performance evaluation of local descriptors. IEEE Trans. Pattern Anal. Mach. Intell, 2005, 27(10):1615-1630.

[5] Jian B, Vemuri B C. Robust point set registration using gaussian mixture models. IEEE Transactions on Pattern Analysis and Machine Intelligence, 2011, 33(8): 1633-1645.

[6] Makni N , Toumi I , Puech P ,et al.A non rigid registration and deformation algorithm for ultrasound & MR images to guide prostate cancer therapies. Engineering in Medicine and Biology Society (EMBC), 2010 Annual International Conference of the IEEE, 2010, 9: 3711-3714.

[7] Chui H , Win L, Schultz R ,et al.A unified non-rigid feature registration method for brain maing. Med. Image Anal, 2003, 7(2):113-130.

[8] Wolfgang W, Shelby B, Ali K,et al.Callstrom, Nassir Navab.Automatic CT-ultrasound registration for diagnostic imaging and image-guided intervention. Med. Image Anal, 2008, (12): 577-585.

[9] Kutter O , Shams R , Navab N. Visualization and GPU-accelerated simulation of medical ultrasound from CT images. Computer methods and programs in biomedicine, 2009, 94(3): 250-266.

[10] Li J X，Jun L，Wei W Z,et al.A novel algorithm for CT-ultrasound registration. Point-of-Care Healthcare Technologies(PHT), 2013 IEEE, 2013, 101(104): 16-18.

[11] Bernhard F, Wolfgang W, Markus M,et al.Automatic ultrasound–MRI registration for neurosurgery using the 2D and 3D LC2 Metric.Med, Image Anal, 2014, 18(8): 1312-1319.

[12] Young K H , Young T O, Jung B K,et al.One Click 3D Ultrasound to MR Registration, 2014 IEEE International Conference on Consumer Electronics.

[13] Christian W, Nassir N. Entropy and Laplacian images: Structural representations for multi-modal registration. Medical Image Analysis, 2012, 16(1): 1-17.

[14] Yu H, Xiang C F, George B. Local joint entropy based non-rigid multimodality image registration, Pattern Recognition Letters, 2013,34(12): 1405-1415.

[15] Mattias P H, Mark J. Towards Real time Multimodal Fusion for Image-Guided Interventions Using Self-similarities. Medical Image Computing and Computer-Assisted Intervention – MICCAI, 2013,(8149): 187-194.

[16] Mattias P H, Mark J, Manav B,et al.MIND: Modality independent neighbourhood descriptor for multi-modal deformable registration.Medical Image Analysis, 2012, 16(7): 1423-1435.

[17] Dante D N, Louis C D , Tal A. Fast and Robust Registration Based on Gradient Orientations: Case Study Matching Intra-operative Ultrasound to Pre-operative MRI in Neurosurgery. Proceedings of the Third international conference on Information Processing in Computer-Assisted Interventions (ICPIA), 2012(7330): 125-134.

[18] Rivaz H, Karimaghaloo Z, Fonov V S,et al.Nonrigid Registration of Ultrasound and MRI Using Contextual Conditioned Mutual Information. IEEE Transactions on Medical Imaging, 2014, 33(3):708-725.

[19] Hassan R, Zahra K, Louis C D. Self-similarity weighted mutual information: A new nonrigid image registration metric. Medical Image Analysis, 2014, 18(2): 343-358.

[20] Guetter C, Xu C , Sauer F ,et al.Learning basednon-rigid ultimodal image registration using Kullback-Leibler divergence. Proc.Int. Conf. Medical mage Comput. Comput.-Assist.Intervent, 2005:255-262.

[21] Wein W, Brunke S , Khamene A,et al.Automatic CT-ultrasound registration for diagnostic imaging and image-guided intervention.Med. Image Anal,2008: 577-585.

[22] Staring M, Heide V D, Klein S A U,et al.Registration of Cervical MRI Using Multi-feature Mutual Information. Medical Imaging,IEEE Transactions, 2009, 28(9): 1412-1421.

[23] Rubner Y T, Carlo G, Leonidas J. The Earth Mover's Distance As a Metric for Image Retrieval. Int. J. Comput. Vision, 2000, 40(2): 99-121.

[24] Klein S, Staring M, Pluim J P W. Evaluation of Optimization Methods for Nonrigid Medical Image Registration Using Mutual Information and B-Splines. Image Processing, IEEE Transactions,2007, 16(12): 2879-2890.

[25] Roche A , Malandain G , Pennec X ,et al.The correlation ratio as a new similarity measure for multimodal image registration. Medical Image Computing Computer Assisted Intervention (MICCAI), 1998:1115-1124.

[26] Hassan R, Sean J S C, Louis C D. Automatic Deformable MRUltrasound Registration for Image-Guided Neurosurgery. Medical Imaging, IEEE Transactions, 2014.

[27] Wolfgang W, Shelby B, Ali K,et al.Automatic CT-ultrasound registration for diagnostic imaging and image-guided intervention.Med. Image Anal, 2008, 12(5): 577-585.

[28] Bernhard F, Wolfgang W, Markus M. Automatic ultrasound–MRI registration for neurosurgery using the 2D and 3D LC2 Metric. Med,Image Anal, 2014, 18(8): 1312-1319.

[29] Roche A, Malandain G, Pennec X,et al.Multimodal Image Registration by Maximization of the Correlation Ratio. Technical Report 3378, INRIA, 1998, 3.

[30] Augustin C. Diffeomorphic surface-based registration for MRUS fusion in prostate brachytherapy. Electrotechnical Conference(MELECON) 16th IEEE Mediterranean, 2012: 903-907.

[30] Martin S, Baumann M, Daanen V,et al.MR prior based automatic segmentation of the prostate in TRUS images for MR/TRUS data fusion. Biomedical Imaging: 2010 IEEE International Symposium,2010, 4: 640-643.

[31] Mattias P H, Mark J, Bartlomiej W P,et al.Towards Realtime Multimodal Fusion for Image-Guided Interventions Using Selfsimilarities. Lecture Notes in Computer Science ofMICCAI, 2013,8149: 187-194.

[32] Roujol S, Ries M, Quesson B,et al.Real-time MR-thermometry and dosimetry for interventional guidance on abdominal organs.Magnetic Resonance in Medicine, 2010, 63: 1080-1087.

[33] Tavard F, Simon A, Hernandez A I,et al.Dynamic registration of cardiac US and CT data using Fourier descriptors and Dynamic Time Warping. Image Processing Theory, Tools and Alications(IPTA), 2012 3rd International Conference on, 2012, 198(203): 15-18.

[34] Denis D S B, Mougenot C, Moonen C T W. Real time adaptive methods for treatment of mobile organs by MRI controlled high intensity focused ultrasound. Magn Reson Med, 2007, 57: 319-330.