基于引導濾波與稀疏表示的醫(yī)學圖像融合

王兆濱，馬一鯤，崔子婧

(蘭州大學信息科學與工程學院 蘭州 730000)

在醫(yī)學影像學中，常用的成像技術(shù)主要有計算機斷層掃描(computed tomography, CT)、磁共振成像(magnetic resonance imaging, MRI)、正電子發(fā)射計算機斷層掃描(positron emission computed tomography, PET)、單光子發(fā)射計算機斷層掃描(single photon emission computed tomography,SPECT)等。其中，通過CT 和MRI 獲得的圖像可以提供解剖信息，具有較高的空間分辨率；前者可區(qū)分骨骼、血管等有密度差異的組織，后者可區(qū)分軟組織但卻無法提供骨骼信息。通過SPECT 和PET獲得的圖像能夠反映人體的功能結(jié)構(gòu)和代謝信息，但其空間分辨率一般較低。因此，若能將人體同一部位的不同醫(yī)學圖像融合，就可用一幅圖像表示多種結(jié)構(gòu)和功能信息[1]。醫(yī)生就能根據(jù)融合圖像做出更精確的診斷。基于以上原因，多模態(tài)醫(yī)學圖像融合已經(jīng)成為圖像處理領(lǐng)域的研究熱點之一。

多模態(tài)醫(yī)學圖像融合是指將不同醫(yī)學成像技術(shù)采集到的同一部位圖像進行融合的過程[2]。獲得的融合圖像與單一成像技術(shù)提供的圖像相比，其信息更加全面。然而，盡管不同成像技術(shù)的圖像具有顯著差異，這些圖像間仍存有大量冗余信息。因此，在有效提取互補信息的同時盡可能剔除冗余信息是獲得高質(zhì)量融合圖像的關(guān)鍵。

根據(jù)融合方法劃分，圖像融合算法可分為空間域和變換域算法。空間域算法根據(jù)融合規(guī)則直接處理圖像每個像素，是邏輯較為簡單的算法。變換域算法則是通過數(shù)學變換方法將圖像轉(zhuǎn)換為頻域表示方式并得到相應的頻域系數(shù)，隨后采用合適的融合規(guī)則對系數(shù)進行融合，最后用相應逆變換得到融合圖像。對于變換域算法來說，變換方法和融合規(guī)則的選擇至關(guān)重要。

目前，已有多種融合算法應用于多模態(tài)醫(yī)學圖像融合領(lǐng)域，如稀疏表示、多尺度變換、神經(jīng)網(wǎng)絡等。然而，這些算法在圖像質(zhì)量和細節(jié)保存方面仍有較多缺陷。因此，本文結(jié)合引導濾波器與自適應稀疏表示，提出了一種新的多模態(tài)醫(yī)學圖像融合算法，用于提升融合圖像的清晰度。

1 自適應稀疏表示模型

由于不同圖像塊的結(jié)構(gòu)差異較大，傳統(tǒng)稀疏表示算法通常需要學習一個高度冗余的字典來滿足圖像重建的要求。然而，這種高冗余字典不僅會使生成的圖像產(chǎn)生視覺偽影，而且還會花費較長的運算時間。與其相比，自適應稀疏表示模型不僅構(gòu)建了一個冗余字典，而且還構(gòu)建了一組緊湊的子字典并在融合過程中自適應地選擇這些子字典進行運算。該模型由文獻[3]首次提出并應用于圖像融合和去噪。此外，自適應聯(lián)合稀疏模型(adaptive joint sparse model, JSM)[4]也被提出，其同樣是基于自適應稀疏表示的概念。隨后，文獻[5]根據(jù)圖像的相似結(jié)構(gòu)將圖像分解為相似模型、平滑模型和細節(jié)模型。其中，對細節(jié)模型使用自適應稀疏表示進行融合，提高了融合效率。文獻[6]將自適應稀疏表示與修正后的空間頻率結(jié)合用于多模態(tài)醫(yī)學圖像融合。

隨后，自適應稀疏表示模型廣泛應用于多尺度變換融合算法中。文獻[7]為更好地提取圖像邊緣和方向信息，利用多個濾波器對輸入圖像進行分解，采用自適應稀疏表示模型對子帶進行融合。文獻[8] 使用非下采樣剪切波變換(non-subsampled shear wave transform, NSST)分解圖像，利用自適應稀疏表示模型融合高頻子帶，得到了較好的融合結(jié)果。文獻[9]也使用了NSST 對輸入圖像進行分解，但自適應稀疏表示模型被用于低頻子帶的融合；而高頻子帶的空間頻率被用于神經(jīng)元的反饋輸入，并且激勵一個脈沖耦合神經(jīng)網(wǎng)絡(pulse coupled neural network, PCNN)進行融合，也得到了較好的結(jié)果。文獻[10] 則利用自適應的卷積稀疏表示模型融合遙感圖像。文獻[11] 使用拉普拉斯金字塔分解輸入圖像，使用自適應稀疏表示融合了每一層子帶，通過實驗證明了該方法應用于醫(yī)學圖像融合的有效性。

在自適應字典的構(gòu)建中，使用的訓練集是由從多幅高質(zhì)量圖像中隨機采樣出的大量圖像塊構(gòu)成的。同時，為保證訓練集中的圖像塊都具有豐富的邊緣結(jié)構(gòu)信息，訓練集需要刪除強度方差小于給定閾值的圖像塊[12-14]。因此假設最終的訓練集為P={p1,p2,···,pM}，即有M個滿足以上條件的圖像塊。然后，對集合P中的圖像塊根據(jù)其梯度主導方向進行分類，構(gòu)建出一組子詞典訓練集。子字典訓練集的構(gòu)建采用方向分配方法[15]，通過計算每個圖像塊的梯度方向直方圖來決定主導方向。如對于圖像塊pm(i,j)∈P，其水平梯度Gx(i,j)和垂直梯度Gy(i,j)通過Sobel 算子計算得出：

隨后，根據(jù)pm中所有像素的梯度幅值和方向來繪制pm的方向直方圖。假設將360°平均劃分為K份，用于構(gòu)建K個子字典，如圖1a 所示為K=4。其次，根據(jù)pm中每個像素的梯度方向?qū)⑺邢袼貏澐譃镵個方向。同時，將屬于同一方向的所有像素的梯度幅值相加作為該方向在直方圖中的統(tǒng)計量，由此得到pm的方向直方圖。最后，選擇直方圖中統(tǒng)計量最大的方向作為pm的主導方向。由此，所有圖像塊根據(jù)其主導方向劃分為K個圖像塊子集{Pk|k=1,2,···,K}∈P。每個子集的圖像塊都有相同的主導方向因而都具有相似的結(jié)構(gòu)。根據(jù)這些子集訓練出的子字典將會更緊湊。

圖1 自適應稀疏表示模型子詞典圖解

此外，模型最終需要構(gòu)建出K+1 個字典D={D0,D1,···,DK}，其中字典D0如圖1b 所示，該字典是由集合P中所有圖像塊訓練得到的，其余子字典{Dk|k=1,2,···,K}則是通過其相應子集學習得出的，如圖1c～圖1f 所示。

在實際融合過程中，將待融合的輸入圖像切分成與pm大小一致的圖像塊，并自適應地選擇合適子字典進行融合。子字典的選擇過程如圖2 所示。設km表示為切分后的某一輸入圖像塊應被分配到第km個子集，km按以下公式得到：

圖2 自適應詞典的選擇過程

式中，θmax=max{θ1,θ2,···,θK}， θK表示該輸入圖像塊的方向直方圖中第K類的統(tǒng)計量，θmax為K個統(tǒng)計量的峰值；k*=argmax{θk|k=1,2,···,K}是θmax的指標，若km=0，則表示輸入的是不規(guī)則圖像塊，將采用字典D0進行融合。

本文應用K-SVD(K-means singular value decomposition)算法對字典進行訓練，訓練集是通過對30 組多模態(tài)醫(yī)學圖像進行采樣得到的10 萬個大小為8×8 的圖像塊。通過實驗確定了參數(shù)K的最佳值為6，最終本文訓練出7 個字典，每張字典的大小均為64×256。

2 引導濾波模型

引導濾波器可較好地保持邊緣信息，同時具有較低的復雜度。由于其良好的保邊平滑性，大量研究將其應用于圖像處理領(lǐng)域。文獻[16] 為在融合中挖掘出更多有效信息，提出基于目標提取和引導濾波增強的圖像融合算法。文獻[17] 引入感知因子對引導濾波器進行改進，對原融合權(quán)重進行修正，明顯改善了融合圖像的邊緣特性。此外，為更好保存有效信息以及提高空間連續(xù)性，文獻[18]通過結(jié)合復剪切波變換(complex shearlet transform,CST)與引導濾波器的優(yōu)點提出了新的圖像融合方案。文獻[19]使用改進的PCNN 模型與引導濾波結(jié)合，更好地匹配了人眼視覺感知。文獻[20]將引導濾波器與小波變換結(jié)合，克服了引導濾波器對噪聲圖像處理的缺陷，并提出了新穎的權(quán)值映射融合規(guī)則。文獻[21] 用引導濾波器和顯著性映射分別計算輸入圖像分解后的權(quán)重圖，得到較好的融合結(jié)果。

此外，引導濾波器也可以看作是一種圖像的多尺度變換法，并且可以結(jié)合圖像的像素值和空間信息來構(gòu)造權(quán)重圖。文獻[22] 基于這種理念結(jié)合隨機游走模型用于融合多聚焦圖像，其客觀評價結(jié)果表明該方法可以得出很好的融合效果。因此，本文算法利用引導濾波器對輸入圖像進行分解，從而得到細節(jié)層和基礎層。下面將介紹本文使用的引導濾波器模型。

引導濾波器可看作是一種具有局部窗口的線性濾波器，定義為：

式中，G為 引導圖像，用于指導輸入圖像I進行濾波，得到輸出結(jié)果O；Wij(G)為I在(i,j)處的像素點權(quán)重值，其根據(jù)G在(i,j)處的局部窗口wk中的所有像素計算得出。

Wij(G)的計算過程如下：

式中， μk和 σk分別為局部窗口wk中的像素平均值和方差；Gi與Gj為兩個相鄰的像素點的像素值； ε為懲罰值。通過上式可得出，當像素點Gi與Gj在邊界時，(Gi-μk)(Gj-μk)為異號，否則為同號。異號時，其權(quán)重較低，遠小于同號時的權(quán)重，因此圖像平坦區(qū)域的像素將會被添加一個較大的權(quán)重值，從而削弱了濾波器平滑模糊的效果，起到保持邊緣細節(jié)的作用。此外，懲罰值ε同樣顯著影響該濾波器的濾波效果。當ε較小時，其濾波效果與上述相同，起到保持邊緣細節(jié)的作用；而當ε較大時，引導濾波器就可以看作是一個均值濾波器，其平滑效果會更加顯著，但卻起不到保邊平滑的效果。

此外，也可以從線性濾波器的角度來理解引導濾波器：

該式表明引導圖像的梯度變化與結(jié)果圖的梯度變化呈線性關(guān)系。因此，一般來說，引導圖像的作用就是在提示輸出圖像不同輸入?yún)^(qū)域的位置和梯度信息，以便更好地區(qū)分和保存輸入圖像的平滑和邊緣區(qū)域。

3 基于引導濾波與自適應稀疏表示的多模態(tài)醫(yī)學圖像融合

本文結(jié)合自適應稀疏表示模型和引導濾波器，提出一種新的多模態(tài)醫(yī)學圖像融合算法，該算法的流程圖如圖3 所示。

圖3 本文算法流程

1) 圖像分解。使用高斯濾波器將輸入圖像A和B分解為細節(jié)層圖像和基礎層圖像，其過程如式(12)～式(14)所示：

式中，G表 示高斯函數(shù)，使用大小為3×3的高斯核。

通過高斯濾波器得到基礎層B1和B2后，將輸入圖像與其基礎層相減得到細節(jié)層D1和D2：

2) 結(jié)合拉普拉斯濾波器和顯著性區(qū)域構(gòu)建基礎層的權(quán)值圖。首先使用大小為3×3的拉普拉斯濾波器L 作用于輸入圖像A和B，獲得其高通圖像H1和H2：

對得到的高斯圖像H1和H2的絕對值進行顯著性區(qū)域計算，構(gòu)建顯著性圖像S1和S2，并得到其相對應的初始權(quán)值圖P1和P2：

式中，Y是大小為(2τg+1)(2σg+1)的高斯低通濾波器，本實驗中 τg和 σg的值設置為5。

若僅使用P1和P2對圖像進行融合會導致融合結(jié)果的拼接感過強，丟失大量邊緣信息，從而使得融合圖像會不符合人類的視覺系統(tǒng)。因此，將P1和P2作為輸入圖像A和B在引導濾波器中的引導圖就可以計算出基礎層的最終權(quán)值圖：

4) 融合細節(jié)層。利用自適應稀疏表示模型對細節(jié)層D1和D2進行融合，從而構(gòu)建出融合細節(jié)層FD。由于醫(yī)學圖像的大小均為256×256，因此不需要圖像配準過程。通過計算細節(jié)層圖像塊的方向直方圖，然后根據(jù)式(5)自適應地選擇不同字典對細節(jié)層D1和D2的圖像塊進行稀疏表示：

式中，Dk1和Dk2為自適應挑選的字典； α1和 α2分別為細節(jié)層D1和D2的稀疏系數(shù)。隨后，對稀疏系數(shù)進行融合。在本文算法中選擇依據(jù)稀疏系數(shù)的方差進行融合。因此，設稀疏系數(shù) α1和 α2的均值分別為μ1和 μ2，則 α1和 α2的方差為：

根據(jù)稀疏系數(shù)的方差大小和稀疏矩陣的稀疏度d1=‖α1‖0、d2=‖α2‖0來指定稀疏系數(shù)的融合規(guī)則：

最后，根據(jù)融合后的稀疏系數(shù)進行細節(jié)層重構(gòu)，得到融合細節(jié)層：

5) 計算融合結(jié)果圖。將融合后的基礎層圖像與融合細節(jié)層圖像相加，得到融合圖像F：

4 實驗結(jié)果與分析

為驗證提出方法的有效性，選擇3 組多模態(tài)腦部病變醫(yī)學圖像進行對比實驗。這些圖片由網(wǎng)站(http://www.med.harvard.edu/aanlib/home.html)獲取。實驗使用的3 組CT 和MRI 圖像為：肉瘤患者的第17 組腦切片、腦膜瘤患者的第17 組腦切片和急性中風言語停止患者的第18 組腦切片。所有的CT 和MRI 圖像的大小均為256×256。

本文算法與拉普拉斯金字塔(laplacian-pyramid,LP)、離散小波變換(discrete wavelet transform,DWT)[27]、輪廓波變換(contourlet transform, CVT)[28,29]、非下采樣輪廓波變換(nonsubsampled contourlet,NSCT)[30]、稀疏表示(sparse representation, SR)和自適應稀疏表示(adaptive sparse representation, ASR)等方法進行了比較。實驗結(jié)果如圖4～圖6 所示。

圖4 肉瘤患者的融合結(jié)果

圖6 急性中風患者的融合結(jié)果

4.1 融合結(jié)果的指標評價

所有算法的融合結(jié)果通過以下8 個指標進行評價：標準方差(standard deviation, SD)、平均梯度(average gradient, AG)、熵(entropy, E)[23]、空間頻率(spatial frequency, SF)[24]、互信息(mutual information,MI)[25]、融合質(zhì)量QAB/F、融合損失LAB/F和融合偽影NAB/F。其中QAB/F、LAB/F和NAB/F這3 個指標是由文獻[26]提出，具有互補信息，且三者之和為1。

SD、 AG、E、 SF僅由融合圖像F計算得出。因此，設F的圖像大小為M×N，F(xiàn)i,j表示F在(i,j)處的像素值，則標準差 SD反應了F的像素離散程度；SD越大表明像素值分布越分散，越小則表明分布更集中：

式中，Pi為像素值為i的像素在F中出現(xiàn)的概率。

空間頻率 SF反映了灰度值的變化率，可用于評估圖像的細節(jié)及紋理的多少。 SF值越大說明F的紋理邊緣信息更豐富。 SF由橫向頻率、縱向頻率、主對角線頻率和副對角線頻率組成：

式中， RF是F中所有處于同一行的相鄰兩像素差值的平方平均數(shù)；同理， CF、MDF、SDF分別是F中所有處于同一列、同一主對角線、同一副對角線的相鄰兩像素差值的平方平均數(shù)。

此外， MI、QAB/F、LAB/F、NAB/F由輸入圖像A和B以及融合圖像F聯(lián)合計算得出。其中，互信息MI表示F從A和B中獲取的信息總量。 MI越大，表明融合圖像從輸入圖像中獲取的信息量越多，融合效果越好：

式中， |W|為滑動窗口 ω滑動的次數(shù)；s(A|ω)和s(B|ω)為A和B在滑動窗口 ω內(nèi)的顯著性特征；SSIM(A,F|ω)和SSIM(B,F|ω)為F在 ω處與A和B的結(jié)構(gòu)相似度。

融合損失LAB/F是對融合過程中丟失信息的度量，即表示輸出中的有價值信息在融合圖像中丟失的程度。因此，其值越低，丟失的信息就會越少：

融合偽影NAB/F用于評價引入融合圖像的無用信息的數(shù)量。這些信息在任何輸入圖像中都沒有相對應的特征。因此，融合偽影本質(zhì)上就是錯誤的信息，其值越小融合效果越好：

不同算法結(jié)果的評價指標數(shù)值由表1～表3 給出。所有指標中的最優(yōu)數(shù)值加粗表示。從表1 可以看出，本文方法在E、 MI、QAB/F和NAB/F指標上的結(jié)果是最優(yōu)，同時在 AG和 SF指標上取得了次優(yōu)。在表2 中，本文算法在 SD、E、 SF、QAB/F和LAB/F指標上取得最優(yōu)結(jié)果；同時，在 AG和 MI指標上取得次優(yōu)。最后，在表3 中，本文算法在E、S F、QAB/F和LAB/F指標上取得最優(yōu)結(jié)果；同時，在 AG和 MI指標上取得次優(yōu)。

表1 第一組實驗融合結(jié)果的客觀評價指標對比

表2 第二組實驗融合結(jié)果的客觀評價指標對比

表3 第三組實驗融合結(jié)果的客觀評價指標對比

從整體來看，本文算法在E、QAB/F上均取得了最優(yōu)結(jié)果，且在 MI和LAB/F指標上也可獲得最優(yōu)或次優(yōu)結(jié)果。這證明了本文算法相比于其他算法減少了輸入圖像的信息損失，保存了最多的輸入信息，同時具有優(yōu)良的紋理保留和細節(jié)信息的能力，因而在 AG和 SF指標上多次取得了最優(yōu)或次優(yōu)結(jié)果。最后，對于 SD和NAB/F指標，盡管本文算法分別存在一組最優(yōu)結(jié)果，但其整體表現(xiàn)不佳，這意味著該算法的融合結(jié)果存在一定的偽影干擾，且融合圖像的像素分布較為集中。

此外，其他6 種算法在某些指標上的表現(xiàn)也較為優(yōu)異，甚至可以取得超越本文算法的結(jié)果，如表1 和表3 所示，LP 算法在 SD指標上取得了最優(yōu)結(jié)果。表2 中的結(jié)果也表明LP 算法取得了次優(yōu)結(jié)果。而SR 算法和ASR 算法在 AG、 MI、 SF這3 個指標上表現(xiàn)較為優(yōu)異，在3 組實驗中多次取得最優(yōu)和次優(yōu)結(jié)果。最后，DWT 算法在NAB/F指標上的表現(xiàn)較為出色，這表明其融合結(jié)果有較少的偽影。

4.2 融合結(jié)果的視覺評價

為了進一步證明本文所提算法的有效性，對圖4～圖6 中所有算法融合結(jié)果的部分細節(jié)和邊界區(qū)域進行了放大處理，以便在視覺評價方面對比本文算法與其他算法。所有算法的放大結(jié)果如圖7～圖9 所示，其中融合圖像中的方框為選中放大的圖像區(qū)域，而放大后的圖像在原圖像的右側(cè)展示并由對應的箭頭連接。

圖7 第一組實驗融合結(jié)果的部分細節(jié)和邊界區(qū)域放大對比

圖9 第三組實驗融合結(jié)果的部分細節(jié)和邊界區(qū)域放大對比

從圖7a～圖9a 可以看出，LP 算法的融合結(jié)果存在亮度過高的情況，導致該算法會丟失MRI 圖像的部分邊緣細節(jié)。而這種情況同樣出現(xiàn)在CVT算法中，如圖9c 所示。

對于DWT 算法，從表2 可以看出，其在NAB/F指標上取得了最優(yōu)結(jié)果。然而，從圖8b 可以看出，其相對應融合圖像的紋理信息和圖像對比度與其他算法相比明顯較差。這意味著該算法保存的信息較少，從表2 和表3 中也可看出，DWT 算法在與信息保存相關(guān)的QAB/F指標上的數(shù)值結(jié)果是最差的。

圖8 第二組實驗融合結(jié)果的部分細節(jié)和邊界區(qū)域放大對比

對于SR 算法和ASR 算法，盡管在指標評價分析中，兩者的融合結(jié)果在與評價細節(jié)和紋理信息保存相關(guān)的 AG、 MI、 SF指標上都能取得優(yōu)異的數(shù)值，如表1～表3 所示；但在視覺分析中，從圖7e～圖9e 以及圖7f～圖9f 中可以看出，這兩種算法融合結(jié)果的放大區(qū)域與其他算法相比更加模糊，這是由于圖像重建的圖像塊效應導致的。而對于本文算法，盡管也會受到影響，但與SR 算法和ASR 算法相比有明顯的提升。

最后，從整體來看，NSCT 算法和本文所提算法在紋理和邊界信息都取得了較優(yōu)異的視覺效果。但是通過比較圖5b 以及圖8a～圖8f 可以看出，6 種經(jīng)典算法都存在不同程度的邊界損失情況，而通過圖8g 可以看出，本文算法與其他6 種算法相比保存了最多的邊界信息。因而，與NSCT 算法相比，本文算法更具魯棒性。

圖5 腦膜瘤患者的融合結(jié)果

4.3 算法的時間復雜度評價

在算法的時間復雜度方面，與傳統(tǒng)的變換域融合算法(LP、DWT、CVT、NSCT) 相比，基于稀疏表示的算法(SR、ASR、本文算法) 更為復雜。對于ASR 算法來說，每個待融合的圖像塊都要根據(jù)其梯度特征自適應地選擇最適合的子字典，因而該算法會比SR 算法更加耗時。而在本文提出的算法中，ASR 模型僅用于細節(jié)層的融合。因此，為評價SR、ASR 和本文算法的時間復雜度，在Inter Core i7-10875H@2.30 GHz 八核處理器、16.0 GB內(nèi)存及NVIDIA GeForce RTX 2060 6 GB 顯卡的硬件環(huán)境下進行實驗。每個算法的程序以前文所述的3 組實驗圖像為輸入重復運行10 次，計算每次的運行時間，并取均值作為最終結(jié)果來評價該算法的時間復雜度。

3 種算法的平均運行時間如表4 所示。其中，ASR 算法的運行時間在3 組實驗中均為最長，而本文算法的運行時間在所有實驗中均是最短。除此之外，通過對比SR 和ASR 算法的3 組實驗可以看出，不同的輸入圖像對其運行時間影響較大，尤其是在第三組實驗中，SR 和ASR 算法的運行時間均高于前兩組實驗。而對于本文算法來說，其運行時間的變化較為穩(wěn)定，受輸入圖像影響較小。

表4 SR、ASR 和本文算法的平均運行時間s

綜上所述，與SR 和ASR 算法相比，本文算法的時間復雜度更低，且不受輸入圖像影響；然而，與LP、DWT、CVT、NSCT 算法相比，該算法的時間復雜度仍有較大改進空間。

5 結(jié) 束 語

本文提出了一種基于自適應稀疏表示和引導濾波的多模態(tài)圖像融合算法。該算法將待融合圖像分解為細節(jié)層和基礎層；利用基于引導濾波器的加權(quán)融合方法融合基礎層；同時，利用自適應稀疏表示模型融合細節(jié)層。實驗結(jié)果表明，本文算法在紋理和邊界信息保存上均取得了優(yōu)異結(jié)果；同時，時間復雜度也優(yōu)于基于稀疏表示的方法。然而，本文算法在融合圖像的對比度、偽影以及塊效應上存在一定缺陷，需要進一步改進；算法的時間復雜度也需要進一步降低。