基于經驗模態(tài)分解提取紋理的圖像融合算法

張寶華，劉 鶴，張傳亭

(內蒙古科技大學信息學院，包頭014010)

引 言

隨著成像技術的突飛猛進，各類精密成像設備推動了醫(yī)學影像學的發(fā)展，為臨床提供豐富的人體醫(yī)學影像，但成像設備種類眾多，其成像機理不同，反映醫(yī)學信息各有側重。為了全面分析醫(yī)學圖像中包含的解剖信息和功能信息，需要對多模態(tài)醫(yī)學圖像進行融合。

醫(yī)學圖像融合技術面向多模態(tài)醫(yī)學圖像，把各種醫(yī)學圖像的信息有機地結合，完成各類醫(yī)學信息融合，不僅有效利用已有醫(yī)學影像，而且還有助于發(fā)掘潛在醫(yī)學信息，輔助醫(yī)院診療。

歐美國家在20世紀80年代已開始在圖像融合領域取得豐碩的成果，處于領先地位;MITIANOUDIS等人［1］通過將稀疏表達應用于系數(shù)選擇中，拓展了圖像融合的視野和手段;國內YAN，LI等人［2-3］提出了基于脈沖耦合神經網絡(pulse coupled neural net-work，PCNN)的圖像融合方法，應用于不同類型的圖像，將神經網絡作為系數(shù)優(yōu)選的方法，豐富了融合系數(shù)選擇方法，取得了較好的融合效果。

目前圖像融合方法包括基于空域變換和頻域變換的方法，其中以小波變換和各類超小波變換為代表的基于多分辨分析的圖像融合方法應用最廣。但小波變換和其改進方法依賴于預先定義的濾波器或基函數(shù)，小波變換會有下采樣操作，變換后圖像會引進偽吉布斯現(xiàn)象，降低融合圖像質量。HUANG等人為達到對非線性和非穩(wěn)態(tài)的數(shù)據(jù)進行自適應和多尺度分析，提出了經驗模態(tài)分解(empirical mode decomposition，EMD)［4］。經驗模式分解作為一種新的多尺度圖像分解方法［5-6］，具有比小波分析［7］更直觀的特征表示方式和更靈活的頻率特性，避免了分解中引進冗余信息，同時EMD對于圖像細節(jié)保護和圖像紋理的提取等方面具有優(yōu)勢［8-11］，適合于對安全性要求較高的醫(yī)學圖像進行多分辨分析。

本文中將2維經驗模態(tài)分解(bidimensional empirical mode decomposition，BEMD)運用到醫(yī)學圖像特征提取中，通過將BEMD分解后的子圖像和趨勢圖輸入神經網絡獲取它們的點火映射圖，提取不同分解層對應醫(yī)學圖像特征。之后將對應于圖像紋理信息和背景信息的系數(shù)分別通過PCNN和雙通道PCNN(daul channel pulse coupled neural network，DCPCNN)選取融合系數(shù)。由于區(qū)別對待代表圖像紋理和背景信息的像素，既保護了圖像中的特征，又有效改善了PCNN在醫(yī)學圖像系數(shù)選擇中的效果。

1 基于BEMD和點火映射圖的特征分類

1.1 BEMD

EMD分解具有優(yōu)越的空間和頻率特性。其易拓展性也推動了BEMD方法的出現(xiàn)，BEMD同樣具有數(shù)據(jù)驅動和良好的自適應性等特點，而且具有多尺度特性。將BEMD用于醫(yī)學圖像處理，可得到源圖像在不同頻率的2維內蘊模函數(shù)(bidimensional intrinsic mode functions，BIMF)和殘差項 R。

BEMD是數(shù)據(jù)驅動的，在BEMD分解過程中，當篩分終止條件一定時，圖像分解得到的BIMF個數(shù)與圖像數(shù)據(jù)本身的特征相關［12］，在紋理分析等領域中得到廣泛的應用。

BEMD分解的算法步驟如下:

(1)初始化，設源圖像為 I(x，y)(其中 x=1，2，…，P;y=1，2，…，Q)，殘差項為 Rj(x，y)，j為分解層數(shù)，當 j=0 時，R(x，y)=I(x，y)。

(2)如果殘差項R(x，y)單調或達到圖像的分解層數(shù)，則算法停止;否則，令:Hk-1(x，y)=Rj-1(x，y)。若 k=1，即 H0(x，y)=Rj-1(x，y)，Hk(x，y)為第k級分解系數(shù)去掉Rj(x，y)的殘余系數(shù)，k為分解層數(shù)，進入篩選過程。

(3)求圖像I(x，y)極值點:將解空間點集分為區(qū)域極大值點集和極小值點集。

(4)平面插值區(qū)域極值點集，得到圖像的上、下包絡面 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)是否為內蘊模函數(shù)(intrinsic mode functions，IMF)，其依據(jù)是:若 S ＜ ξ(ξ為閾值，本文中 ξ=0.2)，則 Dj(x，y)=Hk(x，y)，否則令 k=k+1，轉到步驟(2)。

(7)求殘差項Rj=Rj-1-Dj。若Rj中間部分極值點數(shù)大于2或分解得到的IMF數(shù)目小于指定值，將 Rj轉到步驟(2)，j=j+1。

(8)得到2維分解表達式為:

以上步驟中，Dj是第j個2維內蘊模函數(shù)，Rj是經過j層分解后的趨勢圖像。

1.2 復合型脈沖耦合神經網絡

PCNN是一種新型神經網絡，是依據(jù)動物的大腦視覺皮層上同步脈沖振蕩現(xiàn)象提出的［12］。PCNN模型中，相鄰的一群神經元可通過發(fā)放同步脈沖傳送激勵。當一個或數(shù)個神經元點火，輸出的同步脈沖將被傳送到相鄰的神經元，使之迅速點火，神經元激發(fā)的過程稱為點火，將一幅圖像的像素映射到神經網絡中，就獲得了圖像的點火映射圖。

通過PCNN神經網絡來尋找各個像素間的聯(lián)系，將存在聯(lián)系的像素進行歸類，進而確定這一類像素的特征。根據(jù)一副圖像在神經網絡中各個像素被激發(fā)點火次數(shù)的不同，進而將同一點火次數(shù)的像素分為一類。這就找到了圖像中像素間的聯(lián)系。圖1、圖2分別為計算機斷層掃描成像圖像(computed tomography，CT)和磁共振成像圖像(magnetic resonance imaging，MRI)頭部源圖像。PCNN根據(jù)點火次數(shù)將像素按灰度值的大小進行分類如圖3所示，圖3a～圖3c分別為圖1的1～3次點火映射圖，通過比較可以看到相同次數(shù)的點火點組成的圖像完整表達了圖像輪廓信息，點火次數(shù)越多對應的細節(jié)越清晰。

Fig.1 CT image

Fig.2 MRI image

Fig.3 Firingmap

由于PCNN對圖像中偏暗區(qū)域特征的篩選效果較差，而DCPCNN見圖4，圖中SA和SB代表輸入激勵;FA和FB表示反饋通道輸入;Lij表示連接輸入項;Wij為突觸聯(lián)接權;Vl和Vt為歸一化常數(shù);Uij表示神經元的內部活動項;βA和βB表示連接強度;Yij表示神經元的脈沖輸出，它的值為0或者1;Tij是動態(tài)閾值;αl和αt為調節(jié)對應式子的常量，maximum是比較器。對圖像中偏暗或偏亮區(qū)域特征的處理效果較好［13-14］，將圖像根據(jù)特征分布情況進行分類后分別通過PCNN和DCPCNN可以較好地提取圖像中各個范圍的特征，本文中將PCNN和DCPCNN的復合結構定義為復合型脈沖耦合神經網絡。

Fig.4 Dual-channel PCNN model

1.3 基于點火映射圖的特征提取

圖像中的像素可類比成一個神經元，整個圖像就是一個復雜的神經網絡系統(tǒng)，將圖像的灰度值輸入網絡對應的神經元，當一個像素被激發(fā)時，與它相鄰的像素點就會受到影響，被激發(fā)或者保持熄滅狀態(tài)。通過多次點火激發(fā)就可以根據(jù)各個像素點火次數(shù)，來確定圖像中每個點與附近的點或者偏離它較遠的點有沒有相互間的脈沖聯(lián)系。依據(jù)一副圖像中各個像素點的點火次數(shù)將一副圖像的像素進行分類，點火次數(shù)相同的像素歸為一類，從而將圖像依據(jù)神經網絡中各神經元激發(fā)過程中聯(lián)系強弱的特征映射到圖像中對應的像素，得到了各個像素間的相互聯(lián)系。找到不同分類的像素的灰度極值就可以確定各類像素灰度分布范圍。

BEMD的分解過程實現(xiàn)了圖像從低頻到高頻多尺度分離過程。首先分解出來的固有模態(tài)分量BIMF是圖像所含有的最高頻率分量，該分量的各處頻率都對應著圖像在各處的局部最高頻率。各個模態(tài)分別從不同頻率反映了一幅圖像的特征，再通過PCNN得到圖像各頻段的點火映射圖，根據(jù)BEMD逆變換獲得一副圖像總的點火映射圖。圖5a～圖5h表示圖像經過BEMD分解得到的分解圖像及其點火映射圖，圖5a～圖5d是圖1對應的BIMF和殘差項，圖5e～圖5h分別是圖5a～圖5d對應的點火映射圖，由于表示圖像特征信息的像素點在點火時會受到多次激發(fā)，因而通過將不同點火次數(shù)的像素分類即可提取圖像的特征信息。

Fig.5 Firingmaps of BIMF and R

2 基于特征分類的圖像融合方法

利用CT和MRI多模頭部醫(yī)學圖像進行融合，通過利用BEMD特征提取將圖像區(qū)域分為紋理和背景兩部分，將兩區(qū)域分別建立融合規(guī)則選擇融合系數(shù)，由于輪廓、紋理等信息被較好保護，本方法具備了BEMD和PCNN兩者的優(yōu)勢，改善了融合圖像質量。

2.1 融合步驟

本文中提出的基于特征提取和復合型脈沖耦合神經網絡的醫(yī)學圖像融合方法如下:

(1)將待融合的2幅圖像分別通過BEMD得到BIMF和殘差項。

(2)BIMF和殘差項分別輸入到PCNN中，分別得到各自點火映射圖，并將各自的點火映射圖相加得到總的點火映射圖。

(3)將原始圖像中點火次數(shù)相同的像素點歸為一類，若最大點火次數(shù)為N，則像素點共分為N類(N為自然數(shù))，由于點火次數(shù)高的像素一般對應于圖像的紋理，本文中將點火次數(shù)為n-N的分類定義為紋理類，剩余類定義為背景類，其中:

Fig.6 Algorithm flowchart

式中，k屬于正整數(shù)。

(4)根據(jù)各分類集合的像素極值確定各個分類灰度分布范圍(極小值～極大值)。

(5)將兩幅圖像中紋理類的灰度分布范圍求交集，即若源圖像A和B的最大點火次數(shù)為6，其點火次數(shù)為4～6的分類為紋理類，若集合中像素點灰度分布范圍分別為(25～190)和(60～210)，則取兩者交集(60～190)范圍對應的像素，通過PCNN進行選擇融合。

(6)其余像素對應圖像的偏暗和偏亮部分，通過DCPCNN進行融合。

(7)通過PCNN和DCPCNN得到的融合系數(shù)重構得融合圖像。

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

2.2 融合系數(shù)選擇

2.2.1 紋理區(qū)域融合系數(shù)選擇 通過PCNN選擇紋理區(qū)域的融合系數(shù)，能較好提取圖像中紋理信息。根據(jù)像素點映射的神經元激發(fā)所產生的的點火次數(shù)的大小作為像素點優(yōu)選的指標，來選擇兩幅圖像中對應位置的融合系數(shù)。

2.2.2 背景區(qū)域融合系數(shù)選擇 采用DCPCNN選擇背景區(qū)域融合系數(shù)，DCPCNN可改善PCNN在醫(yī)學圖像中偏暗區(qū)域特征選擇的效果，與傳統(tǒng)的單通道PCNN相比，DCPCNN是由兩個簡化PCNN并行組成，首先通過計算以像素點O(i，j)為中心位置的3×3鄰域中任意3個點的和與其它任意3個點的差值，得到其中最小值和最大值，將最大值和最小值的差 W，經過 1/［1+e-γ×W(i-1，j-1)］(γ 是調節(jié)量)運算得到O(i-1，j-1)的β值。通過選擇兩個通道中神經元的內部活動項Uij來控制像素點的點火狀態(tài)。從而根據(jù)Uij選擇兩幅圖中像素點Uij大的作為融合圖像的像素點。

3 實驗結果及分析

為了比較不同融合方法的融合效果，選擇兩組反映不同醫(yī)學信息的源圖像進行融合實驗，融合前已利用基于光流場的方法對源圖像進行配準，配準誤差小于1個像素，將融合圖像與其它5種傳統(tǒng)融合方法進行比較:基于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種融合方法都實現(xiàn)了源圖像的信息融合，豐富了圖像信息，但基于DWT和contrast pyramid方法的融合圖像部分區(qū)域受偽影干擾，邊緣不清晰，基于PCA方法的融合圖像局部失真明顯。由本文中方法得到的融合圖像更加接近于源圖像，細節(jié)清晰完整，視覺效果更好。

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(xiàn)(A，B 代表源圖像，F(xiàn) 代表融合圖像)、熵和平均梯度?；バ畔⒑饬咳诤蠄D像中包含源圖像的信息的程度，QA，B，F(xiàn)，表達的是融合圖像含有源圖像邊緣信息的豐富程度，信息熵可以衡量圖像攜帶的信息量，值越大，說明融合圖像包含信息量越豐富;平均梯度反映圖像中的特征細微變化，平均梯度大的圖像清晰程度較好，以上指標越大說明融合效果越好。表1、表2中分別顯示兩組實驗融合圖像的客觀評價指標比較結果，可以看到，本文中的方法評價指標均優(yōu)于其它方法，通過主觀評價和客觀指標的比較，說明運用本文中方法得到的融合圖像含有更多源圖像信息，圖像紋理較豐富，細節(jié)突出，融合效果更好。

Table 1 The first group of image fusion evaluation objective comparison

對表1中的實驗結果進行比較可以得出，本文中的方法與經典Laplace、DWT相比，互信息分別提高了 2.7 倍和 2.3 倍，QA，B，F(xiàn)分別提高了 14% 和30.8%，信息熵分別提高了14.9%和9.5%，平均梯度分別提高了0.7%和5%。與其它方法相比，本文中方法也有明顯優(yōu)勢，表2和圖10對第2組實驗的客觀指標也進行了比較，可以得到和第1組相同的結論，即本文中的方法在客觀評價中有較好的效果。

Table 2 The second group of image fusion evaluation objective comparison

4 結論

利用BEMD分析方法將醫(yī)學圖像分解為不同頻率子圖像，通過提取特征建立基于復合型脈沖耦合神經網絡的融合規(guī)則，由于克服了PCNN在偏暗圖像中細節(jié)無法有效提取的問題，本算法能夠從醫(yī)學圖像中提取更多的紋理信息。從主觀和客觀方面與其它圖像融合方法進行了實驗結果比較，證明了該方法在保留邊緣、紋理、細節(jié)信息的有效性。

［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.