王文杰，喬志偉，牛 蕾，席雅睿

(山西大學(xué) 計(jì)算機(jī)與信息技術(shù)學(xué)院，太原030006)

0 引言

計(jì)算機(jī)斷層成像(Computed Tomography, CT)在醫(yī)學(xué)和工業(yè)生產(chǎn)中有非常重要和廣泛的應(yīng)用。絕大部分的商用CT機(jī)使用傳統(tǒng)的濾波反投影(Filtered Back Projection, FBP)算法來(lái)重建圖像已經(jīng)有將近30年的歷史[1]；然而，大劑量的輻射會(huì)增加患者患癌的風(fēng)險(xiǎn)，低劑量的掃描，即數(shù)據(jù)不完備的條件下，F(xiàn)BP算法難以重建出精確的圖像。為了解決這個(gè)問(wèn)題，有兩種常用思路：一是通過(guò)插值方法對(duì)不完備的投影數(shù)據(jù)進(jìn)行補(bǔ)足；二是通過(guò)在迭代過(guò)程中對(duì)重建圖像施加某些約束條件，即對(duì)模型增加正則化項(xiàng)的手段來(lái)解決。一些基于迭代的算法諸如代數(shù)迭代重建(Algebraic Reconstruction Technique, ART)[2]和期望最大化(Expectation Maximization, EM)算法[3]已經(jīng)被應(yīng)用于圖像重建中[4]?；趦?yōu)化的迭代算法可以利用壓縮感知、低秩矩陣等稀疏優(yōu)化技術(shù)對(duì)圖像進(jìn)行精確重建[5]。其中，全變分(Total Variation, TV)最小化約束的迭代重建算法已經(jīng)在各種CT中廣泛的應(yīng)用。

但是，TV模型的主要缺點(diǎn)是它不能恢復(fù)和保留圖像中的精細(xì)結(jié)構(gòu)和紋理，并傾向于產(chǎn)生階梯效應(yīng)(staircase artifacts)。因?yàn)閭鹘y(tǒng)的TV模型僅考慮圖像的局部像素信息，且假設(shè)圖像在單個(gè)對(duì)象內(nèi)是平滑的，只在不同對(duì)象的邊界處具有不連續(xù)的灰度跳躍[6]。

Buades等[7]利用兼顧了圖像中更多全局信息的非局部均值(NonLocal Means, NLM)算法來(lái)進(jìn)行圖像去噪。該算法通過(guò)對(duì)一個(gè)像素點(diǎn)周?chē)袼氐募訖?quán)平均值來(lái)恢復(fù)該像素，因此，相較于傳統(tǒng)的局部算子能更好地恢復(fù)和保留圖像的細(xì)節(jié)。Gilboa等[8]將傳統(tǒng)的變分框架和非局部算子進(jìn)行了有機(jī)的結(jié)合，提出了非局部變分(NonLocal TV, NLTV)框架。相較于局部變分框架,非局部變分框架在基于優(yōu)化的CT重建算法中能得到更加精確的結(jié)果是可以合理預(yù)期的。

Sidky等[9]提出了自適應(yīng)最速下降投影到凸集(Adaptive Steepest Descent Projection Onto Convex Sets, ASD-POCS)算法來(lái)解決圖像重建中基于約束的TV最小化問(wèn)題。其中，自適應(yīng)最速下降(ASD)即算法采用梯度下降法來(lái)使圖像的全變分最小化，同時(shí)步長(zhǎng)因子是自適應(yīng)的。投影到凸集(POCS)即為ART與正約束的合稱，ART保證重建結(jié)果滿足投影數(shù)據(jù)保真項(xiàng)的約束，正約束是因?yàn)槟J(rèn)圖像不存在負(fù)值像素點(diǎn)。ASD-POCS算法可以有效和魯棒地從不完備數(shù)據(jù)中重建出相當(dāng)精確的結(jié)果。它交替地運(yùn)行POCS過(guò)程和TV下降過(guò)程，從而使結(jié)果逐步地接近TV最小化解，同時(shí)滿足數(shù)據(jù)保真項(xiàng)和正則項(xiàng)的要求。

然而，ASD-POCS算法使用梯度下降算法來(lái)求解TV最小化，直接使用梯度下降法解決NLTV最小化的約束問(wèn)題是非常復(fù)雜和不實(shí)用的。Goldstein等[10]提出了分離Bregman (Split Bregman, SB)算法，可以有效求解很多的L1正則最優(yōu)化問(wèn)題。Zhang等[11]把這個(gè)方法推廣到了非局部最優(yōu)化問(wèn)題中。Lou等[12]將非局部算子應(yīng)用到了圖像的反卷積和重建中。實(shí)際上，ASD-POCS算法中的TV下降過(guò)程可以看作是對(duì)POCS步驟后得到的圖像的去噪過(guò)程，因此，本文把ASD-POCS算法作為基礎(chǔ)的重建框架，以期利用其參數(shù)自適應(yīng)的特點(diǎn)，同時(shí)，結(jié)合SB算法求解L1正則化問(wèn)題的高效性來(lái)解決模型中NLTV最小化的問(wèn)題，進(jìn)而提出自適應(yīng)非局部全變分最小化投影到凸集(Adaptive Non-local TV Minimization Projection Onto Convex Sets, ANLTVM-POCS)算法。

1 相關(guān)工作

1.1 CT模型

非約束的、基于優(yōu)化的CT數(shù)學(xué)模型可以表示為如下的形式：

(1)

其中：向量g是大小為M的離散投影數(shù)據(jù)；大小為N的列向量u表示離散圖像；A是一個(gè)大小為M*N的系統(tǒng)矩陣，這里表示平行束CT的二維Radon變換。系統(tǒng)矩陣A可以通過(guò)射線驅(qū)動(dòng)法(Siddon算法[13])和距離驅(qū)動(dòng)算法等來(lái)求取[14]。J(u)是正則項(xiàng)，λ是正則化參數(shù)。向量u*是最優(yōu)化模型的解，即被重建圖像。

那么，核心問(wèn)題便是如何選擇J(u)，在基于TV的重建型中，J(u)是梯度圖像的L1范數(shù)：

J(u)=‖u‖TV=

(2)

其中us,t表示圖像單個(gè)像素點(diǎn)。顯然，傳統(tǒng)的TV模型僅僅考慮了一個(gè)像素毗鄰的另外兩個(gè)像素，相較于非局部算法框架來(lái)說(shuō)，難以保留圖像的細(xì)節(jié)和紋理。

1.2 非局部變分框架

接下來(lái)，簡(jiǎn)要回顧在文獻(xiàn)[8]中介紹的非局部算子。令Ω?R2，x∈Ω，并且u(x):Ω→R是實(shí)函數(shù)。設(shè)w:Ω×Ω→R是關(guān)于參考圖像的非負(fù)對(duì)稱實(shí)函數(shù)，即：w(x,y)=w(y,x)。關(guān)于一像素點(diǎn)的非局部導(dǎo)數(shù)可以定義為：

(3)

則關(guān)于像素點(diǎn)x的非局部梯度算子▽NLu(x):Ω×Ω→Ω定義為其所有偏導(dǎo)數(shù)▽NLu(x,·)組成的向量：

(4)

對(duì)于向量u在點(diǎn)x∈Ω的非局部梯度L2范數(shù)：

(5)

所有像素點(diǎn)的非局部梯度L2范數(shù)也可以組成一幅圖像，即加權(quán)非局部梯度圖像。

對(duì)于向量v1=v1(x,y)，v2=v2(x,y)，向量?jī)?nèi)積定義為：

(6)

根據(jù)向量?jī)?nèi)積的定義及散度和梯度的共軛關(guān)系，

〈▽NLu,v〉=-〈u,divNLv〉

(7)

非局部散度divNLv(x):Ω×Ω→Ω定義為：

(8)

拉普拉斯算子可以寫(xiě)為：

(9)

權(quán)函數(shù)w(x,y)：

(10)

其中：Ga是標(biāo)準(zhǔn)差為a的高斯核函數(shù)，參數(shù)h是濾波參數(shù)。式(10)顯示了權(quán)函數(shù)只有在關(guān)于點(diǎn)x和點(diǎn)y的比較像素塊相似度大時(shí)具有較大的取值。

為了便于計(jì)算機(jī)進(jìn)行離散運(yùn)算，分別給出非局部梯度、非局部梯度L2范數(shù)、非局部散度和拉普拉斯算子的離散形式：

(11)

(12)

(13)

(14)

綜上所述，非局部全變分可以看作加權(quán)非局部梯度圖像的L1范數(shù)：

(15)

其中：wi, j是w(x,y)的離散形式，i和j是圖像的像素索引。

對(duì)于非局部框架，其核心問(wèn)題是對(duì)三個(gè)參數(shù)的選取：比較塊的尺寸(patch size)、搜索窗的尺寸(search window size)和式(10)中濾波參數(shù)h的取值；然而，這三個(gè)參數(shù)之間的耦合影響非常嚴(yán)重，作出一個(gè)恰當(dāng)?shù)倪x擇是一項(xiàng)比較艱巨的任務(wù)。本文在大量實(shí)驗(yàn)的基礎(chǔ)上，嘗試了幾種參數(shù)組合，然后給出一些可以得到相對(duì)精確結(jié)果的經(jīng)驗(yàn)參數(shù)選擇。

理論上，搜索窗的尺寸應(yīng)該是整幅圖像，然而實(shí)際上，一個(gè)過(guò)大的搜索窗會(huì)導(dǎo)致計(jì)算開(kāi)銷的幾何倍增加。許多文獻(xiàn)給出了從11×11到21×21之間的不同選擇，一個(gè)好的搜索窗尺寸應(yīng)該是在兼顧精確結(jié)果和計(jì)算開(kāi)銷的中間選擇。在本次工作中，發(fā)現(xiàn)11×11大小的搜索窗已經(jīng)足夠得出相對(duì)精確的結(jié)果且僅需要相對(duì)較少的計(jì)算開(kāi)銷。

在一些關(guān)于NLTV圖像去噪的研究[15]中提到，大的比較塊尺寸可以增強(qiáng)算法對(duì)于噪聲的魯棒性，較小的比較塊尺寸可以保留更多的低對(duì)比度圖像細(xì)節(jié)。在式(10)中，可以發(fā)現(xiàn)，結(jié)果像素點(diǎn)的取值還和鄰域像素點(diǎn)的高斯均值相關(guān)。對(duì)于不同的成像模型以及不同的噪聲類型，應(yīng)該根據(jù)情況作出恰當(dāng)?shù)倪x擇。在本文的實(shí)驗(yàn)中，大小為3×3的比較塊已經(jīng)可以得到較好的視覺(jué)結(jié)果圖像。

濾波參數(shù)h可以說(shuō)在非局部全變分模型中扮演了靈魂角色。該參數(shù)控制著比較塊的相似度大小，而相似度的準(zhǔn)確性又是非局部框架能得到精確結(jié)果的基本保證。盡管如此，很少有文獻(xiàn)提到在圖像重建中關(guān)于參數(shù)h應(yīng)該如何選擇。在許多關(guān)于非局部框架圖像去噪的文獻(xiàn)中參數(shù)h一般取圖像噪聲方差σ大小的若干倍。在文獻(xiàn)[12]中，對(duì)預(yù)處理圖像進(jìn)行小波變換；然后使用小波系數(shù)的絕對(duì)中位差來(lái)確定參數(shù)h的大小。許多圖像的噪聲估計(jì)算法可以在這里用來(lái)對(duì)參數(shù)h的選擇作出參考，本文不再詳細(xì)展開(kāi)。

1.3 分離Bregman算法

在引言中曾提到，ASD-POCS算法交替運(yùn)行POCS步驟和TV最小化步驟，而TV最小化可以看作是對(duì)POCS步驟后得到的圖像的去噪過(guò)程，因此，可以把TV最小化表述為如下最優(yōu)化模型：

(16)

其中u0是POCS步驟后得到的圖像。

標(biāo)準(zhǔn)ASD-POCS算法使用梯度下降法來(lái)最小化TV。對(duì)于NLTV，使用分離Bregman算法來(lái)求解這個(gè)最優(yōu)化問(wèn)題。引入一個(gè)輔助變量d:Ω×Ω→R來(lái)替換式(16)中的▽NLu，然后可以得出如下的約束最優(yōu)化問(wèn)題：

(17)

s.t.d=▽NLu

式(17)的約束問(wèn)題可以轉(zhuǎn)化為如下無(wú)約束形式：

(18)

通過(guò)引入分離Bregman算法框架，式(18)可以轉(zhuǎn)化為如下Bregman迭代形式：

(20)

式(19)可以分離為分別對(duì)u和d的遞歸迭代：

(21)

(22)

利用Euler-Lagrange公式，可以得出式(21)的最優(yōu)條件：

λ(u-u0)-γdivNL(dk-▽NLu-bk)=0

(23)

根據(jù)式(23)，使用Gauss-Seidel迭代算法來(lái)求取u的近似解(給出關(guān)于一個(gè)像素點(diǎn)在第k+1次迭代的離散形式，同d和b)：

(24)

關(guān)于子問(wèn)題d，可以通過(guò)收縮算子(shrinkage operators)[10]來(lái)求?。?/p>

(25)

最后，b可以通過(guò)直接迭代求解：

(26)

2 本文方法

本文方法主要受到文獻(xiàn)[9,11-12]的啟發(fā)。2008年Sidky等[9]提出的ASD-POCS算法，應(yīng)用正則化方法將目標(biāo)函數(shù)定義為優(yōu)化準(zhǔn)則及懲罰項(xiàng)二者之和的形式，以罰函數(shù)的形式增強(qiáng)解的穩(wěn)定性，以罰函數(shù)的形式增強(qiáng)了解的穩(wěn)定性[16]。2008年的方法改進(jìn)自他們2006年提出的重建算法，2006的重建算法首次利用醫(yī)學(xué)圖像梯度變換圖像具有稀疏性這一先驗(yàn)知識(shí)，結(jié)合凸集投影，無(wú)論是缺失角度投影數(shù)據(jù)(e.g. 0～45°)，還是大間隔角度投影數(shù)據(jù)(e.g. 0～180°間隔6°采樣)，相較于傳統(tǒng)的FBP解析類重建算法，ART和EM等迭代類算法都表現(xiàn)出了極其精確高效的稀疏重建能力。2008年的ASD-POCS算法考慮了投影數(shù)據(jù)中存在多種不一致性的問(wèn)題，新增了數(shù)據(jù)容差控制因子，大幅度提高了算法的魯棒性，改善了2006算法重建輪廓附近存在灰度漸變偽影的缺點(diǎn)，并且該算法在迭代過(guò)程中步長(zhǎng)因子自適應(yīng)的特點(diǎn)，也改善了其他迭代方法容易陷入局部最優(yōu)解的問(wèn)題，能在較少的迭代次數(shù)下得到更加精確的解。

本文方法主要基于ASD-POCS算法，并且針對(duì)ASD-POCS算法中TV模型的缺點(diǎn)，改用NLTV模型來(lái)約束最終解，提高算法保留重建結(jié)果中精細(xì)結(jié)構(gòu)和紋理的特性。

且針對(duì)NLTV傳統(tǒng)迭代方法(如梯度下降法)難以求解的問(wèn)題，利用分離Bregman方法求解L1正則化問(wèn)題的優(yōu)勢(shì)來(lái)求取NLTV最小化。

下面給出提出算法的偽代碼。

1)

Initialization:β=1.0;βred=0.995;α=0.2;αred=0.95;n=2;rmax=0.95;u=0;;

2)

Repeat main loop

3)

up=u;

4)

fori=1 toNd

5)

u=u+βAi(g-Ai·u)/(Ai·Ai)

6)

end for

7)

fori=1 toNi

8)

ifui<0 thenui=0

9)

end for

10)

dd=‖g-g0‖2

11)

dp=‖u-up‖2

12)

if (first iteration)

13)

dtvg=α*dp

14)

end if

15)

up=u;

16)

Initialization:u0=u;d0=0;b0=0;λ=1;γ=dtvg;

17)

fork=1 ton

18)

updateuk+1by equation (24)

19)

updatedk+1by equation (25)

20)

updatebk+1by equation (26)

21)

end for

22)

dg=‖u-up‖2

23)

ifdg>rmax*dpanddd>

24)

dtvg=dtvg*αred

25)

end if

26)

β=β*βred

27)

Until (stopping criteria)

行4)的Nd是投影的總個(gè)數(shù)，行7)的Ni是總像素個(gè)數(shù)。算法中有6個(gè)參數(shù)(β;βred;α;αred;n;rmax)控制整個(gè)算法的迭代過(guò)程，在文獻(xiàn)[9]中有詳細(xì)的敘述，本文不再展開(kāi)。u是算法的初始化圖像，在這里設(shè)定為零。行3)和15)的up是臨時(shí)變量，用來(lái)存儲(chǔ)算法當(dāng)前步驟的結(jié)果并參與步長(zhǎng)自適應(yīng)調(diào)整的計(jì)算。dtvg是算法最小化子問(wèn)題17)到21)行的步長(zhǎng)因子。算法的另一個(gè)重要參數(shù)是數(shù)據(jù)容差，是多種數(shù)據(jù)不一致的簡(jiǎn)化，包括簡(jiǎn)化模型帶來(lái)的誤差、投影數(shù)據(jù)中的噪聲以及X射線的散射等，這使得不可能總是重建出與數(shù)據(jù)完全一致的結(jié)果圖像[9]。參數(shù)確定了重建生成圖像的投影與實(shí)際投影數(shù)據(jù)的歐氏距離(L2范數(shù))。在結(jié)果部分討論了關(guān)于的選擇對(duì)于重建結(jié)果的影響。

把標(biāo)準(zhǔn)ASD-POCS算法的TV梯度下降過(guò)程替換為基于NLTV的分離Bregman迭代過(guò)程。這里有兩個(gè)重要參數(shù)(λ,γ)在分離Bregman的迭代過(guò)程中占據(jù)主導(dǎo)性地位，進(jìn)而影響最終重建結(jié)果的精度。在下一章詳細(xì)展開(kāi)討論。

3 實(shí)驗(yàn)結(jié)果與討論

為了評(píng)估和檢驗(yàn)本文的算法，設(shè)計(jì)了兩組實(shí)驗(yàn)。實(shí)驗(yàn)采用大小為128×128的Sheep-Logan和Forbild頭部仿真模體(http://www.imp.uni-erlangen.de/forbild/english/results/index.htm)，成像采樣旋轉(zhuǎn)中心為第64×64號(hào)像素點(diǎn)。探測(cè)器個(gè)數(shù)為128，單個(gè)探元大小為1。采樣角度范圍為0到180°，所有采樣均勻間隔。系統(tǒng)矩陣采用Siddon算法生成。圖像第一組實(shí)驗(yàn)是通過(guò)理想的不含噪聲的欠采樣數(shù)據(jù)來(lái)進(jìn)行運(yùn)算。使用兩組不同投影角度個(gè)數(shù)(20和30)生成的投影數(shù)據(jù)來(lái)檢驗(yàn)本文算法的稀疏重建能力。第二組實(shí)驗(yàn)，在模擬的投影數(shù)據(jù)中加入高斯噪聲(強(qiáng)度通過(guò)高斯噪聲方差衡量)來(lái)評(píng)估算法的抗噪能力。如圖1所示。

本文通過(guò)重建圖像剖線比較的方式來(lái)評(píng)估算法的重建特性。通過(guò)均方根誤差(Root Mean Squared Error, RMSE)來(lái)定量地評(píng)估重建結(jié)果的精度：

其中：u和f是分別是重建結(jié)果圖像和真值圖像；N是圖像的總像素個(gè)數(shù)。

圖1 仿真模體Fig. 1 Simulation phantoms

3.1 不含噪聲數(shù)據(jù)重建

至于算法的停止條件(stopping criteria)，可以當(dāng)計(jì)算結(jié)果達(dá)到一個(gè)特定的條件(如‖uk-uk+1‖/‖uk‖<10-3)時(shí)結(jié)束循環(huán)，但是為了更直觀地描述本文的算法，在本組實(shí)驗(yàn)中使算法直接循環(huán)500次，以便更全面地通過(guò)算法迭代過(guò)程的RMSE曲線來(lái)分析算法的收斂特性。

參數(shù)γ的恰當(dāng)選擇能使分離Bregman算法的子問(wèn)題u和d收斂的速度更快。在標(biāo)準(zhǔn)分離Bregman算法中，為了使算法更好地收斂，γ既不能過(guò)大也不能過(guò)小。文獻(xiàn)[10]中提到，當(dāng)γ=2λ時(shí)，分離Bregman算法通?？梢匀〉幂^好的收斂結(jié)果。在本文的算法中，參數(shù)γ可以說(shuō)是連接數(shù)據(jù)保真項(xiàng)和正則化項(xiàng)之間的橋梁。本文直接設(shè)定γ等于標(biāo)準(zhǔn)ASD-POCS算法框架中的自適應(yīng)步長(zhǎng)參數(shù)dtvg，以使ASD-POCS算法步長(zhǎng)參數(shù)自適應(yīng)的特點(diǎn)與分離Bregman迭代過(guò)程有機(jī)地結(jié)合起來(lái)。

對(duì)于參數(shù)λ，它影響了算法在保留更多原始投影數(shù)據(jù)內(nèi)容和去除噪聲之間的平衡。如果λ過(guò)大，數(shù)據(jù)保真項(xiàng)的權(quán)重會(huì)增加，換句話說(shuō)，結(jié)果可能將會(huì)過(guò)擬合給定的投影數(shù)據(jù)，相應(yīng)的數(shù)據(jù)中的噪聲也會(huì)更多地被保留；相反，較小的λ將會(huì)增大正則化項(xiàng)的權(quán)重，也就是說(shuō)，算法不僅會(huì)去除更多的噪聲信息，同時(shí)也會(huì)平滑模糊掉圖像中更多的細(xì)節(jié)。一個(gè)好的選擇應(yīng)該是在去噪和盡可能保留圖像細(xì)節(jié)之間的某種折中。因?yàn)閰?shù)γ在本文的算法中是自適應(yīng)的，所以直接設(shè)定λ為1。在文獻(xiàn)[17]中，對(duì)這兩個(gè)參數(shù)(γ和λ)的選擇進(jìn)行了更加詳細(xì)的討論。

關(guān)于濾波參數(shù)h的選擇，經(jīng)過(guò)大量實(shí)驗(yàn)后，對(duì)于Sheep-Logan模體和Forbild頭部模體，分別取0.02和0.03可以得出較為精確的結(jié)果。

圖2是ASD-POCS算法和本文算法在20和30個(gè)投影角度下的Sheep-Logan模體的重建結(jié)果圖像。圖2的實(shí)驗(yàn)結(jié)果表明：傳統(tǒng)TV算法重建的結(jié)果圖像中含有更多的噪點(diǎn)，并且圖像細(xì)節(jié)也被平滑掉；相對(duì)來(lái)說(shuō)，NLTV重建的結(jié)果圖像保留了更多的細(xì)節(jié)，圖像結(jié)構(gòu)邊緣清晰。

為了更清晰地定量比較圖2四幅重建結(jié)果的圖像，在圖3中給出了四幅圖像的垂直于圖像中心的切面垂直中心剖線。對(duì)于20個(gè)角度的重建條件，NLTV的垂直中心剖線擁有更少的波動(dòng)。在30個(gè)投影角度的情況下，NLTV的垂直中心剖線與真值圖像的垂直中心剖線基本重合。結(jié)果表明，相對(duì)于原TV算法，NLTV算法可以獲得更高的精度。

圖2 Sheep-Logan模體不含噪聲投影數(shù)據(jù)重建圖像Fig. 2 Reconstruction results of Sheep-Logan phantom by noise-free projection data

圖3 Sheep-Logan模體不含噪投影數(shù)據(jù)重建結(jié)果圖像垂直中心剖線Fig. 3 Vertical profile of Sheep-Logan phantom reconstructed by noise-free projection data

圖4給出了算法迭代過(guò)程的RMSE趨勢(shì)曲線。顯然，NLTV模型可以在更少的迭代次數(shù)下重建出更好的結(jié)果圖像，且在收斂點(diǎn)處擁有更高的精度。在20個(gè)角度投影數(shù)據(jù)重建的條件下，TV模型最終的收斂均方誤差是0.011，NLTV模型的最終收斂精度是0.003；在30個(gè)角度的條件下，TV模型最終精度是0.002，NLTV模型則得到了十萬(wàn)分之一精度(5.3×10-5)。

圖4 Sheep-Logan模體不含噪重建RMSE收斂趨勢(shì)曲線Fig. 4 RMSE convergence curve of reconstructing Sheep-Logan phantom by noise-free projection data

圖5顯示了30個(gè)和40個(gè)投影角度下Forbild頭部模體的重建結(jié)果。Forbild模體的右側(cè)有蜂窩狀結(jié)構(gòu)，在較小的區(qū)域有密集的灰度變化，對(duì)于重建算法的要求較高?？梢钥吹?，因?yàn)門(mén)V模型對(duì)于單一結(jié)構(gòu)內(nèi)灰度均勻這一先驗(yàn)條件的假設(shè)，TV模型對(duì)于蜂窩狀結(jié)構(gòu)的重建結(jié)果較差，有明顯的放射狀偽影，且TV模型在增加投影角度數(shù)量后，對(duì)于蜂窩狀結(jié)構(gòu)的重建結(jié)果并沒(méi)有明顯的改善；而NLTV模型對(duì)于蜂窩狀結(jié)構(gòu)的恢復(fù)更加精確。

圖5 Forbild頭部模體不含噪聲投影數(shù)據(jù)重建結(jié)果Fig. 5 Reconstruction results of Forbild head phantom by noise-free projection data

表1顯示了6組實(shí)驗(yàn)的均方根誤差結(jié)果。在Sheep-Logan模體的仿真實(shí)驗(yàn)中，TV模型在50個(gè)投影角度下的重建精度與NLTV模型20個(gè)角度下的精度在一個(gè)數(shù)量級(jí)。

表1 6組實(shí)驗(yàn)結(jié)果RMSE對(duì)比 Tab. 1 RMSE comparison of 6 sets of experimental results

3.2 含噪聲數(shù)據(jù)重建

本節(jié)比較了TV和NLTV模型在含噪聲Sheep-Logan模體投影數(shù)據(jù)的重建情況。首先，生成了50個(gè)角度的投影數(shù)據(jù)；然后，在投影數(shù)據(jù)中加入了方差為0.01的高斯白噪聲。

重建結(jié)果如圖6所示。可以看出，TV模型和NLTV模型都顯示出了對(duì)于噪聲的抑制特性，但是，NLTV模型圖像結(jié)構(gòu)之間有著更清晰的邊緣間隔，且整幅圖像看起來(lái)比TV模型更加銳利清晰。

圖6 Sheep-Logan模體含噪聲投影數(shù)據(jù)重建結(jié)果Fig. 6 Construction results of Sheep-Logan phantom by noise-added projection data

圖7顯示的垂直中心剖線結(jié)果表明，NLTV模型比TV模型擁有更高的精度，且在像素灰度值跳躍較大的地方更好地與真值重合。這也進(jìn)一步證實(shí)了重建結(jié)果圖像中顯示的，NLTV模型在重建圖像的不同結(jié)構(gòu)之間有著更銳利的邊緣區(qū)隔。

圖7 Sheep-Logan模體含噪重建結(jié)果圖像垂直中心剖線Fig. 7 Vertical profile of Sheep-Logan phantom reconstructed by noise-added projection data

圖8顯示的含噪聲重建過(guò)程的RMSE收斂曲線進(jìn)一步表明了，NLTV相比傳統(tǒng)TV模型能在更少的迭代次數(shù)下達(dá)到更高的收斂精度。NLTV模型在迭代到30次時(shí)已經(jīng)基本接近最終的收斂精度0.002 2，而TV模型則需要迭代到150次時(shí)，才基本接近最終的收斂精度0.005 5。

圖8 Sheep-Logan 模體含噪重建RMSE 收斂趨勢(shì)曲線Fig. 8 RMSE convergence curve of reconstructing Sheep-Logan phantom by noise-added projection data

3.3 參數(shù)對(duì)重建結(jié)果的影響

在這組實(shí)驗(yàn)中，并不計(jì)劃去尋找能使重建結(jié)果最優(yōu)的參數(shù)組合，僅僅通過(guò)選定測(cè)試參數(shù)的不同值來(lái)討論其對(duì)于重建結(jié)果的影響。

(28)

圖9 數(shù)據(jù)容差對(duì)重建結(jié)果的影響Fig. 9 Effect of data tolerance on reconstruction result

3.3.2 參數(shù)λ

接下來(lái)，簡(jiǎn)單討論平衡參數(shù)λ對(duì)于重建結(jié)果的影響。為了使對(duì)比更加明顯，把噪聲方差增加到0.03。

圖10的結(jié)果顯示，平衡參數(shù)λ對(duì)于重建結(jié)果的影響與數(shù)據(jù)容差有某種相似性：較大的λ可以增加數(shù)據(jù)保真項(xiàng)的權(quán)重，但會(huì)造成重建圖像中含有較多的噪聲，這是因?yàn)闆](méi)有很好地抑制在原始投影數(shù)據(jù)中所含的噪聲內(nèi)容；相反地，較小的λ可以對(duì)原始投影數(shù)據(jù)中的噪聲內(nèi)容進(jìn)行很好的過(guò)濾，但是會(huì)使結(jié)果丟失細(xì)節(jié)。

圖10 平衡參數(shù)λ對(duì)重建結(jié)果的影響Fig. 10 Effect of balance parameter λ on reconstruction result

4 結(jié)語(yǔ)

本文以ASD-POCS算法為基礎(chǔ)，結(jié)合非局部全變分模型的優(yōu)點(diǎn)，采用分離Bregman算法求解算法中的L1約束問(wèn)題，提出了ANLTVM-POCS算法。該算法結(jié)合了ASD-POCS算法參數(shù)自適應(yīng)的優(yōu)點(diǎn)和非局部全變分模型能更好保存圖像精細(xì)結(jié)構(gòu)及紋理的特點(diǎn)，在一定程度上改善了ASD-POCS算法采用TV模型不能很好恢復(fù)圖像細(xì)節(jié)的弱點(diǎn)，同時(shí)該算法也可以簡(jiǎn)單地?cái)U(kuò)展到諸如三維束重建和螺旋重建中。實(shí)驗(yàn)結(jié)果表明，該算法可以在較少的重建次數(shù)下達(dá)到較高的重建精度，且相比ASD-POCS算法達(dá)到收斂點(diǎn)所需迭代次數(shù)更少。由于非局部框架的引入，對(duì)于平衡參數(shù)λ和濾波參數(shù)h的選擇，本文只是在大量實(shí)驗(yàn)結(jié)果的基礎(chǔ)上給出的經(jīng)驗(yàn)參數(shù)。合適的參數(shù)是能重建出精確結(jié)果的保證，能根據(jù)模型的特點(diǎn)開(kāi)發(fā)出更加自動(dòng)精確的參數(shù)選擇算法是未來(lái)的工作重點(diǎn)之一。