CT 圖像重建算法對(duì)比研究

彭 慧，孫維昆，冀東江

（天津職業(yè)技術(shù)師范大學(xué)理學(xué)院，天津 300222）

CT 圖像重建是指通過對(duì)物體進(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 模型可以很好地保持圖像的邊界，對(duì)于分片常數(shù)的圖像重建效果很好，但是該算法常會(huì)在分片光滑區(qū)城產(chǎn)生“階梯效應(yīng)”。SART-TV 算法是一種基于最小化總變差，結(jié)合SART的重建算法。本文分別從稀疏角度和有限角度情況下，采用上述幾種算法對(duì)Shepp-Logan 模型進(jìn)行仿真試驗(yàn)，并比較其重建效果。

1 算法介紹及圖像數(shù)值評(píng)價(jià)指標(biāo)

1.1 FBP算法

FBP 算法的基本思想[9-10]：掃描被測物體，先把獲得的投影函數(shù)做濾波處理，得到修正的投影函數(shù)，再對(duì)修正后的投影函數(shù)進(jìn)行反投影，最后得到重建圖像。FBP 算法的實(shí)現(xiàn)步驟：①在角度θ 下對(duì)原圖像進(jìn)行掃描，得到投影數(shù)據(jù)f（s，θ）。②求投影數(shù)據(jù)f（s，θ）的傅里葉變換，得到f（ω，θ），用濾波器的傳遞函數(shù)|ω|對(duì)f（ω，θ）進(jìn)行濾波，得到Q（ω，θ）。③求Q（ω，θ）的以ω 為變量的一維傅里葉逆變換，得到q（ω，θ）。④用濾波后的數(shù)據(jù)q（ω，θ）進(jìn)行反投影得到重建圖像f（x，y）[10-11]。

1.2 ART算法

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 條射線；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]：①為被重建對(duì)象賦初值，給定ART 的迭代次數(shù)，松弛因子fj(k)=0，k=0，（j=1，2，…，N）。②計(jì)算估計(jì)投影值（i=1，…，M）。③計(jì)算由②得到的估計(jì)投影值和實(shí)際投影值的誤差④利用誤差對(duì)所有射線經(jīng)過的像素點(diǎn)進(jìn)行修正。⑤將上面的結(jié)果作為第一步的初始值，重復(fù)所有步驟，直到滿足收斂條件或者達(dá)到迭代次數(shù)。

1.3 SART算法

SART 算法的基本思想[10，14]：在修正某個(gè)像素值前，需要計(jì)算經(jīng)過像素的所有射線的誤差來對(duì)該像素進(jìn)行修正，并進(jìn)行加權(quán)及歸一化處理之后，再把上述結(jié)果更新到該像素上去，重復(fù)上述過程，直到滿足收斂條件或迭代次數(shù)滿足人為約定的值。

SART 算法的公式

式中：pi為第條射線的實(shí)際投影數(shù)據(jù)；pφ為同一投影角度下的實(shí)際投影數(shù)據(jù)集合為理論投影值；λ為松弛因子[10，14]。

SART 算法的實(shí)現(xiàn)步驟：①為待重建對(duì)象賦初值，給定SART 的迭代次數(shù)，松弛因子fj(k)=0，k=0（j=1，2，…，N）。②計(jì)算理論投影值（i=1，…，M）。③計(jì)算由②得到的理論投影值與實(shí)際投影值的誤差④對(duì)其余射線，重復(fù)上述步驟，計(jì)算其余射線的誤差值。⑤用步驟④得到的所有射線誤差值修正待重建圖像，直到滿足收斂條件或者迭代次數(shù)滿足人為約定的值[10，14]。

1.4 SART-TV算法

SART-TV 算法的基本思想：基于TV 正則化，結(jié)合SART 算法。圖像的TV 表示為[13-15]

式中：f 為圖像函數(shù)，其中每一個(gè)像素由fj，j=1，2，…，N 表示；P 為系統(tǒng)采集的真實(shí)投影數(shù)據(jù)；A 矩陣為射線穿過圖像時(shí)和每一像素相交的長度的集合。要得到圖像的總變差最小，要求圖像總變差的梯度遞減。圖像像素對(duì)應(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）對(duì)fi進(jìn)行TV 正則化約束。④輸出結(jié)果圖像。

1.5 圖像數(shù)值評(píng)價(jià)指標(biāo)

本文分別采用均方誤差（MSE）、峰值信噪比（PSNR）[17]來評(píng)價(jià)圖像質(zhì)量。

式中：ti，j、ri，j分別為原始圖像和重建圖像在像素（i，j）處的像素值；M×N 為圖像的大??；MAX1為原圖的最大像素值。

2 實(shí)驗(yàn)結(jié)果及分析

本文對(duì)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è)。

2.1 稀疏角度投影數(shù)據(jù)重建結(jié)果比較

本研究在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 種評(píng)價(jià)圖像質(zhì)量的數(shù)值指標(biāo)，如表1所示。

表1 4 種算法在2 種稀疏角度下的客觀評(píng)價(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è)稀疏投影情況下均方誤差小。

2.2 有限角度投影數(shù)據(jù)重建結(jié)果比較

本文在[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 種算法在不同有限角下的客觀評(píng)價(jià)指標(biāo)如表2 所示。

圖5 有限角度的重建結(jié)果

表2 4 種算法在不同有限角下的客觀評(píng)價(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é)果

3 結(jié)語

本文利用4 種算法對(duì)Shepp-Logan 頭模型進(jìn)行了實(shí)驗(yàn)，對(duì)稀疏角度和有限角度投影數(shù)據(jù)進(jìn)行圖像重建，并分析每種算法的優(yōu)缺點(diǎn)及適用情況。通過對(duì)重建圖像效果分析可知，用Shepp-Logan 頭模型進(jìn)行實(shí)驗(yàn)時(shí)，迭代類重建算法的重建效果比解析類好，重建效果的優(yōu)劣依次為：SART-TV、SART、ART、FBP。但是，每種算法都有其局限性，并不能在所有情況下都能重建出完美的圖像，同時(shí)迭代類的算法還會(huì)受迭代次數(shù)以及松弛因子的影響且算法耗時(shí)長。此問題還需要進(jìn)一步研究，今后針對(duì)不同的模型進(jìn)行實(shí)驗(yàn)，提出更適用的改進(jìn)算法。