基于Moreau-包絡(luò)的近似平滑迭代磁共振圖像重建算法

劉曉暉，路利軍，馮前進(jìn)，陳武凡

(南方醫(yī)科大學(xué) 生物醫(yī)學(xué)工程學(xué)院，廣州 510515)(*通信作者電子郵箱chenwf@fimmu.com)

0 引言

磁共振成像(Magnetic Resonance Imaging, MRI)因其成像中不使用電離輻射、可以更靈活地選擇成像參數(shù)、對(duì)軟組織結(jié)構(gòu)及代謝功能的完美表達(dá)以及對(duì)許多特定參數(shù)更加敏感等優(yōu)點(diǎn)成為臨床應(yīng)用中性能最卓越的醫(yī)學(xué)成像方式之一[1-2]。然而，MRI特殊的成像機(jī)制導(dǎo)致的相對(duì)較慢的成像速度阻礙了MRI進(jìn)一步的發(fā)展和應(yīng)用：磁共振(MR)圖像對(duì)于運(yùn)動(dòng)非常敏感，因而成像時(shí)間長會(huì)產(chǎn)生更多的運(yùn)動(dòng)偽影，這種情況常常導(dǎo)致需要對(duì)病人進(jìn)行鎮(zhèn)靜或麻醉；較低的空間分辨率限制了MR對(duì)某些和呼吸相關(guān)的部位的成像，如腹部成像和心臟成像；較長的掃描時(shí)間在增加成像成本的同時(shí)也會(huì)限制MRI的適用病人范圍[3]。出于這些原因，許多研究人員在尋求在減少采樣數(shù)據(jù)的同時(shí)又不降低圖像質(zhì)量的成像方法。

從20世紀(jì)70年代MR還未被廣泛地應(yīng)用于臨床以來，MRI的加速問題就被廣泛的探討和研究。MRI中成像數(shù)據(jù)為Fourier變換空間，通常稱為k空間。MRI的采樣時(shí)間主要由k空間的分段采樣模式所決定，在給定一個(gè)采樣序列以后，k空間中的分段數(shù)就決定了MR圖像的信號(hào)采集時(shí)間?；谶@個(gè)原因，MRI的加速通常是采用減少k空間中的分段數(shù)，即對(duì)k空間進(jìn)行降采樣。這種方法利用了k空間的冗余特性，即k空間中的各個(gè)點(diǎn)是由位于成像物體不同空間位置的信號(hào)疊加產(chǎn)生的[4]，主要代表技術(shù)是并行MRI(parallel MRI, pMRI)技術(shù)[5-8]。pMRI技術(shù)利用多組線圈同時(shí)接收MR信號(hào)，利用線圈的空間敏感度差異來代替部分梯度編碼，從而減少梯度編碼步數(shù)，節(jié)省掃描時(shí)間。pMRI無需修改成像序列，在保持掃描視野和空間分辨率不變的情況下減少成像時(shí)間，因此為現(xiàn)代商業(yè)MRI系統(tǒng)廣泛采用，成為目前的常規(guī)掃描方法。然而，由于欠采樣導(dǎo)致的混疊偽影，以及重建算法對(duì)噪聲的放大效果，pMRI方法在實(shí)際臨床應(yīng)用中只能使用小加速倍數(shù)以保證理想的重建質(zhì)量。

近年來在信號(hào)處理領(lǐng)域發(fā)展起來的壓縮感知(Compressed Sensing, CS)[9-12]理論為快速M(fèi)RI技術(shù)提供了一種新的可能。CS理論認(rèn)為，一個(gè)本身稀疏的或者可壓縮的未知信號(hào)，選擇一組合適的感知信號(hào)(即感知矩陣)，在感知信號(hào)個(gè)數(shù)遠(yuǎn)小于原始未知信號(hào)維數(shù)的情況下，通過求解一個(gè)非線性最優(yōu)化問題，可以精確恢復(fù)出原始信號(hào)。在MRI中：1) MR圖像被認(rèn)為是稀疏(如血管MR圖像、三維MR圖像相鄰圖像層之間以及動(dòng)態(tài)MR圖像在時(shí)域上的圖像序列)或是可壓縮的(一幅MR圖像在某個(gè)變換域中的系數(shù)矩陣具有可壓縮性)；2) MRI的感知信號(hào)空間是Fourier空間，其成像方式?jīng)Q定了可以有選擇的掃描k空間中的一部分?jǐn)?shù)據(jù)。上述兩個(gè)條件決定了CS可以應(yīng)用于快速M(fèi)RI中，這種快速M(fèi)RI技術(shù)一般稱為CS-MRI。

CS-MRI通常需要求解一個(gè)包含稀疏正則項(xiàng)的非線性凸優(yōu)化問題。為了更好地表達(dá)MR圖像的稀疏性，CS-MRI的數(shù)學(xué)模型中往往包含兩個(gè)稀疏正則項(xiàng)，即優(yōu)化問題中包含三個(gè)函數(shù)：數(shù)據(jù)保真項(xiàng)、1范數(shù)稀疏正則項(xiàng)和MR圖像的總變分(Total Variation, TV)正則項(xiàng)；而經(jīng)典的近似迭代類算法求解的是包含兩個(gè)函數(shù)(通常為數(shù)據(jù)保真項(xiàng)和一個(gè)正則項(xiàng))的最優(yōu)化問題，所以在使用近似點(diǎn)迭代類算法的CS-MRI方法中，通常的做法是將初始優(yōu)化問題根據(jù)正則項(xiàng)分解為兩個(gè)優(yōu)化子問題，分別對(duì)這兩個(gè)子問題用近似迭代算法進(jìn)行求解得到相應(yīng)的解，再對(duì)求得的兩個(gè)解取平均值得到最終解。這類算法由于每個(gè)子問題的輸入變量相同可以稱為多正則項(xiàng)優(yōu)化問題的并行求解方法。

近似迭代算法主要是針對(duì)優(yōu)化問題中存在的非平滑函數(shù)，如1范數(shù)正則項(xiàng)、核范數(shù)以及TV項(xiàng)等，這些非平滑函數(shù)不能直接求導(dǎo)，因此在優(yōu)化算法中需要求取它們的次梯度，在近似迭代類算法中往往使用近似算子來求解這個(gè)次梯度，近似算子本身即為一個(gè)最優(yōu)化問題的求解過程，通過對(duì)輸入信號(hào)x0迭代地施加近似算子即可求解原始優(yōu)化問題的最優(yōu)解。在前面提到的多正則項(xiàng)優(yōu)化問題的并行求解算法中，兩個(gè)子問題的解在同一個(gè)迭代步驟中是互相獨(dú)立的，即互相不以對(duì)方為輸入，這類算法在兩個(gè)子問題的求解過程中沒有任何優(yōu)化步驟。針對(duì)這個(gè)問題，本文提出一種基于Moreau包絡(luò)的近似平滑迭代算法(Proximal Smoothing Iterative Algorithm, PSIA)，該算法首先將原始優(yōu)化問題中的一個(gè)非平滑正則項(xiàng)的梯度用它的平滑函數(shù)的梯度近似表示，然后把數(shù)據(jù)保真項(xiàng)和平滑后的正則項(xiàng)的線性組合看作一個(gè)平滑可導(dǎo)函數(shù)，最后利用經(jīng)典的近似迭代算法對(duì)優(yōu)化問題進(jìn)行求解。

1 磁共振圖像稀疏重建相關(guān)算法

在文獻(xiàn)[4]中，對(duì)CS-MRI的理論基礎(chǔ)及數(shù)學(xué)模型都作了詳細(xì)的探討，并證明了在附加一個(gè)1范數(shù)小波變換約束項(xiàng)和TV約束項(xiàng)的前提下，可以從降采樣的k空間數(shù)據(jù)中準(zhǔn)確地重建出MR圖像，并提出如式(1)中的非約束凸優(yōu)化問題：

α‖Wx‖1+β‖x‖TV

(1)

其中：x∈RN×N為待重建的MR圖像，N×N為MR圖像尺寸；Fu:RN×N→CM×N(M≤N)，為降采樣Fourier變換，是傳感矩陣作用于Fourier矩陣的結(jié)果，即Fu=SF，其中S∈RM×N，F(xiàn):RN×N→CN×N分別為傳感矩陣和Fourier變換；b∈RM×N為k空間降采樣數(shù)據(jù)；W:RN×N→RN×N為小波變換；‖·‖TV:RN×N→R為MR圖像的TV值；α，β>0為正則項(xiàng)權(quán)重系數(shù)。式(1)中的幾個(gè)范數(shù)定義分別為式(2)～(4)：

(2)

(3)

(4)

式中D1、D2為分離梯度算子，分別表示對(duì)圖像x作水平、垂直兩個(gè)方向滿足Neumann邊界條件的前向有限差分：

(5)

并滿足下面的Neumann邊界條件：

(6)

式中n1、n2滿足n1,n2∈[1,N]。對(duì)于問題式(1)的求解，最主要的困難來自于第二項(xiàng)小波變換1范數(shù)和第三項(xiàng)TV項(xiàng)，這兩項(xiàng)都不能直接求取微分或?qū)?shù)，因而不能對(duì)問題式(1)直接使用梯度法求最小值。針對(duì)這個(gè)問題，文獻(xiàn)[4]提出了共軛梯度(Conjugate Gradient, CG)下降算法來求解問題式(1)，但是這種算法的缺點(diǎn)是收斂速度慢，對(duì)于一幅512×512像素大小的圖像，CG重建出視覺效果比較好的圖像所需要的時(shí)間是2 min～3 min，即使重建一幅256×256像素大小的圖像，也需要十幾秒，當(dāng)對(duì)三維物體成像時(shí)，所需時(shí)間就要以小時(shí)計(jì)，這顯然不能滿足實(shí)際應(yīng)用的需要。

為了提高重建速度，考慮到問題式(1)的正則項(xiàng)是兩個(gè)截然不同的范數(shù)表達(dá)式，數(shù)學(xué)優(yōu)化算法中的分離算法被用來解決式(1)這樣的多個(gè)正則項(xiàng)優(yōu)化問題。分離算法主要包括算子分離算法和變量分離算法。算子分離算法的核心問題是尋找最小化問題所對(duì)應(yīng)的最大單調(diào)算子為零的解，經(jīng)典的算子分離算法有forward-backward算法[13]、Douglas-Rachford算法[14]和Peaceman-Rachford算法[15]，它們都是求解目標(biāo)函數(shù)為兩個(gè)函數(shù)和的最小化問題的算法。與之相對(duì)應(yīng)的，變量分離算法則是使用變方向法(Alternating Direction Method, ADM)求解最小化問題的增廣拉格朗日形式。變量分離算法是Douglas-Rachford算法和Peaceman-Rachford算法的優(yōu)化變形，最早在20世紀(jì)70年代被提出用來解決偏微分方程(Partial Differential Equation, PDE)問題[16]。近來，在文獻(xiàn)[17]中提出了變量分離算法可以用來解1范數(shù)和TV正則優(yōu)化問題。文獻(xiàn)[18]和文獻(xiàn)[19]分別提出算子分離類方法TV1范數(shù)壓縮MRI(TVl1Compressed MR Imaging, TVCMRI)算法與變量分離類方法部分k空間重建算法(Reconstruction from Partial k space algorithm, RecPF)來求解問題式(1)。這兩種算法都比算法CG要快，重建一幅像素為256×256的磁共振圖像所需時(shí)間為2 s～3 s，在保證重建準(zhǔn)確度的前提下，大幅提高了重建速度。然而由于分離算法的特點(diǎn)，這兩種算法人為添加了一定數(shù)量的輔助變量，使得算法代碼的復(fù)雜性較大，從而算法的重現(xiàn)具有一定的難度；另外，這兩種算法中，對(duì)于兩個(gè)正則項(xiàng)求導(dǎo)過程中的近似處理都是簡(jiǎn)單地加上一個(gè)非負(fù)二次項(xiàng)，對(duì)于算法的簡(jiǎn)化沒有貢獻(xiàn)。

針對(duì)以上問題，近似迭代算法被引進(jìn)非平滑正則約束的優(yōu)化問題中。近似迭代算法的基本原理是給出非平滑函數(shù)的近似算子[20]，這個(gè)過程就是一個(gè)求取子問題的最小化的過程。對(duì)給定一個(gè)初始點(diǎn)的原問題式(1)，迭代地算出它的近似點(diǎn)，即可求解得出理論上的最優(yōu)解。這類算法和前述幾類算法不同的地方在于：直接對(duì)非平滑項(xiàng)的近似算子進(jìn)行操作，而不是對(duì)非平滑項(xiàng)求取梯度或是微分[21]。兩種最著名的近似迭代算法是迭代收縮閾值算法(Iterative Shrinkage Thresholding Algorithm, ISTA)[22]和它的加速形式——快速迭代收縮閾值算法(Fast Iterative Shrinkage Thresholding Algorithm, FISTA)[23]，這兩種算法被提出用來求解兩個(gè)函數(shù)和的最小化問題。在快速復(fù)合分離算法(Fast Composite Splitting Algorithm, FCSA)[24-25]中將問題式(1)分解為兩個(gè)子問題，即數(shù)據(jù)保真項(xiàng)與小波變換1范數(shù)的和的最小化問題以及數(shù)據(jù)保真項(xiàng)與TV項(xiàng)的和的最小化問題，分別對(duì)這兩個(gè)問題使用快速迭代收縮閾值算法，將求得的解取平均，即認(rèn)為所得的值為問題式(1)的解。在FCSA中，兩個(gè)最小化子問題的權(quán)重被簡(jiǎn)單地取值為0.5，然而在實(shí)際應(yīng)用中，這兩個(gè)子問題的權(quán)重則根據(jù)實(shí)際所要重建的磁共振圖像的內(nèi)容而有所不同，這就向算法的魯棒性提出了挑戰(zhàn)。本文提出的重建算法則避免了這個(gè)困擾，即不對(duì)問題式(1)進(jìn)行子問題分離，利用近似平滑方法的數(shù)學(xué)特性求解。

2 快速近似平滑迭代算法

2.1 相關(guān)快速重建算法

在給出本文算法之前，先討論目標(biāo)函數(shù)為兩個(gè)凸函數(shù)線性組合的優(yōu)化問題的求解。近似迭代類算法FISTA被用來求解式(15)所示兩個(gè)函數(shù)和的最小化問題：

Γ(x)=p(x)+q(x):x∈RN×N

(7)

其中:p(x)為一個(gè)凸的Lp平滑函數(shù);q(x)為一個(gè)不可微的連續(xù)凸函數(shù)。FISTA提出通過迭代地求解函數(shù)q(x)的近似點(diǎn)便可得到問題式(7)的最優(yōu)解，即：

xk+1=proxq(xk-η▽p(xk));

(8)

其中：xk+1、xk分別為第k+1與第k次迭代求得的解；▽p(x)為函數(shù)p(x)的導(dǎo)數(shù)；proxq(x)為函數(shù)q(x)在x處的近似值函數(shù)。

FISTA算法偽代碼為：

步驟1 參數(shù)及變量初始化：

參數(shù)：正則化參數(shù)α,β，迭代計(jì)數(shù)k=1，近似算子參數(shù)

η=1/Lp。

變量：圖像變量x0=0，y1=0，加速步長變量θ1=1。

步驟2 變量迭代：

Fork=1,2,… do

xk=proxq(yk-η▽p(yk))

yk+1=xk+((θk-1)/θk+1)(xk-xk-1)

End for

算法FCSA將問題式(1)分解為兩個(gè)子問題，分別是小波1范數(shù)正則約束問題和TV正則約束問題，利用FISTA對(duì)這兩個(gè)問題求解后所得解為x1、x2，再取平均即為問題式(1)的解。本文將提出另外一種利用算法FISTA求解問題式(1)的近似平滑迭代算法(PSIA)。

2.2 近似平滑迭代算法

本文需要求解的問題式(1)是由三個(gè)函數(shù)線性組合而成，將問題式(1)簡(jiǎn)化表達(dá)為：

Γ(x)=f(x)+g(x)+h(x):x∈RN×N

(9)

其中f(x)、g(x)和h(x)的表達(dá)式分別為：

(10)

如果要對(duì)問題式(9)使用算法FISTA進(jìn)行求解，則對(duì)照式(8)可得其變量近似迭代步驟為：

xk+1=proxh(xk-η▽(f+g)(xk));

(11)

式中的proxh(x)在已知x時(shí)很容易求得，而導(dǎo)數(shù)▽(f+g)(x)由于函數(shù)g(x)的非平滑性不能直接求取。為了得到▽(f+g)(x)的表達(dá)式，首先要對(duì)函數(shù)(f+g)(x)進(jìn)行平滑近似，由于函數(shù)f(x)是平滑函數(shù)，因此平滑近似主要是對(duì)函數(shù)g(x)進(jìn)行。在凸優(yōu)化理論[26]中，可以通過求解非平滑函數(shù)的Moreau包絡(luò)來對(duì)其進(jìn)行平滑近似，其定義為：

定義1 Moreau包絡(luò)[26]。對(duì)于一個(gè)非平滑不可微函數(shù)g(x):RN×N→R∪{+∞}，它的Moreau包絡(luò)為：

(12)

定理1[27]給定函數(shù)g(x)為一個(gè)下連續(xù)的非平滑凸函數(shù)且μ>0，則由定義1給出的Moreau包絡(luò)為一個(gè)1/μ平滑函數(shù)，且其導(dǎo)數(shù)可由g(x)的近似算子得到：

▽(gμ)(x)=(x-proxg(x))/μ

(13)

其中proxg(x,μ)的表達(dá)式為：

(14)

此時(shí)式(11)中的導(dǎo)數(shù)▽(f+g)(x)中g(shù)(x)部分可由式(13)給出。

本文提出求解優(yōu)化問題式(9)的近似迭代算法為：首先對(duì)非平滑函數(shù)g(x)使用Moreau包絡(luò)進(jìn)行平滑近似，則得到一個(gè)新的優(yōu)化問題：

Γμ(x)=f(x)+gμ(x)+h(x):x∈RN×N

(15)

在優(yōu)化問題式(15)中，由于函數(shù)f(x)、gμ(x)都是定義域?yàn)閤∈RN×N的平滑凸函數(shù)，因此可以將這兩個(gè)函數(shù)的和看成一個(gè)新的函數(shù)K(x)，則式(15)可以表達(dá)為：

K(x)=f(x)+gμ(x):x∈RN×N

(16)

對(duì)問題式(16)，使用算法FISTA進(jìn)行求解，其變量迭代步驟為式(11)。

本文算法偽代碼為：

步驟1 參數(shù)及變量初始化：

參數(shù)：正則化參數(shù)α,β，迭代計(jì)數(shù)k=1，近似算子參數(shù)μ=1/Lf,η=1/(Lf+1/μ)。

變量：所求取圖像變量x0=0，中間變量y1=0，加速步長變量：θ1=1。

步驟2 變量迭代求解：

Fork=1，2，… do

xk+1←proxh(yk-η▽(f+gμ)(yk));

yk+1=xk+((θk-1)/θk+1)(xk-xk-1)

End for

迭代步驟1中▽(f+gμ)的計(jì)算如下：

(17)

▽(gμ)(x)=(x-proxg(x))/μ=

(sign(Wx).×max{0,Wx-μ/Lf})/μ

(18)

(19)

由此可知上述算法PSIA中的所有計(jì)算單元都已知，給定一個(gè)初始點(diǎn)，即可迭代地計(jì)算出重建圖像。

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

3.1 實(shí)驗(yàn)設(shè)置

為了驗(yàn)證算法的有效性，將本文提出的近似平滑迭代算法(PSIA)與第1章中提到的4種經(jīng)典稀疏MR重建算法的重建結(jié)果進(jìn)行對(duì)比，這4種重建算法分別為：CG、TVCMRI、RecPF和FCSA，選取這些算法進(jìn)行對(duì)比是因?yàn)樗鼈兦蠼獾哪繕?biāo)函數(shù)和本文一致，且它們都是經(jīng)典的磁共振稀疏重建算法。本文算法與對(duì)比算法都是在Inter Core i5- 2400、3.10 GHz CPU的PC上進(jìn)行，編程環(huán)境使用Matlab 2015b，每個(gè)算法的迭代次數(shù)都是50。使用4種圖像數(shù)據(jù)進(jìn)行重建結(jié)果對(duì)比，分別是：經(jīng)典的Shepp-Logan體模圖像、人腦MR圖像、手部MR圖像，以及腦部血管MR圖像；為了方便對(duì)比，所有圖像尺寸都設(shè)為256×256像素大小，灰度取值范圍為0～255；實(shí)驗(yàn)降采樣率約為20%，采樣矩陣與實(shí)驗(yàn)圖像數(shù)據(jù)如圖1所示；在實(shí)驗(yàn)中，對(duì)采樣數(shù)據(jù)添加標(biāo)準(zhǔn)差為σ=0.01的高斯白噪聲以模擬采樣中產(chǎn)生的噪聲。

圖1 采樣矩陣與四種參考實(shí)驗(yàn)圖像

算法中的參數(shù)設(shè)置為：正則化參數(shù)α、β分別設(shè)為0.035和0.001，這個(gè)取值與4種對(duì)比算法中相同，以便于算法效果的比較；近似算子參數(shù)的設(shè)置根據(jù)Lipschitz常數(shù)分別設(shè)置為：1，0.5。

3.2 實(shí)驗(yàn)結(jié)果評(píng)估標(biāo)準(zhǔn)

本文主要給出兩大類實(shí)驗(yàn)結(jié)果評(píng)估標(biāo)準(zhǔn)，分別為定性評(píng)估即視覺效果評(píng)估，以及定量評(píng)估包括信噪比(Signal-to-Noise Ratio, SNR)、相對(duì)誤差(Relative Error, RE)和結(jié)構(gòu)相似度指數(shù)(Structural SIMilarity, SSIM)[28]，另外還給出了算法復(fù)雜度對(duì)比(CPU時(shí)間)。其中SNR、RE和SSIM的表達(dá)式為：

(20)

3.3 實(shí)驗(yàn)結(jié)果及討論

四種實(shí)驗(yàn)圖像在五種不同算法(CG、TVCMRI、RecPF、FCSA和PSIA)中的重建結(jié)果與重建結(jié)果對(duì)應(yīng)的差值圖像如圖2～5所示。在差值圖像的表示中，為了比較的公平性，在每個(gè)圖像的結(jié)果表示中，將重建結(jié)果最差的算法對(duì)應(yīng)的差值圖像灰度值作為灰度直方圖范圍。

圖2 Shepp-Logan體模的重建結(jié)果

圖3 腦部MR圖像的重建結(jié)果

圖4 手部MR圖像的重建結(jié)果

從圖2～5中可以看出：本文提出的重建算法PSIA在重建的視覺效果上比之前的幾種經(jīng)典算法都要好，具體體現(xiàn)在差值圖像中的邊緣細(xì)節(jié)的重建以及對(duì)噪聲的消除上。算法CG的重建效果最差，不僅沒有很好地消除引進(jìn)的高斯白噪聲，在圖像的邊緣和解剖學(xué)結(jié)構(gòu)上都丟失了很多細(xì)節(jié)；算法TVCMRI和RecPF的重建結(jié)果在視覺效果上相近，且同時(shí)在細(xì)節(jié)的恢復(fù)和噪聲的消除上優(yōu)于CG；算法FCSA比TVCMRI和RecPF得到更好的重建結(jié)果，具體體現(xiàn)在消除了更多的噪聲，使重建圖像看起來更“干凈”；本文提出的算法PSIA比起FCSA在去噪方面表現(xiàn)更好，具體表現(xiàn)在差值圖中圖像背景部分的噪聲更小，在圖像細(xì)節(jié)的恢復(fù)上也更好，在視覺上與參考圖像以及真實(shí)MR圖像更接近。

四幅圖像在五種重建算法中的SNR與迭代次數(shù)的關(guān)系如圖6所示，用來比較算法的收斂性。

從圖6可以看出，對(duì)于不同的重建圖像，算法CG的效果是最差的，體現(xiàn)為SNR值最低；TVCMRI和RecPF的重建效果相近，這可以從算子分離和變量分離這兩類算法的關(guān)系推測(cè)出原因：變量分離類算法是算子分離類算法的優(yōu)化變形；FCSA和本文中提出的PSIA算法的重建效果比前三種算法要好，而本文提出的算法PSIA在四個(gè)重建圖像上的重建結(jié)果都好過FCSA，這驗(yàn)證了第1章中的討論，即：簡(jiǎn)單地將兩個(gè)子問題的權(quán)重設(shè)為0.5帶來魯棒性問題，導(dǎo)致了雖然都是使用快速迭代收縮閾值算法(FISTA)，PSIA比FCSA的重建效果要更好。

表1～4給出了四種圖像不同算法50次迭代后的重建性能，包括SNR、RE以及SSIM；同時(shí)為了比較算法復(fù)雜度，給出了每個(gè)算法運(yùn)行的CPU時(shí)間。SNR給出了重建圖像與參考圖像的信噪比，SNR的值越大，認(rèn)為重建算法的重建效果越好；RE給出了重建圖像和參考圖像間的差值，RE的值越小，認(rèn)為重建算法的重建效果越好；SSIM給出了重建圖像和參考圖像間的相似度測(cè)量，SSIM的值越大，認(rèn)為重建效果越好。

表1 Sheep-Logan體模的重建性能

從表1～4可以看出，在這三個(gè)評(píng)估參數(shù)方面，PSIA都是最好的，同時(shí)也可以看出算法CG的結(jié)果是最差的，算法TVCMRI和RecPF表現(xiàn)相近，F(xiàn)CSA優(yōu)于前三種算法，比PSIA表現(xiàn)差，這與圖2～5中的視覺評(píng)估以及圖6中的SNR隨迭代次數(shù)變化的結(jié)果一致。在CPU時(shí)間上，算法CG運(yùn)行時(shí)間最長，從TVCMRI開始算法運(yùn)行時(shí)間即加速到十倍數(shù)量級(jí)，RecPF的運(yùn)行時(shí)間和TVCMRI相當(dāng)，F(xiàn)CSA比以上兩個(gè)算法都要快，這得益于快速迭代收縮閾值算法的使用，PSIA在時(shí)間上和FCSA相當(dāng)，這也與兩者使用的核心算法都是同一種算法所得出的推斷一致。

圖5 腦部血管MR圖像的重建結(jié)果

圖6 四種MR圖像重建結(jié)果的SNR

算法SNR/dBRE/%SSIMCPU時(shí)間/sCG10.1826.910.594311.97TVCMRI14.4213.890.62551.83RecPF15.3713.230.64381.81FCSA18.099.100.78221.36PSIA19.617.640.85281.37

表3 手部MR圖像的重建性能

表4 腦部血管MR圖像的重建性能

4 結(jié)語

針對(duì)現(xiàn)有的MR圖像稀疏重建算法中不能在一個(gè)迭代計(jì)算步驟中同時(shí)對(duì)兩個(gè)正則約束項(xiàng)同時(shí)施加近似算子的問題，提出了一種基于Moreau-包絡(luò)的近似平滑迭代算法，該算法具有以下特點(diǎn)：

1)使用基于近似平滑迭代的算法，使算法結(jié)構(gòu)簡(jiǎn)單、易于重現(xiàn)；

2)與同樣基于近似平滑迭代的算法FCSA相比，解決了將原始最小化問題分解為兩個(gè)子問題時(shí)面臨的權(quán)重設(shè)置問題；

3)此算法可以應(yīng)用在不同正則約束稀疏化重建問題，在算法PSIA中，只要給出相對(duì)應(yīng)的正則項(xiàng)的近似算子，即可直接使用該算法。

對(duì)仿真體模以及真實(shí)MR圖像的實(shí)驗(yàn)結(jié)果表明，本文提出的近似迭代算法在重建精度以及算法復(fù)雜性方面都有很好的表現(xiàn)，重建效果比當(dāng)前經(jīng)典的CG、TVCMRI、RecPF和FCSA都要好，在算法復(fù)雜度方面和最快的算法FCSA相當(dāng)。

下一步的工作主要有兩方面：一是對(duì)算法中參數(shù)的優(yōu)化，在本文算法PSIA中，近似算子參數(shù)μ、η取的是固定的常數(shù)，在優(yōu)化算法文獻(xiàn)中，提出可以對(duì)這類參數(shù)進(jìn)行優(yōu)化從而得到更好的重建結(jié)果；二是提出新的正則約束項(xiàng)并利用算法PSIA，尤其是針對(duì)TV約束項(xiàng)，目前已有多種對(duì)于總變分的變形以及新的定義，將這些新的約束用在重建中預(yù)期將得到更好的重建結(jié)果。