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

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

（湖南醫藥學院 醫學院，湖南懷化 418000）

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

1 算法介紹

1.1 直接反投影算法

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

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

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

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

圖1 中心切片定理示意圖

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

式中，F(u,v)表示投影數據的二維頻譜圖。

1.2 濾波反投影算法

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

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

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

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

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

對應的空域離散形式為：

S-L 濾波函數重建圖像的波動幅值較低，抗噪聲性能強和細節分辨率高。S-L 濾波函數對低頻段數據響應效果好，而對于高頻段數據響應就不如R-L 濾波函數。

1.3 迭代重建算法

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

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

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

2 CT圖像的重建結果

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

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

3 圖像質量的評價

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

3.1 圖像重建時間

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

圖4 三種算法重建CT圖像所需的運行時間

3.2 基于峰值信噪比的圖像質量評價

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

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

3.3 基于結構相似性的圖像質量評價

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

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

圖6 比較三種算法重建得CT圖像的結構相似性

3.4 基于平均梯度的圖像質量評價

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

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

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

4 結束語

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