基于自適應聯合字典學習的腦部多模態圖像融合方法

王麗芳,董 俠,秦品樂,高 媛

0 引言

多模態醫學影像在臨床醫療診斷中發揮越來越重要的作用。醫學圖像可分為解剖圖像和功能圖像,計算機斷層成像(Computed Tomography, CT)和磁共振成像(Magnetic Resonance Imaging,MRI)均為解剖圖像;單光子發射計算機斷層成像(Single Photon Emission Computed Tomography,SPECT)和正電子發射斷層成像(Positron Emission Tomography,PET)均為功能圖像。解剖圖像能夠反映器官的解剖結構,分辨率較高;功能圖像能夠顯示人體新陳代謝的功能變化,分辨率較低。單一類型的圖像往往無法為病情診斷提供全面充足的信息[1]。多模態醫學圖像融合技術能夠將互補的多源信息融合到一幅圖像中以改善圖像質量并減少隨機性和冗余性,有助于提高臨床醫療診斷和分析的準確性[2]。

近年來,隨著壓縮感知理論[3]的提出,基于稀疏表示的方法[4-5]被廣泛用于圖像融合領域。Mitianoudis等[6]首次使用滑窗技術將源圖像劃分成圖像塊并用于圖像融合,實驗結果表明該方法能夠很好地提取源圖像的局部顯著特征和保持平移不變性。由于在稀疏域中使用活躍度測量具有可靠性和空間滑窗技術的平移不變性,基于稀疏表示的圖像融合方法一般優于傳統多尺度變換的方法。Yang等[7]首次將稀疏表示的方法用于多聚焦圖像融合,該方法使用冗余的離散余弦變換(Discrete Cosine Transform,DCT)字典表示源圖像,并使用“L1范數最大”的規則融合稀疏表示系數。DCT字典屬于解析字典,具有快速實現的優點,但這種字典僅限于特定類型的信號表示,并且在選擇解析基時需要先驗知識。Aharon等[8]提出K奇異值分解(K-means-based Singular Value Decomposition,K-SVD)算法用于從訓練圖像中學習字典。相對于解析字典,學習的字典包含更豐富的圖像特征信息,表示能力更強,但是時間效率較低。Liu等[9]從40幅自然圖像中隨機采樣10萬個圖像塊作為字典訓練數據,低頻子帶對應的稀疏表示系數使用“L1范數最大”的規則進行融合。Yin等[10]使用K-SVD算法在USC-SIPI數據集[11]上學習字典。上述方法學習的字典均屬于全局訓練字典,全局訓練字典使用的訓練數據來自一些公開的圖像數據庫,是由外部預先收集的圖像組成。在實際應用中,收集合適的圖像數據集并不總是可行的。此外,不同類型圖像的圖像內容千差萬別,全局訓練字典對于不同類型的圖像表現出來的性能差異較大。

醫學圖像具有對比度低、灰度分布不均勻等特點,全局訓練字典對于醫學圖像的自適應性不強,難以適應多模態醫學圖像融合的需要。因此,在基于稀疏表示的多模態醫學圖像融合方法中,如何構建一個自適應于源圖像的過完備字典是關鍵。宗靜靜等[12]將兩幅已配準的源圖像組成一個整體進行聯合稀疏表示,獲得一個共同稀疏部分和不同稀疏部分,其中對不同稀疏部分同樣使用“L1范數最大”的規則進行融合。該方法直接從已配準的源圖像本身學習字典可以提高字典對于源圖像的自適應性。但是,從圖像融合規則來看,傳統的使用稀疏表示系數的L1范數取極大的方式在融合多模態醫學圖像時,易造成融合結果的灰度不連續效應。

針對上述問題,本文提出一種基于自適應聯合字典學習的腦部多模態圖像融合方法。首先,為實現功能圖像和解剖圖像的融合,該方法使用亮度-色調-飽和度(Intensity-Hue-Saturation,IHS)變換將已與解剖圖像配準的功能圖像變換到IHS空間,再將其對應的I分量(代表灰度圖像)和解剖圖像作為待融合的源圖像;其次,不同于解析字典和全局訓練字典,提出的自適應聯合字典的學習不需要任何先驗知識,也不需要外部預先收集的圖像數據,而是將直接從待融合的源圖像中學習得到的子字典組合成自適應聯合字典來提高字典對于每一對輸入圖像的自適應性。同時,在融合階段,提出稀疏表示系數的“多范數”綜合地衡量源圖像塊的活躍度,并使用一種結合“自適應加權平均”與“選擇最大”的無偏規則融合稀疏表示系數。對腦部多模態圖像進行的大量實驗結果表明,本文方法比對比方法可以保留更多的圖像細節信息,取得較好的圖像融合效果。

1 稀疏表示和字典學習

稀疏表示的基本假設是一個信號x∈Rn可以近似地表示為一個超完備字典D∈Rn×m(n

(1)

其中:ε≥0表示允許偏差的精度,問題(1)的求解過程稱為“稀疏編碼”。

(2)

s.t. ?i‖αi‖0≤τ

2 自適應聯合字典

為提高字典對于每一對輸入醫學圖像的自適應性,本文提出使用改進的K-SVD算法[14]分別對兩幅已配準的源圖像進行字典學習得到兩個子字典再組合成自適應聯合字典,使每一幅源圖像均可以用字典中相同的原子表示。自適應聯合字典的學習過程如圖1所示。

在字典學習代價函數(2)中加入支撐完整的先驗信息[17],以源圖像I1為例,將X1訓練的優化問題表示如下:

s.t. ‖A‖0≤τ,A⊙M=0

(3)

其中:字典D1∈Rw×S,稀疏表示系數A∈RS×L。⊙代表點乘,掩膜矩陣M由元素0和1組成,定義為M={|A|=0},等價于若A(i,j)=0則M(i,j)=1,否則為0,因此A⊙M=0能使A中所有零項保持完備。式(3)的求解過程分為稀疏編碼和字典更新兩個階段：

圖1 自適應聯合字典的學習過程

s.t.A⊙M=0

(4)

分別對系數矩陣A中每一列的非零元素進行處理,而保持零元素完備,則式(4)可以轉換為式(5):

(5)

其中:D1i是D1對應A的非零支集的子矩陣；αi是A第i列的非零部分。問題(5)由系數重用正交匹配追蹤(Coefficient Reuse Orthogonal Matching Pursuit,CoefROMP)算法[14]求解,由此可得到更新的稀疏系數矩陣A。

其次,在字典更新階段,式(3)的優化問題可以轉化為:

s.t.A⊙M=0

(6)

則式(6)的補償項可以寫為:

(7)

最后,循環執行上述兩個步驟,直至達到預設的迭代次數為止,輸出字典D1,字典D2的訓練過程同字典D1,然后將得到的子字典D1和D2組合成自適應聯合字典D=[D1,D2],D∈Rw×2S。由于字典更新階段同時更新字典和稀疏表示系數的非零元素,使得字典的表示誤差更小且字典的收斂速度更快。在稀疏編碼階段,考慮到每次迭代時都忽略前一次迭代的表示,本文使用CoefROMP算法根據上次迭代的稀疏表示殘差信息進行系數更新,從而更快地得到所要求問題的解[17]。

3 多模態醫學圖像融合過程

圖像融合方式在很大程度上決定了融合圖像的質量。傳統基于稀疏表示的醫學圖像融合方法采用稀疏表示系數的L1范數取極大的融合方式易造成融合圖像的灰度不連續效應。這是因為對于多模態醫學圖像的同一區域,有可能在一幅源圖像中灰度較高,而在另一幅源圖像中灰度較低,但在兩幅源圖像中都是很“平坦”的,即區域內部灰度變化很小,但這兩幅源圖像在此區域仍會存在細節豐富程度的差異。由于差異很小,稀疏表示系數L1范數取極大的融合方式很難準確地判斷具有較多細節的源圖像塊,得到錯誤判定結果的可能性增大,同時對有噪聲的情況非常敏感,從而導致對此種類型區域的不同圖像塊的判定結果不一致而最終造成融合圖像的灰度不連續效應。

為解決此問題,考慮以下兩個方面：一是如何衡量源圖像塊的活躍度,另一個是建立什么樣的規則融合稀疏表示系數。對于第一個問題,僅僅使用L1范數難以全面地計算源圖像塊的活躍度,本文提出將稀疏表示系數的“多范數”作為源圖像塊的活躍度測量以保留更多的圖像細節信息;對于第二個問題,尋求能夠使融合圖像保留源圖像所有視覺信息的理想的融合規則是困難的[18]。考慮到“加權平均”的融合規則可以有效抑制噪聲,“選擇最大”的融合規則可以保留源圖像最重要的信息,本文提出“自適應加權平均”與“選擇最大”兩者相結合的無偏規則。根據稀疏表示系數的“多范數”的相似度來選擇融合規則,對于相似度大于閾值,即圖像塊細節豐富程度差異較小的情況,使用“自適應加權平均”的規則；反之則使用“選擇最大”的規則融合稀疏表示系數。

3.1 稀疏編碼

(8)

(9)

其中:D是第2章提出的自適應聯合字典。

圖2 多模態醫學圖像融合過程

3.2 多范數無偏融合規則

稀疏系數向量α的L1范數表示向量中各個元素的絕對值之和,可以反映細節信息的重要性,對應的值越大,圖像塊的顯著性越好;稀疏系數向量α的L0范數表示向量中非零項的個數,可以反映細節信息的集中度,對應的值越大,圖像塊包含的細節信息越多。為全面地計算源圖像塊的活躍度以保留更多的圖像細節信息,對于每一個源圖像塊,使用加權指數分別對L1范數和L0范數進行約束,再將兩者的乘積作為兩種范數對源圖像塊的綜合度量結果,公式表示如下:

(10)

(11)

(12)

當稀疏表示系數的“多范數”的相似度大于等于閾值(Si≥T)時,使用“自適應加權平均”的融合規則:

(13)

(14)

ω2=1-ω1

(15)

當稀疏表示系數的“多范數”的相似度小于閾值(Si

(16)

(17)

4 實驗結果與分析

為驗證本文方法的性能,本文進行了大量的仿真實驗,現選取三組已配準的腦部多模態醫學圖像(大小均為256×256像素)作為待融合的源圖像,腦部圖像數據均來源于哈佛大學醫學院[19]。圖3為腦部腫瘤轉移、腦部Ⅱ級星形細胞腫瘤以及腦部阿爾茨海默病的圖像。實驗的硬件環境為Intel Core i3-2350M @ 2.30 GHz處理器、4 GB運行內存,軟件環境為32位Windows 7操作系統、Matlab R2013a。

實驗從主客觀兩個方面將本文方法與其他八種方法進行比較。選取的對比方法有:1)基于多尺度變換的方法，包括:離散小波變換(Discrete Wavelet Transform, DWT)[20]、拉普拉斯金字塔(Laplacian Pyramid, LP)[21]、非下采樣輪廓波變換(Non-Subsampled Contourlet Transform, NSCT)[22];2)基于稀疏表示的方法，包括:基于DCT字典與多范數無偏規則的方法(method based on DCT dictionary and Multi-Norm-unbiased rule,DCT-MN)、基于K-SVD算法的全局訓練字典與多范數無偏規則的方法(method based on the global trained dictionary usingK-SVD algorithm and Multi-Norm-unbiased rule, KSVD-MN)、基于多尺度變換和稀疏表示的方法(method based on Multi-Scale Transform and Sparse Representation,MST-SR)[9]、基于聯合稀疏表示的方法(method based on Joint Sparse Representation, JSR)[12]、基于自適應聯合字典與L1范數最大規則的方法(method based on Adaptive Joint Dictionary and max-L1-norm rule, AJD-L1)。本文提出的基于自適應聯合字典與多范數無偏規則的方法記為AJD-MN(method based on Adaptive Joint Dictionary and Multi-Norm-unbiased rule)。

圖3 三對腦部多模態醫學圖像

圖4 圖像融合結果

對于DWT和LP方法,分解水平均設為3,DWT方法的小波基設為“db6”。NSCT方法使用“9-7”金字塔濾波器和“c-d”方向濾波器,分解水平設為{22,22,23,24}。基于稀疏表示的方法中滑動步長均為1,圖像塊大小均為8×8,字典大小均為64×512,誤差ε=0.01,稀疏度τ=8。KSVD-MN和MST-SR方法的字典訓練樣本是從40幅自然圖像[23]中隨機采樣10萬個圖像塊,MST-SR方法的多尺度變換選擇LP,分解水平設為3。AJD-MN方法使用改進的K-SVD算法,執行6個字典更新周期(Dictionary Update Cycles, DUC)和30次迭代;提出的多范數使用相等的加權指數β=γ=1;經多次實驗驗證,無偏融合規則的閾值T=0.9時融合效果較好。

圖4為不同的融合方法在三組實驗條件下的運行結果,從主觀視覺效果分析,DWT方法的融合圖像的顏色、亮度和對比度明顯失真,存在塊效應和偽影;LP和NSCT方法的融合質量相對較好,但是融合圖像的對比度較低;DCT-MN、KSVD-MN、MST-SR、JSR和AJD-L1方法與上述三種基于多尺度變換的方法(DWT、LP和NSCT)相比,融合圖像的對比度有很大的提升,但在組織和病灶區域存在一些偽影;AJD-MN方法的融合圖像的偽影較少,并且能夠保留更多的圖像細節信息。綜上所述,本文提出的AJD-MN方法的融合圖像在顏色、亮度、對比度、清晰度以及細節的保留度上均優于其他八種方法,病灶清晰,有助于臨床診斷。

為定量地對融合圖像的質量進行客觀評價,本文使用六種指標:標準差(Standard Deviation,SD)[24]、空間頻率(Spatial Frequency,SF)[25]、互信息(Mutual Information,MI)[26]、基于梯度的指標QAB/F[27]、基于通用圖像質量指標Qw[28]和平均結構相似指標(Mean Structural SIMilarity,MSSIM)[29]。其中:SD衡量圖像信息的豐富程度,反映融合圖像的對比度;SF反映圖像的清晰度;MI反映融合圖像從源圖像提取的信息量的多少;QAB/F反映融合圖像對源圖像邊緣特性的保留情況;Qw反映融合圖像在系數相關性、光照和對比度方面與源圖像的關聯;MSSIM指標是融合圖像分別與兩幅源圖像的結構相似指標SSIM的均值,SSIM反映融合圖像與源圖像在亮度、對比度和結構方面的相似度。上述六種指標均是對應的值越大融合圖像質量越高。

表1為圖4～6的三組實驗中不同融合方法的客觀評價指標。從表1可看出,DWT方法在SD、SF、MI、QAB/F、Qw和MSSIM六項指標上均不理想,主要是因為離散小波變換具有平移變化性、走樣和缺乏方向性等缺陷;由于LP方法具有較好的多尺度和多分辨率特性,NSCT具有多尺度、多方向和平移不變性,因此,LP和NSCT方法在上述六項指標上均明顯優于DWT方法;上述三種多尺度變換的方法在SD和SF指標上均低于其他方法,表明融合圖像的對比度可能受到了損失;DCT-MN、KSVD-MN、MST-SR方法在MI指標上優于多尺度變換的方法,而在QAB/F、Qw和MSSIM指標上存在一些低于多尺度變換方法的情況,可能是由于字典對于醫學圖像的自適應性不強,其表示能力不足以提取足夠多的細節信息或者不能精確地表示醫學圖像的特征;JSR和AJD-L1方法在QAB/F和MSSIM指標上存在一些低于多尺度變換方法的情況,可能是由于L1范數取極大的融合規則導致圖像的灰度不連續效應;本文提出的AJD-MN方法存在一些次優值,但是均接近于最優值。經計算, AJD-MN方法的SD、SF、MI、QAB/F、Qw和MSSIM六項指標在上述三組實驗條件下的均值分別為71.078 3、21.970 8、3.679 0、0.660 3、0.735 2和0.733 9。與其他八種方法相比,本文提出的AJD-MN方法在客觀評價方面也具有優勢。

表1 九種融合方法的客觀指標

注:加粗數據為最優值。

5 結語

針對目前全局訓練字典對于腦部醫學圖像的自適應性不強的問題,本文提出使用改進的K-SVD算法自適應地從輸入源圖像中學習子字典并組合成自適應聯合字典的方法；同時,針對稀疏表示系數的L1范數取極大的融合方式易造成圖像的灰度不連續效應的問題,提出用稀疏表示系數的“多范數”全面地衡量源圖像塊的活躍度,并使用“自適應加權平均”與“選擇最大”相結合的無偏規則融合稀疏表示系數。在實驗階段,分別對腦部CT/MRI圖像和腦部MRI/PET圖像進行多次實驗,結果表明,與其他八種方法相比,本文方法能夠取得最佳的圖像融合效果,可以有效輔助臨床醫療診斷。

但仍存在需要進一步研究的問題:本文提出的字典學習方法是將源圖像中所有圖像塊用于字典學習過程,這樣往往會帶來很多冗余信息。如何構造一個非冗余但信息豐富的字典來更好地提高圖像融合質量和時間效率,將是下一步的研究重點。

參考文獻(References)

[1] LI S T, KANG X D, FANG L Y, et al. Pixel-level image fusion: a survey of the state of the art[J]. Information Fusion, 2017, 33: 100-112.

[2] JAMES A P, DASARATHY B V. Medical image fusion: a survey of the state of the art[J]. Information Fusion, 2014, 19(3): 4-19.

[3] DONOHO D L. Compressed sensing[J]. IEEE Transactions on Information Theory, 2006, 52(4): 1289-1306.

[4] LIU Y, WANG Z F. Simultaneous image fusion and denoising with adaptive sparse representation[J]. IET Image Processing, 2015, 9(5): 347-357.

[5] YIN H T. Sparse representation with learned multiscale dictionary for image fusion[J]. Neurocomputing, 2015, 148(148): 600-610.

[6] MITIANOUDIS N, STATHAKI T. Pixel-based and region-based image fusion schemes using ICA bases[J]. Information Fusion, 2007, 8(2): 131-142.

[7] YANG B, LI S T. Multifocus image fusion and restoration with sparse representation[J]. IEEE Transactions on Instrumentation & Measurement, 2010, 59(4): 884-892.

[8] AHARON M, ELAD M, BRUCKSTEIN A.K-SVD: an algorithm for designing overcomplete dictionaries for sparse representation[J]. IEEE Transactions on Signal Processing, 2006, 54(11): 4311-4322.

[9] LIU Y, LIU S P, WANG Z F. A general framework for image fusion based on multi-scale transform and sparse representation[J]. Information Fusion, 2015, 24: 147-164.

[10] YIN H T, LI S T. Multimodal image fusion with joint sparsity model[J]. Optical Engineering, 2011, 50(6): 7007.

[11] WEBER A. The USC-SIPI image database[EB/OL]. [2017-09-20]. http: //sipi.usc.edu/database/.

[12] 宗靜靜, 邱天爽, 郭冬梅. 聯合稀疏表示的醫學圖像融合及同步去噪[J]. 中國生物醫學工程學報, 2016, 35(2): 133-140.(ZONG J J, QIU T S, GUO D M. Simultaneous medical image fusion and de-noising with joint sparse representation[J]. Chinese Journal of Biomedical Engineering, 2016, 35(2): 133-140.)

[13] 董俠, 王麗芳, 秦品樂, 等. 改進耦合字典學習的腦部CT/MR圖像融合方法[J]. 計算機應用, 2017, 37(6): 1722-1727.(DONG X, WANG L F, QIN P L, et al. CT/MR brain image fusion method via improved coupled dictionary learning[J]. Journal of Computer Applications, 2017, 37(6): 1722-1727.)

[14] SMITH L N, ELAD M. Improving dictionary learning: multiple dictionary updates and coefficient reuse[J]. IEEE Signal Processing Letters, 2013, 20(1): 79-82.

[15] TU T M, SU S C, SHYU H C, et al. A new look at IHS-like image fusion methods[J]. Information Fusion, 2001, 2(3): 177-186.

[16] DANESHVAR S, GHASSEMIAN H. MRI and PET image fusion by combining IHS and retina-inspired models[J]. Information Fusion, 2010, 11(2): 114-123.

[17] 練秋生, 石保順, 陳書貞. 字典學習模型、算法及其應用研究進展[J]. 自動化學報, 2015, 41(2): 240-260.(LIAN Q S, SHI B S, CHEN S Z. Research advances on dictionary learning models, algorithms and applications[J]. Acta Automatica Sinica, 2015, 41(2): 240-260.)

[18] ZHANG Q, LIU Y, BLUM R S, et al. Sparse representation based multi-sensor image fusion for multi-focus and multi-modality images: a review[J]. Information Fusion, 2017,40:57-75.

[19] JOHNSON K A, BECKER J A. The whole brain Atlas[DB/OL]. [2017-09-20]. http: //www.med.harvard.edu/aanlib/home.html.

[20] PAJARES G, CRUZ J M D L. A wavelet-based image fusion tutorial[J]. Pattern Recognition, 2004, 37(9): 1855-1872.

[21] BURT P J, ADELSON E H. The Laplacian pyramid as a compact image code[J]. IEEE Transactions on Communications, 1983, 31(4): 532-540.

[22] ZHANG Q, GUO B L. Multifocus image fusion using the nonsubsampled contourlet transform[J]. Signal Processing, 2009, 89(7): 1334-1346.

[23] The imagefusion image database [EB/OL]. [2017-09-20]. http: //www.imagefusion.org.

[24] SHI W Z, ZHU C Q, TIAN Y, et al. Wavelet-based image fusion and quality assessment[J]. International Journal of Applied Earth Observation & Geoinformation, 2005, 6(3/4): 241-251.

[25] ESKICIOGLU A M, FISHER P S. Image quality measures and their performance[J]. IEEE Transactions on Communications, 1996, 43(12): 2959-2965.

[26] QU G H, ZHANG D L, YAN P F. Information measure for performance of image fusion[J]. Electronics Letters, 2002, 38(7): 313-315.

[27] XYDEAS C S, PETROVIC V. Objective image fusion performance measure[J]. Military Technical Courier, 2000, 36(4): 308-309.

[28] PIELLA G, HEIJMANS H. A new quality metric for image fusion[C]// ICIP 2003: Proceedings of the 2003 International Conference on Image Processing. Washington, DC: IEEE Computer Society, 2003: 173-176.

[29] WANG Z, BOVIK A C, SHEIKH H R, et al. Image quality assessment: from error visibility to structural similarity[J]. IEEE Transactions on Image Processing, 2004, 13(4): 600-612.

This work is partially supported by the Natural Science Foundation of Shanxi Province (2015011045).