基于卷積反投影的CT系統(tǒng)參數(shù)標(biāo)定和成像模型

王譽(yù)儒，楊玉聰，官 莊，王志勇

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

電子計(jì)算機(jī)斷層掃描(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個(gè)等距單元的探測(cè)器,探測(cè)器單元視作等距排列的接收點(diǎn)；X射線的發(fā)射器和探測(cè)器相對(duì)位置固定不變,整個(gè)發(fā)射-接收系統(tǒng)繞某固定的旋轉(zhuǎn)中心逆時(shí)針旋轉(zhuǎn)180次,以測(cè)量射線能量,并經(jīng)過增益等處理后得到180組接收信息。

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

求解得到:

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

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

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

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

求解得到:

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

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

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

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

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

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

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

圖2中O為坐標(biāo)原點(diǎn),坐標(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è)角度下的接收信息G,即吸收率沿180個(gè)不同角度,對(duì)于每個(gè)角度下512個(gè)探測(cè)器不同路徑的積分值,要直接根據(jù)該積分值反演出吸收率是困難的。故根據(jù)參考文獻(xiàn)[11],可以通過傅立葉變換和卷積反投影法間接得到吸收率r(x,y)的值。首先,根據(jù)Lambert-bees定理[11]得到:

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

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

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

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

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

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

式中,d為采樣間隔。

3 模型運(yùn)用

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

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

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

從表1中可以看出,探測(cè)器的初始旋轉(zhuǎn)角度為29.666 7°,間隔為1.000 0°,則180個(gè)角度范圍為29.666 7°～208.666 7°；探測(cè)器單元之間的間隔為0.276 7 mm,則512個(gè)單元組成的陣列長(zhǎng)度為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進(jìn)行求解,得到兩個(gè)未知介質(zhì),即介質(zhì)1和介質(zhì)2的吸收率、形狀和位置。

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

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

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

圖4 10個(gè)具體位置的編號(hào)示意圖

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

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

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

依據(jù)某位置介質(zhì)1的吸收率數(shù)據(jù),求得的介質(zhì)2的形狀圖像和幾何參數(shù)標(biāo)注圖像,像素為256×256,各個(gè)位置吸收率的范圍為0.0～8.0,吸收率為0的部分為無介質(zhì)部分,如圖5(a)和圖5(b)所示。從圖5(a)可以看出介質(zhì)形狀很不規(guī)則,為存在很多空洞的網(wǎng)狀圖形,長(zhǎ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中一致的編號(hào)順序,對(duì)10個(gè)具體位置的吸收率求解,得到結(jié)果如表3所示。

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

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

4 結(jié)束語

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