基于優(yōu)化的CT系統(tǒng)標定與濾波反投影圖像重建

李江珊，鄒義華，姚 竹，王志勇

(1.電子科技大學 信息與通信工程學院，四川 成都 611731；2.電子科技大學 數(shù)學科學學院，四川 成都 611731)

本文是對2017年全國大學生數(shù)學建模競賽題目的求解，題目和數(shù)據(jù)由組委會提供[1]?，F(xiàn)有某類型的CT系統(tǒng)，平行的X射線入射到樣本上，經(jīng)過衰減后在另一端被探測器接收，探測器有512個等距的探測單元，可測量接收到的X射線的能量。將整個發(fā)射探測系統(tǒng)繞一個固定的旋轉(zhuǎn)中心旋轉(zhuǎn)180個方向，經(jīng)增益處理后可得到180組數(shù)據(jù)[2]。為減小CT系統(tǒng)的誤差，論文建立了基于幾何分析和優(yōu)化的標定模型。在此基礎上，通過濾波反投影法對圖像進行了重建。最后，利用計算機模擬分析了標定模板的精度和穩(wěn)定性，并進行了新模板設計。

1 模型的建立與求解

1.1 基于直線方程和曲線方程的CT系統(tǒng)標定模型

CT系統(tǒng)的標定隨著CT系統(tǒng)的發(fā)展不斷發(fā)展，其大致可分為兩類，一類是自標定方法，另一類則是基于標定模板的標定方式[3]。自標定方式是指對CT系統(tǒng)參數(shù)進行標定時僅根據(jù)被檢測物體本身的投影。而不依靠其他數(shù)據(jù)。該方法增加了計算量，無法獲得實際分辨率并且具有局限性，所以使用較少。而基于標定模板的標定方式是通過在特定的運動軌跡下采集投影，然后與CT系統(tǒng)參數(shù)建立聯(lián)系，從而得到CT系統(tǒng)的幾何參數(shù)。本文就是建立在第一代CT系統(tǒng)的基礎上，利用基于標定模板的標定方式對系統(tǒng)進行標定。標定的具體方法有基于統(tǒng)計特征的提取法，最小二乘法，自適應遺傳算法等[4]，其中自適應遺傳算法涉及的參數(shù)較多，計算量較大[5]，并且遺傳算法不能保證得到的解為最優(yōu)解；最小二乘法則是需要建立多個目標規(guī)劃模型[6]，較為復雜。本文將問題分割為確定增益系數(shù)和探測器單元間距、確定CT系統(tǒng)X射線的180個方向、確定旋轉(zhuǎn)中心位置等多個模塊，相較于自適應遺傳算法減少了參數(shù)，降低了計算難度，并且在計算CT系統(tǒng)方向部分直接利用幾何代數(shù)進行求解，相較于最小二乘法，也更為簡單。

為了更加直觀地觀察題目給出接收信息，畫出接收信息隨組數(shù)、探測器單元的變化曲線，如圖1所示。顯然，中間類似花瓶形狀的凸起是由模板中的橢圓導致的，彎曲曲線是由模板中的圓導致的?；谝陨险J識，下面給出基于直線方程和曲線方程CT系統(tǒng)標定模型。

1.1.1 確定增益系數(shù)和探測器單元間距

根據(jù)文獻[7]得到比賽數(shù)據(jù)附件2[1]中的接收信息η與相應X射線通過介質(zhì)的長度l成正比：

圖1 接收信息隨組數(shù)、探測器單元的變化曲線

η=Al

(1)

取出位于圓心同一側的經(jīng)過了圓介質(zhì)的3條相鄰X射線的接收信息，假設圓心到其中一條射線的距離為d0，d表示相鄰探測器單元之間的距離，R表示圓的半徑，可得：

(2)

(3)

(4)

1.1.2 確定CT系統(tǒng)X射線的180個方向

分析得到射線的方向和探測器的初始位置大致如圖2所示，以正方形托盤的中心為原點建立如圖2所示的坐標系。

推知第k個方向上編號為n的射線方程式為：

ykn=(tanθk)x+bkn

(5)

圖2 CT系統(tǒng)示意圖

為了求出各組數(shù)據(jù)對應的平行射線與x軸正方向的夾角，取比賽數(shù)據(jù)附件2中兩條經(jīng)過橢圓的射線對應的接收信息值ηk4、ηk5，根據(jù)式(1)和橢圓弦長公式可以得到以下兩組方程：

(6)

(7)

其中，

ik、jk分別表示第k組數(shù)據(jù)中經(jīng)過橢圓的兩條射線的編號。

1.1.3 確定旋轉(zhuǎn)中心位置

根據(jù)直線繞點旋轉(zhuǎn)的性質(zhì)可知，旋轉(zhuǎn)中心到旋轉(zhuǎn)過后直線的距離與旋轉(zhuǎn)中心到未旋轉(zhuǎn)直線的距離相等。

(8)

式中，(x0,y0)表示旋轉(zhuǎn)中心坐標，ki、bi分別表示第i條射線方程的斜率和截距。

設定旋轉(zhuǎn)中心位于正方形托盤內(nèi)，旋轉(zhuǎn)中心到直線的理論距離值必不大于正方形托盤的內(nèi)最長線段(對角線)的長度，即得到優(yōu)化模型的約束條件：

(9)

式中，(x0,y0)表示旋轉(zhuǎn)中心坐標，r表示旋轉(zhuǎn)中心到直線的理論距離值。求解優(yōu)化模型即可得到旋轉(zhuǎn)中心的位置。

1.1.4 模型結果與分析

將增益系數(shù)值和探測器單元間距代入模型求解，求解出在180個方向上對應的180組平行射線與水平方向夾角結果，將其繪制成柱狀圖形式如圖3所示。平行射線與圖2所示坐標系的x軸正方向的夾角范圍為(-60.385 3°,118.695 9°)，相鄰的兩個方向的角度差集中在1°附近。

利用遺傳算法得到旋轉(zhuǎn)中心坐標(在圖2所示的坐標系中)的10組數(shù)據(jù)，將這10組數(shù)據(jù)的均值作為最終的旋轉(zhuǎn)中心坐標，則最后確定的旋轉(zhuǎn)中心坐標為(-9.268 9,6.274 3)，位于第二象限，較為靠近原點。

圖3 CT系統(tǒng)X射線的180個方向示意圖

1.2 基于濾波反投影法的CT圖形重建模型

根據(jù)文獻[2]采用濾波反投影模型重建附件3[1]的CT圖形。

論文的濾波反投影法分為4個部分：濾波器選擇、線性內(nèi)插，反投影重建和確定吸收率。

1.2.1 濾波器選擇

論文分別采用R-L濾波器[8]、S-L濾波器[9]、R-L和S-L混合濾波器、R-L和NEW混合濾波器進行濾波處理，比較濾波效果。

1.2.2 線性內(nèi)插

根據(jù)文獻[10-11]，本文利用每個單元格集合中心的坐標來表征整個單元格，記為(i,j)，但(i,j)可能在n0射線和n0+1射線之間，因此需要進行線性內(nèi)插預測經(jīng)過點(i,j)的射線的接收信息，即：

p1(xij,θm)=

(10)

其中

(11)

Δ=xij-(n0+N0)d

(12)

式中，n0表示512條射線中，位于點(i,j)左側且距離點(i,j)最近的射線編號；Δ表示(i,j)到n0射線的距離。由此即可求出在任意方向上虛擬的探測器單元接收到的經(jīng)過點(i,j)的射線信息。

1.2.3 反投影重建

(13)

式中，p1(xij,θm)=p(xij,θm)f(xij)，p(xij,θm)表示接收信息，f(xij)表示濾波沖激響應函數(shù)。根據(jù)上述分析即可計算得到規(guī)模為256×256的單元的具體衰減系數(shù)。

1.2.4 確定吸收率

論文認為衰減系數(shù)與吸收率成正比[17]，通過建立優(yōu)化模型來求得最優(yōu)的系數(shù)值。

利用濾波反投影法對附件2[1]數(shù)據(jù)進行重建，得到重建后的吸收率衰減系數(shù)矩陣與附件1中的實際數(shù)據(jù)進行比較。采用實際吸收率與理論吸收率的之差的絕對值之和為目標函數(shù)，則優(yōu)化模型的目標是：

(11)

式中，k0(i,j)表示單元格(i,j)的實際吸收率。該優(yōu)化模型是無約束優(yōu)化模型。

結果如表1所示，R-L和NEW混合濾波器的濾波效果最好，因此采用R-L和NEW混合濾波器來重構附件3[1]和附件5[1]的圖像。

表1 不同濾波器的濾波效果

利用濾波反投影模型得到附件3[1]和附件5[1]數(shù)據(jù)重建后樣本的幾何形狀和在正方形托盤中的位置分別如圖4所示，其中白色部分表示有介質(zhì)。而附件4[1]中10個位置處的吸收率具體值如表2所示。

圖4 附件3和附件5樣本形狀與位置

表2 10個位置的吸收率

1.3 標定模型精度、穩(wěn)定性分析與模板改進

1.3.1 標定模型的精度分析

模型的精度是指計算值與真實值的接近程度。論文通過模擬已知的實際情況，設定發(fā)射-接收系統(tǒng)的旋轉(zhuǎn)中心，相鄰探測器單元間距和180個方向角為θk，利用計算機模擬出接收信息。再根據(jù)模型2.1利用模擬得到的接收信息進行CT系統(tǒng)標定，將計算結果與實際結果進行對比，得到探測器單元之間的間距的相對誤差約為3×10-6，X射線與水平方向夾角與實際值的誤差平方和值為0.020 9。根據(jù)上述結果，可以看出參數(shù)標定模型的精度很高。

1.3.2 標定模型的穩(wěn)定性分析

因為接收信息的獲得過程是未知的，所以接收信息具有不確定性，需要檢驗外界擾動對標定結果的影響，即分析標定模型的穩(wěn)定性。為了描述接收信息的不確定性，給每個接收信息加上一個在區(qū)間[0,0.01]之間的隨機擾動。根據(jù)模型2.1利用新的接收信息來計算增益系數(shù)、間距等參數(shù)。通過未擾動數(shù)據(jù)計算的結果對比，旋轉(zhuǎn)中心的橫縱坐標的絕對誤差分別為-0.128 8、0.057 5，相對誤差為1.41%、0.92%。探測器單元之間的間距的相對誤差為-0.83%，增益系數(shù)A′的相對誤差為0.24%。計算得到X射線與水平方向夾角與實際值的誤差平方和值為0.006 3。根據(jù)上述結果，可以看出參數(shù)標定模型的穩(wěn)定性較好。

1.3.3 標定模板與標定模型的改進

模板的尺寸和位置如圖5所示。

(a)長方-圓模板

(b)長方-橢圓-圓模板

利用模型2.1確定增益系數(shù)、探測器單元間距、CT系統(tǒng)X射線的180個方向和旋轉(zhuǎn)中心位置的方法，重新建立相關方程，并對新標定模板重新分析精度和穩(wěn)定性。

對于長方形-圓形模板，其計算增益系數(shù)和探測器單元間距的方法與原模板相同，則與原模板相比，增益系數(shù)和探測器單元間距的精度和穩(wěn)定性不變。而其算出的X射線與水平方向夾角與實際值的誤差平方和值為0.229 9，比原模板的誤差平方和大，這是因為采用新模板計算得到在90°附近的夾角值與實際值相差很大，但在其他位置計算值與實際值相差十分小，有些甚至沒有誤差。因此，可以認為新模板的精度有所提升。給模擬的接收信息增加一個在區(qū)間[0,0.01]之間的隨機擾動，利用擾動后的數(shù)據(jù)對CT系統(tǒng)重新進行標定，得到其與未擾動結果之間的誤差平方和為0.002 8，比原模板小，即模型的穩(wěn)定性更高。

對于長方形-橢圓-圓形模板，其計算增益系數(shù)和探測器單元間距的方法與原模板相同，則與原模板相比，增益系數(shù)和探測器單元間距的精度和穩(wěn)定性不變。其計算出的X射線與水平方向夾角與實際值的誤差平方和值為0.006 8，比原模板和長方形-圓模板的誤差平方和小，說明新模板的精度比原模板和長方形-圓模板的高。加上一個在區(qū)間[0,0.01]之間的隨機擾動，利用擾動后的數(shù)據(jù)對CT系統(tǒng)重新進行標定，得到其與未擾動結果之間的誤差平方和為0.038 4，比原模板大，即模型的穩(wěn)定性更低。

2 結束語

模型2.1先通過圓確定增益和相鄰探測單元的間隔，減少優(yōu)化模型的決策變量，簡化模型，便于計算。模型2.2利用多種濾波器對附件二求解吸收系數(shù)，通過對比選擇了較優(yōu)的濾波器處理附件三和附件五。

模型2.2中預測未知點的接收信息時，可以考慮多項式插值、樣條插值等方法，提高圖形精度。