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

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

(1.河海大學 物聯網工程學院，江蘇 常州 213022;2.南通河海大學 海洋與近海工程研究院，江蘇 南通 226300)

0 引 言

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

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

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

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

基于小波樹和互補分解的CS-MRI重建算法系統框架如圖1所示。首先欠采樣輸入數據得到測量值b，將其傅里葉逆變換預處理生成初始重構圖像，再對其互補分解，生成初始殘差分量和平滑分量，分別進行TV和1范數處理，同時利用結構化稀疏理論[13]，對整幅圖像進行小波樹結構稀疏處理，最后與最小二乘擬合共同構建目標函數。同時，提出DualWaTMRI重構算法進行求解，分別得到重構的殘差分量和平滑分量，由兩分量疊加組成最終的重建圖像。

圖1 算法系統框架

1.1 基于小波樹和互補分解的目標函數

1.1.1 小波樹

1.1.2 互補分解

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

x=L+S

(1)

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

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

1.1.3 目標函數

綜合互補分解和小波樹結構稀疏的優勢，提出基于小波樹和互補分解的CS-MRI重建算法，構建目標函數如下：

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

(2)

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

1.2 DualWaTMRI重構算法

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

1.2.1L子問題求解

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

(3)

為簡便求解小波樹結構稀疏正則項，令zp=GφL，其中G是ζ的二值矩陣[15]，式3改為：

(4)

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

(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時假設L為固定值，求解公式如下：

(8)

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

(9)

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

i=1,2,…,s

(10)

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

(11)

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

1.2.3 重建算法總結

針對目標函數，利用交替最小化將函數分為兩個子問題求解，結合上述重建步驟構建算法DualWaTMRI，該算法重建步驟總結如下：

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

步驟1：初始化。

步驟1.1：對b進行傅里葉逆變換預處理，生成初始圖像x0；

步驟1.2：對x0取合適的閾值μ互補分解得到初始殘差分量S0；

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

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

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

步驟2.2：求解S子問題：給定Ln，結合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 實驗結果與分析

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

表1 四幅MR圖像進行重建時的PSNR結果 dB

2.1 實驗參數：閾值比μ的選取

對DualWaTMRI算法進行仿真，閾值比μ取值從0到1，結果如圖2所示。

從圖2可看出，在μ取值為0.05～0.20時，SNR值急劇上升，在μ大于0.20時，SNR值稍有減少，最后趨于平緩。可以得出，在DualWaTMRI算法互補分解模型中μ取值0.16。

圖2 仿真結果對比

2.2 重構性能的比較

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

圖3 四幅MR圖像的SNR比較

圖4 四幅MR圖像的t比較

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

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

3 結束語

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