王譽儒，楊玉聰，官 莊，王志勇

(1.電子科技大學(xué) 英才實驗學(xué)院,四川 成都 611731；2.電子科技大學(xué) 數(shù)學(xué)科學(xué)學(xué)院,四川 成都 611731)

電子計算機斷層掃描(CT)系統(tǒng)成像如今在醫(yī)療領(lǐng)域應(yīng)用廣泛,例如鼻骨骨折診斷[1],早期股骨頭缺血壞死診斷[2],術(shù)后輔助治療[3]等,CT掃描在電子、材料和航空航天等領(lǐng)域也有著大量的應(yīng)用[4]。而在CT系統(tǒng)成像中,參數(shù)標(biāo)定和圖像重建是主要的問題[5-8]。

本文主要討論一種典型的二維CT系統(tǒng),平行的X射線垂直入射于具有512個等距單元的探測器,探測器單元視作等距排列的接收點；X射線的發(fā)射器和探測器相對位置固定不變,整個發(fā)射-接收系統(tǒng)繞某固定的旋轉(zhuǎn)中心逆時針旋轉(zhuǎn)180次,以測量射線能量,并經(jīng)過增益等處理后得到180組接收信息。

CT系統(tǒng)安裝時往往存在誤差,因此需要借助于已知結(jié)構(gòu)的樣品(稱為模板)標(biāo)定CT系統(tǒng)的參數(shù),并據(jù)此對未知結(jié)構(gòu)的樣品進行成像。要求依據(jù)接收信息,確定CT系統(tǒng)旋轉(zhuǎn)中心在正方形托盤中的位置、探測器單元之間的距離以及該CT系統(tǒng)使用的X射線的180個方向；在標(biāo)定之后,再依據(jù)得到的標(biāo)定參數(shù),確定兩種未知介質(zhì)在正方形托盤中的位置、幾何形狀和吸收率等信息,并給出10個具體位置的吸收率。

1 基于接收信息的參數(shù)標(biāo)定模型

要實現(xiàn)對CT系統(tǒng)參數(shù)標(biāo)定,就是要對旋轉(zhuǎn)角度、探測器間距和旋轉(zhuǎn)中心的標(biāo)定。

1.1 旋轉(zhuǎn)角度的標(biāo)定模型

對于旋轉(zhuǎn)角度的標(biāo)定,模板在探測器上的投影長度與旋轉(zhuǎn)中心無關(guān)。從接收信息中,提取出不同角度時模板的實際投影長度,并由幾何關(guān)系得到模擬投影長度,若實際與模擬的投影長度在某個角度區(qū)間相似程度較高,則認為該角度區(qū)間即為旋轉(zhuǎn)角度區(qū)間。基于此,建立優(yōu)化模型,將初始旋轉(zhuǎn)角度和旋轉(zhuǎn)間隔角度作為決策變量,將實際與模擬的投影長度間的差異大小作為優(yōu)化目標(biāo),得到投影長度相似度最高的初始角度和間隔角度,從而標(biāo)定旋轉(zhuǎn)角度。

根據(jù)參考文獻[9],可以得到X射線透射過介質(zhì)時的規(guī)律為:

式中,A和A0分別為透射和入射的X射線強度,l為材料厚度,r為吸收系數(shù)[9]。我們認為該吸收系數(shù)的定義與該問題中的 “吸收率”為等價概念。

對于具有不同吸收系數(shù)的非均勻物體,可將式(1)推廣為:

式中,L為X射線穿過物體的路徑,r(x)表示L上某點x處的吸收率。對式(2)變換后,可以轉(zhuǎn)換為:

式中,g為X射線的吸收系數(shù)r(x)在某一方向上的積分值,即為X射線成像得到的數(shù)值。由于該問題中,接收信息為探測介質(zhì)吸收衰減后的射線能量在經(jīng)過增益處理后得到的,故認為探測器的接收信息為:

式中,G為接收信息,c為其增益常數(shù),下文驗證了該定義的合理性。

假設(shè)有一束平行光,從某個θ角度對模板進行照射,可以得到模板的模擬投影長度Sm,若改變?nèi)肷浣嵌葎t可以得到模板投影長度隨角度的變化曲線。給定初始角度,對入射方向連續(xù)逆時針旋轉(zhuǎn)360°,可以得到模板的投影長度變化曲線Sm(θ),θ∈[0°,360°],進而得到所有入射角度下的模擬模板投影長度。

若采用X光照射模板,根據(jù)式(2),若接收信息為0,則說明X射線沒有透過介質(zhì),從而可以提取出類似于光線投射時的模板投影長度,但是由于未知探測器單元的間隔,故用接收信息非0的單元數(shù)量表征投影長度。對于該問題,可以得到180個角度下的接收信息非0單元間隔數(shù)量Nm(θ0＋iβ)-1,i＝0,1, …,179, θ0為初始旋轉(zhuǎn)角度,β為180次等間隔旋轉(zhuǎn)的間隔角度,角度的單位均采用角度制。另外,若初始的角度相同,且認為投影接收平板足夠大(板長大于該角度下模板的投影長度),則在某角度下的投影長度與旋轉(zhuǎn)中心的位置無關(guān)。故可以建立目標(biāo)函數(shù):

若 θ0＋iβ 的值超過 360°則減去 360°后, 再代入模型計算。故優(yōu)化目標(biāo)函數(shù)可以表示為:

通過以上分析可知,初始旋轉(zhuǎn)角度θ0和旋轉(zhuǎn)間隔角度β為決策變量。對于初始旋轉(zhuǎn)角度θ0,由于沒有多余的信息可以約束,僅考慮旋轉(zhuǎn)一周內(nèi)的角度,故對θ0限定如下:

考慮到若179β＞360°時,可能會在同一角度重復(fù)投射,且間隔過大不利于對介質(zhì)的檢測,故對β限定如下:

基于以上分析,對于旋轉(zhuǎn)角度的確定問題,建立單目標(biāo)優(yōu)化模型如下:

式中,Sm(θ0＋iβ) 為平行光線投射模板在 θ0＋iβ角度時的投影值,Nm(θ0＋iβ)為從X射線投射模板在θ0＋iβ角度時接收信息非0的探測器單元數(shù)量,Fi(θ0,β)為優(yōu)化目標(biāo)函數(shù)。通過搜索模型的最優(yōu)解可以得到CT系統(tǒng)的180個旋轉(zhuǎn)角度。

1.2 探測器間距標(biāo)定模型

由1.1節(jié)模型可知,在旋轉(zhuǎn)角度為θ0＋iβ時的投影長度為 Sm(θ0＋iβ),i＝0,1, …,179,而對應(yīng)的接收信息非0探測器單元數(shù)量為Nm(θ0＋iβ),而探測器又是等間距并排的,故可以近似認為探測器的間距為:

式中,sd(i)為在旋轉(zhuǎn)角度為θ0＋iβ時探測器單元間距的近似值,i＝0,1,…,179。

對于180個角度的接收信息,每個角度均能得到一個近似的探測器間距值sd,但由于離散化帶來的誤差,會導(dǎo)致不同入射角度時sd的值不同,但這些值都在一定范圍內(nèi)波動,為了減小誤差,得到較為精確的探測器間距值,這里采用最小二乘法對精確值進行估計[10]。

求解得到:

將Sd作為探測器間距值的求解結(jié)果。

1.3 旋轉(zhuǎn)中心標(biāo)定模型

因為探測器陣列的垂直平分線必然都交匯于旋轉(zhuǎn)中心。對1.2節(jié)模型求解,可以得到探測器陣列長度,結(jié)合1.1節(jié)模型的求解結(jié)果,可得到探測器陣列的位置。根據(jù)探測器近似線段的位置和長度,可以確定線段垂直平分線。由于問題中存在180個旋轉(zhuǎn)角度,每個角度均有一條近似線段的垂直平分線,任意兩條垂直平分線的交點pi均可以作為旋轉(zhuǎn)中心的近似位置,故存在C＝161 10個交點,即i＝1,2,…,161 00。

以托盤幾何中心為原點建立如圖1所示坐標(biāo)系,x軸和y軸分別與托盤上邊緣和左邊緣平行,故近似的旋轉(zhuǎn)中心可以表示為pi＝(xi,yi),i＝1,2,…,161 00。與1.2節(jié)模型類似,仍采用相似的最小二乘法對旋轉(zhuǎn)中心的準確位置進行估計[3],假設(shè)旋轉(zhuǎn)中心的準確位置為P＝(X,Y)。

求解得到:

將P＝(X,Y)作為旋轉(zhuǎn)中心位置的求解結(jié)果。

圖1 旋轉(zhuǎn)中心坐標(biāo)系示意圖

2 未知介質(zhì)的確定模型

對于介質(zhì)吸收率的確定:先用傅立葉變換將問題轉(zhuǎn)換到頻域上進行討論[11],通過尋找接收信息和吸收率在頻域上的關(guān)系,用頻域下的接收信息反演出頻域下的吸收率,從而由傅立葉逆變換反演出吸收率的結(jié)果。

對于介質(zhì)形狀和位置的確定:若某處存在介質(zhì),則該處吸收率必為正值,故根據(jù)吸收率是否為0判斷該處是否存在介質(zhì)。在已經(jīng)求出吸收率的情況下,依據(jù)吸收率的大小和分布實現(xiàn)對介質(zhì)形狀和位置的重建。

對于任意介質(zhì),以介質(zhì)上某一點為原點可以建立如圖2所示坐標(biāo)系。

圖2 確定介質(zhì)吸收率模型的坐標(biāo)系示意圖

圖2中O為坐標(biāo)原點,坐標(biāo)系xOy為介質(zhì)坐標(biāo)系,x軸和y軸分別與托盤下邊沿和左邊沿平行,坐標(biāo)系xtOyt為投射坐標(biāo)系,yt軸方向與X射線入射方向一致。要確定介質(zhì)的吸收率就是要確定介質(zhì)坐標(biāo)系下,任意位置(x,y)處的吸收率r(x,y)。

現(xiàn)已知180個角度下的接收信息G,即吸收率沿180個不同角度,對于每個角度下512個探測器不同路徑的積分值,要直接根據(jù)該積分值反演出吸收率是困難的。故根據(jù)參考文獻[11],可以通過傅立葉變換和卷積反投影法間接得到吸收率r(x,y)的值。首先,根據(jù)Lambert-bees定理[11]得到:

式中,Gφ(xt)是在X射線照射角度為φ的情況下,在其方向下的接收信息。對r(x,y)做二維傅立葉變換:

式中,F(u,v)為r(x,y)經(jīng)過傅立葉變換后在頻域的結(jié)果,u和v分別為對應(yīng)的頻域坐標(biāo)。再對Gφ(xt)做一維傅立葉變換:

式中,ρ為頻域極坐標(biāo)系下的極軸,根據(jù)中心切片定理[9]:

進而可以通過對其二維傅立葉變換進行逆變換得到r(x,y)為:

式中,Φ為頻域極坐標(biāo)系下的極角,將式(12)改寫為卷積形式:

根據(jù)Paley-Wiener定理[12],ei2πρxt是 不可積的,但只要介質(zhì)的吸收率無大范圍劇烈變化,則可以保證其吸收率的傅立葉變換中高頻分量幅度不大,所以在實際構(gòu)造卷積函數(shù)時,可以通過對卷積核H(ρ)＝｜ρ｜進行限定。根據(jù)參考文獻[9]可采用如下限定:

式中,d為采樣間隔。

3 模型運用

3.1 參數(shù)標(biāo)定結(jié)果

將所給的模板在托盤中的位置等參數(shù)代入模型中,對模型求解得到模板標(biāo)定的結(jié)果如表1所示。

表1 模板標(biāo)定結(jié)果匯總表

從表1中可以看出,探測器的初始旋轉(zhuǎn)角度為29.666 7°,間隔為1.000 0°,則180個角度范圍為29.666 7°～208.666 7°；探測器單元之間的間隔為0.276 7 mm,則512個單元組成的陣列長度為0.1414 m；在如圖2所示的坐標(biāo)系中,系統(tǒng)的旋轉(zhuǎn)中心為(-9.209 4,6.115 6)mm,即旋轉(zhuǎn)中心距正方形托盤的左邊沿40.790 6 mm,距正方形托盤的下邊沿56.115 6 mm。

3.2 介質(zhì)確定結(jié)果

根據(jù)上文的模型,并用MATLAB進行求解,得到兩個未知介質(zhì),即介質(zhì)1和介質(zhì)2的吸收率、形狀和位置。

圖3(a)和圖3(b)分別為依據(jù)某位置介質(zhì)1吸收率數(shù)據(jù),求得的介質(zhì)1的形狀圖像和幾何參數(shù)標(biāo)注圖像,像素為256×256,各個位置吸收率的范圍為0.0～1.4,吸收率為0的部分為無介質(zhì)部分。由圖3可以看出介質(zhì)形狀大致為橢圓,中間有兩個空洞的橢圓,有兩個吸收率較高的橢圓,且二者存在重疊,重疊部分的吸收率最高,以及一個吸收率不同的小圓。圖3(a)中存在部分斜線,這是由于模型固有缺陷導(dǎo)致的偽跡；從圖3(b)中可以看出,橢圓大致位于托盤中心位置,橢圓介質(zhì)距托盤下邊沿約10.546 9 mm,距右邊沿約26.562 6 mm。

根據(jù)所給數(shù)據(jù),對需具體求解吸收率的10個位置進行編號,如圖4所示。

圖3 介質(zhì)1的求解結(jié)果圖

圖4 10個具體位置的編號示意圖

通過對這10個位置的吸收率求解,可以得到結(jié)果如表2所示。

表2 介質(zhì)1在10個具體位置吸收率的結(jié)果

從表2中數(shù)據(jù)可以看出,位置1、3、8、9和10處的吸收率為0,說明這些位置不存在介質(zhì)。

依據(jù)某位置介質(zhì)1的吸收率數(shù)據(jù),求得的介質(zhì)2的形狀圖像和幾何參數(shù)標(biāo)注圖像,像素為256×256,各個位置吸收率的范圍為0.0～8.0,吸收率為0的部分為無介質(zhì)部分,如圖5(a)和圖5(b)所示。從圖5(a)可以看出介質(zhì)形狀很不規(guī)則,為存在很多空洞的網(wǎng)狀圖形,長為94.921 9 mm,寬為86.328 2 mm。從圖5(b)中可以看出,圖形大致位于托盤中心位置,橢圓介質(zhì)距托盤下邊沿5.468 8 mm,距右邊沿4.687 5 mm。

圖5 介質(zhì)2的求解結(jié)果圖

采用與圖4中一致的編號順序,對10個具體位置的吸收率求解,得到結(jié)果如表3所示。

表3 介質(zhì)2在10個具體位置吸收率的結(jié)果

從表3中數(shù)據(jù)可以看出,位置1、4、5和8處的吸收率為0,說明這些位置不存在介質(zhì)。

4 結(jié)束語

本文針對CT系統(tǒng)的參數(shù)標(biāo)定和未知樣品的成像問題,基于模板的接收信息,建立優(yōu)化模型解決了參數(shù)標(biāo)定問題；采用傅立葉變換和卷積反投影法建立了介質(zhì)吸收率重建模型,并根據(jù)接收信息求解出了未知介質(zhì)的吸收率、形狀和位置,求解效果較好。但是本文模型在物體對X射線的吸收率變化劇烈時會產(chǎn)生誤差,需進一步改進。