(浙江理工大學(xué) 信息學(xué)院,浙江 杭州 310018)
磁共振成像(Magnetic Resonance Imaging,MRI)是一種對(duì)人體無(wú)輻射傷害的診斷方法,具有多方位、多參數(shù)、多模態(tài)等優(yōu)點(diǎn),能夠提供豐富的診斷信息[1],這些優(yōu)點(diǎn)使得MRI成為醫(yī)學(xué)研究的重點(diǎn)領(lǐng)域之一。然而,成像時(shí)間長(zhǎng)是其不可避免的缺點(diǎn),如何重構(gòu)高質(zhì)量的磁共振圖像以便醫(yī)護(hù)人員更好地對(duì)病灶進(jìn)行診斷是研究的重點(diǎn)。壓縮感知[2](Compressed Sensing,CS)打破了傳統(tǒng)的奈奎斯特采樣定理,使得僅采集少量數(shù)據(jù)就能恢復(fù)出高質(zhì)量的信號(hào)成為可能?;趬嚎s感知的磁共振成像[3](CS-MRI)是MRI 的一個(gè)重大突破,利用CS 原理,不僅磁共振數(shù)據(jù)采集量大大減少,而且圖像重構(gòu)質(zhì)量也相應(yīng)提高。
全變差分[4](Total Variation,TV)能夠有效地保留圖像邊緣信息,是一種經(jīng)典的圖像稀疏方法。2016 年,馬杰等[5]提出一種基于全變分正則化低秩稀疏分解的動(dòng)態(tài)MRI重建方法,通過在低秩加稀疏分解模型基礎(chǔ)上,結(jié)合TV 正則化項(xiàng),最終使用非精確增廣拉格朗日算法對(duì)模型求解。該算法在重構(gòu)精度及速度上都有較大提高。Varghees 等[6]提出一種基于TV 最小化模型和本地噪聲估計(jì)技術(shù)以自適應(yīng)地去除圖像中的噪聲,實(shí)驗(yàn)表明,該方法既能保留原始圖像的邊緣信息,又可成功減少Ricaian 噪聲。非局部均值[7](Non-Local Means,NLM)去噪的基本思想是當(dāng)前像素的估計(jì)值由圖像中與它具有相似領(lǐng)域結(jié)構(gòu)的像素加權(quán)平均得到,能夠最大程度地保持圖像的細(xì)節(jié)特征。2015 年,陳創(chuàng)泉[8]對(duì)一些改進(jìn)的NLM 方法加以總結(jié),實(shí)現(xiàn)MRI 圖像去噪,分析利用NLM 技術(shù)進(jìn)行磁共振圖像去噪的關(guān)鍵技術(shù)以供研究人員參考。BM3D(Block Matching of 3D Fil?tering,BM3D)是一種基于NLM 和稀疏表示的性能優(yōu)越的圖像去噪算法[9],該方法通過相似性準(zhǔn)則找到與參考?jí)K相似的二維圖像塊,按照相似度大小堆疊成三維相似組,將其變換到三維空間域后進(jìn)行協(xié)同濾波處理,將最終得到的近似圖像塊進(jìn)行加權(quán)平均后返回到原始圖像塊位置以得到最終恢復(fù)圖像。
BM3D 是公認(rèn)的去噪性能優(yōu)越的算法,基于BM3D 的磁共振圖像重構(gòu)算法陸續(xù)被提出。Eksioglu 等[10]將BM3D算法作為基于噪聲去除的迭代閾值算法(Denoising-based Iterative Thresholding,D-IT)和去噪近似消息傳遞算法(DAMP)的去噪先驗(yàn)進(jìn)行磁共振圖像重構(gòu)[11]。前者的主要思想是將迭代閾值算法中的閾值處理看成是噪聲去除操作,從而使用BM3D 算法實(shí)現(xiàn)噪聲去除,這種方法簡(jiǎn)稱為BM3D-IT-MRI;后者的主要思想是利用現(xiàn)有噪聲去除算法實(shí)現(xiàn)去噪近似消息傳遞算法中的濾波操作。由于磁共振圖像的數(shù)據(jù)具有復(fù)值特性,Eksioglu 等[10]將數(shù)據(jù)分為實(shí)和虛兩個(gè)通道分別加以處理,采用BM3D 算法去除噪聲,這種方法簡(jiǎn)稱為BM3D-AMP-MRI。2015 年,Qu 等[12]提出一種基于塊的非局部算子(Patch-based Nonlocal Operator,PANO),該算法通過PANO 算子盡可能地對(duì)相似塊組進(jìn)行稀疏化,最終實(shí)現(xiàn)磁共振圖像重構(gòu);2017 年,Zha 等[13]提出組稀疏殘差約束的圖像去噪算法,該算法將圖像去噪問題轉(zhuǎn)化成最小化組稀疏殘差問題,首先將BM3D 算法應(yīng)用于噪聲圖像,再對(duì)降噪后的圖像使用組稀疏殘差法估計(jì)出組稀疏集合,相應(yīng)地得到估計(jì)圖像。這種算法是在BM3D 去噪基礎(chǔ)上實(shí)現(xiàn)組稀疏系數(shù)最小化,實(shí)驗(yàn)表明,該算法去噪性能優(yōu)于BM3D 算法。
BM3D-AMP-MRI 是近年來(lái)提出的一種重構(gòu)性能較高的重建算法,其主要思想是將BM3D 作為去噪近似消息傳遞算法的去噪先驗(yàn)。由于GSRC 是一種去噪性能較BM3D優(yōu)越的去噪算法,為獲得更高的磁共振圖像重建精度,本文將GSRC 作為去噪近似消息傳遞的去噪先驗(yàn)。將DAMP 迭代重構(gòu)的圖像估計(jì)值和噪聲標(biāo)準(zhǔn)差作為GSRC 算法的初始噪聲圖像和噪聲標(biāo)準(zhǔn)差,經(jīng)GSRC 去噪后得到的圖像以進(jìn)一步實(shí)現(xiàn)D-AMP 算法的殘差更新。最終實(shí)驗(yàn)結(jié)果表明,基于GSRC 的磁共振圖像重建算法不僅有效保留了原始圖像的細(xì)節(jié)信息,而且提高了重建質(zhì)量,并具有較強(qiáng)的魯棒性。
傳統(tǒng)的磁共振圖像重構(gòu)模型如式(1)所示。
其中,X為待重構(gòu)的MR 圖像,F(xiàn)是傅里葉操作算子,U為欠采樣矩陣,y是測(cè)量的K 空間數(shù)據(jù),Ψ代表稀疏變換。式(1)中前半部分表示數(shù)據(jù)的保真項(xiàng),后半部分λ‖ΨX‖1代表數(shù)據(jù)的懲罰項(xiàng)或者稀疏項(xiàng)目,參數(shù)λ用來(lái)平衡數(shù)據(jù)的保真項(xiàng)和稀疏項(xiàng),使得求解出的X更加接近原始圖像。
組稀疏模型表示如式(2)所示。
組稀疏基本思想如下:X∈RN是一幅清晰的圖像,將X分割成m個(gè)尺寸固定的圖像塊,將這m個(gè)重疊的圖像塊表示為xi,i=1,2,...,m。在固定的窗口中(大小為W×W)尋找n個(gè)相似塊,這些相似塊組表示為Xi,定義為Xi={xi,1,xi,2,...,xi,n}。xi,n表示第i組中第n個(gè)相似塊,Ai表示相似塊組Xi的組稀疏系數(shù)。Di表示稀疏字典,本文使用的是PCA 子字典,使得?!ぁ現(xiàn)表示Frobenius 范數(shù),‖Ai‖p為模型的稀疏項(xiàng),‖·‖p表示lp范數(shù),p的取值為0 或者1,λi是正則化參數(shù)。式(2)表明整個(gè)圖像能夠通過組稀疏編碼的集合{Ai}進(jìn)行表示。
對(duì)于磁共振成像而言,式(1)中原始圖像X未知,其研究目的在于重構(gòu)出與X近似的。本文基本思想是使用迭代磁共振重構(gòu)算法,每次迭代都會(huì)得到一個(gè)與原始磁共振圖像近似的含有噪聲的圖像,假設(shè)為R,從R中提取n個(gè)噪聲圖像塊ri,i=1,2,...m,找到噪聲圖像塊ri的相似塊,并組成相似塊組Ri={ri,1,ri,2,...,ri,n}。此時(shí)磁共振圖像去噪轉(zhuǎn)化為利用組稀疏編碼從Ri中恢復(fù)Xi。
當(dāng)獲得所有組稀疏系數(shù){A_Ei},通過Xi重建原始磁共振圖像X。
對(duì)比式(2)和式(3),將求解問題轉(zhuǎn)化為縮小估計(jì)的組稀疏系數(shù)A_E與原始圖像的組稀疏系數(shù)A之間的誤差。因此,可以將問題轉(zhuǎn)化為通過定義殘差進(jìn)行相應(yīng)的求解計(jì)算,即組稀疏系數(shù)A_E與組稀疏系數(shù)A的差值,用Red表示,如式(4)所示。最終將求解問題轉(zhuǎn)換為如何減小組稀疏系數(shù)殘差Red問題,其模型表示如式(5)所示。
真實(shí)的組稀疏編碼A是未知的,本文首先對(duì)迭代過程中噪聲圖像R使用BM3D 方法去噪,將去噪后的結(jié)果表示為Z,用Z的組稀疏系數(shù)取代未知A。因此式(5)中,只有A_Ei是待求解的變量。
文獻(xiàn)[13]證明式(5)在第t+1 次迭代中,計(jì)算方式如式(6)所示。
其中,Sλ(?)是軟閾值操作算子,表示第i次重構(gòu)相似塊組。相似塊組Redi對(duì)應(yīng)參數(shù)λi設(shè)置如式(7)所示。
其中,σi表示Red的估計(jì)方差,c 是一個(gè)非常小的常量。在第t+1 次迭代過程中,噪聲標(biāo)準(zhǔn)差計(jì)算如式(8)所示。
其中,γ表示一個(gè)小常數(shù)。當(dāng)求出A_Ei后,重構(gòu)的圖像塊組,最終通過聚合所有圖像塊組得到去噪后的圖像。
D-AMP 算法是一種基于噪聲去除的算法,本文使用該算法對(duì)磁共振圖像進(jìn)行重構(gòu),并將GSRC 去噪算法作為該重構(gòu)算法的先驗(yàn)知識(shí)。最終結(jié)果表明,該方法能夠重構(gòu)出質(zhì)量相對(duì)較高的圖像?;诮M稀疏殘差去噪的近似消息傳遞磁共振圖像重構(gòu)流程如下:①輸入零填充的磁共振圖像X0=(FΩ表示欠采樣傅里葉算子),迭代殘差z0=y;②計(jì)算初步得到的估計(jì)值;③計(jì)算噪聲標(biāo)準(zhǔn)差,N為y的長(zhǎng)度;④使用GSRC 算法對(duì)Riter去噪,;⑤計(jì)算Onsager校正項(xiàng);⑥更新殘差ziter=y-FΩ Xiter+qiter;⑦判斷迭代次數(shù)是否達(dá)到最大次數(shù),滿足則輸出最終重構(gòu)值,否則返回步驟②。
本文使用的實(shí)驗(yàn)數(shù)據(jù)如圖1(a)—圖1(c)所示,實(shí)驗(yàn)使用3 組數(shù)據(jù)比較幾種重建算法的重建表現(xiàn)。這3 組數(shù)據(jù)分別是人的腦部[14](Brain)、心臟[15](Heart)及上半身圖[14](Bust)。其中,Heart 的數(shù)據(jù)來(lái)源于文獻(xiàn)[14]中三維心臟的第十幀數(shù)據(jù)。采樣模式不同會(huì)影響圖像重構(gòu)質(zhì)量,本次實(shí)驗(yàn)將選用偽徑向[16]和笛卡爾采樣兩種模式對(duì)原始磁共振數(shù)據(jù)進(jìn)行欠采樣,圖1(d)和圖1(e)展示了采樣率為20%的偽徑向采樣和30% 的笛卡爾采樣mask。為了驗(yàn)證算法的有效性和準(zhǔn)確性,本文算法均運(yùn)行在Matalab R2016a 軟件上,主機(jī)CPU 的配置為Intel Core i5,內(nèi)存為8G,操作系統(tǒng)為Windows 10。
對(duì)比算法有基于BM3D 去噪的迭代閾值磁共振圖像重構(gòu)算法[10](BM3D-IT-MRI)、基于塊的非局部算子PA?NO[12]以及基于BM3D 的近似消息傳遞算法[10](BM3DAMP-MRI)算法。為了衡量重構(gòu)效果和算法性能,本實(shí)驗(yàn)計(jì)算峰值信噪比(Peak to Signal Ratio,PSNR)和相對(duì)L2范數(shù)誤差(RelativeL2Norm Error,RLNE)。PSNR 是一種常用的評(píng)價(jià)圖像重構(gòu)質(zhì)量好壞的指標(biāo),PSNR 值越高表明圖像重構(gòu)性能越好。RLNE 用來(lái)表示重構(gòu)圖像與原始圖像的偏差,RLNE 值越小,表明重構(gòu)圖像越接近原始圖像。除使用PSNR 和RLNE 兩個(gè)指標(biāo)外,將會(huì)對(duì)重構(gòu)圖像的細(xì)節(jié)進(jìn)行放大分析,觀察不同重構(gòu)算法保留原始圖像細(xì)節(jié)能力上的差異。
Fig.1 Experimental source data and two sampling masks圖1 實(shí)驗(yàn)源數(shù)據(jù)以及兩種采樣mask
表1、表2 展示了在采樣率為20% 和30% 的條件下,磁共振圖像重構(gòu)的PSNR 和RLNE 值。其中,加粗?jǐn)?shù)據(jù)表示相同條件下得到的最高PSNR 和最低的RLNE 值。根據(jù)表1 計(jì)算得出,在20% 采樣率下,本文算法重構(gòu)的PSNR值較BM3D-IT-MRI、PANO 及BM3D-AMP-MRI 算法平均增加3.8dB、1.3dB 和1dB,RLNE 平均降低0.060 7、0.016 5、0.017 7。因此,可以判定本文重構(gòu)性能高于其它3 種算法。對(duì)比發(fā)現(xiàn),相同采樣率下,偽徑向采樣模式下重構(gòu)圖像的PSNR 均高于笛卡爾采樣模式下的PSNR 值、RLNE值均低于笛卡爾采樣模式下的RLNE 值,可知偽徑向采樣是一種優(yōu)于笛卡爾采樣的模式。比較表1 和表2 可知,對(duì)于同一源數(shù)據(jù),相同采樣模式下,20% 欠采樣的圖像重構(gòu)質(zhì)量比30% 欠采樣重構(gòu)的圖像質(zhì)量差,說明采樣率越高,獲取到源數(shù)據(jù)的有用信息越多,最終重構(gòu)質(zhì)量也越高。
Table 1 The PSNR/RLNE of reconstructed image under 20% under-sampling表1 20% 欠采樣條件下圖像重構(gòu)的PSNR/RLNE
Table 2 The PSNR/RLNE of reconstructed image under 30% under-sampling表2 30% 欠采樣條件下圖像重構(gòu)的PSNR/RLNE
20% 偽徑向采樣模式下,Bust 圖像的重構(gòu)細(xì)節(jié)放大圖及放大10 倍后的重構(gòu)誤差如圖2 所示(彩圖掃OSID 碼可見),第一行從左到右分別為原圖、BM3D-IT-MRI 重構(gòu)細(xì)節(jié)圖、PANO 重構(gòu)細(xì)節(jié)圖、BM3D-AMP-MRI 重構(gòu)細(xì)節(jié)圖以及本文算法的重構(gòu)細(xì)節(jié)圖;第二行從左到右分別為20% 偽徑向采樣mask、BM3D-IT-MRI 重構(gòu)誤差圖、PANO 重構(gòu)誤差圖、BM3D-AMP-MRI 重構(gòu)誤差圖以及本文算法重構(gòu)誤差圖。圖2 的第一行內(nèi)容中,較小的紅色框內(nèi)表示待放大重構(gòu)區(qū)域,較大的紅色框內(nèi)表示經(jīng)放大后的待觀察重構(gòu)區(qū)域,對(duì)比發(fā)現(xiàn),本文算法對(duì)原始細(xì)節(jié)信息保留得更好,更加接近原始圖像。并且,由圖2 的第二行可以看出,本文算法的重構(gòu)誤差最小。圖3 展示了30% 笛卡爾采樣模式下,Brain 圖像重構(gòu)的細(xì)節(jié)放大圖及放大10 倍后的重構(gòu)誤差圖(彩圖掃OSID 碼可見),第一行從左到右分別為:原圖、BM3D-IT-MRI 重構(gòu)細(xì)節(jié)圖、PANO 重構(gòu)細(xì)節(jié)圖、BM3DAMP-MRI 重構(gòu)細(xì)節(jié)圖以及本文算法的重構(gòu)細(xì)節(jié)圖;第二行從左到右分別為:30% 笛卡爾采樣mask、BM3D-IT 重構(gòu)誤差圖、PANO 重構(gòu)誤差圖、BM3D-AMP-MRI 重構(gòu)誤差圖以及本文算法重構(gòu)誤差圖。通過對(duì)比圖2 和圖3 結(jié)果可知,本文提出的重建算法能夠取得令人滿意的結(jié)果,重構(gòu)后的磁共振圖像邊緣較清晰,視覺效果較對(duì)比算法佳。
Fig.2 Experimental results(1)圖2 實(shí)驗(yàn)結(jié)果展現(xiàn)(一)
Fig.3 Experimental results(2)圖3 實(shí)驗(yàn)結(jié)果展現(xiàn)(二)
采樣因子大小對(duì)重構(gòu)算法會(huì)產(chǎn)生影響,采樣因子越大,圖像欠采樣率越低,表明欠采樣數(shù)據(jù)越低,采集到的原圖像有效信息越少。本文分別研究在采樣因子為2~10 的條件下,圖像重構(gòu)質(zhì)量變化情況。圖4(a)和圖4(b)展示了偽徑向采樣下,采樣因子對(duì)4 種算法重構(gòu)質(zhì)量的影響。觀察可知,當(dāng)采樣因子增加時(shí),重構(gòu)的PSNR 值呈下降趨勢(shì),同時(shí)RLNE 值呈增長(zhǎng)趨勢(shì),說明圖像重構(gòu)質(zhì)量在下降。
由圖4(a)和圖4(b)可知,本文算法與BM3D-AMPMRI 算法重構(gòu)的PSNR 和RLNE 非常接近,兩者重構(gòu)質(zhì)量相當(dāng)。但本文算法的重構(gòu)性能明顯優(yōu)于BM3D-IT-MRI 和PANO 算法。圖4(c)和圖4(d)展示了在笛卡爾采樣模式下,采樣因子對(duì)圖像重構(gòu)的影響。同樣可以得出,采樣因子大小與圖像重構(gòu)質(zhì)量成反比。但是可以清晰地看出,本文算法重構(gòu)的PSNR 值和RLNE 值明顯優(yōu)于其它3 種算法。不同于偽徑向采樣,笛卡爾采樣條件下本文算法的重構(gòu)性能明顯優(yōu)于BM3D-AMP-MRI 的重建性能。一方面,表明本文算法在笛卡爾采樣模式下具有更加優(yōu)異的重建性能;另一方面,說明不同的采樣模式對(duì)圖像重構(gòu)結(jié)果有著至關(guān)重要的影響。
Fig.4 The PSNR and RLNE of reconstructed Brain image圖4 Brain 圖像重構(gòu)的PSNR 與RLNE
本文提出基于組稀疏殘差去噪的磁共振圖像重構(gòu)算法,將組稀疏殘差作為去噪近似消息傳遞算法的去噪先驗(yàn),實(shí)現(xiàn)磁共振圖像重建。在Matlab 仿真平臺(tái)上,對(duì)比了3種磁共振圖像重建算法,研究不同采樣模式和采樣因子條件下,各算法的重建性能。實(shí)驗(yàn)表明,本文提出的算法具有更高的PSNR 值和更低的RLNE 值,同時(shí)保留了原始圖像更多的細(xì)節(jié)信息,具有一定的應(yīng)用價(jià)值。由于實(shí)際應(yīng)用中,采集到的磁共振更多的是三維K 空間數(shù)據(jù),因此,如何將本文算法拓展到三維動(dòng)態(tài)磁共振成像中,是下一步研究的重點(diǎn)。