基于radon算法的平行束形狀CT系統參數標定及成像

趙吉祥，姚興貴，李國強，許鵬飛

（安徽農業大學 工學院，安徽 合肥 230036）

CT(Computed Tomography)可以在不破壞樣品的情況下，利用樣品對射線能量的吸收特性對生物組織和工程材料的樣品進行斷層成像，由此獲取樣品內部的結構信息.CT系統安裝時往往存在誤差，造成偽影，從而影響成像質量，因此需要對安裝好的CT系統進行參數標定，并據此對未知結構的樣品進行成像.因為CT系統工作過程中由于誤差易造成成像質量下降，所以需求一種準確合理的標定方式顯得尤為重要.

例如，如下圖1所示的正方形托盤上放置兩個均勻固體介質，組成的標定模板，模板的幾何信息如圖中所示.

圖1 標定模板的幾何信息

需要根據正方形托盤上放置的兩個均勻固體介質組成的標定模板及其接收信息等數據，確定CT系統旋轉中心在正方形拖盤中的位置、探測器單元之間的距離以及該CT系統使用的X射線的180個方向,并進行以下其他問題的研究.

本文主要符號說明如下：

序號 符號 符號說明1 I0 X射線的入射強度2 I1 穿過介質后射線強度3 λ線衰減系數4射線穿過的對象厚度5 R小圓半徑X 6橢圓長半軸7 b橢圓短半軸a 8探測器單元間間距9 ρ介質內部的吸收強度d 10 （X，Y） 轉轉中心坐標

首先根據模板及其接收信息，確定CT系統旋轉中心在正方形托盤中的位置、探測器單元之間的距離以及該CT系統使用的X射線的180個方向.所有數據只有0和1，其中每一點的數值反映了該點的吸收強度.經分析可知，介質內部每點的吸收強度均為1，即介質內每一點處對X射線的吸收強度相同，為均勻的.所以射線能量在介質內部每點處的衰減能量相同，射線在介質內部所損失的總的能量僅與射線照射時在介質內所通過的距離有關.整個發射-接收系統繞某固定的旋轉中心逆時針旋轉180次，對每一個X射線方向，在具有512個等距單元的探測器上測量經位置固定不動的二維待檢測介質吸收衰減后的射線能量，并經過增益等處理后得到180組接收信息.通過分析可知，數據為0點處表示X射線從發射到接收過程中并未通過介質，而數據不為0點處表示X射線在發射-接收系統繞固定旋轉中心旋轉過程中通過介質，有部分能量被介質吸收.所以可根據excel表示出所有數據不為0的點，即表示整個發射-接收系統在旋轉過程中所通過的所有點.

1 投影模型

1.1 模型的準備

首先根據問題數據分析可知所有的吸收能力為1的點構成的圖形即為圖1（模板示意圖）；通過MATLAB編程畫出圖形如下圖2所示.通過上述的分析可知，探測器上所接收到的射線能量信息是經介質吸收衰減后的相應的射線能量信息.而由X射線成像的原理[1]，當入射強度為I0時，X射線通過厚度為x的物體輸出的X射線強度I1與該物體的衰減系數λ有關，即

圖2 投影分布圖

結合數據分析可知，本文中X射線通過厚度為x的物體輸出的X射線強度I1應為

其中I1為穿過介質后射線強度，I0為X射線的入射強度，λ為線衰減系數（吸收效率），x為射線穿過的對象厚度.

1.2 模型的建立

因為X射線能量的損失量與射線通過介質內部距離x有關，因此應求出整個發射-接收系統在旋轉過程中每條射線穿過介質的距離，做出幾何關系圖[2]如下圖3所示，建立投影模型如下，以正方形托盤的左下角為坐標原點B（0，0）建立直角坐標系.

圖3 投影幾何關系圖

其中T為接收信息數據，右上角直角坐標的坐標原點為 A（G，H），橢圓長軸 40，短軸 15，故有：

其中ρ為介質內部的吸收強度，θ為以（G，H）為坐標原點的坐標系下射線與橫軸的夾角，即為初始角度，α為旋轉中心坐標系下旋轉中心向射線做切線時，此切線與橫軸之間的夾角.

1.3 模型的求解

通過MATLAB編程求解上述模型，任意選取五次旋轉的探測器接收的數值，然后求解此非線性方程組，方程組如下所示：

通過求解由公式（5）～（9）構成的非線性方程組，結果僅為選取一部分得出的結果，因此需仿照上述算法求出多組的值，取平均值作為求解X射線的初始入射角度和探測器單元之間的間距聯合幾個不同入射角度下的投影數據求解出X射線的旋轉角度和旋轉中心.通過求解多個相鄰入射角度之間差值的均值，進一步得到該CT系統使用的X射線的180個方向.

2 驗證模型

2.1 探測器單元之間距離的檢驗

通過Excel編程將相關數據進行分類標定，得到如下圖4所示的圖形.然后結合整個發射-接收系統在旋轉過程中規律及數據分析方法求解問題.

通過對圖形及數據的分析可知，當小圓沒被大圓遮擋時說明小圓一直完全暴露在平行射線之間，而對應的接收器上面的探測器個數F可由數據分析得出.根據小圓直徑R，可通過得出探測器單元之間的距離.

圖4 Excel表示數據圖

2.2 旋轉角度間距的檢驗

切線轉角指整個發射-接收系統繞某固定的旋轉中心逆時針旋轉過程中，大圓剛好開始遮擋小圓時的切線轉至大圓剛好完全遮擋小圓切線的轉角.

通過對圖4以及相關數據分析可知，可通過建立旋轉方向模型較簡單求解該CT系統使用的X射線的180個方向.在此模型中需進行假設如下，整個發射-接收系統在旋轉過程每次都是轉過均勻的角度，因為整體接收器長度較小，又分為526個小的單位等距探測器，所以每兩個探測器之間的距離均較小，所以可不考慮因小圓兩邊射線剛好與小圓相切，因此由此所造成的誤差可忽略不計.通過分析旋轉規律可做出如下圖5所示的兩條切線.

圖5 切線圖

上圖中L2直線表示平行直線的最邊上一條射線剛好與小圓相切，小圓剛開始被大圓遮擋，L1直線則表示轉過幾次之后完全從小圓離開，小圓已經完全被大圓遮擋.易知Ω=11°，由數據可知旋轉次數K范圍大致在10到11之間，即10≤K≤11，取平均值K=10.5，因此

與模型一求解結果相近，同時驗證模型一的合理性.

2.2 反投影重建的檢驗

通過模型一的求解結果，反投影重建出相應的圖像，并將其與標準模板示意圖進行比較.根據模型一的求解結果，得到重建圖像如圖6（b），通過比較重建前后的圖像，二者基本一致，說明模型一的求解結果具有較好的精度.

圖6 標準模板重建前后的圖像

3 平行束濾波反投影重建算法模型

3.1 模型的準備

平行束濾波反投影重建算法[3]的具體過程是先把探測器上獲得的投影數據與濾波器函數進行卷積運算，得到各方向上卷積濾波后的投影數據，然后再把它們進行反投影，得到每一矩陣單元的CT值，進而得到相關被掃描物體的數據，形象信息.

基于平行束濾波反投影重建算法模型，首先本文利用MATLAB編程對數據進行圖像處理，圖像如下圖7所示.

圖7 原始數據圖

然后分析求解該未知介質的吸收率，算法步驟如下：對式（1）兩邊取對數得，然后求解可得：

即表示每一點處的吸收率，其中I1為穿過介質后射線強度，I0為X射線的入射強度，λ為線衰減系數，X為射線穿過的對象厚度.

3.2 模型的求解

通過以上模型的建立，利用MATLAB編程求解，得到重建后的圖像（圖8）如下圖所示.

圖8 反投影重建圖

4 當投影數據存在噪聲情況時

根據數據分析，首先利用MATLAB軟件編程使數據信息圖像化.做出的圖像分別如下圖9所示，圖為MATLAB編程求得的相應的整個發射-接收系統繞某固定的旋轉中心逆時針旋轉過程中介質內部有射線經過后的圖像.

由圖9分析可知，數據有噪聲存在，應該首先進行濾波去噪聲處理，采用卷積濾波及反投影重建算法進行編程求解.

圖9 原始數據圖

具體平行束濾波反投影重建算法如下[3]：

（1）選擇合適的濾波器.經查相關的文獻資料可知，本文選用Hamming濾波器進行濾波反投影重建算法的濾波器.

（2）把在固定角度ω下測得的投影P(Xr,ω)與R-L濾波器進行卷積濾波，得到濾波后的投影f(Xr,ω).

（3）對于每一個ω，把經R-L濾波器濾波后的投影f(Xr,ω)反投影于滿 Xr=rcos(Ψ-ω)的射線上的所有各點（r,Ψ）.

（4）將步驟三中的反投影值對所有0≤ω≤2π進行累加，得到重建后的圖像數據.

得到的該位置介質的相關信息，得到重建圖像如下所示.

圖10 反投影重建圖

5 模型的誤差分析

5.1 建立圖像分析模型求解探測器單元間距d時，直接根據圖形及數據進行分析求解，可能會由于數據分析不夠準確，導致結果誤差較大.

5.2 建立旋轉方向模型求解該CT系統使用的X射線的180個方向時，也是僅根據相關圖像及相關數據進行分析求解，同時模型假設整個發射-接收系統繞某固定的旋轉中心逆時針旋轉過程中每次轉過均勻角度，這與實際旋轉方式相差也可能存在一定誤差，導致結果存在誤差.同時僅考慮從小圓剛開始被大圓遮擋到小圓完全被大圓遮擋這段短暫時間的旋轉過程，對最終的旋轉結果也存在一定的影響.

6 模型的評價

6.1 本文緊緊圍繞著整個發射-接收系統繞某固定的旋轉中心逆時針旋轉過程中的特點，緊密圍繞圖像特點，結合接收器接收到的數據特點以及X射線衰減方程規律.由于介質內部每點處對射線的吸收能力即吸收率相同，所以僅需考慮距離對X射線衰減的影響，即射線穿過介質的距離越長射線能量損失越多.于是建立投影模型，使問題得到簡化，模型簡單易懂，具有較好的代表性.

6.2 在利用投影重建算法時[4]，反投影圖像會引入星狀偽跡，即原來圖像中密度為零的點，重建后不一定為零，使圖像失真.在利用濾波反投影重建算法時，我們可以在投影重建以前把投影數據先行修正，再把修正后的投影數據進行反投影運算而求出無偽跡的圖像.其中，圖像的重建精度依賴于具體的濾波函數.

6.3 雖然濾波反投影算法[5]能獲得較好的重建質量，但是CT重建的數據量龐大，使得圖像重建的計算非常耗時，尤其是對于高分辨率圖像重建，因此提高該算法的重建速度是一個重要的研究方向.