摘 要 介绍了一种基于最大互信息原理的图像配准技术。并针对基于最大互信息图像配准的不足，研究了基于Harris角点算子的多模态医学图像配准。在计算互信息的时候，采用部分体积插值法计算联合灰度直方图。在优化互信息函数的时候采用了改进的遗传算法将配准参数收敛到最优值附近。实验结果表明本方法具有较高的配准精度和稳定性。
1 引言
互信息是信息论的一个基本概念，是两个随机变量统计相关性的测度。Woods用测试图像的条件熵作为配准的测度，用于PET 到MR 图像的配准。Collignon 、Wells[1] 等人用互信息作为多模态医学图像的配准测度。以互信息作为两幅图像的相似性测度进行配准时，如果两幅基于共同解剖结构的图像达到最佳配准时，它们对应的图像特征互信息应为最大。最大互信息法几乎可以用在任何不同模式图像的配准中，特别是当其中一个图像的数据部分缺损时，所以这种方法广泛用于多模态图像的配准中。但是，当待匹配图像是低分辨率、图像包含的信息不够充分或两幅待匹配图像的重叠部分较少时，基于互信息的配准目标函数就会极不光滑，出现较多局部最优解，为目标函数最优解的搜索带来较大的难度。但由于该测度不需要对不同成像模式下图像灰度间的关系作任何假设，也不需要对图像进行分割或任何预处理，因此，该测度可以被广泛地应用于CT-MR，PET-MR 等多种图像的配准工作。
基于最大互信息的图像配准取得了很大的成功，但是它也存在一些缺陷，例如任何图像和一幅只有一种灰度值的图像（例如灰度值为255 的全黑图像）配准，无论几何变换怎样，他们的联合灰度直方图都是一样的。因此在实际情况中我们就会碰到一个问题，一幅图像可能包括大量的同类区域（例如天空、大海等），那么这样的图像就不太适合用最大互信息的方法进行配准，实际上基于图像灰度的配准方法都不是太适合这样的情况。此外，基于最大互信息的图像配准因为要进行全局参数优化搜索，配准时间也比较长。Rangarajan 等提出了一种利用互信息匹配形状特征点进行配准的策略. 该策略针对待配准的两幅图像，首先分别提取出形状特征点的集合，并定义这两个集合它们的互信息，然后使之最大化，以达到配准。
由于角点是景物轮廓线上曲率的局部极大点，对掌握景物的轮廓特征具有决定作用。一旦找到了景物的轮廓特征点也就大致掌握了景物的形状。虽然角点相对于其他的特征来说比较少，但是使用大量特征实现配准必会导致算法复杂度的提高。出于提高速度而又不降低配准精度的考虑，角点是一个很好的配准特征。因此，在本文中我们将要对基于角点特征的图像配准做一些初步的探讨。
2 配准方法
2.1 变换和插值模型
我们将研究的范围限制在二维脑断层图像的配准. 因为脑组织受到颅骨的严密保护，所以脑部运动可以近似为刚体运动，即内部无相对运动. 同时，假设待配准的图像经过预处理后具有相同的空间比例.我们的目标是寻求空间变换T，使MI(T) 最大.。针对前面所做假设，令T = T1*T2；其中，T1 为平移矩阵，T2为旋转矩阵。
最近邻插值法的精确度很低，而双线性插值法会产生新的灰度值。这对于联合直方图的统计是不利的。因为新加入的灰度值使得随着Tα的一些小变动，联合直方图中就会增加新的象素对，或者减少象素对，从而互信息值变化比较大，也就是互信息函数曲线会不光滑，这样不利于优化。因此为了消除新产生的灰度值的不利影响，我们在配准过程中引入了另外一种插值法：PV（Partial Volume）插值法。从产生插值图像这个方面来说，PV 插值法不能算是一种插值方法，它是专门针对联合直方图的更新而设计的。和双线性插值法一样，PV 插值法也是利用点Tα(X)的四个最近邻点和权值 ，可是，不同于双线性插值法的是，PV 插值法不是根据最近邻点的加权平均所得到的灰度值，从而更新联合直方图，而是根据权值 使周围四个象素点都贡献于联合直方图的统计，如图1所示。
可用公式表示为：
(2-1)
(2-2)
图1 pv插值法
2.2 特征点的提取
由于角点是景物轮廓线上曲率的局部极大点，对掌握景物的轮廓特征具有决定作用。一旦找到了景物的轮廓特征点也就大致掌握了景物的形状。直观的讲，角点就是图像上所显示的物体边缘拐角所在的位置点。
Harris角点检测法是一种基于图像灰度的检测方法，这类方法主要通过计算点的曲率及梯度来检测角点。该方法是由Harris和Stephen提出来的，也叫Plessey角点检测法。其基本思想与Moravec角点算子相似，但对其作了许多改进。
Moravec角点算子计算各象素沿小同方向的平均灰度变化，选取最小值作为对应象素点的角点响应函数。定义在一定范围内具有最大角点响应的象素点为角点。假设图像的灰度定义为I那么平移(x，y)所得到的灰度变化的计算公式为：
(2-3)
这里W表示图像窗口，平移(x，y)表示了四个方向：水平、垂直、对角线和反对角线，即(0，1)，(1，0)，(1，1)，(-1，1)。
Moravec角点算子简单快速，但是它存在一些缺点，Harris角点算子正是针对这些缺点做了很大的改进。
首先，Moravec角点算子是各向异性的，因为它的角点响应只计算了四个方向。故为了包含所有的方向，Harris角点算子对式(2-3)进行了展开：
(2-4)
这里一阶微分可以由下面的式子近似：
(2-5)
因此，E可以表示
(2-6)
这里
(2-7)
为了避免噪声的影响，这里w采用高斯平滑窗口：
(2-8)
其次，Moravec角点算子对强边界敏感，这是因为它的响应值只考虑了E的最小值。Harris角点算子则利用了E在平移方向上的变化。
在平移方向(x，y)上的E可以表示如下
（2-9）
这里2×2的矩阵M为
(2-10)
可以看出，E和局部自相关函数联系非常紧密。设α，β为矩阵M的特征值，则α，β与局部自相关函数的主曲率成比例。当两个曲率都低时，局部自相关函数是平坦的，那么窗口图像区域的灰度值近似为常量；当只有一个曲率高而另一个曲率低时，局部自相关函数呈脊形，那么E只有当沿山脊移动时变化小，这就表示是边缘；当两个曲率都高时，局部自相关函数是尖峰，那么E在任意方向上移动都会增加，这就表示是角点。因此我们可以由α，β的值判断是否是角点。为了不对M进行分解求特征值，可以采用Tr(M)和Det(M)来代替α，β，其中
(2-11)
从而形成对矩阵M与旋转无关的描述：
(2-12)
其中：且，K是随高斯函数和微分模板变化的变常量，一般推荐取为0.04。
只有当图像中象素的R值大于一定的门限，且在周围的八个方向上是局部极大值时才认为该点是角点。
2.3 多元互信息
互信息可用熵来描述，熵表达的是一个系统的复杂性或者是不确定性。一幅图像的熵反映了该图像中像素灰度的分布情况，灰度级别越多，灰度越分散，熵就越大。反映在直方图上就是灰度动态范围应用充分，且平坦。
图像的熵是对图像概率分布的一种表述。熵的定义为：
其中，H(A)，H(B)分别为A和B的信息熵，H(A|B)和H(B|A)分别为给定B的条件下A 的条件熵以及给定A的条件下B的条件熵，H(A，B)为A和B的联合熵。a∈A，b∈B， pA(a)、pB(b) 分别表示图像A和B的概率分布，pAB(a，b)表示2幅图像的联合分布。互信息可以用信息熵来表示，其关系式为：
I(A，B)=H(A)-H(A|B)=H(B)-H(B|A)=H(A)+H(B)-H(A，B)
2.4 优化算法
这里采用遗传算法搜索变换参数，它的并行性可以避免搜索陷入局部极值。由对标准遗传算法的分析可以知道，影响遗传算法应用的主要因素有两个：一是容易早熟；二是收敛速度慢。这里针对医学图像配准对标准遗传算法进行了改进，改进后的遗传算法
1）编码方式
在针对2D人脑图像配准的算法中，有三个待寻优参数：旋转角a、一幅图像相对于另一幅图像沿X轴方向的平移rX、一幅图像相对于另一幅图像沿Y轴方向的平移rY。对这三个参数采用实数编码方式，减少编码解码所耗费的时间，改善算法的搜索效率，提高配准的速度。
2）适应度表示
这里用两幅图像的角点互信息作为图像配准的相似性测度，这使得适应度函数的描述非常简单且易于实现。在某一变换T下，个体的适应度为：
3）轮盘赌法和最优保存策略
在用轮盘赌法进行个体选择时有可能产生随机误差，导致当前种群中适应度最高的个体没被选中，使其在下一代中得不到繁衍。为了避免这种现象，这里采样最优保存策略，把每一代种群中适应度最高的个体直接复制到下一代，对剩下的N-1个个体采用轮盘赌法进行选择。直接保存最优个体可以保证最优个体在下一代中出现，改善局部搜索能力，提高收敛速度。
适应度高的个体变异概率小，在小范围内搜索；相反，适应度低的个体变异概率大，在较大范围内搜索。
3 实验结果
在这一部分我们使用本文提出的配准方法，对CT和PET不同设备采集的图像做刚性配准．以CT 图像为配准参考图像(refernce image)，以PET图像为可变动配准图像(floating image)。结合Harris角点和互信息技术，使用改进的遗传算法做配准实验。
图2、图3分别是提取了Harris角点后的CT图和PET图。我们应用角点特征的互信息结合遗传算法，计算最佳配准参数（旋转角度α，和平移参数x，y)。
图2 提取Harris 图3 提取Harris
角点后的CT图 角点后的PET图
本程序应用遗传算法做变换参数的搜索，根据提取的角点处的局部子块计算互信息（如提取了角点的CT图的角点局部子块与原未提取角点的PET图计算互信息），找出变换后两幅图像互信息最大的变换参数。图4为配准过程及结果的效果图。
图4
4 结束语
本文提出了基于角点特征的多模态医学图像配准方法。本文采用了Harris角点提取算子，应用互信息作为相似性度量。应用改进的遗传算法做最佳变换参数的搜索。实验结果表明该方法配准速度较快，精度好，是一种有效的自动配准方法。下一步工作中我们要对提取的角点特征做相应处理，去掉部分特征不明显的角点信息，进一步提高配准的速度。
参考文献
[1] 罗欣等.多元互信息在超光谱图像自动配准中的应用. 计算机工程与应用，2006.27(3)：3-8
强赞霞，彭嘉雄，王洪群．基于互信息的分层遥感图像配准方法【J】．计算机工程与应用，2004.40(13)：31—33
白继伟，赵永超，张兵等．基于包络线消除的超光谱图像分类方法研究【J】.计算机工程与应用，2003；39(13)：88—90
Neville R A，Sun L X，Staenz K．Detection of Keystone in Image Spectrometer Data[C]．In：Shen S S，Lewis P E eds．Algorithms and Technologies for Multispectral，Hypempectral and Uhraspectral Imagery X，Bellingham：SPIE，2004：208—217
周永新，罗述谦.【基于形状特征点最大互信息的医学图像配准】.计算机辅助设计与图形学学报2002.14(7)654- 658
张 煜，刘哲星，李树祥，陈武凡.基于互信息量的医学图像自动弹性配准. 中国生物医学工程学报，2004.23(3)：279-281
张二虎，卞正中. 基于最大熵和互信息最大化的特征点配准算法. 计算机研究与发展，2004.41(7)：1194-1198
C. Harris and M. Stephens. A combined corner and edge detector. Proc. 4th Alvey Vision Conf.，1998：147~151.
