非參數(shù)貝葉斯分類字典學(xué)習(xí)的MRI重建方法

朱 路，曹賽男，劉 松，劉媛媛，李康康

(華東交通大學(xué) 信息工程學(xué)院，江西 南昌 330013)

0 引 言

為減少磁共振成像(magnetic resonance imaging，MRI)的數(shù)據(jù)采集時間，相關(guān)研究者提出了許多圖像重構(gòu)方法。其中，基于數(shù)據(jù)或圖像本身學(xué)習(xí)出的字典訓(xùn)練方法在MRI重構(gòu)中得到廣泛應(yīng)用[1-5]。文獻(xiàn)[3]中基于自適應(yīng)學(xué)習(xí)稀疏變換的DLMRI算法，通過數(shù)據(jù)本身學(xué)習(xí)獲得字典，使信號具有更稀疏的表示。文獻(xiàn)[4]提到的GradDLRec算法采用K-SVD訓(xùn)練字典，需要預(yù)先設(shè)定好字典的大小、信號的稀疏度、噪聲方差等參數(shù)。Zonoobi等[5]利用圖像塊之間的相似性進(jìn)行圖像塊聚類，對動態(tài)MR圖像進(jìn)行分類字典學(xué)習(xí)，提高了字典的結(jié)構(gòu)表示能力。然而，該方法采用K-means執(zhí)行圖像塊分組，需憑經(jīng)驗預(yù)設(shè)聚類個數(shù)，缺乏自適應(yīng)性。Zhou等[6]提出基于BPFA(beta process factor analysis)的字典學(xué)習(xí)算法，解決了字典學(xué)習(xí)中的參數(shù)選擇問題。由于Zhou等通過貝塔-伯努利先驗來構(gòu)造稀疏，未利用圖像的稀疏先驗和非局部信息，所以重構(gòu)質(zhì)量仍可進(jìn)一步提升。

本文利用梯度域稀疏和非局部相似性，提出基于非參數(shù)貝葉斯分類字典學(xué)習(xí)的MRI重構(gòu)方法。通過在梯度域中進(jìn)行字典學(xué)習(xí)，使字典的稀疏表示能力更強(qiáng)。采用無限高斯混合模型，為圖像塊自動確定聚類個數(shù)。引入基于BPFA的非參數(shù)字典學(xué)習(xí)方法，分別為每一類圖像塊訓(xùn)練字典，降低對參數(shù)選擇的依賴性，提高圖像的結(jié)構(gòu)化表達(dá)。

1 基于梯度域和分類字典的MRI重構(gòu)模型

1.1 壓縮感知與字典學(xué)習(xí)

壓縮感知理論指出，若信號在特定域上具有稀疏性或可壓縮性，則通過少量測量即可重建整個原始信號。壓縮感知理論主要包括信號的稀疏性、測量矩陣和非線性重建方法3個方面。MR圖像本身不具有稀疏性，但可利用小波、差分、字典等變換方式來表示MR圖像的稀疏性。其中，通過構(gòu)造自適應(yīng)字典的變換方法得到了廣泛應(yīng)用[5-7]。對于一幅MR圖像，可通過構(gòu)造一個字典D和相應(yīng)的稀疏表示系數(shù)矩陣來近似表示。字典D∈Rn×K， 由K列dj向量構(gòu)成，每一列dj稱為一個原子。其目標(biāo)函數(shù)為

(1)

1.2 梯度域MR圖像的非局部相似先驗信息

MR圖像具有變換域稀疏、局部平滑、非局部相關(guān)等多種特性，充分利用MR圖像的不同先驗信息受到研究者的廣泛重視。近年來的研究結(jié)果表明[8,9]，數(shù)據(jù)越稀疏，自適應(yīng)學(xué)習(xí)字典的稀疏表示能力越強(qiáng)。因此，可考慮先進(jìn)行圖像的稀疏變換，在稀疏的變換域中再進(jìn)行自適應(yīng)字典學(xué)習(xí)。同時，圖像的非局部相似性在圖像去噪、重建領(lǐng)域中得到了極為廣泛地應(yīng)用[5,10,11]。在一幅MR圖像中取一小塊，則可以在該圖像的其它多處區(qū)域找到與此塊相似的小塊，這種相似特性即為非局部相似性。由于差分變換是線性的，因此在梯度域中同樣具有這種相似性。圖1顯示了MR腦圖像的非局部相似性，將原始像素做差分后，在橫向差分域和縱向差分域中依然具有非局部相似特性。

圖1 MR原始圖像和梯度域圖像中的非局部相似性

1.3 梯度域和分類字典的MRI重構(gòu)模型

基于上述兩種先驗，本文提出基于梯度域和分類字典學(xué)習(xí)的MRI重構(gòu)模型，如圖2所示。

圖2 基于梯度域和分類字典學(xué)習(xí)的MRI重構(gòu)框架

首先將原始圖像轉(zhuǎn)換到梯度域，使圖像的表示更稀疏；同時利用圖像的非局部相似性將圖像塊聚類，分別為每一類圖像塊進(jìn)行字典訓(xùn)練，具有保持圖像局部結(jié)構(gòu)和邊緣紋理的優(yōu)點(diǎn)。該重構(gòu)模型的數(shù)學(xué)表述為

(2)

2 基于非參數(shù)貝葉斯分類字典學(xué)習(xí)的MRI重構(gòu)方法

2.1 非局部相似塊聚類

在式(2)中，字典Dj根據(jù)每一類圖像塊的結(jié)構(gòu)特點(diǎn)訓(xùn)練而來。在對圖像塊聚類時，有別于Zonoobi等直接對圖像塊做非局部相似性處理和利用K-means方法確定聚類個數(shù)，本文對梯度域圖像塊做均值削減并利用無限高斯混合模型自動聚類。本文對梯度域圖像塊的分類過程[12]如下：

(3)

第n個PG變化向量的條件分布為

(4)

當(dāng)混合系數(shù)向量π服從狄利克雷分布時，式(4)則為無限高斯混合模型。所有變量的聯(lián)合概率分布為

(5)

此聯(lián)合概率的后驗參數(shù)，可由EM(expectation maximization)算法[13]估計。推斷出所有的后驗參數(shù)后，從而獲得所有PG變化向量的聚類個數(shù)Ng。

2.2 非參數(shù)貝葉斯字典學(xué)習(xí)與優(yōu)化方法

(6)

(7)

這里采用交替迭代方法(alternating direction method，ADM)[15]對式(6)進(jìn)行變量分離，得到P1、P2、P3這3個子問題，并對3個子問題分別求解

(8)

(2)P2求解：對差分圖像變量ω進(jìn)行更新，同樣采用最小二乘法進(jìn)行求解，并令β=P， 結(jié)果為

(9)

(3)P3求解：更新每一組圖像塊的字典Dj和對應(yīng)的稀疏表示系數(shù)矩陣Γj。大部分字典學(xué)習(xí)方法需要預(yù)先設(shè)置字典大小、稀疏度、噪聲方差等參數(shù)，當(dāng)參數(shù)設(shè)置與圖像實(shí)際偏離較大時，這些字典學(xué)習(xí)方法的重構(gòu)性能將大大下降。因此，本文提出基于BPFA的非參數(shù)分類字典學(xué)習(xí)方法，該模型表示如下

(10)

其中，dk表示第j組字典矩陣D的第k個原子，εl表示噪聲。Ip表示p×p大小的單位矩陣。稀疏系數(shù)αl被分解成兩個向量的哈達(dá)瑪乘積，sl表示圖像塊稀疏表示的權(quán)重向量；zl表示圖像塊稀疏表示的二值向量；⊙為哈達(dá)瑪向量乘積，表示逐個元素對應(yīng)相乘。字典中各個原子的使用情況由zlk∈{0,1}決定，表示此原子是否參與對信號的表示。每個zlk均服從參數(shù)為πk的Bernoulli分布，參數(shù)πk又由Beta分布生成。這種概率分布設(shè)置，使得當(dāng)K為一個足夠大的值時，有相當(dāng)數(shù)量的原子都未參與對信號的表示。由于該算法的稀疏先驗用Truncated Bernoulli-Beta的先驗假設(shè)來代替，因此不需要預(yù)先設(shè)置稀疏度、字典大小等參數(shù)，是一種非參數(shù)字典學(xué)習(xí)的方法。

為求解后驗概率，將BPFA模型與優(yōu)化方法相結(jié)合[16,17]，通過吉布斯采樣(Gibbs sampling)更新BPFA的后驗參數(shù)，其更新方程為

(11)

算法1：基于BPFA的分類字典學(xué)習(xí)算法

輸入：梯度域圖像塊Plωk，超參數(shù)a，b，c，d，e，f。

輸出：更新的梯度域圖像塊Plωk+1。

從N個梯度圖像塊Plωk中提取N個PG，計算其均值μi和PG變化向量。通過貝葉斯變分推斷估計無限高斯混合模型中的參數(shù)，對每一類PG變化向量執(zhí)行以下操作：

(1)構(gòu)造字典矩陣Dj=[dj1,…,djK]，djk～N(0,P-1Ip)；

(3)權(quán)重向量和噪聲向量的精度值采樣γs～Gamma(c,d),γε～Gamma(e,f)；

(4)對第j組第i個圖像塊進(jìn)行以下操作：

構(gòu)造稀疏系數(shù)：αi=zi⊙si；

根據(jù)P1、P2和P3子問題的求解過程，本文重構(gòu)方法的具體步驟歸納為算法2。

算法2：基于非參數(shù)貝葉斯分類字典學(xué)習(xí)的MRI重構(gòu)算法

輸入：下采樣K空間數(shù)據(jù)v，循環(huán)次數(shù)K，k=1。

輸出：重構(gòu)圖像Ik+1。

重復(fù)以下步驟直到終止條件：

(1)根據(jù)式(8)，更新圖像I；

(2)根據(jù)式(9)，更新差分圖像變量ω；

(3)根據(jù)式(10)，通過Gibbs采樣更新D和αl；

(5)令k=k+1， 若k

3 實(shí)驗及分析

為驗證本文提出方法的重建性能，仿真實(shí)驗選用4張典型的MR測試圖像，如圖3所示。選用當(dāng)前兩種經(jīng)典的MRI重建算法進(jìn)行性能比較，分別為DLMRI算法和GradDLRec算法。

圖3 測試圖像

在所有的實(shí)驗中，均采用二維隨機(jī)壓縮采樣方式，并設(shè)置25%、38%、46%、57%、64%這5種從低到高不等間隔的采樣率。采用峰值信噪比(PSNR)和高頻誤差范數(shù)(HFEN)評價算法的重建質(zhì)量。其中，PSNR從像素之間的灰度差別整體進(jìn)行評價，HFEN利用高頻信息從圖像的邊緣和細(xì)節(jié)特征進(jìn)行評價。PSNR值越高且HFEN值越低，說明重建質(zhì)量越好。

首先對MR圖像在無噪聲干擾情況下的重建性能進(jìn)行測試。然后，對MR圖像在噪聲標(biāo)準(zhǔn)方差為15、20、25、30、35情況下的重建性能進(jìn)行測試，用于比較3種重建算法的魯棒性能。實(shí)驗仿真軟件為MATLAB，實(shí)現(xiàn)平臺是1.70 GHz Intel?64 GB RAM的計算機(jī)。

3.1 參數(shù)設(shè)置

在測試實(shí)驗中，本文所提算法的參數(shù)設(shè)置如下：抽取的圖像塊大小為P1/2=8，圖像塊的重疊抽取步長為r=1。非局部相似塊的搜索窗口大小為w×w=31×31。圖像塊的個數(shù)N為MR圖像的大小。N個PG變化的初始聚類個數(shù)為Ng=32。當(dāng)某類PG變化的個數(shù)大于104時，字典的原子大小設(shè)為K=512，反之，字典的原子大小設(shè)為K=256；基于BPFA的分類字典學(xué)習(xí)算法采用10次循環(huán)迭代。將超參數(shù)a，b，c，d，e，f經(jīng)驗性地設(shè)置為a=b=1，c=d=e=f=10-6。對比算法DLMRI算法和GradDLRec算法的代碼均取自原作者主頁。其中，DLMRI算法的參數(shù)設(shè)置如其論文中所建議的：字典大小K=36，圖像塊大小P1/2=6，圖像塊的重疊抽取步長為r=1，β=36，λ=140，圖像塊的個數(shù)為200×P，稀疏閾值T0=5；GradDLRec算法除增加參數(shù)λ1=60，v2=3外，其它設(shè)置均與DLMRI算法保持一致。

3.2 稀疏采樣性能比較

首先對無噪聲情況下的4幅MR圖像進(jìn)行測試。圖4是Brain1圖像在采樣率為25%時的重建圖像(第一行)及重建殘差圖像(第二行)?？梢杂^察，DLMRI算法和GradDLRec算法的重構(gòu)殘差圖中的外部邊緣和內(nèi)部紋理仍有部分保留，左側(cè)和右側(cè)的居中位置處均呈現(xiàn)出白色短線條。其中，DLMRI算法呈現(xiàn)的短線條更明顯，且外部邊緣的上側(cè)也出現(xiàn)了白色短線條。而在本文算法的重構(gòu)殘差圖中，外部邊緣的任意位置都未出現(xiàn)邊緣短線。這說明，使用DLMRI算法和GradDLRec算法重建出的圖像丟失了Brain1圖像的部分邊緣和紋理信息；而使用本文算法重構(gòu)的圖像誤差最小，邊緣結(jié)構(gòu)和局部紋理均保持得很好。這是因為本文算法利用了圖像的稀疏性和非局部相似性兩種先驗信息，增加了圖像的結(jié)構(gòu)性表達(dá)，故獲得了細(xì)節(jié)上更加豐富的圖像。

圖4 Brain1圖像在3種算法下的重建效果

表1詳細(xì)列出了4幅測試圖像在不同重建算法下的重構(gòu)結(jié)果。表中加粗的地方表示3種重建算法的重構(gòu)結(jié)果中最大的PSNR值和最小的HFEN值。由表1可以看出，在各幅圖像中和各個采樣率下，DLMRI算法的PSNR值最低、HFEN值最高；GradDLRec算法的PSNR值、HFEN值優(yōu)于DLMRI算法；本文算法的PSNR值、HFEN值均優(yōu)于前兩種算法。相比GradDLRec算法，平均每幅圖像、每個采樣率，本文所提算法的PSNR值平均提高了2.9 dB，HFEN值平均降低了0.45。由此可知，本文所提算法與DLMRI算法和GradDLRec算法相比，重構(gòu)質(zhì)量有顯著提高。

表1 不同算法重建圖像的PSNR/HFEN值

3.3 有噪聲稀疏采樣性能比較

為測試在有噪聲情況下的稀疏采樣性能，在K空間加入不同量級的高斯白噪聲。圖5是在不同的噪聲標(biāo)準(zhǔn)方差下，3種算法在采樣率為25%時對Shoulder1圖像的重建性能曲線。從圖中可以看出，無論在何種噪聲量級中，本文算法的PSNR值都高于GradDLRec算法和DLMRI算法，特別是在高噪聲情況下，PSNR值有顯著提高。

圖5 Shoulder1圖像在不同量級噪聲情況下的重建性能比較

圖6給出了當(dāng)噪聲標(biāo)準(zhǔn)方差為20、采樣率為25%時，3種方法對Shoulder1圖像的重建效果?？梢钥闯觯杀疚乃惴ㄖ貥?gòu)的圖像更加清晰，對Shoulder1圖像中的線條和紋理結(jié)構(gòu)保持的最好。

圖6 Shoulder1圖像在3種算法下的重建圖像

例如，在圖中的白色方框位置，由DLMRI算法和GradDLRec算法重構(gòu)出的圖像由于受到噪聲的干擾，出現(xiàn)了一些龜裂線條。而由本文算法重構(gòu)出的圖像，方框位置內(nèi)的線條和細(xì)小結(jié)構(gòu)均保持得很好。這是因為GradDLRec算法和DLMRI算法都采用KSVD訓(xùn)練字典，這種方法需要經(jīng)驗性地設(shè)置字典大小、稀疏度等參數(shù)。而在實(shí)際應(yīng)用中，預(yù)先設(shè)置的參數(shù)往往會因具體每幅MR圖像的差異而偏離真實(shí)值?；贐PFA模型的分類字典學(xué)習(xí)方法直接跳過這些設(shè)置環(huán)節(jié)，根據(jù)每幅MR圖像數(shù)據(jù)自動推斷參數(shù)，所以學(xué)習(xí)出的字典具有更強(qiáng)的表示能力。

3.4 收斂性分析

為測試3種重建算法的收斂性能，觀察3種重建算法在每次迭代時的PSNR值和HFEN值。圖7是Brain2圖像在采樣率為46%時的PSNR值隨迭代次數(shù)變化的曲線。從中可以看出，3種算法都具有良好的收斂性能。其中，DLMRI算法的收斂速度最快，但達(dá)到收斂時的PSNR值最低，為45.45 dB；GradDLRec算法的收斂速度相對較慢，達(dá)到收斂時的PSNR值相對較高，為49.21 dB；本文算法的收斂速度相對較慢一些，但達(dá)到收斂時的PSNR值最高，為51.39 dB，相對提高了5.84 dB和2.17 dB。

圖7 3種算法的PSNR值與迭代次數(shù)的關(guān)系曲線

圖8是Brain2圖像在采樣率為46%時的HFEN值隨迭代次數(shù)變化的曲線，可以看出，3種算法都有良好的收斂性能。這是因為3種重建算法都是在字典學(xué)習(xí)和稀疏表示之間交替迭代求解，使其非負(fù)的目標(biāo)函數(shù)單調(diào)遞減。本文算法由于在梯度域中執(zhí)行圖像塊的無限高混合聚類和利用BPFA進(jìn)行字典學(xué)習(xí)，導(dǎo)致收斂速度有所降低，但重構(gòu)效果最好。

圖8 3種算法的HFEN值與迭代次數(shù)的關(guān)系曲線

3.5 計算復(fù)雜度分析

假設(shè)圖像塊個數(shù)為N，圖像塊大小為P，字典大小為K，則稀疏閾值為T0=0.15P。在每次迭代中，DLMRI算法在稀疏表示和圖像重構(gòu)兩個階段的計算復(fù)雜度分別為O(PKNT0) 和O(Nlog(N))； GradDLRec算法在梯度變換、稀疏表示和圖像重構(gòu)3個階段的計算復(fù)雜度分別為O(Nlog(N))、 O(PKNT0)、 O(Nlog(N))。 在本文算法中，假設(shè)每一類圖像塊的個數(shù)為Nj，則在稀疏表示階段的計算復(fù)雜度為O(P(N1+N2+…+Ng))。 本文算法由于利用了非局部相似性，故又增加非局部相似塊的搜索計算時間。雖然本文算法的計算復(fù)雜度有所增加，但能獲得質(zhì)量更好的MR圖像。在獲得相同質(zhì)量的MR圖像時，可以采用更低的采樣率，這對醫(yī)學(xué)成像而言具有重要意義。見表1，當(dāng)要求Brain1重構(gòu)圖像的PSNR值達(dá)到47 dB以上時，GradDLRec算法和DLMRI算法分別需要38%和46%以上的采樣率，而本文算法只需要使用低于25%的采樣率即可滿足要求。

4 結(jié)束語

本文利用差分變換稀疏和非局部相似性這兩種先驗信息，提出基于非參數(shù)貝葉斯分類字典學(xué)習(xí)的MR圖像重建模型。實(shí)驗結(jié)果相比經(jīng)典的MRI重建方法而言，本文所提算法在相同的采樣率下的重建質(zhì)量更優(yōu)，尤其在含有噪聲情況下圖像質(zhì)量有顯著提高，驗證了該重構(gòu)方法的有效性。由于采用了分類訓(xùn)練字典的策略以及對非局部相似塊的搜索計算，本文算法運(yùn)行時間較長，后續(xù)工作將在分類字典學(xué)習(xí)和非局部相似塊搜索方面做進(jìn)一步優(yōu)化。