基于小波变换医学图像融合算法的对比分析

胡俊峰 唐鹤云 钱建生（徐州医学院医学影像学院，徐州 006）

2（中国矿业大学信息与电气工程学院，徐州 221008）

引言

不同模态的医学图像都有各自的优缺点，如果通过图像融合技术将不同成像方式的互补信息综合在一起，就能为医学诊断和治疗提供更加充分有效的信息依据［1-3］。如何将这些表现细节信息不同的图像数据转化为更丰富和清晰的信息，就成为医学图像融合所要解决的关键问题。医学图像融合不仅广泛用于疾病诊断，而且在外科手术和放射治疗等的计划设计、方案实施以及疗效评估方面发挥着重要作用。

近年来，小波变换由于其多分辨率、时频局部等特性，被广泛应用于图像处理领域，并且成为图像融合领域的重要算法。小波变换［4］能够聚焦到分析对象的任意细节（即“变聚焦”特性），因此被人们誉为“数学显微镜”，并在信号处理、图像处理、模式识别、语音识别等众多非线性科学领域得到了成功的应用［5-6］。

基于小波变换的融合算法研究已经取得了很大进展，但每一种算法都只针对个别问题，普遍适用性算法的研究很少。现在大多数算法都是基于多聚焦图像提出的，没有考虑到医学图像灰度分布与多聚焦图像的不同。由于医学图像灰度分布不均匀，主要信息集中在低频部分，即使病变部位突变信息也不明显，需要重点考虑低频融合规则。文中介绍了小波变换的基本理论，以CT和MRI灰度图像为主要研究对象，通过调整低频和高频融合规则，在像素级上深入对比分析了各种融合规则对融合结果的影响。提出了低频能量取大、高频方差取大相结合的融合算法，比基于传统小波融合规则的融合质量及各项指标都有明显提高；在此基础上，提出了低频能量取大、高频系数绝对值取大相结合的融合改进算法，在各种算法比较中最优，并且验证了方法的有效性。理论分析和实验结果证明：选取合适的融合规则对融合结果影响很大，本研究提出的算法简单有效，对小波变换的医学图像融合算法的深入研究具有较大的指导意义。

1 基于小波变换的图像融合机理与过程

1986年，从事信号处理与研究工作的Mallat受到多尺度思想的启发，将其引入到小波函数的构造方法中，形成了统一的小波函数构造理论［6］，为小波函数的应用奠定了基础。同时，Mallat又研究了小波变换（wavelet translation，WT）的离散化问题，给出了与快速傅里叶变换（fast Fourier transform，FFT）相对应的快速小波算法—Mallat算法［5］，还把它应用到图像处理中。小波变换是一种类似于金字塔变换的图像处理方法，继承了金字塔变换的优点，也可以说，它是一种广义的金字塔变换方法。之后，小波变换理论在图像处理领域得到了较为广泛的应用。

1.1 小波变换的图像融合机理

对图像进行小波变换，实质上是利用小波变换所构造的高通、低通滤波器对图像进行高通、低通滤波，将其分解到不同频率下的不同特征域上。采用Mallat金字塔算法，只由系数就可以实现图像的分解与重构。

二维小波分解的Mallat快速算法，其分解公式为［4-6］

式中，h、v、d分别表示水平、垂直和对角分量，L为低通滤波（low-pass）算子，H为高通滤波（high-pass）算子，L*和H*分别为L和H的共扼转置矩阵，j为分解层数。

因此，基于小波分解将源图像分解成低频成分和高频成分，得到4个不同的频带，即iLL、iHL、iLH、iHH。低频成分iLL代表了图像的轮廓信息，高频成分包括iHL、iLH、iHH3个频带，代表了图像的细节信息，分别保留了原图水平、垂直和对角方向的高频信息，图1下标中1、2表示分解层数。

图1 二维图像的小波分解Fig.1 Decompose of Wavelet Transform

基于小波变换的重构，就是采用不同的融合规则和融合策略，进行小波系数融合和重构。对于低频部分融合策略，主要是提取原图的轮廓信息，由于图像的能量主要集中在图像的轮廓部分，经提取获得的子带图像iLL部分较清晰。对iLL继续进行小波分解，j层小波分解后可得到（3j+1）个频带。对于高频iHL、iLH、iHH部分，融合策略主要是提取图像的细节信息，包括边缘、纹理等特征，所以较黑暗。

相应的小波重构算法为

1.2 小波变换图像融合的一般过程

小波变换是空间和频率的局部转换，能有效地从图像中提取信息，达到高频处空间细分、低频处频率细分的目的，其融合过程如图2所示，基本步骤如下［7］。

步骤1：取两幅已配准的图像，分别标记为图像A和图像B。

步骤2：选取合适的小波基函数。

步骤3：将待融合的原始图像进行j层小波分解，得到3j个位于不同尺度、具有不同空间分辨率和频率特性的高频子图像系数和1个位于最高层（j层）的低频子图像系数。

步骤4：选取合适的融合算法。对两幅图像所对应分解层的不同频率分量进行融合，得到融合后的3j个高频子图像和1个低频子图像。由于低频区域和高频区域所代表的意义不同，相应地所采用的融合算法（融合规则）也可不同。

步骤5：将融合后的子图像重构（小波逆变换），得到融合后的融合结果图像。

图2 小波分解融合过程Fig.2 Schematic diagram of wavelet fusion

2 小波变换常用的融合规则

基于小波变换的图像融合方法的关键技术之一在于融合规则的选取。近年来，人们进行了大量小波图像融合算法研究，获得了多种小波系数的融合规则。由于小波分解的本质是采用不同的滤波器，将源图像分解成低频成分和高频成分，低频成分代表了图像的概貌信息，高频成分代表了图像的细节信息。因此，应该根据其系数矩阵的各自特点，采用不同的融合规则和融合策略，以期在融合图像最大限度地获取源图像的信息。图像重构实际上是对融合后的小波系数矩阵进行相应的小波逆变换，即利用相应的小波重构滤波器组（高通、低通滤波器）对融合后的小波系数进行滤波，这样就可以得到融合后的图像。下面介绍几种常见的小波系数融合规则。

2.1 低频系数融合规则

通过小波分解得到的低频成分代表了图像的概貌信息，所以低频系数都是正的变换值。常用的小波低频系数的融合规则有均值法、加权平均法、极大值法等。由于低频系数都是正的变换值，在变换时可以根据具体的图像和目的对该变换值进行取均值，或者取极大值。

1）均值法（average，AVE）［8］，即融合系数是各源图像对应系数的均值。

2）加权平均（weighted average，WAV）准则，即将两幅输入源图像A、B的低频子带系数各自乘上一个权重系数，融合而成新的图像F，该算法实现简单。

3）选取系数绝对值（coefficient absolute value，CAV）最大的准则［9］。

2.2 高频系数融合规则

通过小波分解得到的高频成分代表了图像的细节信息，也就是图像中的显著特征点，如边缘、亮线及纹理，其中包含了3个高频子带。3个高频子带的系数是一些在零附近的变换值，在这些子带中，亮度急剧变化的点的对应的变换值较大。为了获取尽可能多的细节信息，需要对高频系数采用不同的融合规则和融合策略，常用的高频系数融合规则有如下

1）小波分解系数绝对值极大法（coefficient absolute value，CAV）。

2）基于区域的最大值法（region coefficient absolute value，RCAV）。基于区域的方法是在选定区域根据特定融合规则取得该区域像素最大值的算法，是小波分解系数绝对值极大法的变形。由于图像中某一局部区域内的各像素间往往具有较强的相关性，在该区域选取像素最大值，就能更好地获取源图像的细节信息。相比在整幅图像取最大值算法，区域内像素点的相关性更强，提取出的边缘、纹理特征更加准确可信，可有效减少虚假信息和伪轮廓，从而获得视觉效果更佳、细节更丰富的融合效果。

3）基于区域能量的图像融合算法（region energy，REN）［10］。基于区域能量融合依然是以图像像素值来衡量，所不同的是引入了区域能量的概念。首先计算选定区域高频成分的区域能量及匹配度，然后比较匹配度与给定阈值的大小，进而比较匹配度与区域能量的大小。

4）基于边缘强度的自适应融合法（edge intensity，EDI）［11］。基于边缘强度的自适应融合法的基本思想：选定M×N窗口，计算像素点的边缘，并以边缘强度为权重，对两幅源图像的高频系数进行求和。

5）区域方差（region variance，RVA）最大准则［12］。

3 改进的基于小波变换的融合规则

通过对基于小波变换融合规则进行深入研究，发现目前各种文献中常用小波变换融合规则在医学图像融合中效果并不十分理想。在理论分析、实验研究并与文献对比的基础上，由于图像的能量主要集中在图像的低频部分，本研究提出低频能量取大、高频方差取大相结合的融合算法，进而又提出了低频能量取大、高频系数绝对值取大相结合的融合算法，以验证融合规则选取的重要性。

3.1 低频系数融合规则

选取基于区域能量取大的融合规则（REN）。

Ek，j（x，y）表示低频系数矩阵中以（x，y）为中心的区域大小为M×N的能量，定义为

区域均值的定义为

定义两幅图像对应区域能量的匹配度为

若Matk，j≤Thr，说明两图像区域能量相差较大，则融合后的低频系数为

若Matk，j＞Thr，说明两幅图像的区域能量相差较小，则融合后的低频系数为

式中，Thr为图像匹配度阈值（本研究选取0.7，由实验获取），Matk，j为两幅图像对应区域能量的匹配度；k代表图像A、B、F，分别为源图像A、B和融合图像F；J为分界解层数，区域大小选为3×3，下同。

3.2 高频系数融合规则

分别采用基于方差取大（RVA）的小波变换融合方法和小波分解系数绝对值极大法（CAV）。

3.2.1 区域方差（RVA）最大准则［12］

3.2.2 小波分解系数绝对值极大法（CAV）

如下式，有

4 实验研究

4.1 融合实验步骤

步骤1：选择已配准后的CT、MRI灰度图像（见图3中（a）和（b））作为融合源图像，PNG格式［7］，大小为256像素×256像素，实验平台为Mobile AMD SempronTMProcessor 3200+1.6GHz，512MB内存，Matlab7.01。

步骤2：选取小波基和分解层数。融合实验的目的是要对比分析各种融合规则的融合性能，以期提出与检验新的改进算法。因此，在实验中需要首先设定参与融合的小波基和分解层数，以保证融合条件的一致性。实验和文献表明，Db3小波基在图像的分解中效果最优［13-14］，本研究选取该小波基参与融合。医学灰度图像信息主要集中在低频部分，高频细节信息不丰富，所以分解层数选取为2层；多聚焦灰度图像信息分布一般较为均匀，高低频信息较平均，所以选择5层分解。

步骤3：讨论融合规则对融合性能的作用，从而提出新的有效改进融合算法。选取目前常用的图像融合规则进行组合，获得融合结果，并与本研究提出的融合算法比较，从中选出最优算法。

步骤4：为了验证上述实验结果的正确性和普遍有效性，选取Lena标准图像，大小为256像素×256像素，采用左右分别聚焦的方式获得多聚焦源图像（图4中（b）和（c）），重点步骤3进行各种融合规则的融合性能对比实验，以验证改进融合算法的有效性。

步骤5：采用临床图像CT、MRI灰度图像（图7中（a）和（b））作为融合源图像，重做步骤3进行各种融合规则的融合性能对比实验，以验证改进融合算法的普遍有效性。

步骤6：对融合结果进行综合评价。

4.2 融合性能评价方法

融合评价采用主观视觉分析与客观标准共同进行，客观标准［15］选取均值、标准差、信息熵、平均梯度、空间分辨率，分别对每组进行融合结果评价；为了对融合结果边缘信息进行有效检测，增加了边缘信息评价因子［16］。

边缘信息评价因子是一种新的基于边缘信息的性能评价因子，反映了源图像与融合图像边缘信息的传递量，大小越接近1说明边缘传递越好，融合效果也就越好。其公式为

式中，QABF表示融合图像F相对于源图像A、B的整体边缘保留量，QAF和QBF分别表示融合图像F相对于源图像A、B的边缘保留量，gA和gB分别表示源图像A，B的边缘强度。

5 实验结果与分析

表1给出了不同融合规则相互组合融合结果的客观性能评价，图3为相应的融合结果图像。使用的客观评价标准分别是均值、标准差、信息熵、平均梯度、空间频率和边缘信息评价因子。其中：图像均值反映融合图像像素的灰度平均值，即平均亮度；标准差反映了图像灰度相对于灰度平均值的离散情况，值越大表明图像灰度级分布分散，可以看出更多的信息；信息熵表示图像中所包含的信息，值越大越好；互信息表明融合图像从源图像中所获取的信息，值越大融合效果越好；平均梯度在整体上反映图像细节反差程度和纹理变化特征的能力；空间频率反映图像空间域的总体活跃程度越高影像越清晰；边缘信息评价因子反映源图像与融合图像边缘信息的传递量，大小越接近1说明边缘传递越好，融合效果也就越好。

表1 不同融合规则对图像融合性能指标影响比较Tab.1 Comparison the image fusion of different fusion rules

图3 不同规则图像融合效果。（a）CT图像；（b）MRI图像；（c）AVE＆CAV融合图像；（d）CAV＆CAV融合图像；（e）CAV＆REN融合图像；（f）CAV＆EDI融合图像；（g）CAV＆RVA融合图像；（h）REN＆CAV融合图像；（i）REN＆RVA融合图像；（j）RVA＆REN融合图像Fig.3 Image fusion results of different fusion rules.（a）CT image；（b）MRI image；（c）AVE＆CAV rule fusion image；（d）CAV＆CAV rule fusion image；（e）CAV＆REN rule fusion image；（f）CAV＆EDI rule fusion image；（g）CAV＆RVA rule fusion image；（h）REN＆CAV rule fusion image；（i）REN＆RVA rule fusion image；（j）RVA＆REN rule fusion image

由图3和表1中的数据可知，选取不同的融合规则，对融合性能有较大的影响。低频平均法（AVE），不论高频采用绝对值取大（CAV）、能量取大（REN）［9］、边缘强度取大（EDI）还是方差取大（RVA）融合性能表现为均值、标准差、熵、平均梯度、空间频率和边缘评价因子较低，融合图像表现为较暗，灰度及边缘保持较差（见图3（c））。低频绝对值取大（CAV）、高频采用上述融合规则组合有一定优势，融合性能的各项指标都有提高，图像细节被加强，高频表现突出，图像显得比较尖锐。其中，低频和高频都采用最大值融合规则的，在本组表现最优（见图3（d）），分别优于文献［10-12］中提出的算法（见图3中的（e）、（f）和（g））。但是，低频和高频都采用最大值融合规则只是对图像像素值直接取大，方法太过简单，造成图像亮度和对比度过大，当两幅图像灰度值差别较大时容易忽略其中灰度值小的一副图像信息。低频能量取大（REN）与高频上述融合规则组合同样有一定优势，融合性能的各项指标比低频平均法（AVE）都有提高，图像细节被加强，高频表现突出，图像显得比较尖锐；在与边缘强度取大（EDI）融合规则组合时会产生边界伪影。其中，低频能量取大（REN）与高频采用最大值（CAV）融合规则的在本组表现最优（见图3（h））。相比低频绝对值取大（CAV），均值略微减小、平均梯度稍微升高，变现为亮度稍微降低，图像细节更加清晰。而低频方差取大（RVA）与上述各种高频规则组合时，融合结果虽然数据有较好的表现，但是融合图像表现了不同程度的方块效应（见图3（j））。产生这些不同的原因首先因为图像的能量主要集中在低频部分，因此对低频部分采用基于区域能量取大的融合规则，可以尽可能多地获取图像的主要信息。医学图像更是如此，从图4中可以看出，CT、MRI图像直方图分布集中、层次感不强，而LENA图像直方图分布分散、层次感较强。这样，导致了医学图像分解后高频部分信息不多，因此对于高频部分可直接比较像素大小，采用直接取大的算法，就可以获得足够的高频信息。另外，由于小波分解合成都要进行边界延拓，这些损失的信息量均是小波逆变换不能恢复的损失，如果对低频采用方差取大的算法，将导致对边界的过度开发，虽然边缘因子增加，但是融合图像会出现方块效应，放大后尤为明显。

实验结果说明3点：一是在对融合图像对比度要求不高时，低频平均法（AVE）、高频绝对值取大（CAV）可以满足基本需要；二是低频方差取大效果可产生方块效应，总的来说低频部分采用能量算法优于最大值法和方差最大法；三是高频方差取大优于高频能量取大。

综合以上3点得出如下结论：低频能量取大（REN）、高频方差取大（RVA）的融合规则组合能有效地保证融合质量，提高融合效果（见图3（i））。但是，通过对比实验及主客观评价结果，可以看到低频能量取大（REN）、高频绝对值取大（CAV）的融合规则组合（见图3（h））优于低频能量取大（REN）、高频方差取大（RVA）的融合规则组合（见图3（i）），说明高频选择较简单算法仍有优势。其原因是医学图像融合时分解层数可以稍低，Db3小波可以做1～2层分解，基于方差取大等复杂算法的优势在高频部分无法充分发挥。

图4 不同图像的灰度直方图分布。（a）CT图像；（b）MRI图像；（c）Lena图像Fig.4 Histogram of the images.（a）CT image；（b）MRI image；（c）Lena image

图5 不同规则图像融合效果。（a）Lena图像；（b）左聚焦图像；（c）右聚焦图像；（d）CAV＆CAV规则融合图像；（e）CAV＆REN融合图像；（f）CAV＆RVA融合图像；（g）REN＆RVA规则融合图像；（h）REN＆CAV规则融合图像Fig.5 Image fusion results of different fusion rules.（a）Lena image；（b）right focus image；（c）left focus image；（d）CAV＆CAV rule fusion image；（e）CAV＆REN rule fusion image；（f）CAV＆RVA rule fusion image；（g）REN＆RVA rule fusion image；（h）REN＆CAV rule fusion image

图5给出了图3和表1中各种算法规则验证实验的融合结果，表2给出了客观评价指标。可以看出，低频能量取大（REN）、高频绝对值取大（CAV）的融合规则组合（见图5（g））在各种算法组合中性能最优，并优于低频能量取大（REN）、高频方差取大（RVA）的融合规则组合（见图5（h））。实验结果同样表明上述结论的正确性。为了能够清晰地表明实验结果，把图5中各图的眼部细节分别放大（见图6），直观显示了融合结果，从而进一步验证了本算法的有效性。

图7给出了临床CT、MRI图像各种算法规则验证实验的融合结果，表3为客观评价指标参数。可以看出，低频方差取大（CAV）、高频绝对值取大（CAV）的融合规则组合（见图7（c））和低频绝对值取大（CAV）、高频方差取大（RVA）的融合规则组合（见图7（d））在各种算法组合中性能较优。但是，主观视觉评价在临床无法使用，图7（c）像素值取大造成了MRI图像信息基本消失，而图7（d）高频方差取大进一步强化了边界信息。低频能量取大（REN）、高频绝对值取大（CAV）的融合规则组合（见图7（f））在各种算法组合中性能最优，并优于低频能量取大（REN）、高频方差取大（RVA）的融合规则组合（见图7（e））。实验结果同样表明上述结论的正确性和算法的有效性。

表2 不同融合规则对图像融合性能指标影响比较（Lena）Tab.2 Comparison the image fusion of different fusion rules（Lena）

图6 不同规则图像融合效果。（a）Lena图像；（b）左聚焦图像；（c）右聚焦图像；（d）CAV＆CAV规则融合图像；（e）CAV＆REN融合图像；（f）CAV＆RVA融合图像；（g）REN＆RVA规则融合图像；（h）REN＆CAV规则融合图像Fig.6 Image fusion results of different fusion rules.（a）Lena image；（b）right focus image；（c）left focus image；（d）CAV＆CAV rule fusion image；（e）CAV＆REN rule fusion image；（f）CAV＆RVA rule fusion image；（g）REN＆RVA rule fusion image；（h）REN＆CAV rule fusion image

6 讨论与结论

与传统的基于小波融合算法和相关文献报道的融合算法相比，本研究提出的小波改进算法考虑了医学图像灰度分布与多聚焦图像的不同，分析了不同的高低频融合规则的特性，根据具体图像特点，选取不同的高低频组合融合规则。由于医学图像灰度分布不均匀，主要信息集中在低频部分，即使病变部位突变信息也不明显，需要重点考虑低频融合规则。从各种低频和高频组合融合规则的仿真实验结果可知，对低频部分采用基于区域能量取大的融合规则，可以尽可能多地获取图像的主要信息，而基于区域方差取大的融合规则对图像纹理特征比较敏感。所以，本算法可以取得明显的效果。但是，由于本研究所选用的图像无论是标准的CT/MRI图像，还是Lena多聚焦图像和临床实际获取的CT/MRI图像都是灰度图像，如果选用彩色PET/CT图像或遥感彩色图像，本方法是否有效还需要进一步实验验证。

小波变换能把图像各自携带的不同特征与细节分解到不同尺度下，分别作为源图像的近似信息和细节信息，在多个分解层、多个频带上分别以不同算子进行融合。通过调整小波变换低频和高频融合规则，深入对比分析了各种融合规则对医学图像融合性能的影响。提出了基于小波变换的融合改进算法，在各种算法比较中最优，并且验证了方法的有效性。理论分析和实验结果证明：选取合适的融合规则对融合结果影响很大，所提出的算法简单有效。实验结果对小波变换的医学图像融合算法的深入研究具有较大的指导意义。

图7 不同规则图像融合效果。（a）CT源图像；（b）MRI源图像；（c）CAV＆CAV融合图像；（d）CAV＆RVA融合图像；（e）REN＆RVA融合图像；（f）REN＆CAV融合图像Fig.7 Image fusion results of different fusion rules.（a）CT image；（b）MRI image；（c）CAV＆CAV rule fusion image；（d）CAV＆RVA rule fusion image；（e）REN＆RVA rule fusion image；（f）REN＆CAV rule fusion image

表3 不同融合规则对图像融合性能指标影响比较Tab.3 Comparison the image fusion of different fusion rules

［1］Matsopouios GK，Marshall S，Brunt J.Multi-resolution morphological fusion of MR and CT images of the human brain［J］.IEEE Transactions on Image and Signal Processing，1994，141（3）：137-142.

［2］Wang Anna，Sun Haijing，Guan Yueyang.The application of wavelet transform to multi-modality medical image fusion//Proceedings of IEEE international Conference on Networking，Sensing and Control.Chicago：IEEE，2006：270-274.

［3］张泾周，李婷，吴疆.医学图像的小波变换融合算法研究［J］.中国生物医学工程学报，2008，27（4）：521-525.

［4］Mallat S.G.，A theory for multiresolution signal decomposition：the wavelet representation［J］.IEEE Trans.on Pattern Analysis and Machine Intelligence，1989，11（7）：674-693.

［5］Mallat S G，Zhong S.Characterization of signal from mutiscale edges［J］.IEEE Trans on Pattern Analysis and Machine Intelligence，1992，14（7）：710-732.

［6］Jian Muwei，Dong Junyu，Zhang Yang.Image fusion based on wavelet transform//Proceedings of the 8th ACIS International Conference on Software Engineering，Artificial Intelligence，Networking，and Parallel/Distributed Computing.Chicago：IEEE，2007：713-718.

［7］胡俊峰，唐鹤云，钱建生.不同图像格式在医学图像融合中的性能研究［J］.中国生物医学工程学报，2010，29（2）：310-313.

［8］Li H，Manjunath B S，Mitra S K.Multisensor image fusion using the wavelet transform［J］.Graphical Models and Image Processing，1995，27（3）：235-244.

［9］刘贵喜，杨万海.基于小波分解的图像融合方法及性能评价［J］.自动化学报，2002，28（6）：65-72.

［10］Zhang Z，Blum R.S.A Categorization of multiscaledecomposition-based image fusion schemes with a performance study for a digital camera application［J］.Proceedings of the IEEE，1999，87（8）：1315-1326.

［11］郭志强.基于区域特征小波变换图像融合方法［J］.武汉理工大学学报，2005，27（2）：61-63.

［12］黄彩霞，陈家新.基于小波系数区域相似度的医学图像融合［J］.计算机应用研究，2008，25（1）：274-276.

［13］刘凌湘，陈武凡，周凌宏，等.基于小波变换的CT/PET图像融合最佳参数研究［J］.中国生物医学工程学报，2010，29（4）：498-503.

［14］胡俊峰，唐彩银，巩萍.基于小波变换的CT/SPECT图像融合最佳层数选取［J］.中国医疗设备，2009，24（3）：10-12.

［15］Wang Zhijun，Djemel Z，Costas A，et al.A comparative analysis of image fusion methods［J］.IEEE Transactions on Geoscience and Remote Sensing，2009，9（18）：2137-2143.

［16］Xydeas CS，Petrovic V.Objcetive image fusion performance measure［J］.Electronics Letters，2000，36（4）：305-309.