結(jié)合稀疏理論與非下采樣剪切波變換的多模態(tài)醫(yī)學(xué)圖像融合算法

，，

(1.浙江理工大學(xué)機械與自動控制學(xué)院，杭州 310018；2.浙江工商大學(xué)信息與電子工程學(xué)院，杭州 310018)

0 引 言

近年來，高精度成像技術(shù)的不斷發(fā)展，給臨床醫(yī)學(xué)的診斷與治療帶來了巨大便利。然而，由于臨床醫(yī)學(xué)領(lǐng)域的特殊性，要求圖像不但能夠呈現(xiàn)直觀的器官輪廓，而且要有清晰的紋理效果，依靠單模態(tài)圖像很難達到這些要求[1-2]。比如，CT圖像能夠反映密度較大的骨骼結(jié)構(gòu)，但對密度較小的組織缺乏清晰成像的能力；相反，MRI圖像只適用于顯示密度較小的組織，而對于大密度的骨骼結(jié)構(gòu)的顯示明顯不適用。因此，為便于臨床醫(yī)學(xué)診斷，需要對單一模態(tài)圖像取長補短，進行多模態(tài)的醫(yī)學(xué)圖像融合。

醫(yī)學(xué)圖像融合算法種類眾多，就分解工具而言，主要有小波變換、輪廓波變換以及剪切波變換等。與剪切波變換相比，小波變換[3]線奇異性較差，對二維圖像信息的捕捉能力不足，容易產(chǎn)生偽吉布斯現(xiàn)象；而輪廓波變換[4]算法相對復(fù)雜，運行效率不高。這兩種變換在計算復(fù)雜度、運行效率、捕捉結(jié)構(gòu)信息能力等方面均遜于剪切波變換[5]。由于傳統(tǒng)的剪切波變換存在下采樣操作，為防止出現(xiàn)偽吉布斯現(xiàn)象，人們對其進行改進，提出了非下采樣剪切波變換[6](Non-subsampled shearlet transform , NSST)。

NSST優(yōu)勢明顯，應(yīng)用簡便，近年來已逐漸成為研究的熱點。如Cao等[7]和Padma等[8]利用NSST進行醫(yī)學(xué)圖像融合，與傳統(tǒng)方式相比，這類方法在融合圖像在紋理、清晰度、亮度以及對比度等方面都獲得了一定的改進，但通過NSST獲得的低頻系數(shù)近似為零項極其有限，即對于圖像的低頻系數(shù)信息無法稀疏的表示，若對其進行直接融合處理，會對圖像重要特征的提取帶來不便，結(jié)果將會造成低頻子帶在融合處理過程中無法最佳的獲取圖像中理想位置的信息。

在圖像融合處理中，研究人員希望盡可能地提高圖像子帶系數(shù)的稀疏度并利用較少的信息獲得最優(yōu)的融合結(jié)果，因此研究者關(guān)注到了稀疏理論。Yu等[9]、Shahdoosti等[10]以及Yin等[11]成功將稀疏理論與字典學(xué)習(xí)思想應(yīng)用到圖像融合處理過程中，有效提高了子帶系數(shù)的稀疏度，但稀疏融合處理是在單尺度基礎(chǔ)上對源圖像進行直接訓(xùn)練，雖然能夠?qū)υ磮D像中的重要特征進行提取，但對于多尺度數(shù)據(jù)的分析能力有限。稀疏理論與多尺度分析的有效結(jié)合，既能滿足提取圖像顯著特征的要求，又可以提高數(shù)據(jù)分析能力，使得融合圖像達到醫(yī)療診斷對清晰度的要求。

基于以上分析，為提高圖像系數(shù)的稀疏度，增強視覺效果，進一步凸顯圖像的細節(jié)信息，本文結(jié)合NSST和稀疏理論對源圖像進行處理。首先，通過NSST濾波器分解源圖像得到高低頻子帶系數(shù)，其中高頻子帶包含源圖像的紋理和細節(jié)等重要的信息，低頻圖像作為源圖像的近似圖像；其次，低頻區(qū)域融合通過稀疏理論的方法來改善其較差的稀疏性，高頻系數(shù)處理采用相對標準差比較法(Relative standard deviation comparison method，RSDC)，通過對相對標準差以及能量值的求取來選擇最終的融合高頻系數(shù)；最后，將低高頻的融合系數(shù)通過NSST逆變換得到最終的融合圖像，并通過空間頻率(Spatial frequency, SF)、標準差(Standard deviation, SD)、平均梯度(Average gradient, AG)、邊緣信息評價因子 (Edge based similarity measure , QABF)、結(jié)構(gòu)相似(Structure similarity, SSIM)等客觀指標對融合后的圖像進行評價。以上幾個客觀評價指標值越大，表明圖像融合效果越好，指標定義見文獻[12]。

1 剪切波及稀疏理論

1.1 剪切波理論

剪切波以合成小波理論為基礎(chǔ)，結(jié)合仿射系統(tǒng)和多尺度分析構(gòu)建模型。當維數(shù)n=2時，其定義式為:

(1)

(2)

NSST操作主要分為兩個步驟，即多尺度剖分以及方向局部化。多尺度剖分主要是源圖像經(jīng)過濾波器以及金字塔組分解成j+1個子帶，其中一個低頻子帶和j個高頻子帶;方向局部化：標準化的剪切波分解器在經(jīng)過從偽極網(wǎng)格系統(tǒng)到笛卡爾系統(tǒng)變換以及傅里葉變換，卷積等相關(guān)運算即可實現(xiàn)非下采樣剪切波的方向局部化。圖1是NSST的3層分解圖，源圖像經(jīng)過第一層分解就得到兩個分支，其中一個分支為低頻子帶，另一分支為高頻子帶，要獲得圖像的多層分解，需要對上一層分解得到的低頻子帶繼續(xù)分解即可。NSST的支撐基是形狀為梯形的一對區(qū)間，如圖2所示，每個梯形區(qū)間寬度為22j，高度為2j。

圖1 NSST分解

圖2 NSST梯形支撐區(qū)間

1.2 稀疏理論

近年來稀疏理論日益完善，在圖像處理領(lǐng)域得到了廣泛的應(yīng)用。稀疏表示主要是將圖像領(lǐng)域明顯的結(jié)構(gòu)信息聚集到字典中較少的原子圖像，而原子圖像的集成可以構(gòu)成一個過完備的字典，由于該字典具有冗余性，信號可以憑借字典中少量原子的線性組合來表示，該表示中原子數(shù)目最少的表示即為稀疏表示，其數(shù)學(xué)模型可以表示為：

(3)

上述稀疏模型中最關(guān)鍵的是字典選擇，而常規(guī)方法主要包括學(xué)習(xí)方式和解析方式兩種，其中學(xué)習(xí)方式主要通過ODL算法、K-SVD算法等算法來訓(xùn)練過完備字典，使其形態(tài)豐富，匹配圖像結(jié)構(gòu)良好。本文通過K-SVD算法來訓(xùn)練字典，其目標優(yōu)化函數(shù)如下

(4)

其中:Y=[y1,y2,…,yN]∈Rn×N為樣本訓(xùn)練矩陣;X∈Rn×N為稀疏系數(shù)矩陣,Xi為X的第i列的稀疏系數(shù)矩陣;T為稀疏度。

由于圖像高頻領(lǐng)域信息較少，冗余較多，為改善圖像融合后的效果，本文在圖像低頻領(lǐng)域進行字典學(xué)習(xí)，式(4)可等價為:

(5)

其中：Ws表示分析變換因子；I為圖像訓(xùn)練矩陣。

通過上述方式構(gòu)建出的學(xué)習(xí)字典能夠?qū)ο∈瓒炔桓叩膱D像低頻領(lǐng)域進行稀疏表示，使得學(xué)習(xí)字典具有NSST的分析能力，同時，對NSST低頻系數(shù)的再次稀疏，可以實現(xiàn)雙重稀疏，進一步提高低頻領(lǐng)域的稀疏度，以更好地提取圖像的相關(guān)特征。

2 NSST與稀疏理論相結(jié)合的融合算法

本文融合策略主要分為三個步驟。首先，利用NSST分解待融合的圖像，得到不同頻率的系數(shù)子帶；其次，根據(jù)不同頻率系數(shù)子帶自身的不同特性設(shè)計相應(yīng)的融合規(guī)則。由于低頻子帶反映源圖像的輪廓，包含了源圖像的大量能量，但其稀疏性不佳，融合后的圖像效果較差，因此為提高其稀疏度，本文引入稀疏理論對低頻區(qū)域進行處理；而圖像的高頻區(qū)域則反映源圖像的細節(jié)信息，因此為凸顯其局部細節(jié)，本文通過相對標準差比較法進行判斷選取。最后，將低高頻區(qū)域各自的融合結(jié)果通過NSST重構(gòu)，獲得最終的融合結(jié)果。具體融合框架圖如圖3所示，源圖像經(jīng)過NSST分解成低高頻子帶，低頻子帶融合通過訓(xùn)練字典，稀疏表示等步驟進行處理，高頻子帶則通過相應(yīng)規(guī)則進行融合，處理后的低高頻子帶系數(shù)經(jīng)過NSST逆變換得到融合圖像。

圖3 融合框架圖

2.2 融合策略

2.2.1 基于稀疏理論的低頻系數(shù)融合算法

圖像的輪廓信息主要集中在低頻子帶系數(shù)中，傳統(tǒng)低頻的融合方案如方差取大法、梯度加權(quán)法、空間頻率加權(quán)法等方法基本能夠較好保留源圖像的輪廓信息，但都無法有效改善低頻系數(shù)較差的稀疏性，導(dǎo)致最終的融合效果不夠理想，因此，為提高其稀疏度，改善融合效果，本文將稀疏理論應(yīng)用到圖像低頻子帶融合處理中。相對于傳統(tǒng)方法，本文稀疏理論的算法不但顧及了源圖像的輪廓，而且使其稀疏度得到了有效提升，此外，該算法也使融合圖像的清晰度和亮度得到了一定的改善，融合的效果較為理想。具體方法為:

a) 利用NSST分解配準后尺寸為M×N的源圖像A和B得到低頻子帶系數(shù)gLA(i,j)和gLB(i,j),利用大小為n×n的滑動窗口分別對兩圖的低頻子帶進行分塊處理,得到數(shù)量為(N+n-1)×(M+n-1)的圖像子塊，將所獲得的圖像子塊轉(zhuǎn)換為列向量,從而得到樣本訓(xùn)練集矩陣VA和VB;

c) 依據(jù)式(4)理論采用K-SVD算法對VAS和VBS進行訓(xùn)練，得到過完備字典D;通過OMP算法預(yù)估兩樣本稀疏矩陣VAS和VBS的稀疏系數(shù)，得到矩陣αA和αB;

(6)

(7)

(8)

(9)

(10)

e) 將式(10)中VF的列向量重新轉(zhuǎn)換為子塊，重構(gòu)該子塊即可得到融合后的低頻子帶系數(shù)。

2.2.2 高頻子帶系數(shù)融合

醫(yī)學(xué)診斷領(lǐng)域為便于對微小細節(jié)的觀察，減少誤判，要求圖像的邊緣紋理盡量清晰，而傳統(tǒng)處理方式如區(qū)域能量法、邊緣強度取大法、區(qū)域?qū)Ρ榷确ǖ却嬖诩毠?jié)處理模糊不清和邊緣帶撕裂等問題，因此，為改善視覺效果，便于醫(yī)學(xué)診斷，本文提出了一種相對標準差比較法。相對于傳統(tǒng)方式，本文高頻融合方法不但有效彌補了這些方法的不足，而且使得高頻區(qū)域的邊緣點更加突顯，邊緣信息的豐富性更優(yōu)，融合圖像的邊緣帶更加明了，對于醫(yī)學(xué)判斷更加有利。具體步驟如下：

相對標準差是衡量圖像融合效果的重要指標，其值越小，表明融合圖像越清晰，圖像的融合效果越好。此處為便于計算，區(qū)域窗口大小設(shè)為3×3 ,相對標準差定義式如下：

(11)

其中:M∈{1,2}。若M值為1，則I(i,j)表示源圖像在(i,j)處的高頻系數(shù)；若M值為2，則I(i,j)表示融合圖像在(i,j)處的高頻系數(shù)。

此處能量值是用于衡量圖像邊緣化程度高低的重要指標，若能量值越大，圖像就越清晰，融合效果也就越好，具體公式為：

E(x,y)=∑(i,j)∈Ω(x,y)P2(i,j)T(i,j)

(12)

其中：Ω(x,y)是以(x,y)為中心的區(qū)域窗口，P(i,j)為區(qū)域差值絕對值系數(shù),T(i,j)則是對應(yīng)的窗口矩陣。

為了突出系數(shù)變化劇烈的程度，增強對比度，提高融合質(zhì)量，此處窗口矩陣T(i,j)設(shè)為

(13)

T中元素數(shù)值相差相對較大，可以處理系數(shù)變化劇烈的區(qū)域。

(14)

α(i,j)=

2.3 圖像融合步驟

本文圖像融合具體過程有以下幾個步驟：

b) 采用K-SVD算法對低頻子帶進行樣本訓(xùn)練。

c) 利用大小為n×n的滑動窗口分別對低頻區(qū)域進行分塊處理，對獲得的圖像子塊序列進行優(yōu)化，并結(jié)合OMP算法進行預(yù)估，得到稀疏系數(shù)矩陣αA和αB;通過稀疏矩陣αA和αB,樣本訓(xùn)練字典矩陣D以及樣本均值矩陣進行處理并將融合矩陣重新轉(zhuǎn)化為圖像塊即可得到融合后的低頻子帶系數(shù)GL(i,j)。

e) 將上述融合結(jié)果進行NSST重構(gòu)逆變換，得到融合結(jié)果圖像F;選擇相關(guān)評價指標對融合后的圖像從定性和定量角度進行分析。

3 圖像結(jié)果分析

為驗證本文算法的有效性，本文選取了灰度圖像與彩色圖像進行實驗仿真，并與文獻[13—17]中的算法進行對比，通過SF指標、SD指標、AG指標、QABF指標以及SSIM指標對融合后的圖像進行評價。

為了使本文算法的有效性更有說服力，本文所有實驗均在windows 10操作系統(tǒng)上的MATLAB 2013a環(huán)境下進行仿真。仿真實驗過程中，滑動窗口大小設(shè)為8, NSST分解層數(shù)均設(shè)為3層。

3.1 灰度圖像仿真結(jié)果與分析

實驗選擇了三組像素尺寸大小為256×256 , 顏色深度為8位灰度的圖像作為源圖像進行仿真，如圖4-圖6所示。其中，圖4(a)-(b)分別為急性腦卒患者腦部的CT圖以及MRI圖；圖5(a)-(b)則是多發(fā)性腦梗塞圖MR-T1/MR-T2 ;圖6(a)-(b)分別為正常腦部的CT圖以及MRI圖。

為了對比分析，本文選擇了5組算法進行對比。其中圖4-圖6中的(c)圖為應(yīng)用李亦等[13]提出的香農(nóng)熵加權(quán)稀疏表示算法進行融合，該算法將圖像子帶分為低頻與低頻、低頻與高頻、高頻三部分進行自適應(yīng)加權(quán)處理；圖4-圖6中的(d)圖是利用陳軼鳴等[14]提出的算法進行處理，低頻系數(shù)融合引入稀疏理論，通過空間頻率激勵脈沖耦合神經(jīng)網(wǎng)絡(luò)，根據(jù)點火次數(shù)來選擇高頻系數(shù)；圖4-圖6中的(e)圖是采用王雷等[15]提出的平移不變剪切波變換的融合方法，該算法對低頻系數(shù)融合采用區(qū)域系數(shù)絕對值和權(quán)重系數(shù)規(guī)則，高頻子帶則采用支持向量值激勵的自生成神經(jīng)網(wǎng)絡(luò)的融合方法；圖4-圖6中的(f)圖是利用殷明等[16]提出的非下采樣四元數(shù)剪切波的融合算法，低頻區(qū)域通過改進稀疏表示進行融合，高頻領(lǐng)域結(jié)合改進空間頻率，邊緣能量和區(qū)域相似匹配度法進行處理；圖4-圖6中的(g)圖是采用嚴春滿等[17]提出的多聚焦圖像融合算法，利用過完備稀疏表示理論對圖像塊進行稀疏分解，再經(jīng)過顯著性系數(shù)的融合，圖像塊的重構(gòu)，圖像塊序列的重排得到最終的融合結(jié)果；圖4-圖6中的(h)圖是利用本文提出的算法得到的融合圖像。

圖4 急性腦卒腦部患者的CT/ MRI融合結(jié)果

圖5 多發(fā)性腦梗塞MR-T1/MR-T2融合結(jié)果

圖6 正常腦部CT/MRI融合結(jié)果

從定性角度分析，觀察圖4-圖6可以發(fā)現(xiàn)，本文算法得到的融合圖像圖4-圖6中的(h)圖在圖像亮度，對比度，紋理邊緣清晰度，視覺效果等方面均明顯優(yōu)于其它幾幅圖像，不僅使源圖像的輪廓與紋理得到了有效保留，而且使圖中的微小細節(jié)得到了呈現(xiàn)，對臨床醫(yī)學(xué)診斷與分析提供了強有力的支持。

觀察圖4不難看出，圖4(c)-(e)以及圖4(g)四幅圖邊緣紋理方面有些許模糊，微小細節(jié)沒有得到有效展示，而圖4(f)的紋理方面清晰度尚可，但中心區(qū)域亮度有些不足；觀察圖5 (c)-(d)兩幅圖亮度尚可，但紋理表現(xiàn)欠佳，而圖5(e)-(g)圖紋理尚可，然而對比度有所欠缺；圖6中幾幅圖源圖像的輪廓保留較好，但圖6(c)-(d)圖亮度不強，對比度不夠理想，圖6(e)-(g)圖組織信息邊緣有些許模糊，視覺效果有些不足。

從定量角度分析，本文選取了5種指標對融合后的圖像進行分析，相關(guān)融合數(shù)據(jù)如表1-表3所示。

表1 圖4中不同融合算法得到的相關(guān)融合指標

表2 圖5中不同融合算法得到的相關(guān)融合指標

表3 圖6中不同融合算法得到的相關(guān)融合指標

觀察表1-表3可知，首先利用陳軼鳴等[13]算法得出的融合圖像的評價指標普遍偏低，各項數(shù)值均不突出，這可能由于該算法精度不高所致，同時也說明此算法在圖像融合領(lǐng)域有較大的局限性；其次，其它幾種算法在進行圖像處理時表現(xiàn)尚可，指標值相差不大，部分指標表現(xiàn)突出，比如表2中(g)圖SD指標、表1中(c)圖的SF指標、表3中(f)圖的SSIM指標等更是達到最優(yōu)，然而，幾種算法中某些指標也略微偏低，如表1-表3中(e)圖的SD指標值均偏低，這也說明此種算法得到的融合圖像紋理模糊，在改善圖像細節(jié)方面稍顯不足；最后，觀察各表中本文算法得到的融合圖像指標值，表1-表3中通過本文算法得到指標值大都最優(yōu)，此外，部分未達到最優(yōu)值的指標與最高指標值相差較小，這充分說明了本文算法具有一定的優(yōu)越性。

3.2 彩色圖像仿真結(jié)果與分析

為了進一步驗證本文算法的有效性，本文選取了兩組彩色圖像SPECT/MRI進行實驗仿真。如圖8(a)-(b)為醫(yī)學(xué)領(lǐng)域轉(zhuǎn)移性肺癌SPECT/MRI圖像，而圖9(a)-(b)為醫(yī)學(xué)上正常腦部的SPECT/MRI圖像。同上述灰度圖像仿真一致，圖8(c)-(h)以及圖9(c)-(h)依次通過文獻[13-17]中算法和本文提出的算法融合獲得；此外，表4和表5分別為圖8和圖9中各算法得到的相關(guān)評價指標。

計算機上處理彩色圖像通常采用RGB顏色系統(tǒng)，因此基于圖像融合考慮，需要實現(xiàn)RGB系統(tǒng)與HIS系統(tǒng)之間的轉(zhuǎn)換。此處在進行彩色圖像融合處理時首先要將其從RGB系統(tǒng)映射到HIS系統(tǒng)；其次，提取其中的I(強度)分量與灰度圖像進行融合；再次，將源彩色圖像中H(色度)分量與S(飽和度)分量與上述融合結(jié)果結(jié)合；最后，進行HIS系統(tǒng)與RGB系統(tǒng)之間的轉(zhuǎn)換即可得到融合后的彩色圖像。融合框架如圖7所示。

圖7 彩色圖像融合框架

圖8 轉(zhuǎn)移性肺癌SPECT/MRI圖像融合效果

圖9 正常腦部SPECT/MRI圖像融合效果

評價指標融合圖像(c)(d)(e)(f)(g)(h)SF9.66158.39079.35039.79199.810010.0031SD56.110751.836454.327656.776956.153357.7371AG0.04330.03970.04160.04350.04200.0429QABF0.48140.41010.45060.44610.46010.4588SSIM0.75390.67110.72030.71970.73500.7558

表5 圖9中不同算法得到的相關(guān)融合指標

結(jié)合圖8、圖9不難發(fā)現(xiàn)，應(yīng)用本文算法得到的融合結(jié)果圖在圖像亮度、紋理清晰度、局部對比度、視覺效果方面均優(yōu)于其它圖像，不僅保留了源圖像的輪廓信息，而且凸顯了重要的細節(jié)特征，對臨床醫(yī)學(xué)的診斷能夠帶來較大的便利。

從定性角度分析，圖8(d)-(e)和圖9 (d)-(e)融合效果圖細節(jié)模糊，紋理不清，源圖像的邊緣帶未盡可能的被保留，丟失了大量細節(jié)信息，不適用于醫(yī)學(xué)診斷，同時也表明該算法在彩色圖像領(lǐng)域應(yīng)用有一定的缺陷；而圖8 (c)和圖8(f)-(g)三幅圖以亮度尚可，但紋理仍不夠清晰，而圖9 (c)和圖9(f)-(g)三幅圖紋理清晰度有所改善，但圖像亮度不佳，對比度不強，視覺效果有所欠缺，對醫(yī)學(xué)診斷與治療仍然不便；觀察本文算法得到的彩色圖像，紋理清晰，細節(jié)明了，視覺效果良好，主觀上有力驗證了本文算法的有效性。

從定量角度分析，觀察表4-表5可以發(fā)現(xiàn)，陳軼鳴等[13]算法得到的融合指標全部偏低，指標值均不突出，融合效果較差；而其它指標值雖然相差不大，但某些指標值如表4中(e)的 SF指標略低，說明其信息豐富性不強，同理，表5中(f)的QABF指標值偏低，也表明其圖像邊緣存在一定的模糊；觀察表4-表5中(h)圖的指標值基本占優(yōu)，即使未取得最優(yōu)的指標值也接近最高，充分說明了本文算法表現(xiàn)良好，有一定的優(yōu)越性，再次驗證了本文算法的有效性。

4 結(jié)束語

圖像融合技術(shù)的不斷發(fā)展，極大地改善了人們的生活方式，也促使醫(yī)學(xué)領(lǐng)域產(chǎn)生了巨大的變革。本文在研究非下采樣剪切波技術(shù)的基礎(chǔ)上提出了一種新的融合算法，即在低頻領(lǐng)域采用引入稀疏理論對其進行融合處理，提高其稀疏度，而對于高頻子帶處理，本文采用相對標準差比較法進行融合，以提高其細節(jié)紋理的清晰度。仿真實驗表明，本文算法得到的融合圖像細節(jié)清晰，亮度適中，圖像信息豐富，同相關(guān)算法比較，無論在主觀視覺上還是客觀指標上均具有一定的優(yōu)越性。