1、1计算机学报2009 年 8 期基于等效子午面与互信息量的医学图像配准卢振泰，冯衍秋，冯前进，陈武凡（南方医科大学生物医学工程学院，广州 510515）摘要： 本文提出了一种新的基于等效子午面和互信息量的三维医学图像快速配准算法EMP-MI 算法。传统的互信息量的方法需要考虑整个三维数据的信息，计算复杂度大，无法满足临床需要。而本算法将三维数据的配准转化为二维数据的配准，在保证精度的前提下，减少了配准所需时间。本文创新点在于利用主成分分析计算出图像的等效子午面并将图像转化到标准坐标系下，从而将质心和等效子午面粗配准，精细配准时只需要对浮动图像进行微小的调整计算等效子午面的互信息量，这就大大提高
2、了配准速度，减少了陷入局部极值的可能。实验结果表明这种先整体后局部的方法能准确，快速地处理图像刚性配准问题，特别适用于三维医学图像的配准。关键词：等效子午面 互信息量 图像配准 主成分分析 局部极值中图分类号 TP391 文献标识码 A 文章编号引言医学图像配准 1-2是指对于一幅医学图像寻求一种（或一系列）空间变换，使它与另一幅图像上的对应点达到空间上的一致。这种一致是指人体上的同一解剖点在两张图像上有相同的空间位置，配准的结果应使两幅图像上所有的解剖点匹配。图像配准的主要目的是去除或者抑制待配准图像和参考图像之间几何上的不一致，包括平移、旋转和形变。它是图像分析和处理的关键步骤，是图像对比
3、、数据融合、变化分析和目标识别的必要前提。多模态医学图像配准是目前医学图像处理领域中的一个热点问题，对于临床诊断和治疗有重要意义。医学图像可以提供有病变组织或器官的大小、形状、空间关系等详细信息，比如CT图像可以显示骨骼结构和组织密度分布情况；MR图像和超声(US) 图像提供的则是软组织的信息；PET，SPECT能反映人体的功能和代谢信息。单一模态的图像往往不能提供医生所需要的足够的信息, 通常要将不同模态的图像融合在一起, 得到更丰富的信息以便了解病变组织或器官的综合情况, 从而做出准确的诊断或制订合适的治疗方案。目前已有多种图像配准方法，大体可以分为基于图像特征和基于图像灰度的配准方法。基
4、于图像特征的方法一般通过寻找图像中比较明显的解剖结构特征来计算变换参数，提取的特征包括点、线、边缘2和轮廓等。该方法计算效率比较高，但是配准的精度取决于特征提取的准确性与否。大部分情况下，医学图像特征点的位置比较复杂，很难进行准确地提取。基于图像灰度的配准方法不需要提取特征点，而是利用图像的灰度信息进行配准。由于充分利用了全部的灰度信息，因此一般更容易得到较精确的配准结果。其中互相关法、最小平方差法 1-2和互信息量法 3是几种常用的方法。互相关法一般实现起来简单，但是计算代价庞大；最小平方差法对图像数据的部分缺失和图像的背景噪声不是非常敏感，具有较强的鲁棒性，但由于需要通过迭代搜寻最优参数，
5、故其执行效率也不是很高，而且其解严重依赖于初始值的选择，不恰当的初始值往往会使解陷入局部极值。互信息量技术是最近几年提出的进行多模态图像配准一种有效方法。它是一种自动的，基于像素灰度的方法，它不需要选择标志点或提取图像特征，不需要假设图像中的灰度值存在某种线性关系，因而在世界范围内得到广泛应用，特别是医学图像处理领域。本文考虑到以上问题，提出了一种新的特征和灰度相结合的算法基于等效子午面与互信息量的快速三维医学图像配准算法。该算法综合了基于特征和基于灰度配准算法快和准的优点，将图像的轮廓看作数据点集合，利用主成分分析找出图像的等效子午面并将其转化到标准坐标系下，使等效子午面与平面重合，从而完成
6、了粗配准；精细配准时只需要对浮动图像进行较小的调整计算 平面的互信息XY XY量，这就大大提高了配准速度，减少了陷入局部极值的可能。1. 等效子午面对于三维医学图像而言，由于是无标记配准，所以无法利用标记点实现快速配准。而如果利用整个三维数据的信息，计算复杂度太大，无法满足临床要求。能否找到三维空间中的一个“特殊”平面，将三维数据的配准简化为二维数据的配准，从而在保证精度的前提下，提高配准的速度？我们知道子午面是地理学中的一个重要概念，它是过地球南北两极所确定的轴线(地轴) 的平面，地球的运动可以通过子午面来表示。推而广之，对于任意曲线，绕曲线端点所确定的轴线旋转 360 度都可以形成一个三维
7、体，而过轴线的平面即为该三维物体的子午面，三维物体的运动可以简化为子午面的运动。对于医学图像，由成像设备产生的断层序列图像也可以重建为三维体，但该三维物体一般是非常不规则的，故不能求取一般意义上的子午面。为此本文提出了等效子午面的概念，从而可以利用等效子午面来表示三维物体的运动。定义：对于断层序列图像重建出的三维体，分别求出其第一主轴和与之正交的第二主轴，而由第一主轴和第二主轴所确定的平面称之为等效子午面。对于二维平面上的点列，可以利用主成分分析 4-6的方法找出点列的第一主轴和第二主轴。第一主轴为穿过点集的最好直线，即点集中所有的点到该直线的距离的平方和最小。与之正交的直线即为第二主轴，如图
8、 1(a)所示，红色为第一主轴，绿色为第二主轴。在三维空间中，第一主轴和第二主轴，确定了一个平面，它满足集合中所有点到这个平面的垂直距离的平方和最小，即等效子午面，如图 1(b)所示。由3于三维医学图像一般是非常不规则的，所以等效子午面是唯一确定的。本文利用图像的“特征点”信息，确定一个“虚拟”的平面，这个平面可以表征整个三维刚性物体的运动，从而可以快速、有效地进行物体的粗配准。该方法实质上是利用主成分分析进行降维，因而不可避免地会造成信息的丢失，但丢失的信息对确定等效子午面没有影响，所以对配准结果也没有影响。(a) 2-D (b) 3-D图 1 数据点集及主轴（分别用红色，绿色和蓝色表示）F
9、ig.1 The red, green, and blue lines represent the first, the second, and the third principal axes, respectively.三维图像的配准关键是寻找图像的等效子午面。首先将图像二值化，找出待配准物体所在的区域，则区域内点集可以表示为三维向量 ，这里 是 轴， 轴和 轴(,)|1,TiiXxyzn ,iixyzxyz方向的坐标值。点集的均值 即为图像的中心：u(1)1niiuX协方差矩阵 表示为：C(2)1nTiiCu对协方差矩阵 进行特征值分解，得到特征向量矩阵 和特征值矩阵 。 的第一行为最大
10、特征值对应EVE的特征向量，最后一行为最小特征值对应的特征向量。则有：(3)()YXu该式称为主成分变换（PCA transform） 。它的作用是以点集的中心为原点，以协方差矩阵 的特征向量所C指方向为 轴， 轴和 轴的方向建立一个新的坐标系统，将物体沿着自身的主轴对准，从而等效子午xyz面与 平面重合。该变换实质上是一种旋转变换， 即为旋转矩阵。如果主轴与 轴， 轴和 轴的XYExyz夹角分别为 , 和 ，则有：xyy(4)()()xyzRR从而可以得到：(5)3121arcsin()/osiyx yze4图 2 是用头部三维模型说明主成分变换的过程。图 2(a)是一头部三维图像，(b)为
11、旋转后图像，分别用黄色和红色表示。(c)为(a)和(b)叠加后图像。(d) ，(e) 分别为(a)和(b)进行主成分变换后的结果，(f)为(d)和(e)叠加后图像。从图中可以看出主成分变换是先将图像的中心移到直角坐标系的原点，然后将主轴旋转至轴， 轴和 轴。xyz(a)模型 1 (b)模型 2 (c)变换前叠加情况(d)变换后的模型 1 (e)变换后的模型 2 (f)变换后叠加情况图 2 头部三维模型及其变换后的结果Fig.2 The PCA Transform of head2. 基于等效子午面的互信息量配准2.1 互信息量互信息是信息论中的一个测度, 用来度量两个随机变量之间的相似性 3。
12、在多模态医学图像配准问题中，虽然两幅图像来源于不同的成像设备，但是它们基于人体共同的解剖信息，所以当两幅图像的空间位置完全一致时，它们的对应象素的灰度互信息达到最大值，即：一幅图像表达的关于另外一幅图像的信息最多。这个信息就是互信息量。对于给定图像A和图像B，它们的熵和联合熵定义为： (6)()()log()AaHPa(7)Bbb(8),(,)(,)l(,)AABa其中 ， ， 分别为A，B的概率分布和联合概率分布，可以通过归一化各自的灰度直方图和联合灰AP度直方图而得到。互信息量定义为：(9)(,)()(,)MIABHAB5作为一种相似性测度，互信息量取得了巨大的成功，特别在医学图像领域。互
13、信息配准方法目前被公认为是配准精度和鲁棒性最好的回溯性配准方法之一 7-8。2.2 二值化二值化在本算法中起到很大的作用，因为通过二值化得到的区域的形状直接影响到后面求取变换参数的速度与精度。本文结合 Otsu 算法 9和最大互信息量算法，提出了一种新的改进算法。首先利用 Otsu算法求出全局阈值作为初值，然后在初值附近搜索，计算原图像与二值化图像的互信息量，当互信息量最大时的阈值即为最优值。这种方法不仅考虑到图像的灰度分布，而且考虑到像素的空间关联信息以及分割图像与原图像之间的内在联系。所以对于照明不均匀、有突发噪声，或者背景灰度变化比较大的图像都能够得到较好的结果。2.3 基于等效子午面的
14、配准由于同一病人，不同模态的图像都有相似的外轮廓，因而等效子午面也基本相同。等效子午面是三维数据中一个非常特别的平面，集合中每一点到该平面的距离的平方和最小，该平面也是唯一确定的。从而两个三维图像的配准可以简化为等效子午面的对准，即等效子午面互信息量最大时就说明图像完全对准。本文提出的先整体后局部的方法由于先在整体上对图像进行了粗配准，使得图像已经基本对齐，从而避免了陷入局部极值的可能性。精确配准时等效子午面上的信息量已经完全可以满足统计分析的需要，与计算整个三维数据的互信息量相比，这就大大减少了计算量。一般优化算法的缺点是当两幅图像有较大的平移和旋转时，计算时间比较长且容易陷入局部极值，出现
15、误配。为了解决这一问题，本算法利用主成分分析精确确定等效子午面的位置，并将其移动到三维图像的中间层，即 平面，如图 3 所示。从而方便计算互信息量。XY图 3 经 PCA 变换后等效子午面的位置Fig.3 The EMP after PCA transforms. 由于变换后的图像的中心和等效子午面已经基本对齐，这时只需要将一幅图像固定作为参考图像，另外一幅作为浮动图像在较小的范围内调整图像的平移量和旋转角度，计算调整后图像的中间层与参考图像中间层的互信息量。通过优化算法找出达到最大互信息量时的变换参数。这样就减少了陷入局部极6二值化粗配准： PCA 变换精配准：最大化等效子午面上的互信息量值
16、的可能性。2.4 新算法描述本文算法主要步骤描述如下, 算法流程见图 4 1) 给定待配准图像 A，B ,对其进行二值化，得到向量表示；2) 分别计算 A，B 的质心和协方差矩阵；3) 计算 PCA 变换；4) 利用 PCA 变换将原图像转化到标准坐标系下，完成粗配准；5) 小范围内调整浮动图像的位置，计算中间层（等效子午面）的互信息量。通过优化算法求出变换参数，互信息量最大时图像完全配准。图 4 算法流程图Fig.4 Flowchart of the registration process3 实验结果我们通过四组的实验，将新算法（记作 EMP-MI）与基于互信息量的方法（记作 MI）在配准
17、精度和运行时间上进行比较。优化算法均采用 POWELL 算法 10，初始参数设置为零，搜索顺序为：。利用直方图估计概率密度分布时采用 64 位灰度级，则可在速度和精度上均达到较好(,)xyzxyztt效果，故本文所有实验均采用 64 位灰度级。用 Matlab7.0 在 PC 机上(P4 3.0 GHz CPU, 512 M 内存, WindowsXP SP2 操作系统)进行了实验，分别对不同模态的三维图像进行配准。对于已经配准的两幅图像一幅作为目标图像，另一幅作为源图像。将源图像先绕 x, y 和 z 轴分别进行旋转( 单位：度)，再沿 x, y 和7z 轴分别进行平移(单位：像素 )。旋转
18、角度和平移量是随机产生的。如果计算出的变换参数与真实值相差1 个像素或 1 度，即认为配准成功，每套数据进行 100 次实验，配准成功的次数与实验总数的比值即为成功率。3.1 模拟脑部 MR 图像的配准MR-T1 加权和 MR-PD 加权 MR 图像，大小均为 181217181，像素大小为 1.0 mm3。图像取自BrainWeb (www.bic.mni.mcgill.ca/brainweb)，如图 5 所示：(a) (b) (c) (d)图 5 实验数据 (a) 源图像： T1 横断面 (b) 目标图像：PD 横断面 (c) 源图像冠状面(d) 目标图像冠状面。Fig.5. Images
19、 for registration experiments. (a) Source volume T1 MRI brain, axial view. (b) Target volume PD MRI brain, axial view. (c) Source, coronal view. (d) Target, coronal view.MR-PD 图像作为目标图像，MR-T1 图像作为源图像，本实验中旋转角度和平移量满足-10，10 均匀分布。表 1 为实验结果，分别为运行时间(单位：分钟)和成功率。表 1 配准结果Tab.1 Registration results算法 时间(分钟) 成功
20、率MIEMP-MI10.143.69100%100%从实验结果中可以看出对模拟数据进行较小的变形时，两种算法均可以达到 100%成功率，但 EMP-MI 算法所用时间较少，执行效率为 MI 算法的 2.7 倍。32 同一病人脑部 MR T1 加权与 T2 加权图像的配准。图像取自 BIS 医学图像库，该图像库提供各种疾病的影像资料(http:/nova.nlm.nih.gov/Mayo )，供医生学习参考。本实验中所用图像为一脑部患有肿瘤的病人的 T1 加权和 T2 加权 MR 图像，已经由专家配准好。其中 T1 图像大小为 25625662，像素大小为 0.9370.9373.0 mm，T2
21、 图像大小为 25625662，像素大小为 0.9370.9371.5 mm。将 T1 图像进行平移和旋转，平移量和旋转角度满足-20，20 均匀分布。实验结果见表 2：表 2 配准结果Tab.2 Registration results算法 时间(分钟) 成功率8MIEMP-MI25.167.0494%96%EMP-MI 方法成功率为 96%，而 MI 方法成功率为 94%，运行时间为 EMP-MI 方法的 3.6 倍。图 6 是利用 EMP-MI 算法配准后的结果，由上至下分别为 T1，T2 图像和将 T2 边缘（红色）加到 T1后的图像。从左到右依次为第 24，28，32 层横断面。图
22、6：EMP-MI 算法配准结果Fig.6 the registration result of EMP-MI33 同一病人脑部 MR 图像与 PET 图像的配准。所用图像取自 http:/www.cs.tut.fi/kurssit。已经由专家配准好，见图 7，由上至下分别为 MR 和 PET图像，从左到右依次为横断面、矢状面和冠状面。MR 图像大小为 12812863，像素大小为 2.12.12.4 mm，PET 图像大小为 12812867，像素大小为 2.12.12.4 mm。图 7：实验所用图像：MR(第一行)和 PET(第二行) 图像，从左到右依次为横断面、矢状面和冠状面。9Fig.7
23、. Axial, sagital, and coronal view of the MR (first row) and PET (second row) brain volumes used in registration由于 PET 图像成像质量差，空间分辨率比较低，与前两组实验相比有一定困难。为了避免出现误配的情况，我们计算等效子午面及其它两个与之垂直的平面的互信息量，即计算 平面， 平面和XYZ平面的互信息量。对 PET 图像进行平移和旋转。平移量满足-10，10 均匀分布，旋转角度满足-YZ20，20 均匀分布。实验结果见表 3表 3 配准结果Tab.3 Registration r
24、esults算法 时间(分钟) 成功率MIEMP-MI11.684.6190%92%由配准结果可以看出对于分辨率较低的 PET 图像 EMP-MI 算法仍然能够达到较好的效果，成功率为92%，高于 MI 算法的成功率。而所用时间约为 MI 算法的 1/3。34 同一病人脑部 MR 图像与 SPECT 图像的配准图像取自 BIS 医学图像库 (http:/nova.nlm.nih.gov/Mayo)。见图 8，由上至下分别为 MR 和 SPECT 图像，从左到右依次为横断面、矢状面和冠状面。MR 图像和 SPECT 图像已经由专家配准好。图像大小均为 256256124，像素大小为 0.9370
25、.9371.5 mm。图 8：实验所用图像： MR(第一行)和 SPECT(第二行) 图像，从左到右依次为横断面、矢状面和冠状面。Fig.8. Axial, sagital, and coronal view of the MR (first row) and SPECT (second row) brain volumes used in registration为了避免出现误配的情况，该组实验仍然计算等效子午面及其它两个与之垂直的平面的互信息量。对 SPECT 图像进行平移和旋转。平移量满足-20 ，20均匀分布，旋转角度满足 -10，10 均匀分布。实验结果见表 4：表 4 配准结果Ta
26、b.4 Registration results算法 时间(分钟) 成功率10MIEMP-MI36.608.4375%92%通过比较可以看出 EMP-MI 算法在精度和速度上具有较为明显的优势。利用 MI 算法配准时成功率比较低，仅为 75%，而利用 EMP-MI 算法配准成功率为 92%，所用时间是 MI 算法的 23%。4结论本文首次提出了等效子午面的概念并将其应用到三维医学图像的配准中。首先利用主成分分析的方法确定等效子午面完成粗配准，然后计算等效子午面上的互信息量进行精确配准。从实验结果可以看出，这种先整体后局部的方法能够准确、快速地处理刚性配准问题，特别适用于三维医学图像的配准。但当
27、图像数据存在缺失或有较大变形时，得到的等效子午面往往不能代表原三维物体的运动，因而容易导致误配。对于这一类情况我们还需进一步研究，以改善该方法的通用性。参考文献：1. J B Antoine Maintz , M A Viergever. A survey of medical image registration. Medical Image Analysis，1998, 2 (1) : 1-362. B. Zitova, J. Flusser. Image Registration Methods: A Survey. Image and Vision Computing, 2003, 2
28、1(11): 977-10003. F Maes, A Collignon, Dirk Vandermeulen, et al. Multimodality Image Registration by Maximization of Mutual Information. IEEE Trans Med Imaging, 1997, 16 (2): 189-1984. Alpert NM, Bradshaw J F. Principal axes transformation：a method for image registration. J. Nucl. Med, 1990, 31(10):
29、 1717-17225. H Bulow, L Dooley, D Wermser. Application of principal axes for registration of NMR image. Pattern Recognition Letters, 2000, 21(4): 329-3366. Rafael C. Gonzales, Richard E. Woods, Digital Image Processing. Prentice Hall, 20027. C. Studholme, D. L. Hill, and D. J. Hawkes. Automated thre
30、e-dimensional registration of magnetic resonance and positron emission tomography brain images by multiresolution optimization of voxel similarity measures. Med. Phys, 1997, 24(1): 25-358. A. Carrillo, J. L. Duerk, J. S. Lewin, D. L. Wilson. Semiautomatic 3-D image registration as applied to interve
31、ntional MRI liver cancer treatment. IEEE Trans Med. Imag, 2000, 19(3): 1751859. N. Otsu. A threshold selection method from gray-level histograms. IEEE Trans SMC, 1979, 9(1): 6266.10. Yang-Ming Zhu, Steven M. Cochoff. Likelihood Maximization Approach to Image Registration. IEEE Trans Image Process, 2002, 11(12): 14171426
