基于经验模态分解提取纹理的图像融合算法

张宝华，刘 鹤，张传亭

(内蒙古科技大学信息学院，包头014010)

引 言

随着成像技术的突飞猛进，各类精密成像设备推动了医学影像学的发展，为临床提供丰富的人体医学影像，但成像设备种类众多，其成像机理不同，反映医学信息各有侧重。为了全面分析医学图像中包含的解剖信息和功能信息，需要对多模态医学图像进行融合。

医学图像融合技术面向多模态医学图像，把各种医学图像的信息有机地结合，完成各类医学信息融合，不仅有效利用已有医学影像，而且还有助于发掘潜在医学信息，辅助医院诊疗。

欧美国家在20世纪80年代已开始在图像融合领域取得丰硕的成果，处于领先地位;MITIANOUDIS等人［1］通过将稀疏表达应用于系数选择中，拓展了图像融合的视野和手段;国内YAN，LI等人［2-3］提出了基于脉冲耦合神经网络(pulse coupled neural net-work，PCNN)的图像融合方法，应用于不同类型的图像，将神经网络作为系数优选的方法，丰富了融合系数选择方法，取得了较好的融合效果。

目前图像融合方法包括基于空域变换和频域变换的方法，其中以小波变换和各类超小波变换为代表的基于多分辨分析的图像融合方法应用最广。但小波变换和其改进方法依赖于预先定义的滤波器或基函数，小波变换会有下采样操作，变换后图像会引进伪吉布斯现象，降低融合图像质量。HUANG等人为达到对非线性和非稳态的数据进行自适应和多尺度分析，提出了经验模态分解(empirical mode decomposition，EMD)［4］。经验模式分解作为一种新的多尺度图像分解方法［5-6］，具有比小波分析［7］更直观的特征表示方式和更灵活的频率特性，避免了分解中引进冗余信息，同时EMD对于图像细节保护和图像纹理的提取等方面具有优势［8-11］，适合于对安全性要求较高的医学图像进行多分辨分析。

本文中将2维经验模态分解(bidimensional empirical mode decomposition，BEMD)运用到医学图像特征提取中，通过将BEMD分解后的子图像和趋势图输入神经网络获取它们的点火映射图，提取不同分解层对应医学图像特征。之后将对应于图像纹理信息和背景信息的系数分别通过PCNN和双通道PCNN(daul channel pulse coupled neural network，DCPCNN)选取融合系数。由于区别对待代表图像纹理和背景信息的像素，既保护了图像中的特征，又有效改善了PCNN在医学图像系数选择中的效果。

1 基于BEMD和点火映射图的特征分类

1.1 BEMD

EMD分解具有优越的空间和频率特性。其易拓展性也推动了BEMD方法的出现，BEMD同样具有数据驱动和良好的自适应性等特点，而且具有多尺度特性。将BEMD用于医学图像处理，可得到源图像在不同频率的2维内蕴模函数(bidimensional intrinsic mode functions，BIMF)和残差项 R。

BEMD是数据驱动的，在BEMD分解过程中，当筛分终止条件一定时，图像分解得到的BIMF个数与图像数据本身的特征相关［12］，在纹理分析等领域中得到广泛的应用。

BEMD分解的算法步骤如下:

(1)初始化，设源图像为 I(x，y)(其中 x=1，2，…，P;y=1，2，…，Q)，残差项为 Rj(x，y)，j为分解层数，当 j=0 时，R(x，y)=I(x，y)。

(2)如果残差项R(x，y)单调或达到图像的分解层数，则算法停止;否则，令:Hk-1(x，y)=Rj-1(x，y)。若 k=1，即 H0(x，y)=Rj-1(x，y)，Hk(x，y)为第k级分解系数去掉Rj(x，y)的残余系数，k为分解层数，进入筛选过程。

(3)求图像I(x，y)极值点:将解空间点集分为区域极大值点集和极小值点集。

(4)平面插值区域极值点集，得到图像的上、下包络面 Eup(x，y)和 Edown(x，y)，进一步求出图像Hk(x，y)的均值:

(5)令 Hk(x，y)=Hk-1(x，y)-Emean(x，y)，判断筛选过程是否满足停止条件S，如果不满足则转步骤(3)，其中S为:

(6)判断Hk(x，y)是否为内蕴模函数(intrinsic mode functions，IMF)，其依据是:若 S ＜ ξ(ξ为阈值，本文中 ξ=0.2)，则 Dj(x，y)=Hk(x，y)，否则令 k=k+1，转到步骤(2)。

(7)求残差项Rj=Rj-1-Dj。若Rj中间部分极值点数大于2或分解得到的IMF数目小于指定值，将 Rj转到步骤(2)，j=j+1。

(8)得到2维分解表达式为:

以上步骤中，Dj是第j个2维内蕴模函数，Rj是经过j层分解后的趋势图像。

1.2 复合型脉冲耦合神经网络

PCNN是一种新型神经网络，是依据动物的大脑视觉皮层上同步脉冲振荡现象提出的［12］。PCNN模型中，相邻的一群神经元可通过发放同步脉冲传送激励。当一个或数个神经元点火，输出的同步脉冲将被传送到相邻的神经元，使之迅速点火，神经元激发的过程称为点火，将一幅图像的像素映射到神经网络中，就获得了图像的点火映射图。

通过PCNN神经网络来寻找各个像素间的联系，将存在联系的像素进行归类，进而确定这一类像素的特征。根据一副图像在神经网络中各个像素被激发点火次数的不同，进而将同一点火次数的像素分为一类。这就找到了图像中像素间的联系。图1、图2分别为计算机断层扫描成像图像(computed tomography，CT)和磁共振成像图像(magnetic resonance imaging，MRI)头部源图像。PCNN根据点火次数将像素按灰度值的大小进行分类如图3所示，图3a～图3c分别为图1的1～3次点火映射图，通过比较可以看到相同次数的点火点组成的图像完整表达了图像轮廓信息，点火次数越多对应的细节越清晰。

Fig.1 CT image

Fig.2 MRI image

Fig.3 Firingmap

由于PCNN对图像中偏暗区域特征的筛选效果较差，而DCPCNN见图4，图中SA和SB代表输入激励;FA和FB表示反馈通道输入;Lij表示连接输入项;Wij为突触联接权;Vl和Vt为归一化常数;Uij表示神经元的内部活动项;βA和βB表示连接强度;Yij表示神经元的脉冲输出，它的值为0或者1;Tij是动态阈值;αl和αt为调节对应式子的常量，maximum是比较器。对图像中偏暗或偏亮区域特征的处理效果较好［13-14］，将图像根据特征分布情况进行分类后分别通过PCNN和DCPCNN可以较好地提取图像中各个范围的特征，本文中将PCNN和DCPCNN的复合结构定义为复合型脉冲耦合神经网络。

Fig.4 Dual-channel PCNN model

1.3 基于点火映射图的特征提取

图像中的像素可类比成一个神经元，整个图像就是一个复杂的神经网络系统，将图像的灰度值输入网络对应的神经元，当一个像素被激发时，与它相邻的像素点就会受到影响，被激发或者保持熄灭状态。通过多次点火激发就可以根据各个像素点火次数，来确定图像中每个点与附近的点或者偏离它较远的点有没有相互间的脉冲联系。依据一副图像中各个像素点的点火次数将一副图像的像素进行分类，点火次数相同的像素归为一类，从而将图像依据神经网络中各神经元激发过程中联系强弱的特征映射到图像中对应的像素，得到了各个像素间的相互联系。找到不同分类的像素的灰度极值就可以确定各类像素灰度分布范围。

BEMD的分解过程实现了图像从低频到高频多尺度分离过程。首先分解出来的固有模态分量BIMF是图像所含有的最高频率分量，该分量的各处频率都对应着图像在各处的局部最高频率。各个模态分别从不同频率反映了一幅图像的特征，再通过PCNN得到图像各频段的点火映射图，根据BEMD逆变换获得一副图像总的点火映射图。图5a～图5h表示图像经过BEMD分解得到的分解图像及其点火映射图，图5a～图5d是图1对应的BIMF和残差项，图5e～图5h分别是图5a～图5d对应的点火映射图，由于表示图像特征信息的像素点在点火时会受到多次激发，因而通过将不同点火次数的像素分类即可提取图像的特征信息。

Fig.5 Firingmaps of BIMF and R

2 基于特征分类的图像融合方法

利用CT和MRI多模头部医学图像进行融合，通过利用BEMD特征提取将图像区域分为纹理和背景两部分，将两区域分别建立融合规则选择融合系数，由于轮廓、纹理等信息被较好保护，本方法具备了BEMD和PCNN两者的优势，改善了融合图像质量。

2.1 融合步骤

本文中提出的基于特征提取和复合型脉冲耦合神经网络的医学图像融合方法如下:

(1)将待融合的2幅图像分别通过BEMD得到BIMF和残差项。

(2)BIMF和残差项分别输入到PCNN中，分别得到各自点火映射图，并将各自的点火映射图相加得到总的点火映射图。

(3)将原始图像中点火次数相同的像素点归为一类，若最大点火次数为N，则像素点共分为N类(N为自然数)，由于点火次数高的像素一般对应于图像的纹理，本文中将点火次数为n-N的分类定义为纹理类，剩余类定义为背景类，其中:

Fig.6 Algorithm flowchart

式中，k属于正整数。

(4)根据各分类集合的像素极值确定各个分类灰度分布范围(极小值～极大值)。

(5)将两幅图像中纹理类的灰度分布范围求交集，即若源图像A和B的最大点火次数为6，其点火次数为4～6的分类为纹理类，若集合中像素点灰度分布范围分别为(25～190)和(60～210)，则取两者交集(60～190)范围对应的像素，通过PCNN进行选择融合。

(6)其余像素对应图像的偏暗和偏亮部分，通过DCPCNN进行融合。

(7)通过PCNN和DCPCNN得到的融合系数重构得融合图像。

上述方法的流程图如图6所示。

2.2 融合系数选择

2.2.1 纹理区域融合系数选择 通过PCNN选择纹理区域的融合系数，能较好提取图像中纹理信息。根据像素点映射的神经元激发所产生的的点火次数的大小作为像素点优选的指标，来选择两幅图像中对应位置的融合系数。

2.2.2 背景区域融合系数选择 采用DCPCNN选择背景区域融合系数，DCPCNN可改善PCNN在医学图像中偏暗区域特征选择的效果，与传统的单通道PCNN相比，DCPCNN是由两个简化PCNN并行组成，首先通过计算以像素点O(i，j)为中心位置的3×3邻域中任意3个点的和与其它任意3个点的差值，得到其中最小值和最大值，将最大值和最小值的差 W，经过 1/［1+e-γ×W(i-1，j-1)］(γ 是调节量)运算得到O(i-1，j-1)的β值。通过选择两个通道中神经元的内部活动项Uij来控制像素点的点火状态。从而根据Uij选择两幅图中像素点Uij大的作为融合图像的像素点。

3 实验结果及分析

为了比较不同融合方法的融合效果，选择两组反映不同医学信息的源图像进行融合实验，融合前已利用基于光流场的方法对源图像进行配准，配准误差小于1个像素，将融合图像与其它5种传统融合方法进行比较:基于Laplace、离散小波变换(discrete wavelet transform，DWT)、主成分分析(principal component analysis，PCA)、gradient pyramid 和 contrast pyramid的融合方法，Laplace方法进行3层分解;DWT采用Daubechies小波基进行3层小波分解，高频采用取大低频取平均的方法。第1组实验图像为图1和图2，对其采用6种融合方法的比较结果见图7;第2组实验图像为图8和图9，6种融合方法的比较结果见图10，可以明显看到，6种融合方法都实现了源图像的信息融合，丰富了图像信息，但基于DWT和contrast pyramid方法的融合图像部分区域受伪影干扰，边缘不清晰，基于PCA方法的融合图像局部失真明显。由本文中方法得到的融合图像更加接近于源图像，细节清晰完整，视觉效果更好。

Fig.7 Fusion results comparison of group 1

Fig.8 CT image

Fig.9 MRI image

Fig.10 Fusion results comparison of group 2

作者选取了4个指标对融合图像进行客观评价，分别是互信息(mutual information，MI)、边缘强度 QA，B，F(A，B 代表源图像，F 代表融合图像)、熵和平均梯度。互信息衡量融合图像中包含源图像的信息的程度，QA，B，F，表达的是融合图像含有源图像边缘信息的丰富程度，信息熵可以衡量图像携带的信息量，值越大，说明融合图像包含信息量越丰富;平均梯度反映图像中的特征细微变化，平均梯度大的图像清晰程度较好，以上指标越大说明融合效果越好。表1、表2中分别显示两组实验融合图像的客观评价指标比较结果，可以看到，本文中的方法评价指标均优于其它方法，通过主观评价和客观指标的比较，说明运用本文中方法得到的融合图像含有更多源图像信息，图像纹理较丰富，细节突出，融合效果更好。

Table 1 The first group of image fusion evaluation objective comparison

对表1中的实验结果进行比较可以得出，本文中的方法与经典Laplace、DWT相比，互信息分别提高了 2.7 倍和 2.3 倍，QA，B，F分别提高了 14% 和30.8%，信息熵分别提高了14.9%和9.5%，平均梯度分别提高了0.7%和5%。与其它方法相比，本文中方法也有明显优势，表2和图10对第2组实验的客观指标也进行了比较，可以得到和第1组相同的结论，即本文中的方法在客观评价中有较好的效果。

Table 2 The second group of image fusion evaluation objective comparison

4 结论

利用BEMD分析方法将医学图像分解为不同频率子图像，通过提取特征建立基于复合型脉冲耦合神经网络的融合规则，由于克服了PCNN在偏暗图像中细节无法有效提取的问题，本算法能够从医学图像中提取更多的纹理信息。从主观和客观方面与其它图像融合方法进行了实验结果比较，证明了该方法在保留边缘、纹理、细节信息的有效性。

［1］MITIANOUDIS N，STATHAKIT.Optimal contrast correction for ICA-based fusion of multimodal images［J］.IEEE Sensors Journal，2008，8(12):2016-2025.

［2］YAN Ch M，GUO B L，YI M.Multifocus image fusion method based on improved LP and adaptive PCNN［J］.Control and Decision，2012，27(5):703-708.

［3］LIM L，LIY J，WANGH M.Fusion algorithm of infrared and visible images based on NSCT and PCNN［J］.Opto-Electronic Engineering，2010，37(6):90-95.

［4］HUANG N E，SHEN Z，LONG SR，et al.The empiricalmode decomposition and the Hilbert spectrum for nonlinear and nonstationary time series analysis［J］.Proceeding of Royal Society，1998，A454(6):903-995.

［5］FLANDRIN P，RILLING G.Empiricalmode decomposition as a filter bank［J］.IEEE Signal Processing Letters，2004，11(2):112-114.

［6］CHENG JS，YU D J，YANG Y.Research on the intrinsic mode function(IMF)criterion in EMDmethod［J］.Mechanical Systems and Signal Processing，2006，20(4):817-824.

［7］YUW J，GUGH，YANGW.Fusion algorithm of infrared polarization images based on wavelet transform［J］.Laser Technology，2013，37(3):289-292(in Chinese).

［8］DAMERVAL C，MEIGNEN S，PERRIER V.A fast algorithm for bidimensional EMD［J］.IEEE Signal Processing Letters，2005，12(10):701-704.

［9］NUNES JC，BOUAOUNE Y，DELECHELLE E，et al.Image analysis by bidimensionalempiricalmode decomposition［J］.Image and Vision Computing，2003，21(12):1019-1026.

［10］ZHANG Y，SUN Zh X，LIW H.Texture synthesis based on directional empiricalmode decomposition［J］.Journal of Computer Aided Design ＆ Computer Graphics，2007，19(4):515-520(in Chinese).

［11］LIZh G，ZHANG S J，ZHOU JZh.Research of evaluationmethod of infrared countermeasure jamming effectbased on image feature［J］.Laser Technology，2013，37(3):413-416(in Chinese).

［12］ECKHORN R，REITBOECK H J，ARNDT M，et al.Feature linking via synchronization among distributed assemblies:simulation of results from cat cortex［J］.Neural Computation，1990，2(3):293-307.

［13］ZHANG B H，LÜ X Q，ZHANG Ch T.Multi-focus image fusion algorithm based on composite incentive model in surfacelet domain ［J］.Opto-Electronic Engineering，2013，40(5):88-99(in Chinese).

［14］WANG Zh B，MA Y D.Medical image fusion using m-PCNN［J］.Information Fusion，2008，9(2):176-185.