張寶華，劉 鶴，張傳亭

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

引 言

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

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

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

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

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

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

1.1 BEMD

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

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

BEMD分解的算法步驟如下:

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

(2)如果殘差項(xiàng)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級(jí)分解系數(shù)去掉Rj(x，y)的殘余系數(shù)，k為分解層數(shù)，進(jìn)入篩選過(guò)程。

(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)，判斷篩選過(guò)程是否滿足停止條件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)求殘差項(xiàng)Rj=Rj-1-Dj。若Rj中間部分極值點(diǎn)數(shù)大于2或分解得到的IMF數(shù)目小于指定值，將 Rj轉(zhuǎn)到步驟(2)，j=j+1。

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

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

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

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

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

Fig.1 CT image

Fig.2 MRI image

Fig.3 Firingmap

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

Fig.4 Dual-channel PCNN model

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

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

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

Fig.5 Firingmaps of BIMF and R

2 基于特征分類(lèi)的圖像融合方法

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

2.1 融合步驟

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

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

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

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

Fig.6 Algorithm flowchart

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

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

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

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

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

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

2.2 融合系數(shù)選擇

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

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

3 實(shí)驗(yàn)結(jié)果及分析

為了比較不同融合方法的融合效果，選擇兩組反映不同醫(yī)學(xué)信息的源圖像進(jìn)行融合實(shí)驗(yàn)，融合前已利用基于光流場(chǎng)的方法對(duì)源圖像進(jìn)行配準(zhǔn)，配準(zhǔn)誤差小于1個(gè)像素，將融合圖像與其它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組實(shí)驗(yàn)圖像為圖1和圖2，對(duì)其采用6種融合方法的比較結(jié)果見(jiàn)圖7;第2組實(shí)驗(yàn)圖像為圖8和圖9，6種融合方法的比較結(jié)果見(jiàn)圖10，可以明顯看到，6種融合方法都實(shí)現(xiàn)了源圖像的信息融合，豐富了圖像信息，但基于DWT和contrast pyramid方法的融合圖像部分區(qū)域受偽影干擾，邊緣不清晰，基于PCA方法的融合圖像局部失真明顯。由本文中方法得到的融合圖像更加接近于源圖像，細(xì)節(jié)清晰完整，視覺(jué)效果更好。

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

Table 1 The first group of image fusion evaluation objective comparison

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

Table 2 The second group of image fusion evaluation objective comparison

4 結(jié)論

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