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