基于radon算法的平行束形狀CT系統(tǒng)參數(shù)標(biāo)定及成像

趙吉祥，姚興貴，李國(guó)強(qiáng)，許鵬飛

（安徽農(nóng)業(yè)大學(xué) 工學(xué)院，安徽 合肥 230036）

CT(Computed Tomography)可以在不破壞樣品的情況下，利用樣品對(duì)射線能量的吸收特性對(duì)生物組織和工程材料的樣品進(jìn)行斷層成像，由此獲取樣品內(nèi)部的結(jié)構(gòu)信息.CT系統(tǒng)安裝時(shí)往往存在誤差，造成偽影，從而影響成像質(zhì)量，因此需要對(duì)安裝好的CT系統(tǒng)進(jìn)行參數(shù)標(biāo)定，并據(jù)此對(duì)未知結(jié)構(gòu)的樣品進(jìn)行成像.因?yàn)镃T系統(tǒng)工作過程中由于誤差易造成成像質(zhì)量下降，所以需求一種準(zhǔn)確合理的標(biāo)定方式顯得尤為重要.

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

圖1 標(biāo)定模板的幾何信息

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

本文主要符號(hào)說明如下：

序號(hào) 符號(hào) 符號(hào)說明1 I0 X射線的入射強(qiáng)度2 I1 穿過介質(zhì)后射線強(qiáng)度3 λ線衰減系數(shù)4射線穿過的對(duì)象厚度5 R小圓半徑X 6橢圓長(zhǎng)半軸7 b橢圓短半軸a 8探測(cè)器單元間間距9 ρ介質(zhì)內(nèi)部的吸收強(qiáng)度d 10 （X，Y） 轉(zhuǎn)轉(zhuǎn)中心坐標(biāo)

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

1 投影模型

1.1 模型的準(zhǔn)備

首先根據(jù)問題數(shù)據(jù)分析可知所有的吸收能力為1的點(diǎn)構(gòu)成的圖形即為圖1（模板示意圖）；通過MATLAB編程畫出圖形如下圖2所示.通過上述的分析可知，探測(cè)器上所接收到的射線能量信息是經(jīng)介質(zhì)吸收衰減后的相應(yīng)的射線能量信息.而由X射線成像的原理[1]，當(dāng)入射強(qiáng)度為I0時(shí)，X射線通過厚度為x的物體輸出的X射線強(qiáng)度I1與該物體的衰減系數(shù)λ有關(guān)，即

圖2 投影分布圖

結(jié)合數(shù)據(jù)分析可知，本文中X射線通過厚度為x的物體輸出的X射線強(qiáng)度I1應(yīng)為

其中I1為穿過介質(zhì)后射線強(qiáng)度，I0為X射線的入射強(qiáng)度，λ為線衰減系數(shù)（吸收效率），x為射線穿過的對(duì)象厚度.

1.2 模型的建立

因?yàn)閄射線能量的損失量與射線通過介質(zhì)內(nèi)部距離x有關(guān)，因此應(yīng)求出整個(gè)發(fā)射-接收系統(tǒng)在旋轉(zhuǎn)過程中每條射線穿過介質(zhì)的距離，做出幾何關(guān)系圖[2]如下圖3所示，建立投影模型如下，以正方形托盤的左下角為坐標(biāo)原點(diǎn)B（0，0）建立直角坐標(biāo)系.

圖3 投影幾何關(guān)系圖

其中T為接收信息數(shù)據(jù)，右上角直角坐標(biāo)的坐標(biāo)原點(diǎn)為 A（G，H），橢圓長(zhǎng)軸 40，短軸 15，故有：

其中ρ為介質(zhì)內(nèi)部的吸收強(qiáng)度，θ為以（G，H）為坐標(biāo)原點(diǎn)的坐標(biāo)系下射線與橫軸的夾角，即為初始角度，α為旋轉(zhuǎn)中心坐標(biāo)系下旋轉(zhuǎn)中心向射線做切線時(shí)，此切線與橫軸之間的夾角.

1.3 模型的求解

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

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

2 驗(yàn)證模型

2.1 探測(cè)器單元之間距離的檢驗(yàn)

通過Excel編程將相關(guān)數(shù)據(jù)進(jìn)行分類標(biāo)定，得到如下圖4所示的圖形.然后結(jié)合整個(gè)發(fā)射-接收系統(tǒng)在旋轉(zhuǎn)過程中規(guī)律及數(shù)據(jù)分析方法求解問題.

通過對(duì)圖形及數(shù)據(jù)的分析可知，當(dāng)小圓沒被大圓遮擋時(shí)說明小圓一直完全暴露在平行射線之間，而對(duì)應(yīng)的接收器上面的探測(cè)器個(gè)數(shù)F可由數(shù)據(jù)分析得出.根據(jù)小圓直徑R，可通過得出探測(cè)器單元之間的距離.

圖4 Excel表示數(shù)據(jù)圖

2.2 旋轉(zhuǎn)角度間距的檢驗(yàn)

切線轉(zhuǎn)角指整個(gè)發(fā)射-接收系統(tǒng)繞某固定的旋轉(zhuǎn)中心逆時(shí)針旋轉(zhuǎn)過程中，大圓剛好開始遮擋小圓時(shí)的切線轉(zhuǎn)至大圓剛好完全遮擋小圓切線的轉(zhuǎn)角.

通過對(duì)圖4以及相關(guān)數(shù)據(jù)分析可知，可通過建立旋轉(zhuǎn)方向模型較簡(jiǎn)單求解該CT系統(tǒng)使用的X射線的180個(gè)方向.在此模型中需進(jìn)行假設(shè)如下，整個(gè)發(fā)射-接收系統(tǒng)在旋轉(zhuǎn)過程每次都是轉(zhuǎn)過均勻的角度，因?yàn)檎w接收器長(zhǎng)度較小，又分為526個(gè)小的單位等距探測(cè)器，所以每?jī)蓚€(gè)探測(cè)器之間的距離均較小，所以可不考慮因小圓兩邊射線剛好與小圓相切，因此由此所造成的誤差可忽略不計(jì).通過分析旋轉(zhuǎn)規(guī)律可做出如下圖5所示的兩條切線.

圖5 切線圖

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

與模型一求解結(jié)果相近，同時(shí)驗(yàn)證模型一的合理性.

2.2 反投影重建的檢驗(yàn)

通過模型一的求解結(jié)果，反投影重建出相應(yīng)的圖像，并將其與標(biāo)準(zhǔn)模板示意圖進(jìn)行比較.根據(jù)模型一的求解結(jié)果，得到重建圖像如圖6（b），通過比較重建前后的圖像，二者基本一致，說明模型一的求解結(jié)果具有較好的精度.

圖6 標(biāo)準(zhǔn)模板重建前后的圖像

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

3.1 模型的準(zhǔn)備

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

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

圖7 原始數(shù)據(jù)圖

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

即表示每一點(diǎn)處的吸收率，其中I1為穿過介質(zhì)后射線強(qiáng)度，I0為X射線的入射強(qiáng)度，λ為線衰減系數(shù)，X為射線穿過的對(duì)象厚度.

3.2 模型的求解

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

圖8 反投影重建圖

4 當(dāng)投影數(shù)據(jù)存在噪聲情況時(shí)

根據(jù)數(shù)據(jù)分析，首先利用MATLAB軟件編程使數(shù)據(jù)信息圖像化.做出的圖像分別如下圖9所示，圖為MATLAB編程求得的相應(yīng)的整個(gè)發(fā)射-接收系統(tǒng)繞某固定的旋轉(zhuǎn)中心逆時(shí)針旋轉(zhuǎn)過程中介質(zhì)內(nèi)部有射線經(jīng)過后的圖像.

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

圖9 原始數(shù)據(jù)圖

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

（1）選擇合適的濾波器.經(jīng)查相關(guān)的文獻(xiàn)資料可知，本文選用Hamming濾波器進(jìn)行濾波反投影重建算法的濾波器.

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

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

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

得到的該位置介質(zhì)的相關(guān)信息，得到重建圖像如下所示.

圖10 反投影重建圖

5 模型的誤差分析

5.1 建立圖像分析模型求解探測(cè)器單元間距d時(shí)，直接根據(jù)圖形及數(shù)據(jù)進(jìn)行分析求解，可能會(huì)由于數(shù)據(jù)分析不夠準(zhǔn)確，導(dǎo)致結(jié)果誤差較大.

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

6 模型的評(píng)價(jià)

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

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

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