自適應步長非局部全變分約束迭代圖像重建算法

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

(山西大學 計算機與信息技術學院，太原030006)

0 引言

計算機斷層成像(Computed Tomography, CT)在醫學和工業生產中有非常重要和廣泛的應用。絕大部分的商用CT機使用傳統的濾波反投影(Filtered Back Projection, FBP)算法來重建圖像已經有將近30年的歷史[1]；然而，大劑量的輻射會增加患者患癌的風險，低劑量的掃描，即數據不完備的條件下，FBP算法難以重建出精確的圖像。為了解決這個問題，有兩種常用思路：一是通過插值方法對不完備的投影數據進行補足；二是通過在迭代過程中對重建圖像施加某些約束條件，即對模型增加正則化項的手段來解決。一些基于迭代的算法諸如代數迭代重建(Algebraic Reconstruction Technique, ART)[2]和期望最大化(Expectation Maximization, EM)算法[3]已經被應用于圖像重建中[4]?；趦灮牡惴梢岳脡嚎s感知、低秩矩陣等稀疏優化技術對圖像進行精確重建[5]。其中，全變分(Total Variation, TV)最小化約束的迭代重建算法已經在各種CT中廣泛的應用。

但是，TV模型的主要缺點是它不能恢復和保留圖像中的精細結構和紋理，并傾向于產生階梯效應(staircase artifacts)。因為傳統的TV模型僅考慮圖像的局部像素信息，且假設圖像在單個對象內是平滑的，只在不同對象的邊界處具有不連續的灰度跳躍[6]。

Buades等[7]利用兼顧了圖像中更多全局信息的非局部均值(NonLocal Means, NLM)算法來進行圖像去噪。該算法通過對一個像素點周圍像素的加權平均值來恢復該像素，因此，相較于傳統的局部算子能更好地恢復和保留圖像的細節。Gilboa等[8]將傳統的變分框架和非局部算子進行了有機的結合，提出了非局部變分(NonLocal TV, NLTV)框架。相較于局部變分框架,非局部變分框架在基于優化的CT重建算法中能得到更加精確的結果是可以合理預期的。

Sidky等[9]提出了自適應最速下降投影到凸集(Adaptive Steepest Descent Projection Onto Convex Sets, ASD-POCS)算法來解決圖像重建中基于約束的TV最小化問題。其中，自適應最速下降(ASD)即算法采用梯度下降法來使圖像的全變分最小化，同時步長因子是自適應的。投影到凸集(POCS)即為ART與正約束的合稱，ART保證重建結果滿足投影數據保真項的約束，正約束是因為默認圖像不存在負值像素點。ASD-POCS算法可以有效和魯棒地從不完備數據中重建出相當精確的結果。它交替地運行POCS過程和TV下降過程，從而使結果逐步地接近TV最小化解，同時滿足數據保真項和正則項的要求。

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

1 相關工作

1.1 CT模型

非約束的、基于優化的CT數學模型可以表示為如下的形式：

(1)

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

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

J(u)=‖u‖TV=

(2)

其中us,t表示圖像單個像素點。顯然，傳統的TV模型僅僅考慮了一個像素毗鄰的另外兩個像素，相較于非局部算法框架來說，難以保留圖像的細節和紋理。

1.2 非局部變分框架

接下來，簡要回顧在文獻[8]中介紹的非局部算子。令Ω?R2，x∈Ω，并且u(x):Ω→R是實函數。設w:Ω×Ω→R是關于參考圖像的非負對稱實函數，即：w(x,y)=w(y,x)。關于一像素點的非局部導數可以定義為：

(3)

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

(4)

對于向量u在點x∈Ω的非局部梯度L2范數：

(5)

所有像素點的非局部梯度L2范數也可以組成一幅圖像，即加權非局部梯度圖像。

對于向量v1=v1(x,y)，v2=v2(x,y)，向量內積定義為：

(6)

根據向量內積的定義及散度和梯度的共軛關系，

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

(7)

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

(8)

拉普拉斯算子可以寫為：

(9)

權函數w(x,y)：

(10)

其中：Ga是標準差為a的高斯核函數，參數h是濾波參數。式(10)顯示了權函數只有在關于點x和點y的比較像素塊相似度大時具有較大的取值。

為了便于計算機進行離散運算，分別給出非局部梯度、非局部梯度L2范數、非局部散度和拉普拉斯算子的離散形式：

(11)

(12)

(13)

(14)

綜上所述，非局部全變分可以看作加權非局部梯度圖像的L1范數：

(15)

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

對于非局部框架，其核心問題是對三個參數的選?。罕容^塊的尺寸(patch size)、搜索窗的尺寸(search window size)和式(10)中濾波參數h的取值；然而，這三個參數之間的耦合影響非常嚴重，作出一個恰當的選擇是一項比較艱巨的任務。本文在大量實驗的基礎上，嘗試了幾種參數組合，然后給出一些可以得到相對精確結果的經驗參數選擇。

理論上，搜索窗的尺寸應該是整幅圖像，然而實際上，一個過大的搜索窗會導致計算開銷的幾何倍增加。許多文獻給出了從11×11到21×21之間的不同選擇，一個好的搜索窗尺寸應該是在兼顧精確結果和計算開銷的中間選擇。在本次工作中，發現11×11大小的搜索窗已經足夠得出相對精確的結果且僅需要相對較少的計算開銷。

在一些關于NLTV圖像去噪的研究[15]中提到，大的比較塊尺寸可以增強算法對于噪聲的魯棒性，較小的比較塊尺寸可以保留更多的低對比度圖像細節。在式(10)中，可以發現，結果像素點的取值還和鄰域像素點的高斯均值相關。對于不同的成像模型以及不同的噪聲類型，應該根據情況作出恰當的選擇。在本文的實驗中，大小為3×3的比較塊已經可以得到較好的視覺結果圖像。

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

1.3 分離Bregman算法

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

(16)

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

標準ASD-POCS算法使用梯度下降法來最小化TV。對于NLTV，使用分離Bregman算法來求解這個最優化問題。引入一個輔助變量d:Ω×Ω→R來替換式(16)中的▽NLu，然后可以得出如下的約束最優化問題：

(17)

s.t.d=▽NLu

式(17)的約束問題可以轉化為如下無約束形式：

(18)

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

(20)

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

(21)

(22)

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

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

(23)

根據式(23)，使用Gauss-Seidel迭代算法來求取u的近似解(給出關于一個像素點在第k+1次迭代的離散形式，同d和b)：

(24)

關于子問題d，可以通過收縮算子(shrinkage operators)[10]來求?。?/p>

(25)

最后，b可以通過直接迭代求解：

(26)

2 本文方法

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

本文方法主要基于ASD-POCS算法，并且針對ASD-POCS算法中TV模型的缺點，改用NLTV模型來約束最終解，提高算法保留重建結果中精細結構和紋理的特性。

且針對NLTV傳統迭代方法(如梯度下降法)難以求解的問題，利用分離Bregman方法求解L1正則化問題的優勢來求取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是投影的總個數，行7)的Ni是總像素個數。算法中有6個參數(β;βred;α;αred;n;rmax)控制整個算法的迭代過程，在文獻[9]中有詳細的敘述，本文不再展開。u是算法的初始化圖像，在這里設定為零。行3)和15)的up是臨時變量，用來存儲算法當前步驟的結果并參與步長自適應調整的計算。dtvg是算法最小化子問題17)到21)行的步長因子。算法的另一個重要參數是數據容差，是多種數據不一致的簡化，包括簡化模型帶來的誤差、投影數據中的噪聲以及X射線的散射等，這使得不可能總是重建出與數據完全一致的結果圖像[9]。參數確定了重建生成圖像的投影與實際投影數據的歐氏距離(L2范數)。在結果部分討論了關于的選擇對于重建結果的影響。

把標準ASD-POCS算法的TV梯度下降過程替換為基于NLTV的分離Bregman迭代過程。這里有兩個重要參數(λ,γ)在分離Bregman的迭代過程中占據主導性地位，進而影響最終重建結果的精度。在下一章詳細展開討論。

3 實驗結果與討論

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

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

其中：u和f是分別是重建結果圖像和真值圖像；N是圖像的總像素個數。

圖1 仿真模體Fig. 1 Simulation phantoms

3.1 不含噪聲數據重建

至于算法的停止條件(stopping criteria)，可以當計算結果達到一個特定的條件(如‖uk-uk+1‖/‖uk‖<10-3)時結束循環，但是為了更直觀地描述本文的算法，在本組實驗中使算法直接循環500次，以便更全面地通過算法迭代過程的RMSE曲線來分析算法的收斂特性。

參數γ的恰當選擇能使分離Bregman算法的子問題u和d收斂的速度更快。在標準分離Bregman算法中，為了使算法更好地收斂，γ既不能過大也不能過小。文獻[10]中提到，當γ=2λ時，分離Bregman算法通?？梢匀〉幂^好的收斂結果。在本文的算法中，參數γ可以說是連接數據保真項和正則化項之間的橋梁。本文直接設定γ等于標準ASD-POCS算法框架中的自適應步長參數dtvg，以使ASD-POCS算法步長參數自適應的特點與分離Bregman迭代過程有機地結合起來。

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

關于濾波參數h的選擇，經過大量實驗后，對于Sheep-Logan模體和Forbild頭部模體，分別取0.02和0.03可以得出較為精確的結果。

圖2是ASD-POCS算法和本文算法在20和30個投影角度下的Sheep-Logan模體的重建結果圖像。圖2的實驗結果表明：傳統TV算法重建的結果圖像中含有更多的噪點，并且圖像細節也被平滑掉；相對來說，NLTV重建的結果圖像保留了更多的細節，圖像結構邊緣清晰。

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

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

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

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

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

圖5顯示了30個和40個投影角度下Forbild頭部模體的重建結果。Forbild模體的右側有蜂窩狀結構，在較小的區域有密集的灰度變化，對于重建算法的要求較高?？梢钥吹?，因為TV模型對于單一結構內灰度均勻這一先驗條件的假設，TV模型對于蜂窩狀結構的重建結果較差，有明顯的放射狀偽影，且TV模型在增加投影角度數量后，對于蜂窩狀結構的重建結果并沒有明顯的改善；而NLTV模型對于蜂窩狀結構的恢復更加精確。

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

表1顯示了6組實驗的均方根誤差結果。在Sheep-Logan模體的仿真實驗中，TV模型在50個投影角度下的重建精度與NLTV模型20個角度下的精度在一個數量級。

表1 6組實驗結果RMSE對比 Tab. 1 RMSE comparison of 6 sets of experimental results

3.2 含噪聲數據重建

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

重建結果如圖6所示?？梢钥闯觯琓V模型和NLTV模型都顯示出了對于噪聲的抑制特性，但是，NLTV模型圖像結構之間有著更清晰的邊緣間隔，且整幅圖像看起來比TV模型更加銳利清晰。

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

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

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

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

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

3.3 參數對重建結果的影響

在這組實驗中，并不計劃去尋找能使重建結果最優的參數組合，僅僅通過選定測試參數的不同值來討論其對于重建結果的影響。

(28)

圖9 數據容差對重建結果的影響Fig. 9 Effect of data tolerance on reconstruction result

3.3.2 參數λ

接下來，簡單討論平衡參數λ對于重建結果的影響。為了使對比更加明顯，把噪聲方差增加到0.03。

圖10的結果顯示，平衡參數λ對于重建結果的影響與數據容差有某種相似性：較大的λ可以增加數據保真項的權重，但會造成重建圖像中含有較多的噪聲，這是因為沒有很好地抑制在原始投影數據中所含的噪聲內容；相反地，較小的λ可以對原始投影數據中的噪聲內容進行很好的過濾，但是會使結果丟失細節。

圖10 平衡參數λ對重建結果的影響Fig. 10 Effect of balance parameter λ on reconstruction result

4 結語

本文以ASD-POCS算法為基礎，結合非局部全變分模型的優點，采用分離Bregman算法求解算法中的L1約束問題，提出了ANLTVM-POCS算法。該算法結合了ASD-POCS算法參數自適應的優點和非局部全變分模型能更好保存圖像精細結構及紋理的特點，在一定程度上改善了ASD-POCS算法采用TV模型不能很好恢復圖像細節的弱點，同時該算法也可以簡單地擴展到諸如三維束重建和螺旋重建中。實驗結果表明，該算法可以在較少的重建次數下達到較高的重建精度，且相比ASD-POCS算法達到收斂點所需迭代次數更少。由于非局部框架的引入，對于平衡參數λ和濾波參數h的選擇，本文只是在大量實驗結果的基礎上給出的經驗參數。合適的參數是能重建出精確結果的保證，能根據模型的特點開發出更加自動精確的參數選擇算法是未來的工作重點之一。