基于小波樹和互補(bǔ)分解的CS-MRI重建算法

裴 穎，朱金秀，楊語晨，吳文霞

(1.河海大學(xué) 物聯(lián)網(wǎng)工程學(xué)院，江蘇 常州 213022;2.南通河海大學(xué) 海洋與近海工程研究院，江蘇 南通 226300)

0 引 言

核磁共振成像(magnetic resonance imaging，MRI)是醫(yī)學(xué)成像，應(yīng)用廣泛。目前MRI應(yīng)用的關(guān)鍵在于快速成像。Nyquist采樣定理需兩倍帶寬，不符合實(shí)際應(yīng)用[1]，現(xiàn)常用方法是在壓縮感知(compressed sensing，CS)[2-3]框架下，從欠采樣k空間中重建MRI數(shù)據(jù)，能有效減少采樣時間，達(dá)到快速成像的目的。

關(guān)于CS-MRI重建算法的研究有很多。例如，Lusting等[4]利用全變分(total variation，TV)和小波構(gòu)建目標(biāo)函數(shù)，提出共軛梯度算法(CG)，但重建時間有待提高；FISTA算法[5]通過計算更合適的起始點(diǎn)以加快收斂速度；TVCMRI[6]、RecPF[7]和FCSA[8]算法利用算子分割、變量分割思想求解聯(lián)合正則算子，提高重建速度和質(zhì)量；文獻(xiàn)[9]利用圖像在頻域的特性優(yōu)化測量矩陣并提出迭代加權(quán)算法，提高了重建精度；文獻(xiàn)[2,10]利用MR圖像低秩特性進(jìn)行奇異值分解，但過程過于復(fù)雜；文獻(xiàn)[11-12]提出結(jié)構(gòu)稀疏理論，圖像不僅在小波域有稀疏性，其小波稀疏系數(shù)也有特定的四叉樹結(jié)構(gòu)，在此基礎(chǔ)上文獻(xiàn)[13]提出YALL1算法，利用小波樹結(jié)構(gòu)稀疏代替小波稀疏構(gòu)建目標(biāo)函數(shù)，以提高圖像稀疏度；文獻(xiàn)[14]提出WaTMRI算法，聯(lián)合小波樹結(jié)構(gòu)稀疏和小波稀疏分別擁有的結(jié)構(gòu)稀疏和稀疏性，聯(lián)合改善圖像質(zhì)量；Park等[15]提出互補(bǔ)分解，將完整圖像分為平滑和殘差兩個分量，僅將平滑分量用于TV，殘差分量用于1范數(shù)，以解決全變分導(dǎo)致的細(xì)節(jié)過平滑問題；文獻(xiàn)[16]利用貪婪算法提高重建速度，但需要確定圖像稀疏度，缺乏實(shí)際性；文獻(xiàn)[17]基于p范數(shù)構(gòu)建重建算法，計算較為復(fù)雜。

重建算法常用小波稀疏和TV構(gòu)建目標(biāo)函數(shù)，但TV會在濾除噪聲的同時平滑圖像的邊緣紋理信息。因此，文中提出一種基于小波樹和互補(bǔ)分解的CS-MRI(MRI based on complementary dual decomposition and wavelet tree sparsity，DualWaTMRI)重建算法,利用互補(bǔ)分解抑制TV造成的細(xì)節(jié)過平滑問題，小波樹結(jié)構(gòu)稀疏用于補(bǔ)充小波稀疏的先驗(yàn)信息，并與同類重建算法進(jìn)行了對比。

1 基于小波樹和互補(bǔ)分解的CS-MRI重建算法

基于小波樹和互補(bǔ)分解的CS-MRI重建算法系統(tǒng)框架如圖1所示。首先欠采樣輸入數(shù)據(jù)得到測量值b，將其傅里葉逆變換預(yù)處理生成初始重構(gòu)圖像，再對其互補(bǔ)分解，生成初始?xì)埐罘至亢推交至?，分別進(jìn)行TV和1范數(shù)處理，同時利用結(jié)構(gòu)化稀疏理論[13]，對整幅圖像進(jìn)行小波樹結(jié)構(gòu)稀疏處理，最后與最小二乘擬合共同構(gòu)建目標(biāo)函數(shù)。同時，提出DualWaTMRI重構(gòu)算法進(jìn)行求解，分別得到重構(gòu)的殘差分量和平滑分量，由兩分量疊加組成最終的重建圖像。

圖1 算法系統(tǒng)框架

1.1 基于小波樹和互補(bǔ)分解的目標(biāo)函數(shù)

1.1.1 小波樹

1.1.2 互補(bǔ)分解

考慮到對整幅圖像進(jìn)行TV處理會同時平滑噪聲和圖像的細(xì)節(jié)信息，引入互補(bǔ)分解模型[13]，將MR圖像看成是平滑分量和殘差分量的疊加：

x=L+S

(1)

其中，x表示MR圖像；L表示灰度值緩慢變化的平滑分量，一般代表圖像背景信息；S表示灰度值快速變化的殘差分量，一般代表圖像中前景信息。

為避免TV過平滑圖像，僅對L進(jìn)行TV處理，S進(jìn)行小波稀疏處理。文中利用邊緣檢測算子Sobel[18]對圖像進(jìn)行互補(bǔ)分解，令初始?xì)埐罘至縎0為分解得到的邊緣信息，初始平滑分量L0為0。Sobel算子通過核與像素值做卷積和運(yùn)算，選取合適的閾值提取邊緣信息，閾值計算方法為T=μ*255，其中μ是閾值比。

1.1.3 目標(biāo)函數(shù)

綜合互補(bǔ)分解和小波樹結(jié)構(gòu)稀疏的優(yōu)勢，提出基于小波樹和互補(bǔ)分解的CS-MRI重建算法，構(gòu)建目標(biāo)函數(shù)如下：

α‖L‖TV+β(‖φS‖1+

(2)

其中,L，S分別是平滑和殘差分量；A是測量矩陣；b是測量值；φ是小波稀疏基；ζ是小波樹形分組的集合；g表示其中一組；α和β是調(diào)優(yōu)參數(shù)，用于平衡所占比重。式2中第一項(xiàng)最小二乘擬合項(xiàng)用于保證重建圖像的準(zhǔn)確度，第二項(xiàng)平滑分量的TV項(xiàng)用于抑制噪聲，避免圖像過平滑，第三項(xiàng)殘差分量的1范數(shù)用于保證小波稀疏，第四項(xiàng)小波樹結(jié)構(gòu)稀疏項(xiàng)用于保證圖像的結(jié)構(gòu)稀疏性。這四個正則項(xiàng)相互補(bǔ)充，增加圖像先驗(yàn)信息，提高算法的魯棒性。

1.2 DualWaTMRI重構(gòu)算法

針對式2中的L和S兩個未知量，利用交替最小化方法構(gòu)建DualWaTMRI重構(gòu)算法進(jìn)行求解，將目標(biāo)函數(shù)分為L子問題(此時S為固定值)和S子問題(此時L為固定值)，同時將這兩個子問題交替迭代，最后重建圖像x即是求解得到的L和S值之和。

1.2.1L子問題求解

求解L子問題時，假設(shè)S為固定值，則最小化求解時可省略常數(shù)項(xiàng)。L子問題公式表示如下：

(3)

為簡便求解小波樹結(jié)構(gòu)稀疏正則項(xiàng)，令zp=GφL，其中G是ζ的二值矩陣[15]，式3改為：

(4)

其中，λ是調(diào)優(yōu)參數(shù)；ζ中共有s組g，gi表示第i組g，即第i個具有父子依賴性的小波系數(shù)組，i=1,2,…,s。式4中亦包含zp和L兩個未知量，采用交替最小化進(jìn)一步細(xì)化：

(1)L子問題中zp未知量按gi分組求解：

β‖zpgi‖2)

(5)

上式可采用分組軟閾值法求解，令rpi=(GφL)gi：

i=1,2,…,s

(6)

由上式求解zpgi值，將其按位置線性組合生成zp。

(2)L子問題中L未知量求解公式如下：

其中，zp值由式6解得。式7可采用FISTA[5]算法求解得到L值。

1.2.2S子問題求解

與L子問題類似，求解未知量S時假設(shè)L為固定值，求解公式如下：

(8)

同樣的，令zr=GφS，上式更新為：

(9)

(1)未知量zr可通過分組軟閾值法得到，令rri=(GφS)gi，求解公式為：

i=1,2,…,s

(10)

(2)未知量S求解公式如下：

(11)

上式亦可采用FISTA方法求解。

1.2.3 重建算法總結(jié)

針對目標(biāo)函數(shù)，利用交替最小化將函數(shù)分為兩個子問題求解，結(jié)合上述重建步驟構(gòu)建算法DualWaTMRI，該算法重建步驟總結(jié)如下：

輸入：最大迭代次數(shù)N，k空間測量值b，閾值比μ。

步驟1：初始化。

步驟1.1：對b進(jìn)行傅里葉逆變換預(yù)處理，生成初始圖像x0；

步驟1.2：對x0取合適的閾值μ互補(bǔ)分解得到初始?xì)埐罘至縎0；

步驟1.3：令初始平滑分量L0=0，n為迭代次數(shù)，初始化n=1，ρ=1/Lf。

步驟2：交替最小化求解。

步驟2.1：求解L子問題：給定Sn-1，結(jié)合Ln-1求解式5得到zp值，將zp值代入式7求解Ln；

步驟2.2：求解S子問題：給定Ln，結(jié)合Sn-1求解式10得到zr值，將zr值代入式11求解Sn；

步驟2.3：令n=n+1，判斷n是否等于N，若不等于，回到步驟2.1，否則跳到步驟2.4；

步驟2.4：重建圖像是平滑分量和殘差分量的疊加：x=L+S(此時L=Ln，S=Sn)。

2 實(shí)驗(yàn)結(jié)果與分析

對Chest，Heart，Brain和Shoulder等四幅標(biāo)準(zhǔn)MRI圖像進(jìn)行實(shí)驗(yàn)。將DualWaTMRI與CG[4]、TVCMRI[6]、FCSA[8]和WaTMRI[15]算法進(jìn)行比較。根據(jù)經(jīng)驗(yàn)設(shè)α=0.02，β=3.5e-2，λ=β/2，Lf=1。同時選小波為稀疏基，偽高斯[6]為采樣模板，測量矩陣為采樣模板下的部分傅里葉矩陣。結(jié)果如表1所示。

表1 四幅MR圖像進(jìn)行重建時的PSNR結(jié)果 dB

2.1 實(shí)驗(yàn)參數(shù)：閾值比μ的選取

對DualWaTMRI算法進(jìn)行仿真，閾值比μ取值從0到1，結(jié)果如圖2所示。

從圖2可看出，在μ取值為0.05～0.20時，SNR值急劇上升，在μ大于0.20時，SNR值稍有減少，最后趨于平緩?？梢缘贸觯贒ualWaTMRI算法互補(bǔ)分解模型中μ取值0.16。

圖2 仿真結(jié)果對比

2.2 重構(gòu)性能的比較

將MR圖像分別通過CG、TVCMRI、FCSA、WaTMRI和DualWaTMRI等算法進(jìn)行仿真，結(jié)果分別如圖3、4和表1所示。

圖3 四幅MR圖像的SNR比較

圖4 四幅MR圖像的t比較

圖3和圖4給出了四幅MR圖像通過不同算法得到的仿真結(jié)果折線圖(SNR和t)。表1給出了SNR具體值。從圖3和表1中可以看出，四幅MR圖像通過DualWaTMRI算法得到的信噪比最高，比WaTMRI算法分別高出約1.69 dB、0.63 dB、0.60 dB和1.70 dB，更優(yōu)于其他算法。但圖4中DualWaTMRI算法上升趨勢最慢，表明重建時間長，這是由于目標(biāo)函數(shù)中有兩個未知量，需要交替求解，因此文中算法通過犧牲少量計算時間，達(dá)到增加信噪比的目的。

綜上，DualWaTMRI算法與同類CS-MRI算法相比，雖然由于先驗(yàn)信息豐富導(dǎo)致圖像的重建時間有所增加，但重建圖像的質(zhì)量有一定改善。

3 結(jié)束語

為改善TV造成的圖像細(xì)節(jié)模糊問題，引入互補(bǔ)分解模型，與小波樹結(jié)構(gòu)稀疏相結(jié)合，提出一種基于小波樹和互補(bǔ)分解的CS-MRI重建算法。利用互補(bǔ)分解模型將MR圖像分為平滑和殘差兩個分量，僅將平滑分量用于TV，殘差分量用于1范數(shù)，利用兩個分量各自的特性，很好地保留了圖像的細(xì)節(jié)信息，同時利用樹結(jié)構(gòu)稀疏特性，豐富圖像的先驗(yàn)信息，使得在同等數(shù)目的掃描數(shù)據(jù)下，能得到更好的重建圖像，或者只需更少的掃描數(shù)據(jù)，就能獲得同等質(zhì)量的重建圖像。實(shí)驗(yàn)結(jié)果表明，與現(xiàn)有的基于整幅圖像全變分、基于小波樹的算法相比，該算法以部分運(yùn)算時間為代價很好地改善了MR圖像的重建質(zhì)量。