張寶華，劉 鶴，張傳亭

(內(nèi)蒙古科技大學(xué)信息學(xué)院，包頭014010)

引 言

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

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

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

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

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

1 基于BEMD和點(diǎn)火映射圖的特征分類

1.1 BEMD

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

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

BEMD分解的算法步驟如下:

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

(2)如果殘差項R(x，y)單調(diào)或達(dá)到圖像的分解層數(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ù)，進(jìn)入篩選過程。

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

(4)平面插值區(qū)域極值點(diǎn)集，得到圖像的上、下包絡(luò)面 Eup(x，y)和 Edown(x，y)，進(jìn)一步求出圖像Hk(x，y)的均值:

(5)令 Hk(x，y)=Hk-1(x，y)-Emean(x，y)，判斷篩選過程是否滿足停止條件S，如果不滿足則轉(zhuǎn)步驟(3)，其中S為:

(6)判斷Hk(x，y)是否為內(nèi)蘊(yùn)模函數(shù)(intrinsic mode functions，IMF)，其依據(jù)是:若 S ＜ ξ(ξ為閾值，本文中 ξ=0.2)，則 Dj(x，y)=Hk(x，y)，否則令 k=k+1，轉(zhuǎn)到步驟(2)。

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

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

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

1.2 復(fù)合型脈沖耦合神經(jīng)網(wǎng)絡(luò)

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

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

Fig.1 CT image

Fig.2 MRI image

Fig.3 Firingmap

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

Fig.4 Dual-channel PCNN model

1.3 基于點(diǎn)火映射圖的特征提取

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

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

Fig.5 Firingmaps of BIMF and R

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

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

2.1 融合步驟

本文中提出的基于特征提取和復(fù)合型脈沖耦合神經(jīng)網(wǎng)絡(luò)的醫(yī)學(xué)圖像融合方法如下:

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

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

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

Fig.6 Algorithm flowchart

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

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

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

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

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

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

2.2 融合系數(shù)選擇

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

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

3 實驗結(jié)果及分析

為了比較不同融合方法的融合效果，選擇兩組反映不同醫(yī)學(xué)信息的源圖像進(jìn)行融合實驗，融合前已利用基于光流場的方法對源圖像進(jìn)行配準(zhǔn)，配準(zhǔn)誤差小于1個像素，將融合圖像與其它5種傳統(tǒng)融合方法進(jìn)行比較:基于Laplace、離散小波變換(discrete wavelet transform，DWT)、主成分分析(principal component analysis，PCA)、gradient pyramid 和 contrast pyramid的融合方法，Laplace方法進(jìn)行3層分解;DWT采用Daubechies小波基進(jìn)行3層小波分解，高頻采用取大低頻取平均的方法。第1組實驗圖像為圖1和圖2，對其采用6種融合方法的比較結(jié)果見圖7;第2組實驗圖像為圖8和圖9，6種融合方法的比較結(jié)果見圖10，可以明顯看到，6種融合方法都實現(xiàn)了源圖像的信息融合，豐富了圖像信息，但基于DWT和contrast pyramid方法的融合圖像部分區(qū)域受偽影干擾，邊緣不清晰，基于PCA方法的融合圖像局部失真明顯。由本文中方法得到的融合圖像更加接近于源圖像，細(xì)節(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個指標(biāo)對融合圖像進(jìn)行客觀評價，分別是互信息(mutual information，MI)、邊緣強(qiáng)度 QA，B，F(xiàn)(A，B 代表源圖像，F(xiàn) 代表融合圖像)、熵和平均梯度。互信息衡量融合圖像中包含源圖像的信息的程度，QA，B，F(xiàn)，表達(dá)的是融合圖像含有源圖像邊緣信息的豐富程度，信息熵可以衡量圖像攜帶的信息量，值越大，說明融合圖像包含信息量越豐富;平均梯度反映圖像中的特征細(xì)微變化，平均梯度大的圖像清晰程度較好，以上指標(biāo)越大說明融合效果越好。表1、表2中分別顯示兩組實驗融合圖像的客觀評價指標(biāo)比較結(jié)果，可以看到，本文中的方法評價指標(biāo)均優(yōu)于其它方法，通過主觀評價和客觀指標(biāo)的比較，說明運(yùn)用本文中方法得到的融合圖像含有更多源圖像信息，圖像紋理較豐富，細(xì)節(jié)突出，融合效果更好。

Table 1 The first group of image fusion evaluation objective comparison

對表1中的實驗結(jié)果進(jìn)行比較可以得出，本文中的方法與經(jīng)典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組實驗的客觀指標(biāo)也進(jìn)行了比較，可以得到和第1組相同的結(jié)論，即本文中的方法在客觀評價中有較好的效果。

Table 2 The second group of image fusion evaluation objective comparison

4 結(jié)論

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