基于交叉累计剩余熵的图像配准中插值方法的改进

相 艳，贺建峰，张云春，蔡 莉

（1.云南省计算机技术应用重点实验室，昆明 650500；2.昆明理工大学信息工程与自动化学院，昆明 650500；3.云南大学软件学院，昆明 650091）

基于交叉累计剩余熵的图像配准中插值方法的改进

相 艳1，2，贺建峰2，张云春3，蔡 莉3

（1.云南省计算机技术应用重点实验室，昆明 650500；2.昆明理工大学信息工程与自动化学院，昆明 650500；3.云南大学软件学院，昆明 650091）

交叉累计剩余熵（CCRE）比传统互信息在配准强噪声图像时更具优势，但采用部分体积（PV）插值的CCRE在网格点容易产生局部极值，不利于变换参数的优化。针对该问题，研究基于3阶B样条函数的PV插值（BPV）、哈宁窗sinc函数的PV插值（HPV）和Blackman-Harris窗sinc函数的PV插值（BHPV）方法在CCRE中的应用，提出一种新的插值方法。该方法采用灵活的邻域中心，将插值点对联合直方图贡献的权重分散到临近的9个点上，并使用高斯函数作为PV插值的核函数，避免权重突变。实验结果表明，与BPV，HPV和BHPV插值方法相比，该方法对噪声图像的配准率较高，配准速度较快，更适合应用于CCRE的计算。

图像配准；交叉累积剩余熵；部分体积插值；高斯函数；3阶B样条函数；sinc函数

DO I：10.3969/j.issn.1000-3428.2015.10.037

1 概述

基于互信息（Mutual Information，M I）的图像配准方法具有自动化程度高、配准精度高等优点。但当图像采集过程中受噪声影响严重，或2幅待匹配图像的重叠部分较少或分辨率较低时，互信息曲线会出现局部极值，变得极不光滑。 累积剩余熵（Cumulative Residual Entropy，CRE）是由文献［1］提出的一种新的信息度量方式，该方法将香农熵定义中的概率密度函数用概率分布函数来代替。Wang等学者将累计剩余熵用在图像配准中，提出将2幅图像的交叉累计剩余熵（Cross Cumulative Residual Entropy，CCRE）作为配准的相似性测度［2］。CCRE是统计大于某一灰度值的所有直方图信息，保持了局部区域灰度信息的连续性，从而克服噪声对局部极值的影响，能得到较 M I更光滑的配准目标函数［3-4］。

尽管如此，CCRE仍然会出现局部极值，其原因一方面与待配准的2幅图像存在一定的局部匹配有关，另一方面与浮动图像经空间变换后的插值方法有关。CCRE与M I一样需要计算2幅图像的联合直方图，因此，MI的插值方法可供借鉴。传统的插值方法有最近邻法、双线性法等，但两者精度都不高。目前被广泛采用的是部分体积插值法（Partial Volume，PV）。该方法无需进行浮动图像的灰度插值，而是直接更新联合直方图。然而，PV插值的一个明显缺陷是在网格点产生联合熵的突变，从而引起局部极值［5-6］。 针对这一问题，文献［6］利用均匀B样条基函数拓展了PV插值法的邻域，文献［7-8］提出一种新的基于 Hanning窗的PV插值方法，文献［9］用Blackman-Harris窗sinc函数取代传统 PV插值法的核函数。上述方法都是将某点对联合直方图的影响从4个邻域点扩展到16个，有效抑制了MI函数曲线中的局部极值现象。

本文将PV插值方法以及上述3种改进的PV插值方法应用到CCRE的计算中，验证3种改进方法同样能抑制PV插值的CCRE曲线中网格点的局部极值，但运行时间较长。同时提出一种新的插值方法，使用高斯函数作为PV插值的核函数，将插值点对联合直方图的贡献分散到其临近的9个点上，以缓解网格点权重的突变，得到平滑的CCRE曲线，并且提高运行速度。

2 交叉累计剩余熵的计算原理

CRE是一种新的对随机变量的信息测度。设χ是R中的一个随机变量，并且F（λ）是累计剩余分布函数。把香农熵定义中的概率密度函数用其自身的累计剩余分布函数代替，则χ的CRE定义为：

其中，χ∈R+。

累积剩余熵比香农熵更具有一般性，其定义在连续和离散区域都是合理的［10］。图像 X和图像 Y的CCRE则定义为：

具体用IR（χ）表示参考图像，IT（χ）表示浮动图像，变换矩阵为 g（χ；μ），μ为一系列的变换参数。配准问题即转化为寻找最佳变换参数 μ，使得经过变换矩阵变换后的浮动图像 IT°g（χ；μ）与参考图像IR（χ）的CCRE达到最大值，即：

定义P（l，k；μ）为（IT°g（χ；μ），IR）的联合概率密度，PT（l；μ）和PR（k）分别代表浮动图像和参考图像的边缘概率密度，LT和LR分别代表浮动和参考图像的离散密度点集，则2幅图像的CCRE可具体表示为：

3 不同插值方法及其在CCRE计算中的应用

3.1 PV插值法及其改进

计算式（4）中的联合概率密度时，浮动图像的像素经过空间变换映射到参考图像中时通常为非整数点，需要进行插值处理。常用的PV插值法是根据权重分配原则，将每对像素对联合直方图的贡献按权重分配给相邻的4个像素对。然而，当插值点正好落在网格点处时，则该网格点处权重为1，而其他相邻的3个网格点权重为0。因此，非网格点和网格点的PV插值所更新的像素对比例为4：1，联合直方图的杂散程度的差异明显，导致CCRE容易在整数倍位移处出现局部极值。文献［6-9］的3种插值方法都是将插值点对联合直方图的贡献从4个邻域点扩展到16个，采用不同的PV插值核函数，从而使得插

值点的权重分布更合理。PV插值改进方法示意图如图1所示。

图1 PV插值改进方法示意图

如图1所示，设浮动图像 IT（χ）上某点 P经空间变换后与参考图像 IR（χ）上的点 q对应，点 q的邻域点为qi（i=1，2，…，16），其联合直方图的计算如下：

其中，ωi为每个领域点qi的权重；

在式（6）中，diχ和 diy为点qi与q在 χ和 y方向上的距离，而PV插值核函数f有以下3种形式：

（1）f为3阶 B样条函数［6，11］：

（2）f为Hanning窗sinc函数的近似［7-8］：

（3）f为Blackman-Harris窗sinc函数［9］：

这3种PV插值的改进方法的共同点是：当点q落在非网格点上时，参与联合直方图统计的邻域点数目有16个；而当点q落在网格点q6上时，参与联合直方图统计的邻域点数目是最少的，有9个。能分配到权重的最多邻域点数与最少邻域点数的比例为16：9，比传统 PV插值法的比例4：1要小。这个比例越小，说明联合直方图的概率分布越分散［9］，CCRE出现局部极值的机率降低。

3.2 本文插值方法

上述插值方法存在一定的不足：（1）能分配到权重的最多邻域点数与最少邻域点数的比例没有达到1:1。（2）插值点权重需要分散到邻域16个点，应用到CCRE配准中速度较慢。（3）利用插值点与邻域点的空间距离来决定权重大小，但是固定的邻域范围影响了权重分配的公平性。这类插值方法都以 q的整数点 q6为中心形成16个点的邻域。当q较靠近q6时，这种邻域范围是比较合理的。但是当 q远离q6时，则这一邻域范围对于权重分配不太合理。例如当q与q6在 χ和 y方向的距离都为0.9时，q与q1的欧氏距离为2.687 0，与q2和q5的欧氏距离为2.102 4，与q4和q13的欧式距离为2.195 4。而q与邻域范围外的q19和 q23的欧式距离为2.102 4，与q18和q22的欧氏距离为2.284 7，与q20和q24欧式距离为2.230 7。若按空间距离的远近来决定权重分配的大小，则这些邻域范围外的点也理应被分配到权重，但在这3种插值方法中因为不在16个邻域点之中而被忽略。

根据以上分析，本文提出一种基于9个邻域点的插值方法。

如图2所示，其中，虚线网格表示整数点之间的平分线。设浮动图像上某点 P经空间变换后与参考图像IR（χ）上的点 q对应，q的整数点为 q5。选取 q5，q6，q8，q9中与点 q距离最近的点作为邻域中心，将权重分配给包括该中心在内的 9个邻域点。例如图 2中 q离 q5较近，则权重分配给 qi（i=1，2，…，9）；而另一插值点 q′的整数点同样为q5，但其离 q9较近，则其权重分配给以 q9为中心的9个邻域点，即

图2 本文插值方法示意图

确定了9个邻域点后，各点的联合直方图也采

用诸如式（5）和式（6）进行计算，只是式（5）中的16个点变为9个点。同时，本文提出一种新的核函数f来计算式（6）。

根据权重分配的原理，核函数 f应当满足［6］：（1）qi与 q距离越小，式（6）中的 ωi越大；（2）ωi非负经过合理设计，形如式（10）的高斯函数也能满足以上要求：

式（10）中标准差σ决定了权重ωi的分布，归一化系数α的作用是将9个领域点的权重和调整到约等于1。经分析，当时，权重和最接近1。式（10）可简化为：

式（11）的函数曲线如图 3所示。可以看出，f（χ）在χ=1.5时并非为0。意味着与q相距最远的网格点（两者χ和y方向的距离为1.5）也被分配到权重，即此时参与联合直方图统计的邻域点数目仍有9个。因此，无论q落在何处，能分配到权重的最多邻域点数与最少邻域点数的比例为9：9，比BPV，HPV，BHPV的比例还小。

图3 本文方法采用的高斯函数

将Vanderbilt图像配准数据集［12］中training-001的PET第11断层图像作为参考图像，将其沿χ轴平移-1.5～1.5像素作为浮动图像，同时为参考和浮动图像添加密度为0.1的椒盐噪声。采用基于PV，BPV，HPV，BHPV以及本文的插值方法计算参考和浮动图像的CCRE，平移步长为0.1像素，得到如图4所示的结果。可以看出，PV插值的CCRE在网格点位置出现了扰动，而其他4种改进方法的CCRE曲线都较光滑。随着噪声的增加，PV插值的CCRE在网格点位置将有更为明显的波动，这些局部极值会造成优化算法寻优的困难，甚至陷入局部极值造成配准失败。

图4 浮动图像沿χ轴平移的交叉累计剩余熵曲线

4 实验与结果分析

为进一步验证插值方法的性能，本文设计了2组实验。第1组实验的参考和浮动图像为training-001的PET第11断层图像，并添加了密度为0.3的椒盐噪声。设刚体变换参数T=［tχ，ty，θ］，其中，tχ，ty分别为浮动图像沿χ和y轴的平移，θ为绕图像中心点的旋转角度。两者真实的配准参数为T=［5，5，5］。使用步长为2的Pow ell算法对不同插值方法得到的CCRE进行寻优，结果如表1所示。

表1 不同插值方法对同幅PET图像的配准结果

从表1的配准结果可知，PV插值法的CCRE使得Powell难以收敛，迭代次数达到16次，最终得到错误的配准参数。BPV，HPV，BHPV和本文方法的CCRE都能得到正确的配准结果，且本文方法的χ轴平移误差最小。另外，本文方法的运行时间较BPV缩短了87.83%，较HPV缩短了92.99%，较BHPV缩短了98.25%，配准效率显著提高。

第2组实验的参考图像为BrainWeb数据集［13］中正常脑部磁共振 T1加权第90断层图像，浮动图像为相应断层的T2加权像，并根据文献［14］的方法添加标准差为22的莱斯噪声。两者真实的配准参数为T=［6，-8，5］。使用步长为5的Powell算法对5种不同的插值方法得到的CCRE进行寻优，结果如表2所示。

表2 不同插值方法对T1和T2图像的配准结果

同样，PV插值法得到错误的配准参数，而BPV，BHPV出现了配准失误。HPV和本文方法的CCRE都能得到正确的配准结果，且本文方法的精度最高，运行时间较HPV缩短了95.12%。

5 结束语

本文针对使用PV插值计算CCRE时在网格点容易产生局部极值的问题，提出一种基于高斯函数的改进PV插值方法。通过对高斯函数的方差及归一化系数的合理设计，可以满足权重分布对插值核函数的要求。该方法的优势为：（1）采用灵活的邻域中心来决定邻域范围，保证了权重分配的公平性。（2）分配到权重的最多邻域点数与最少邻域点数的比例达到1:1。（3）插值点权重分散到邻域9个点，相比16个点BPV，HPV，BHPV，大幅提高了配准效率。实验结果表明，将本文插值方法应用到CCRE中时，即使待配准图像受到强噪声影响，也能优化得到正确的参数，且运行速度较快。下一步将研究该插值方法在其他相似性测度中的应用。

［1］ Rao M，Chen Y，Vemuri B C，et al.Cumulative Residual Entropy：A New Measure of Information［J］.IEEE Transactions on Information Theory，2004，50（6）：1220-1228.

［2］ Wang Fei，Vemuri B C.Non-rigid Multi-modal Image Registration Using Cross-cumulative Residual Entropy［J］. International Journal of Computer Vision，2007，74（2）：201-215.

［3］ Hasan M，Pickering M R，Jia Xiuping.Robust Automatic Registration of Multimodal Satellite Images Using CCRE with Partial Volume Interpolation［J］.IEEE Transactions on Geoscience and Remote Sensing，2012，50（10）：4050-4061.

［4］ 江万寿，彭芳媛，岳春宇，等.利用Ratio梯度和交叉累积剩余熵进行多源遥感影像匹配［J］.武汉大学学报：信息科学版，2009，34（9）：1047-1050.

［5］ 冯 林，严 亮，贺明峰，等.图像配准中确定性扰动PV插值算法［J］.计算机辅助设计与图形学学报，2005，17（5）：908-914.

［6］ 彭景林，章 兢，李树涛.基于改进PV插值和混合优化算法的医学图像配准［J］.电子学报，2006，34（5）：963-965.

［7］ Lu Xuesong，Zhang Su，Su He，et al.Mutual Information-based Multimodal Image Registration Using a Novel Joint Histogram Estimation［J］.Computerized Medical Imaging and Graphics，2008，32（3）：202-209.

［8］ Lu Xuesong，Zhang Su，Su He，et al.Non-rigid Medical Image Registration with Joint Histogram Estimation Based on Mutual Information［J］.Transactions of Tianjin University，2007，13（6）：452-455.

［9］ 陈伟卿，华顺刚，欧宗瑛.互信息医学图像配准中 PV插值算法的改进［J］.计算机工程与应用，2010，46（20）：113-115.

［10］ 李 超，陈 钱，钱惟贤.基于交叉累计剩余熵的多光谱图像配准方法［J］.红外与激光工程，2013，42（7）：1866-1870.

［11］ Lehmann T M，Gonner C，Spitzer K.Survey：Interpolation Methods in Medical Image Processing［J］.IEEE Transactions on Medical Imaging，1999，18（11）：1049-1075.

［12］ Fitzpatrick M.The Retrospective Image Registration Evaluation Project Dataset［DB/OL］.［2014-08-01］.http：// www.insight-journal.org/rire/download-training-data.php.

［13］ Evans A C.BrainWeb：Simulated MRI Volumes for Normal Brain［DB/OL］.［2014-08-01］.http：//brainweb. bic.mni.mcgill.ca/brainweb/selection-normal.html.

［14］ Aja-Fernández S，Tristán-Vega A.Influence of Noise Correlation in Multiple-coil Statistical Models with Sum of Squares Reconstruction［J］.Magnetic Resonance in Medicine，2012，67（2）：580-585.

编辑 金胡考

Improvement of Interpolation Method in Cross Cumulative Residual Entropy-based Image Registration

XIANG Yan1，2，HE Jianfeng2，ZHANG Yunchun3，CAI Li3
（1.Key Lab of Computer Technologies Application of Yunnan Province，Kunming 650500，China；2.Institute of Information Engineering and Automation，Kunming University of Science and Technology，Kunming 650500，China；3.School of Software，Yunnan University，Kunming 650091，China）

The key strength of the Cross Cumulative Residual Entropy（CCRE）over the popular Mutual Information（M I）method is that the former has significantly larger noise immunity.But CCRE using conventional Partial Volume（PV）interpolation will result in the emergency of the local extremes on grid points，which may ham per the optimization algorithm from getting transformation parameters.In order to solve this problem，three improved PV interpolation methods are studied，including 3-order B-spline PV interpolation（BPV），Hanning windowed sinc PV interpolation（HPV）and Blackman-Harris windowed sinc PV interpolation（BHPV）.Meanwhile，a new interpolation method is proposed which uses flexible neighborhood center and makes the interpolation point to distribute the weight of the joint histogram to its adjacent 9 points.Moreover，it uses a Gaussian function as the PV interpolation kernel function to overcome the mutation of weight.Experimental result show s that the registration accuracy and speed of the proposed method is higher and faster，compared with BPV，HPV and BHPV method.So it is more suitable for CCRE computing.

image registration；Cross Cumulative Residual Entropy（CCRE）；Partial Volume（PV）interpolation；Gaussian function；3-order B-spline function；sinc function

1000-3428（2015）10-0199-05

A

TP391

相 艳，贺建峰，张云春，等.基于交叉累计剩余熵的图像配准中插值方法的改进［J］.计算机工程，2015，41（10）：199-203.

英文引用格式：Xiang Yan，He Jianfeng，Zhang Yunchun，et al.Improvement of Interpolation Method in Cross Cumulative Residual Entropy-based Image Registration［J］.Computer Engineering，2015，41（10）：199-203.

国家自然科学基金资助项目（11265007）；云南省计算机技术应用重点实验室开放基金资助项目（2403660120）。

相 艳（1979-），女，讲师、硕士，主研方向：医学图像处理；贺建峰（通讯作者），教授、博士；张云春，讲师、博士；蔡 莉，副教授、硕士。

2014-10-23

2014-11-24E-m ail：sharonxiang@126.com