CT 圖像重建算法的比較研究

陳宗桂，董曉軍，張英俊，管海辰

（湖南醫(yī)藥學(xué)院 醫(yī)學(xué)院，湖南懷化 418000）

計(jì)算機(jī)斷層掃描(Computerized Tomography，CT)在肺部疾病和腦出血檢查方面優(yōu)于其他設(shè)備[1-2]，是X 射線影像診斷的重要組成部分。CT 掃描過程中X 射線的生物學(xué)效應(yīng)、圖像質(zhì)量和X 射線劑量三者之間的關(guān)系，一直是人們關(guān)注的問題[3-5]。若患者檢查時(shí)X 射線的輻射劑量盡可能低，則采集圖像信噪比低，分辨率差[6-7]。為此，該研究對(duì)比采用濾波反投影、直接反投影和迭代重建算法得到CT 圖像，以峰值信噪比、結(jié)構(gòu)相似性和平均梯度評(píng)估重建圖像的質(zhì)量，總結(jié)出最合適的算法，以輔助臨床醫(yī)師提高病灶檢出率。

1 算法介紹

1.1 直接反投影算法

設(shè)f(x,y) 為被檢測(cè)物體截面線性吸收系數(shù)的分布函數(shù)，一束強(qiáng)度為I0的X 射線從被檢測(cè)物體透射出后X 射線強(qiáng)度降為I。X 線強(qiáng)度I滿足比爾定理，即：

式中，θ是X 線束的投影方向與直角坐標(biāo)系中X軸的夾角，表示螺旋掃描中X 線束的投影方向；R是單束X 線束在投射方向上的一維函數(shù)。令gθ(R)=ln(I0/I)得到式（2）：

gθ(R)是穿過物體后X 線束的投影值。從式（2）可以看出投影值是對(duì)沿著投照方向?qū)Ρ粰z體截面線性吸收系數(shù)的分布函數(shù)f(x,y)求線積分值，而圖像重建是從每個(gè)方向上獲得不同的投影值，通過反投影重建以求得物體組織密度的分布函數(shù)f(x,y)。

CT 圖像重構(gòu)的理論基礎(chǔ)[8]是傅里葉投影切片定理，它指出f(x,y)在某一方向上的投影函數(shù)gθ(R)的一維傅里葉變換函數(shù)gθ(ρ)是f(x,y)的二維傅里葉變換函數(shù)F(p,θ)(極坐標(biāo)形式)在(p,θ)平面上過原點(diǎn)的一條直線，如圖1 所示。

圖1 中心切片定理示意圖

因此，為了從不同投照方向獲得足夠的數(shù)據(jù)，并對(duì)每一個(gè)方向的數(shù)據(jù)進(jìn)行傅里葉變換，將變換后的數(shù)據(jù)填滿頻域空間。對(duì)完整頻域數(shù)據(jù)進(jìn)行傅里葉逆變換就可得到物體斷面組織密度的圖像f(x,y)，如式（3）：

式中，F(xiàn)(u,v)表示投影數(shù)據(jù)的二維頻譜圖。

1.2 濾波反投影算法

為了抑制CT 圖像上的“星芒狀偽影”，需要對(duì)采集的原始數(shù)據(jù)進(jìn)行濾波處理。因此，令u=ρcosθ，v=ρsinθ，則式（3）可進(jìn)一步變形為：

根據(jù)傅里葉變換的性質(zhì)，頻域范圍內(nèi)的圖像濾波可以通過空域范圍內(nèi)的圖像進(jìn)行卷積操作來完成，因此式（4）又可以寫成式（5）：

式中，h(R) 為濾波函數(shù) |ρ|的空域形式，因而gθ(R)h(R)δ(xcosθ+ysinθ-R)表示對(duì)空域圖像進(jìn)行卷積得到初始圖像即濾波反投影方法。

在濾波反投影算法中，濾波器的設(shè)計(jì)是影響圖像精度至關(guān)重要的因素。理想濾波器是具有無限頻帶和積分結(jié)果不收斂。根據(jù)Paley–Wiener 準(zhǔn)則，理想濾波器是不存在的[9]。但是，在實(shí)際成像過程中，對(duì)采集的數(shù)據(jù)增加一個(gè)窗口函數(shù)就可以得到理想的數(shù)據(jù)精度。因?yàn)镃T 圖像分辨率有限，使得探測(cè)器采集數(shù)據(jù)能量分布都集中在低頻部分。所以，理想濾波函數(shù) |ρ|必須加窗，換言之，僅保留頻域空間上的低頻數(shù)據(jù)。

該研究采用一種比較平滑的窗函數(shù)sinc(ρ/2ρ0)rect(ρ/2ρ0)來約束濾波函數(shù) |ρ|，從而得到了S-L 濾波函數(shù)。其頻域公式為：

對(duì)應(yīng)的空域離散形式為：

S-L 濾波函數(shù)重建圖像的波動(dòng)幅值較低，抗噪聲性能強(qiáng)和細(xì)節(jié)分辨率高。S-L 濾波函數(shù)對(duì)低頻段數(shù)據(jù)響應(yīng)效果好，而對(duì)于高頻段數(shù)據(jù)響應(yīng)就不如R-L 濾波函數(shù)。

1.3 迭代重建算法

迭代重建算法將圖像數(shù)據(jù)視為未知圖像矩陣，將在任何投影角處收集的投影數(shù)據(jù)與圖像初始猜測(cè)解進(jìn)行比較，通過比較結(jié)果更新模擬圖像，直至模擬圖像逼近原始圖像。算法的原理是首先設(shè)置一組模擬圖像矩陣，從不同角度采集數(shù)據(jù)，并與實(shí)際采集的投影數(shù)據(jù)進(jìn)行比較。通過反向投影比較結(jié)果來更新模擬圖像，并且可以重復(fù)獲得重構(gòu)圖像，直到滿足迭代結(jié)束條件。圖2是CT圖像迭代重建的系統(tǒng)框圖。每次迭代使用一個(gè)方程組，即每次迭代都只考慮一個(gè)輻射的影響和貢獻(xiàn)，通過一定次數(shù)的迭代逐次近似所需圖像。

圖2 迭代重建算法的迭代過程

迭代重建(Algebraic Reconstruction Technique，ART)算法先設(shè)定模擬圖像矩陣投影值qi。當(dāng)投影角度為θ時(shí)，計(jì)算模擬圖像矩陣的投影值qi和實(shí)際投影值pi。通過誤差Δp=pi-qi對(duì)投影角度θ下的所有像素點(diǎn)進(jìn)行校正。迭代重建的公式為:

2 CT圖像的重建結(jié)果

在該研究中，采用三種不同的算法重構(gòu)CT圖像，如圖3 所示。圖3（b）是采用同步迭代算法重建的圖像。在實(shí)際應(yīng)用中，由于圖像尺寸大，為了求解圖像上每一個(gè)像素值需要建立多個(gè)方程組，計(jì)算量大，復(fù)雜程度高，所以提出了迭代重建方法。圖3（b）是采用式（8）進(jìn)行迭代重建CT 圖像的。該算法重建圖像時(shí)對(duì)噪聲不敏感，它與其他代數(shù)迭代算法不同的地方在于每次迭代不止對(duì)一條投影方向上的像素進(jìn)行修正，還修正了與該直線平行的所有投影方向上的投影數(shù)據(jù)，這是該算法能有效抑制圖像重建過程中的噪聲的主要原因。隨著迭代次數(shù)的無限增加，該算法得到的圖像無限逼近原始圖像。但是，算法收斂速度慢，在重構(gòu)圖像的邊緣和細(xì)節(jié)存在不同程度的模糊。圖3（c）是采用直接反投影算法得到的圖像。由于在反投影不斷疊加的過程中，中間低頻信號(hào)不斷疊加，而周圍高頻信號(hào)不斷缺失，使得圖像中感興趣區(qū)域的邊緣出現(xiàn)了“星芒狀偽影”。這種情況在投影數(shù)據(jù)少的時(shí)候更明顯。圖3（d）是采用濾波反投影算法得到的圖像。探測(cè)器采集的是0～180°方向的投影信號(hào)，濾波反投影算法是分別對(duì)每一個(gè)方向的投影數(shù)據(jù)進(jìn)行一維傅里葉變換后填滿頻譜圖。所有投影信號(hào)在頻域中濾波，即乘以加權(quán)因子r。通過傅里葉逆變換將濾波后的投影數(shù)據(jù)恢復(fù)到時(shí)域，得到原始圖像。每一次濾波后投影信號(hào)被反投影且最后重疊。該算法采用的濾波器是S-L濾波器。S-L 濾波器不是在頻域中加窗截?cái)嘣肼暫透蓴_信號(hào)，而是通過一些相對(duì)平滑的窗口函數(shù)對(duì)空域圖像卷積。因此，S-L 濾波后圖像振蕩幅值更小，圖像的結(jié)構(gòu)特征顯示更清晰。

圖3 三種不同算法重建得到顱腦CT圖像

3 圖像質(zhì)量的評(píng)價(jià)

該研究采用三種算法對(duì)CT 圖像進(jìn)行重建，如果想知道原始圖像的質(zhì)量是否有所提高，則需要采用一系列的定量指標(biāo)進(jìn)行評(píng)價(jià)。圖像信息越詳細(xì)，圖像邊緣灰度變化越快，圖像越清晰，人類視覺分辨率越高，結(jié)構(gòu)特征越明顯。該研究采用峰值信噪比、結(jié)構(gòu)相似性、平均梯度對(duì)重構(gòu)后的圖像邊緣結(jié)構(gòu)特征、清晰度和圖像相似性進(jìn)行評(píng)估。

3.1 圖像重建時(shí)間

為了提高X 射線的檢查效率，縮短設(shè)備的成像時(shí)間，改進(jìn)圖像重建算法的運(yùn)行速度就顯得尤其重要。為此，選擇的五組圖片的尺寸分別是128 pixel×128 pixel、256 pixel×256 pixel、512 pixel×512 pixel、1 024 pixel×1 024 pixel 和2 048 pixel×2 048 pixel。采用不同的算法重建相同尺寸的圖像，比較這三種算法的運(yùn)行時(shí)間，如圖4 所示。通過圖4 可知，針對(duì)相同尺寸的圖像，直接反投影和濾波反投影的重建時(shí)間基本一致。但是，迭代重建算法的運(yùn)行時(shí)間明顯大于另外兩種算法。另外，圖像矩陣從128 pixel×128 pixel 增至2 048 pixel×2 048 pixel，三種重建算法的運(yùn)行時(shí)間都變長。

圖4 三種算法重建CT圖像所需的運(yùn)行時(shí)間

3.2 基于峰值信噪比的圖像質(zhì)量評(píng)價(jià)

隨著現(xiàn)代醫(yī)療影像設(shè)備的迅速發(fā)展，其無論是在影像成像方面還是在圖像后處理方面，都有著廣泛的應(yīng)用前景，新的成像模式和圖像處理技術(shù)一直在提高和完善。因此，醫(yī)學(xué)圖像質(zhì)量評(píng)價(jià)已成為制定治療方案、評(píng)估外科手術(shù)、預(yù)后療效評(píng)估的重要參考。高質(zhì)圖像能如實(shí)反映人體內(nèi)部的解剖結(jié)構(gòu)，為影像診斷和治療提供豐富的信息。在精準(zhǔn)影像服務(wù)中，圖像質(zhì)量在保證圖像準(zhǔn)確度方面起著重要作用。目前，評(píng)價(jià)圖像質(zhì)量的方法有主觀評(píng)價(jià)方法和客觀評(píng)價(jià)法[10-13]。主觀評(píng)價(jià)方法以觀察者為主體，對(duì)被測(cè)圖像的優(yōu)缺點(diǎn)進(jìn)行評(píng)估。客觀評(píng)價(jià)法是采用相應(yīng)算法計(jì)算圖像的一些指標(biāo)，通過這些指標(biāo)來評(píng)價(jià)圖像優(yōu)劣。該研究采用均方誤差（Mean Squared Error，MSE）和峰值信噪比（Peak Signal to Noise Rate，PSNR）評(píng)價(jià)重構(gòu)圖像[14]。均方誤差法先計(jì)算初始圖像和重建圖像的差值，再除以行與列的乘積，如式（9）。均方值反映了該算法重建得到圖像與初始圖像的逼近程度。該方法的優(yōu)點(diǎn)是簡單、快速，可以將誤差量化表達(dá)且具有一定的參考價(jià)值。根據(jù)式（10）可求出不同算法下重構(gòu)圖像的峰值信噪比。若重建后圖像的PSNR 值越高，就代表該圖的微細(xì)結(jié)構(gòu)特征反映好，失真少。該研究通過比較峰值信噪比來判斷三種算法重建得到的CT 圖像的質(zhì)量，如圖5 所示。針對(duì)五組不同大小的圖像，每一種算法重建得到圖像的PSNR 變化基本一致。但是，和另外兩種算法相比，濾波反投影算法的峰值信噪比最好，失真最小。直接反向投影算法重建圖像的峰值信噪比最小，微細(xì)結(jié)構(gòu)特征反映差。

圖5 比較三種算法重建得到的CT圖像的峰值信噪比

3.3 基于結(jié)構(gòu)相似性的圖像質(zhì)量評(píng)價(jià)

結(jié)構(gòu)相似度(SSIM)是對(duì)亮度相似性、對(duì)比度相似性以及結(jié)構(gòu)相似性這三個(gè)參量的聯(lián)合度量，是評(píng)價(jià)兩幅圖像特征結(jié)構(gòu)相似性的有效度量。通過比較這三部分的信息，從而獲得兩幅圖像的整體度量[15]。最后采用求均值的方法計(jì)算得重建前后兩幅圖像的結(jié)構(gòu)相似信息。其表達(dá)式為：

式中，μx、μy分別表示圖像X 和Y 的均值，σx、σy分別表示圖像X 和Y 的標(biāo)準(zhǔn)差，分別表示圖像X 和Y 的方差。常數(shù)C1、C2用來防止分母可能等于零時(shí)，分式出現(xiàn)不穩(wěn)定現(xiàn)象。該研究采用SSIM比較三種不同算法重建得到圖像的結(jié)構(gòu)相似性以評(píng)價(jià)算法的優(yōu)缺點(diǎn)，如圖6 所示。從圖6 可知，針對(duì)相同尺寸的圖像，濾波反投影重建得到圖像最接近原始圖像。而直接反投影未對(duì)采集數(shù)據(jù)進(jìn)行預(yù)處理，使得重建圖像存在“星芒狀偽影”與原始圖像偏離程度最大。

圖6 比較三種算法重建得CT圖像的結(jié)構(gòu)相似性

3.4 基于平均梯度的圖像質(zhì)量評(píng)價(jià)

平均梯度（Average Gradient）是指圖像邊界點(diǎn)附近的灰度變化差異，可用來表示圖像清晰度[16]。平均梯度反映圖像中感興趣區(qū)域的細(xì)節(jié)差異與紋理特征變化[17-18]。一般來說圖像平均梯度越大，表達(dá)信息越全面，圖像清晰度越高。因此，其值越大越好。圖像平均梯度計(jì)算公式如下：

式中，M×N表示圖像矩陣，ΔIx與ΔIy分別表示采用水平方向與垂直方向的結(jié)構(gòu)元素對(duì)圖像形態(tài)學(xué)處理后得到的梯度信息。通過平均梯度評(píng)價(jià)三種不同算法重建后的CT 圖像的質(zhì)量，如圖7 所示。從圖7 可知，濾波反投影和迭代重建算法在不同圖像矩陣中的平均梯度最好，這說明濾波反投影法重建得到圖像能夠較好反映微小細(xì)節(jié)反差與邊緣特征變化，而直接反投影算法得到的CT 圖像的平均梯度最低，說明該算法對(duì)于邊緣特征變化的反應(yīng)較差。

圖7 比較三種算法重建得CT圖像的平均梯度

4 結(jié)束語

該研究采用了三種不同算法重建CT 圖像，并比較重建圖像的質(zhì)量差異。從運(yùn)行時(shí)間分析可知，對(duì)于相同尺寸的圖像，濾波反投影算法的重建時(shí)間短，重建過程中圖像損失的細(xì)節(jié)和失真都較少，且接近原始圖像。直接反投影算法和濾波反投影算法的運(yùn)行時(shí)間相近，但是前者重建的圖像效果較差。迭代重建算法與其他兩種算法比較，耗時(shí)最長。從算法的復(fù)雜程度和收斂程度進(jìn)行分析可知，三種算法都收斂。但是，迭代重建算法為了使重建圖像收斂至初始圖像，迭代次數(shù)必須接近無窮。這必然會(huì)增加計(jì)算機(jī)工作的復(fù)雜性和運(yùn)行時(shí)間。從抗噪性分析可知，直接反投影算法的峰值信噪比較低，說明該算法重建得到的圖像魯棒性較差。雖然迭代重建算法的峰值信噪比比直接反投影算法高但是其收斂慢，圖像邊緣和細(xì)節(jié)存在不同程度的模糊。從結(jié)構(gòu)相似性和平均梯度分析可知，濾波反投影算法的結(jié)構(gòu)相似性和平均梯度比另外兩種算法好，說明該算法重建得到CT 圖像的細(xì)節(jié)分辨最清楚，重建得到CT 圖像與原始圖像最接近。因而在實(shí)際應(yīng)用中，針對(duì)相同尺寸的CT 圖像重建，濾波反投影算法的運(yùn)行時(shí)間最短，重建圖像失真小，抗噪性強(qiáng)，分辨率高。對(duì)于圖像重建算法的未來研究重點(diǎn)應(yīng)該是提高重建算法的抗噪性，縮短運(yùn)行時(shí)間，同時(shí)降低算法重建過程的復(fù)雜性，加快算法收斂。