袁小君,蔣明峰,楊曉城,李 楊

(浙江理工大學(xué) 信息學(xué)院,杭州 310018)

1 引言

磁共振成像(Magnetic Resonance Imaging,MRI)是以核磁共振(Nuclear Magnetic Resonance,NMR)為物理基礎(chǔ)的一種成像技術(shù),該技術(shù)通過在靜磁場中將一定頻率的射頻脈沖施加到人體上,以刺激人體內(nèi)的氫質(zhì)子引起共振現(xiàn)象.與計(jì)算機(jī)X 射線斷層掃描成像技術(shù)(Computed Tomography,CT)相比,MRI 不僅無需使用造影劑和沒有電離輻射,而且可以直接做出橫斷面、矢狀面和冠狀面等優(yōu)點(diǎn)[1].壓縮感知[2]理論的提出,打破了傳統(tǒng)的奈奎斯特采樣定律的限制.基于壓縮感知的磁共振成像(CS-MRI)[3]技術(shù)大大推動了磁共振成像的發(fā)展.基于CS-MRI 的重構(gòu)問題可以表示如下:

圖像去噪是一個較為成熟的研究領(lǐng)域,實(shí)現(xiàn)去噪的途徑之一是利用圖像的稀疏性.常用的圖像稀疏方法有傅里葉變換、全變差變換、小波變換以及字典學(xué)習(xí)等.近年來,基于非局部自相似塊思想[4]的低秩去噪方法在CS-MRI 領(lǐng)域中得到應(yīng)用.非局部自相似性的原理通過在圖像中尋找與參考塊相似的相似塊組,利用該相似塊組所組成的矩陣具有的低秩屬性進(jìn)行圖像去噪.求解矩陣的秩函數(shù)是NP 難問題,常使用核范數(shù)來近似代替秩函數(shù)的求解.其中最經(jīng)典的方法是核范數(shù)最小化[5]算法.NNM 通過對待修復(fù)的矩陣奇異值分解,對奇異值進(jìn)行閾值處理來求解核范數(shù)最小化問題,該方法稱為奇異值閾值算法(Singular Value Thresholding,SVT).由于NNM 平等的對每個奇異值進(jìn)行收縮,忽略了矩陣奇異值的差異.在此基礎(chǔ)上,Dong 提出了加權(quán)核范數(shù)最小化(Weighted Nuclear Norm Minimization,WNNM)[6]算法,針對不同大小的奇異值,設(shè)定不同的權(quán)值進(jìn)行收縮,得到了更好的秩最小化問題的解,實(shí)驗(yàn)表明較NNM 其保留更多原始圖像的邊緣信息,去噪性能更高.研究表明,Schatten-p 范數(shù)是一種比核范數(shù)更加逼近低秩矩陣的范數(shù),Xie 等[7]在2016年提出加權(quán)Schatten-p 范數(shù)最小化(Weighted Schatten p-norm Minimization,WSNM)去噪問題,該方法比WNNM 具有更好的去噪性能.

CS-MRI 重建算法的研究一直是MRI 領(lǐng)域研究的熱門之一.迭代閾值算法(Iterative Shrinkage Thresholding Algorithm,ISTA)[8]是一種實(shí)現(xiàn)簡單、計(jì)算量小并且性能較為穩(wěn)定的算法.Zhang 等[9]利用小波變換實(shí)現(xiàn)圖像的稀疏化結(jié)合ISTA 算法最終實(shí)現(xiàn)了CS-MRI 重構(gòu).由于ISTA 算法具有收斂速度慢的缺點(diǎn),Beck 等[10]提出了快速迭代軟閾值算法(Fast Iterative Shrinkage Thresholding Algorithm,FISTA),從探究多步級的方法來加速算法收斂速度或者從優(yōu)化小波子帶寬參數(shù)兩個方向來實(shí)現(xiàn)優(yōu)化重構(gòu)問題.在ISTA 和FISTA 算法的基礎(chǔ)上,Huang 等[11]提出了快速混合分裂算法(Fast Composite Splitting Algorithm,FCSA).FCSA 算法結(jié)合了變量分裂和運(yùn)算分裂兩種算法,將原問題分解為最小化全局差分約束問題和最小化L1范數(shù)約束的兩個子問題.通過線性結(jié)合兩個子問題的解并加入加速步驟得到最終解.除此以外,交替方向乘子法(Alternating Direction Method of Multipliers,ADMM)[12]算法是以增廣拉格朗日算法為基礎(chǔ)的一種常用的圖像重構(gòu)算法.ADMM 在迭代過程中將目標(biāo)函數(shù)進(jìn)行變量分離,通過分而治之的思想將原問題分解為多個交替的子問題,使得原問題的求解轉(zhuǎn)換成較為簡單的子問題求解.近似消息傳遞算法(Approximate Message Passing,AMP)[13]是在消息傳遞算法的基礎(chǔ)上,用來求解圖像的結(jié)構(gòu)化稀疏約束的重構(gòu)模型.針對圖像重構(gòu)領(lǐng)域,AMP 通過結(jié)合圖像的結(jié)構(gòu)先驗(yàn)信息,在其迭代過程中濾出噪聲,迭代逼近原始圖像[14].文獻(xiàn)[15]構(gòu)建了基于去噪的近似消息傳遞算法(Denoising-based Approximate Message Passing,DAMP),利用經(jīng)典的降噪模型作為算法的先驗(yàn)知識實(shí)現(xiàn)DAMP 迭代過程中的濾波操作,最終重構(gòu)出圖像.

本文提出了一種基于自適應(yīng)低秩去噪的近似消息傳遞的磁共振圖像重構(gòu)算法.我們使用去噪的近似消息傳遞重構(gòu)算法,在重構(gòu)過程中使用自適應(yīng)的加權(quán)Schattenp 范數(shù)最小化去噪方法,根據(jù)估計(jì)的噪聲標(biāo)準(zhǔn)差來設(shè)定WSNM 的圖像塊大小和相似塊的數(shù)量.使用真實(shí)的全采樣磁共振圖像作為實(shí)驗(yàn)對象,實(shí)驗(yàn)表明本文提出的方法具有較高的重建性能并且能夠保留更多的圖像信息.

2 基于自適應(yīng)低秩去噪的磁共振圖像重構(gòu)算法

2.1 自適應(yīng)Schatten-p 范數(shù)最小化

矩陣的核范數(shù)定義如下:

矩陣的核范數(shù)等于矩陣進(jìn)行SVT 分解后矩陣奇異值的絕對值之和.核范數(shù)最小化(NNM)問題可以表示成如下所示:

以X的核范數(shù)作為約束條件,其中M是含有噪聲的觀測矩陣,核范數(shù)最小化問題就是從已知的M中恢復(fù)矩陣X.其中λ 為正參數(shù).的求解過程如算法1 所示.

算法1.核范數(shù)最小化求解(NNM)輸入:觀測矩陣M 和閾值.∧X輸出:.λ M=U ∑VT;//對M 使用SVT 分解;))S λ(∑iii=max(∑ i i-λ,0;∧X=US λ(∑)VT;

NNM 通過對奇異值進(jìn)行軟閾值操作,進(jìn)而保留了圖像的低秩結(jié)構(gòu)特性和信息.由于,NNM 中軟閾值操作設(shè)定每個奇異值的權(quán)重均相同,這種設(shè)置忽略了奇異值自身帶來的物理意義,并不是低秩求解問題的最優(yōu)解.在此基礎(chǔ)上提出加權(quán)核范數(shù)最小化(WNNM)方法.該方法給予大奇異值大的權(quán)值、小奇異值小權(quán)值,充分考慮了奇異值的物理意義,使得圖像中更多有用的信息被保留下來.其加權(quán)核范數(shù)定義為:

其中,wi表示奇異值權(quán)重.當(dāng)wi=1時,WNNM 問題變成了NNM 問題.因此,NNM 是WNNM 的一個特殊情況.基于WNNM 的圖像去噪模型如下所示:

給定含噪圖像M,利用圖像的非局部自相似性找到對應(yīng)的相似塊組Mj,利用式(5)估計(jì)對應(yīng)的原始清晰圖像塊Xj;其中是圖像的噪聲方差.待求得整幅圖像的相似塊組后,對圖像塊組進(jìn)行聚合就能夠求得去噪后的清晰圖像.

除核范數(shù)以外,Schatten-p 范數(shù)是秩函數(shù)的一種非凸近似函數(shù).文獻(xiàn)[16,17]提出了使用Schatten-p 范數(shù)來實(shí)施低秩正則化的方法,定義矩陣奇異值的Lp范數(shù)為Schatten-p 范數(shù),其中 0 ＜p≤1.從理論上講,Schatten-p范數(shù)能夠更準(zhǔn)確地恢復(fù)信號.在WNNM 的基礎(chǔ)上,提出了加權(quán)Schatten-p 范數(shù)最小化(WSNM).

矩陣X∈Rn×m的加權(quán)Schatten-p 范數(shù)定義如下:

其中,w=[w1,w2,···,wmin{n,m}]表示權(quán)重向量,式(6)的p次方為:

則利用加權(quán)Schatten-p 范數(shù)最小化進(jìn)行圖像去噪描述成如下數(shù)學(xué)模型:

權(quán)重向量wk是WSNM 去噪的關(guān)鍵步驟之一,求解如下:

其中,n表示圖像的相似塊數(shù)量,ε是一個非常小的正數(shù),用來防止除數(shù)為0,c=表示矩陣的第k奇異值.文獻(xiàn)[7]證明WSNM 可以等價地轉(zhuǎn)化為獨(dú)立的非凸權(quán)重Lp范數(shù)子問題:

其中,δi為觀測矩陣M的第i個奇異值.式(10)的求解可以通過廣義軟閾值(Generalized Soft-Thresholding,GST)[18]算法獲得全局最優(yōu)解,最終求得:

本文使用WSNM 去噪時,根據(jù)當(dāng)前的噪聲標(biāo)準(zhǔn)差大小自適應(yīng)地設(shè)置WSNM 的圖像塊大小以及相似塊的個數(shù),該設(shè)置對磁共振圖像重構(gòu)結(jié)果具有微弱的影響,不同噪聲標(biāo)準(zhǔn)差下WSNM 去噪的圖像塊大小以及相似塊個數(shù)設(shè)置如表1所示.

表1 不同噪聲標(biāo)準(zhǔn)差下WSNM 去噪的圖像塊大小及相似塊的設(shè)置

2.2 基于自適應(yīng)WSNM 的磁共振圖像重構(gòu)算法

近似消息傳遞算法是基于迭代軟閾值的信號重建技術(shù).AMP 迭代過程中,殘差的更新使用了Onsager 校正項(xiàng)來求解,Onsager 校正項(xiàng)使用蒙特卡洛方法近似求解.而去噪近似消息傳遞算法(DAMP)是基于噪聲去除的AMP 算法,以圖像去噪來實(shí)現(xiàn)AMP 的濾波操作.本文將改進(jìn)的自適應(yīng)WSNM 方法作為AMP 迭代濾波的去噪函數(shù)以實(shí)現(xiàn)圖像壓縮感知重構(gòu).用Dσt表示整個去噪算法,借助獨(dú)立同分布向量來估計(jì)去噪函數(shù)的散度.MRI 的觀察模型如下:

其中,是 向量形式的重構(gòu)圖像;Fu表示欠采樣傅里葉操作算子;y是觀測數(shù)據(jù),N表示觀測數(shù)據(jù)y的長度大小;η是觀察噪聲.

基于CS-MRI 的DAMP 模型如下:

其中,b是滿足b～N(0,I) 的隨機(jī)向量,τ是一個數(shù)值很小的常數(shù),文獻(xiàn)[15]證明了Onsager 校正項(xiàng)的存在能夠保證DAMP 算法的快速收斂.在本章中,我們使用WSNM 算法實(shí)現(xiàn)迭代過程中的圖像去噪處理Dσ,我們稱此算法為基于自適應(yīng)WSNM 去噪的DAMP 磁共振圖像重構(gòu)算法(簡稱WSNM-AMP-MRI).算法2 展示了基于自適應(yīng)WSNM 的近似消息傳遞算法.

算法2.WSNM-AMP-MRI 重構(gòu)算法Fux0=FTu y輸入:觀測數(shù)據(jù);欠采樣操作算子 ;初始 ;y輸出:重構(gòu)磁共振圖像.x for then rk=xk-1+FTu zk-1σk=‖zk-1‖2k≤K■Nxk=Dσk(rk)WSNM ok=zk-1D′σk(rk)WSNM/Mzk=y-Fuxk+okend k=k+1

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

為了驗(yàn)證本章提出的WSNM-AMP-MRI 的重構(gòu)性能,我們選用兩張MR 圖像進(jìn)行模擬實(shí)驗(yàn),它們分別是心臟(Heart)和大腦(Brain)磁共振圖像,其中大腦數(shù)據(jù)來源于文獻(xiàn)[19],心臟數(shù)據(jù)來源于文獻(xiàn)[20].使用了笛卡爾采樣和偽徑向采樣兩種模式對原始圖像欠采樣.圖1分別展示了使用的實(shí)驗(yàn)磁共振圖像以及采樣mask (采樣率為25%的條件下).

圖1 測試數(shù)據(jù)及采樣mask

對比實(shí)驗(yàn)有平移不變離散小波變換(Shift-Invariant Discrete Wavelet,SIDWT)[21],基于方向性小波變換(Patchbased Directional Wavelets,PBDW)[22]和基于自相似性的非局部算子(Patch-based Nonlocal Operator,PANO)[23]的磁共振圖像重構(gòu)算法,與WSNM-AMP-MRI 方法進(jìn)行比較.對比實(shí)驗(yàn)中,SIDWT、PBDW 和PANO 重構(gòu)方法將設(shè)置最優(yōu)參數(shù)使得重構(gòu)結(jié)果最佳.實(shí)驗(yàn)中,WSNMAMP-MRI 重建算法的參數(shù)設(shè)置如下:正則化參數(shù)λ=1e-6,p=0.7.同時,使用峰值信噪比(Peak Signal to Noise Ratio,PSNR) 和相對L2范數(shù)誤差(RelativeL2Norm Error,RLNE)作為客觀評價指標(biāo)衡量不同算法的重構(gòu)效果.兩種客觀評價指標(biāo)的定義請參考式(19)和式(20).每種方法的計(jì)算結(jié)果是通過重復(fù)10 次相同的重建過程后計(jì)算結(jié)果的平均值.

其中,x是原始清晰圖像,y是重構(gòu)圖像,圖像的大小是M×N,μx和 μy表示x和y的均值;σx和 σy分別是x和y的方差,σxy是x和y的協(xié)方差;C1和C2是常數(shù),作用是保持穩(wěn)定性.

3.1 實(shí)驗(yàn)結(jié)果

圖2和圖3分別展示了在不同的采樣因子和采樣模式下,相應(yīng)算法重構(gòu)出的磁共振圖像的細(xì)節(jié)放大圖以及誤差圖.圖2展示了采樣因子為4 時的偽徑向采樣下的心臟重構(gòu)結(jié)果展示圖.從細(xì)節(jié)放大圖可以看到,本章提出的WSNM-AMP-MRI 算法能夠保留更多的局部細(xì)節(jié)信息.從放大了10 倍的重構(gòu)誤差圖可以看出,本章提出的算法重構(gòu)誤差明顯小于對比的3 種重構(gòu)算法的誤差.圖3展示采樣因子為4 的笛卡爾采樣下的大腦(Brain)重構(gòu)結(jié)果,可以看出笛卡爾采樣重構(gòu)結(jié)果產(chǎn)生了平移偽影,而為徑向采樣重構(gòu)圖像不存在這個問題.

表2和表3記錄了在欠采樣因子為4 和6 的兩種情況下,不同算法重構(gòu)算法的重構(gòu)PSNR和RLNE值.從表中數(shù)據(jù)可以看出,隨著采樣因子的增加,重構(gòu)的PSNR和RLNE分別在降低和增加,說明圖像重構(gòu)的質(zhì)量在降低.越高的采樣因子重構(gòu)出的圖像質(zhì)量相對越差.同時可以看出,偽徑向采樣下圖像的PSNR和RLNE值優(yōu)于笛卡爾采樣下相應(yīng)的重構(gòu)指標(biāo).總之,在笛卡爾和偽徑向兩種采樣模式下.本文提出的WSNMAMP-MRI 可獲得更高的PSNR值和更小的RLNE值,表明了本文提出的方法能夠得到更好的重建效果,證明了算法的有效性.

圖2 采樣因子為4 的偽徑向采樣下的心臟重構(gòu)結(jié)果

3.2 不同加速因子下圖像重構(gòu)性能比較

為研究不同采樣率對重構(gòu)結(jié)果的影響,本實(shí)驗(yàn)研究了在加速因子為2、4、6、8、10 時使用偽徑向采樣和笛卡爾采樣模式下的重構(gòu)結(jié)果.本次實(shí)驗(yàn)選用心臟(heart)磁共振圖像作為實(shí)驗(yàn)對象.圖4和圖5展示了不同加速因子下圖像重構(gòu)的PSNR和RLNE走向圖.從圖中可以看出 當(dāng)加速因子增加時,重構(gòu)的PSNR值在減小,RLNE值在增大,說明圖像的重構(gòu)質(zhì)量在降低.即采樣因子越高,圖像的重構(gòu)效果越差.但是在相同采樣因子的條件下,本文所提算法的PSNR值始終高于對比實(shí)驗(yàn)的PSNR值,同時RLNE值低于對比實(shí)驗(yàn)的RLNE值,說明本方法在不同的采樣模式以及加速因子下,其重構(gòu)質(zhì)量均高于相同條件下對比實(shí)驗(yàn)的重構(gòu)質(zhì)量.同時可以發(fā)現(xiàn),相同采樣率下,笛卡爾采樣的重構(gòu)質(zhì)量要低于偽徑向采樣模式下的圖像重建質(zhì)量.

圖3 采樣因子為4 的笛卡爾采樣下的大腦重構(gòu)結(jié)果

表2 采樣因子為4 時不同算法下的圖像重構(gòu)質(zhì)量(PSNR/RLNE)

表3 采樣因子為6 時不同算法下的圖像重構(gòu)質(zhì)量(PSNR/RLNE)

3.3 Schatten-p 范數(shù)取值對算法重構(gòu)的影響

本節(jié)研究p范數(shù)的取值對磁共振圖像重構(gòu)的影響大小.圖6展示了在偽徑向和笛卡爾采樣兩種模式下,p范數(shù)取不同值時,腦部圖像重構(gòu)的PSNR值.從圖中可以看出,兩種采樣模式下,p的取值對圖像重構(gòu)的影響并不相同.整體上看,p范數(shù)的取值對兩種采樣模式下圖像的重構(gòu)質(zhì)量影響并不是非常大的.相比較笛卡爾采樣模式下的重構(gòu)質(zhì)量而言,p范數(shù)對偽徑向采樣下圖像的重構(gòu)質(zhì)量影響更大.我們知道當(dāng)p=1 時,WSNM就轉(zhuǎn)化為WNNM 問題.從圖中可以看出,基于Schatten-p范數(shù)的加權(quán)范數(shù)最小化去噪的圖像重構(gòu)PSNR 值高于WNNM,說明WSNM 是一種更加優(yōu)秀的逼近低秩矩陣的方法.實(shí)驗(yàn)中,當(dāng)p的取值范圍在[0.5,0.7]時能夠取得相對較好的重構(gòu)性能.

3.4 噪聲大小對重構(gòu)結(jié)果的影響

本小節(jié)通過比較幾種算法在含有噪聲的測量數(shù)據(jù)下的重構(gòu)性能,來驗(yàn)證本章提出的算法的魯棒性.表4和表5給出在采樣率為25%的條件下,測量噪聲分別為10 dB 和20 dB 時磁共振圖像重構(gòu)的PSNR值和RLNE值.從表4中可以看出,在測量噪聲為10 dB 時,本章提出的算法具有更高的PSNR和更低的RLNE值.同理,表5說明在噪聲為20 dB 條件下本章算法仍然具有更高的重構(gòu)質(zhì)量.通過將表4、表5與表2中的重構(gòu)質(zhì)量指標(biāo)值對比,有測量噪聲的干擾條件下圖像重構(gòu)質(zhì)量與沒有測量噪聲時的重構(gòu)質(zhì)量相比普遍降低;并且隨著噪聲的增加,磁共振圖像重構(gòu)的質(zhì)量在降低.因此可以得出結(jié)論,無論在有測量噪聲還是沒有測量噪聲的條件下,本章提出的方法均優(yōu)于對比的3 種圖像重構(gòu)算法.

圖4 不同加速因子的偽徑向采樣下的心臟圖像重構(gòu)質(zhì)量

圖5 不同加速因子的笛卡爾采樣下的心臟圖像重構(gòu)質(zhì)量

圖6 范數(shù) p 對腦部圖像重構(gòu)的影響

表4 采樣因子為4、測量噪聲為10 dB 時的4 種算法重建圖像的PSNR(dB)/RLNE

表5 采樣因子為4、測量噪聲為20 dB 時的4 種算法重建圖像的PSNR(dB)/RLNE

4 結(jié)論與展望

本文提出了一種基于自適應(yīng)低秩去噪的近似消息傳遞磁共振圖像重構(gòu)算法.將加權(quán)Schatten-p 范數(shù)最小化作為圖像的低秩約束,再結(jié)合DAMP 的特性,根據(jù)每次迭代過程中的噪聲標(biāo)準(zhǔn)差設(shè)定自適應(yīng)的圖像塊大小以及相似塊的個數(shù),最終重構(gòu)出質(zhì)量相對較好的磁共振圖像.實(shí)驗(yàn)結(jié)果表明,與近幾年提出的幾種重建算法比較,本文的方法可以獲得更高的峰值信噪比和相對L2范數(shù)誤差.與此同時,在測量數(shù)據(jù)加噪聲的條件下觀察算法的重構(gòu)性能,驗(yàn)證了該算法具有魯棒性.

在本文的研究中,參數(shù)調(diào)優(yōu)是手動進(jìn)行的,在未來的工作中設(shè)計(jì)一個能夠自適應(yīng)調(diào)整參數(shù)的算法顯得迫在眉睫.同時,本文提出的方法是針對二維磁共振圖像的重建,而實(shí)際生活中也會采集動態(tài)磁共振成像進(jìn)行病灶診斷.在考慮到算法重構(gòu)時間和性能的基礎(chǔ)上,如何將本文提出的二維重建方法使用到三維動態(tài)磁共振成像是接下來的研究目標(biāo)之一.