彭 慧,孫維昆,冀東江
(天津職業(yè)技術師范大學理學院,天津 300222)
CT 圖像重建是指通過對物體進行不同角度的射線掃描,根據檢測到的數據來重建物體內部截面圖像的技術[1]。目前主要的圖像重建算法有2 類,分別為解析類重建算法和迭代類重建算法[2]。解析重建算法的特點為算法簡單,計算量小,重建速度快,但在投影數據不完備或投影數據噪聲大的情況下,獲得的圖像效果并不好[2]。解析重建算法中最常用的是濾波反投影算法(FBP)[3]。目前運用較多的迭代類重建算法有代數重建算法(ART)和聯(lián)合代數重建算法(SART)[4],該類算法可以獲得較好的重建圖像,但是耗時比較長[5]。ART 算法是由Gordon 等[6]于1970 年提出,該算法在低劑量重建時不能很好地起到去噪的作用,但是相比于FBP,其重建出來的圖像更好。SART 算法是由Anderson 等[7]于1984 年提出,該算法在ART 的基礎上進行了改進,可以解決ART 算法受噪聲影響的問題,但是由于該算法不具備高速收斂性,進行一次圖像更新需要所有的投影數據,所以在重建效率上比ART 算法低。Sidky 等[8]提出了基于最小化圖像總變差的優(yōu)化準則的迭代重建算法(TV),TV 模型可以很好地保持圖像的邊界,對于分片常數的圖像重建效果很好,但是該算法常會在分片光滑區(qū)城產生“階梯效應”。SART-TV 算法是一種基于最小化總變差,結合SART的重建算法。本文分別從稀疏角度和有限角度情況下,采用上述幾種算法對Shepp-Logan 模型進行仿真試驗,并比較其重建效果。
FBP 算法的基本思想[9-10]:掃描被測物體,先把獲得的投影函數做濾波處理,得到修正的投影函數,再對修正后的投影函數進行反投影,最后得到重建圖像。FBP 算法的實現步驟:①在角度θ 下對原圖像進行掃描,得到投影數據f(s,θ)。②求投影數據f(s,θ)的傅里葉變換,得到f(ω,θ),用濾波器的傳遞函數|ω|對f(ω,θ)進行濾波,得到Q(ω,θ)。③求Q(ω,θ)的以ω 為變量的一維傅里葉逆變換,得到q(ω,θ)。④用濾波后的數據q(ω,θ)進行反投影得到重建圖像f(x,y)[10-11]。
ART 算法的基本思想[10,12-13]:求解方程p=Au,其中u=[u1,u2,…,uN]T為待估計的圖像,p=[p1,p2,…,pM]T為系統(tǒng)采集的真實投影數據,A 為系統(tǒng)矩陣,其維度為M×N。
ART 算法的公式
式中:0 <λ<1 為松弛因子;k 為迭代次數;N 為圖像的大?。籭 為第i 條射線;j 為第j 個像素點;fj(k)為第j個像素點的灰度值;pi為第i 條射線的真實投影數據;wij為第i 條射線與第j 個像素相交的長度[10,12,14]。
ART 算法的實現步驟[10,12,14]:①為被重建對象賦初值,給定ART 的迭代次數,松弛因子fj(k)=0,k=0,(j=1,2,…,N)。②計算估計投影值(i=1,…,M)。③計算由②得到的估計投影值和實際投影值的誤差④利用誤差對所有射線經過的像素點進行修正。⑤將上面的結果作為第一步的初始值,重復所有步驟,直到滿足收斂條件或者達到迭代次數。
SART 算法的基本思想[10,14]:在修正某個像素值前,需要計算經過像素的所有射線的誤差來對該像素進行修正,并進行加權及歸一化處理之后,再把上述結果更新到該像素上去,重復上述過程,直到滿足收斂條件或迭代次數滿足人為約定的值。
SART 算法的公式
式中:pi為第條射線的實際投影數據;pφ為同一投影角度下的實際投影數據集合為理論投影值;λ為松弛因子[10,14]。
SART 算法的實現步驟:①為待重建對象賦初值,給定SART 的迭代次數,松弛因子fj(k)=0,k=0(j=1,2,…,N)。②計算理論投影值(i=1,…,M)。③計算由②得到的理論投影值與實際投影值的誤差④對其余射線,重復上述步驟,計算其余射線的誤差值。⑤用步驟④得到的所有射線誤差值修正待重建圖像,直到滿足收斂條件或者迭代次數滿足人為約定的值[10,14]。
SART-TV 算法的基本思想:基于TV 正則化,結合SART 算法。圖像的TV 表示為[13-15]
式中:f 為圖像函數,其中每一個像素由fj,j=1,2,…,N 表示;P 為系統(tǒng)采集的真實投影數據;A 矩陣為射線穿過圖像時和每一像素相交的長度的集合。要得到圖像的總變差最小,要求圖像總變差的梯度遞減。圖像像素對應的梯度定義為[13-15]
于是,圖像的總變差(TV)寫成下列形式
式中:fs,t為二維圖像第s 行,第t 列像素的灰度值[14-16]。
SART-TV 重建算法步驟:①輸入原始數據P,姿SART,kSART,KTV,for i=1:k,重建結果圖像。②用式(2)計算SART,得到fi。③用式(5)對fi進行TV 正則化約束。④輸出結果圖像。
本文分別采用均方誤差(MSE)、峰值信噪比(PSNR)[17]來評價圖像質量。
式中:ti,j、ri,j分別為原始圖像和重建圖像在像素(i,j)處的像素值;M×N 為圖像的大??;MAX1為原圖的最大像素值。
本文對Shepp-Logan 頭模型進行圖像重建仿真模擬,用Matlab R2018b 軟件在PC 機(4.0 GB 內存,3.10 GHz CPU)上實現。重建物體的大小為256×256,一個掃描角度下探測元個數設置為256 個。
本研究在180°范圍內依次等間隔取樣45 個投影和60 個投影,用FBP、ART、SART、SART-TV 4 種算法重建圖像。在ART 算法中,kART=17,姿ART=0.2;在SART 算法中,kSART=5,姿SART=0.3;在SART-TV 算法中,kSART=10,姿SART=0.1,kTV=25。圖1 和圖2 分別為45 個、60 個投影情況下重建出的圖像。
圖1 45 個投影角度下重建結果
圖2 60 個投影角度下重建結果
從主觀角度來看,ART、SART、SART-TV 算法重建出的圖像質量優(yōu)于FBP。因為這3 種迭代算法重建出的圖像只有較少的偽影,而FBP 算法重建出的圖像偽影較重。為了客觀分析不同算法重建圖像的質量,本文計算了2 種評價圖像質量的數值指標,如表1所示。
表1 4 種算法在2 種稀疏角度下的客觀評價指標
由表1 可以看出,ART、SART、SART-TV 算法的均方誤差小,峰值信噪比大,重建出的圖像要好于FBP 算法。隨著投影數據的增多,4 種算法的均方誤差逐漸變小,峰值信噪比逐漸增大。由此可見,在稀疏采樣的情況下,利用ART、SART、SART-TV 算法重建圖像比較合適,相比之下SART-TV 重建結果更好。
圖3 和圖4 分別為45 個、60 個投影角度下每種迭代算法隨著迭代次數增加時MSE 的結果。
圖3 45 個投影角度下MSE 結果
圖4 60 個投影角度下MSE 結果
從圖3、圖4 中可以看出,隨著迭代次數的增多,3種算法的均方誤差逐漸減小,且逐漸趨于穩(wěn)定,所以這3 種迭代算法是收斂的。并且在60 個稀疏投影下,每種算法的均方誤差都比45 個稀疏投影情況下均方誤差小。
本文在[0°,130°]、[0°,150°]以及[0°,170°]范圍內,每隔1°均勻采樣,采用FBP、ART、SART、SART-TV 算法進行圖像重建。在ART 算法中,kART=17,姿ART=0.2;在SART 算法中,kSART=5,姿SART=0.3;在SART-TV 算法中,kSART=10,姿SART=0.1,kTV=25。
有限角度的重建結果如圖5 所示。從圖5 可以看出,FBP 算法重建出的圖像條狀偽影較多,其他算法雖然也有偽影,但是比FBP 的更少。用FBP 算法進行重建時圖像的細節(jié)恢復程度也沒有迭代類算法的恢復程度好,并且獲得的重建圖像的效果從好到差依次為:SART-TV、SART、ART、FBP。隨著角度的增多,4種算法重建出的圖像逐漸變好,獲得的圖像數據越多,重建出來的圖像也越接近原始圖像。4 種算法在不同有限角下的客觀評價指標如表2 所示。
圖5 有限角度的重建結果
表2 4 種算法在不同有限角下的客觀評價指標
從表2 可以看出,ART、SART、SART-TV 算法的均方誤差小,峰值信噪比大,FBP 的均方誤差比其他3種算法大,峰值信噪比比其他3 種算法小,相比之下SART-TV 算法的均方誤差最小,峰值信噪比最大。所以,在有限角度的情況下,相比其他3 種算法,使用SART-TV 重建得到的圖像更好。
在[0°,130°]、[0°,150°]、[0°,170°]下每種迭代算法隨著迭代次數增加時MSE 的結果如圖6—圖8所示。從圖6—圖8 中可以看出,隨著角度的增大,每種算法的MSE 都是減小的,同時隨著迭代次數的增多,3 種算法的均方誤差逐漸減小,且逐漸趨于穩(wěn)定。于是,在有限角度情況下,這3 種迭代算法是收斂的。
圖6 [0°,130°]下MSE 結果
圖7 [0°,150°]下MSE 結果
圖8 [0°,170°]下MSE 結果
本文利用4 種算法對Shepp-Logan 頭模型進行了實驗,對稀疏角度和有限角度投影數據進行圖像重建,并分析每種算法的優(yōu)缺點及適用情況。通過對重建圖像效果分析可知,用Shepp-Logan 頭模型進行實驗時,迭代類重建算法的重建效果比解析類好,重建效果的優(yōu)劣依次為:SART-TV、SART、ART、FBP。但是,每種算法都有其局限性,并不能在所有情況下都能重建出完美的圖像,同時迭代類的算法還會受迭代次數以及松弛因子的影響且算法耗時長。此問題還需要進一步研究,今后針對不同的模型進行實驗,提出更適用的改進算法。