夏彬珂，李天昊，鄧明興，，覃思義

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

電子計(jì)算機(jī)斷層掃描(CT)[1]可以在不破壞樣品的情況下，利用樣品對(duì)射線能量的吸收特性，對(duì)生物組織和工程材料的樣品進(jìn)行斷層成像，由此獲取樣品內(nèi)部的結(jié)構(gòu)信息。在典型的二維[2]CT系統(tǒng)中，平行入射的X射線垂直于探測(cè)器平面，每個(gè)探測(cè)器單元可看成一個(gè)接收點(diǎn)，且等距排列。X射線的發(fā)射器和探測(cè)器相對(duì)位置固定不變，整個(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)過(guò)增益等處理后得到180組接收信息。

CT系統(tǒng)安裝時(shí)往往存在誤差，從而影響成像質(zhì)量，因此需要對(duì)安裝好的CT系統(tǒng)進(jìn)行參數(shù)標(biāo)定，即借助于已知結(jié)構(gòu)的樣品(稱為模板)標(biāo)定CT系統(tǒng)的參數(shù)，并據(jù)此對(duì)未知結(jié)構(gòu)的樣品進(jìn)行成像[3]。

CT系統(tǒng)參數(shù)的標(biāo)定是利用探測(cè)單元測(cè)得的信息重建精確圖像的基礎(chǔ)?；跇?biāo)定后的CT系統(tǒng)主要參數(shù)，即可根據(jù)探測(cè)單元接收的信息準(zhǔn)確構(gòu)建未知介質(zhì)的圖像并獲取相關(guān)信息。難點(diǎn)在于根據(jù)標(biāo)定模板的吸收率以及模板的接收信息反推出CT系統(tǒng)的相關(guān)參數(shù)。這就需要建立CT系統(tǒng)的標(biāo)定模型，標(biāo)定CT系統(tǒng)的旋轉(zhuǎn)中心點(diǎn)、探測(cè)器單元之間的距離以及該CT系統(tǒng)使用X射線的180個(gè)方向，進(jìn)而利用這些標(biāo)定的參數(shù)進(jìn)行后面CT圖像的重建，獲取未知介質(zhì)的相關(guān)信息。

本文利用X射線在介質(zhì)中的衰減規(guī)律以及CT圖像重建的相關(guān)知識(shí)，建立了CT參數(shù)標(biāo)定的單目標(biāo)擬合優(yōu)化模型以及CT成像的濾波反投影圖像重建模型。首先結(jié)合X射線在介質(zhì)中的衰減規(guī)律以及Radon變換的相關(guān)理論，基于最小二乘原理建立單目標(biāo)擬合優(yōu)化模型。在求解上采用分組優(yōu)化的方法，利用遺傳算法對(duì)旋轉(zhuǎn)中心點(diǎn)以及探測(cè)單元間距進(jìn)行分組優(yōu)化多次求解，求解結(jié)果基本收斂于同一數(shù)值。之后在標(biāo)定的CT系統(tǒng)參數(shù)下，結(jié)合CT圖像重建相關(guān)知識(shí)，運(yùn)用傅里葉變換以及傅里葉切片定理，建立了CT圖像的濾波反投影圖像重建模型，基于探測(cè)器單元的接收信息，求解圖像重建模型得到了正方形托盤上各個(gè)位置的吸收率信息，并構(gòu)建出了介質(zhì)在托盤上的位置和形狀。最后通過(guò)對(duì)原數(shù)據(jù)信息進(jìn)行圖像重建和對(duì)比，驗(yàn)證了模型的可行性。

本文解決的問(wèn)題來(lái)自2017年高教社標(biāo)全國(guó)大學(xué)生教學(xué)競(jìng)賽本科組A題，文中所提及的“問(wèn)題一”“問(wèn)題二”，以及“附件1”“附件2”等信息請(qǐng)參閱該題目信息。

1 CT系統(tǒng)參數(shù)標(biāo)定模型

在本模型中，已知標(biāo)定用模板的幾何形狀和吸收率數(shù)據(jù)，且若已知CT系統(tǒng)的旋轉(zhuǎn)中心和旋轉(zhuǎn)角度以及探測(cè)單元的位置，則可以通過(guò)Lambert-Beers定理計(jì)算接收到的射線能量。然后將探測(cè)器接收到的信息乘上系統(tǒng)的增益因子，即可得到系統(tǒng)最后反映出的接收信息?，F(xiàn)已知CT系統(tǒng)在180個(gè)未知方向上的接收信息，所以將CT系統(tǒng)旋轉(zhuǎn)中心的坐標(biāo)、180個(gè)方向的角度以及探測(cè)單元的間距作為決策變量，利用題目附件2的數(shù)據(jù)，通過(guò)最小二乘法建立損失函數(shù)。將最小化該損失函數(shù)作為優(yōu)化目標(biāo)，利用遺傳算法求解決策變量的值。

1.1 探測(cè)器實(shí)際接收信息的計(jì)算

對(duì)于具有旋轉(zhuǎn)效果的二維平面束CT系統(tǒng)，其系統(tǒng)結(jié)構(gòu)可抽象為如圖1所示。

圖1 平行束CT掃描示意圖

圖1中s軸為探測(cè)單元所在的一維坐標(biāo)軸，可繞旋轉(zhuǎn)中心旋轉(zhuǎn)。其坐標(biāo)原點(diǎn)為過(guò)旋轉(zhuǎn)中心與作s軸的垂線與s軸的交點(diǎn)，易知該點(diǎn)落在第256與257個(gè)探測(cè)單元中間。由此，那么從左至右第i個(gè)探測(cè)單元在s上的坐標(biāo)為[4]：

(1)

式中，d為探測(cè)單元之間的距離。

當(dāng)CT系統(tǒng)的旋轉(zhuǎn)中心與介質(zhì)物體坐標(biāo)系的坐標(biāo)原點(diǎn)重合時(shí)，若CT系統(tǒng)當(dāng)前的旋轉(zhuǎn)角度為θ，則經(jīng)過(guò)s軸上s點(diǎn)的探測(cè)線方程可記為：

s=xcosθ+ysinθ

(2)

由X射線吸收規(guī)律可知，CT系統(tǒng)旋轉(zhuǎn)角度為θ時(shí)，第i個(gè)探測(cè)器接收到的X射線衰減信息可記為[5，6]：

p(si,θ)=?μ(x,y)δ(xcosθ+ysinθ-si)dxdy(3)

式中，μ(x,y)為介質(zhì)物體在點(diǎn)(x0,y0)處的吸收率。而當(dāng)CT系統(tǒng)的旋轉(zhuǎn)中心點(diǎn)與坐標(biāo)系坐標(biāo)原點(diǎn)不重合時(shí)，記旋轉(zhuǎn)中心點(diǎn)在該坐標(biāo)系中的坐標(biāo)為x0,y0),經(jīng)坐標(biāo)變換后可得：

psi,θ=?μx-x0,y-y0δ[(x-x0)cosθ+
(y-y0)sinθ-si]dxdy

(4)

1.2 最小二乘擬合優(yōu)化模型

題目附件2中的數(shù)據(jù)反映的是經(jīng)過(guò)增益處理后X射線衰減信息。假設(shè)增益處理都是線性的。對(duì)于得到的第i條射線在第j次旋轉(zhuǎn)中穿越介質(zhì)物體的總吸收率pij，pij=p(si,θj)引入增益因子λ，使λpij與題目附件2中的數(shù)據(jù)具有相同的含義，反映的結(jié)果為系統(tǒng)最終輸出的接收信息。在此基礎(chǔ)上基于最小二乘原理構(gòu)建損失函數(shù)：

(5)

式中，tij表示第i個(gè)探測(cè)單元在第j次旋轉(zhuǎn)中測(cè)得的經(jīng)過(guò)增益等處理后的信息，即題目附件2中的第i行第j列數(shù)據(jù)。

以使該損失函數(shù)最小為目標(biāo)，建立如下優(yōu)化模型：

(6)

式中，(x0,y0)為CT系統(tǒng)旋轉(zhuǎn)中心點(diǎn)的坐標(biāo)，θj為CT系統(tǒng)第j次旋轉(zhuǎn)的旋轉(zhuǎn)角。優(yōu)化模型的決策變量為x0,y0,λ,d,θj(j=1,2,,180)。式中pij=psi,θj，其計(jì)算公式為：

pij=?μx-x0,y-y0δ[(x-x0)cosθ+
(y-y0)sinθ-si]dxdy

(7)

1.3 標(biāo)定模型求解方法

針對(duì)問(wèn)題的實(shí)際情況對(duì)以上擬合優(yōu)化的模型進(jìn)行了適應(yīng)性的求解預(yù)處理，包括：

1)對(duì)原始圖片通過(guò)雙線性插值[7]將像素從256×256擴(kuò)增為1024×1024，增加求解的精度；

3)對(duì)接收信息數(shù)據(jù)進(jìn)行等距分組，以30個(gè)旋轉(zhuǎn)角為一組的組內(nèi)間距將所有數(shù)據(jù)分為30組，每組含有6個(gè)旋轉(zhuǎn)方向。通過(guò)減少每次優(yōu)化求解的決策變量個(gè)數(shù)，以減少求解運(yùn)算的復(fù)雜度；

4)通過(guò)對(duì)接收信息數(shù)據(jù)的分析，減小旋轉(zhuǎn)角的約束范圍，使求解時(shí)結(jié)果能更快的收斂到最優(yōu)解。

最后，我們針對(duì)問(wèn)題一進(jìn)行適應(yīng)性修正后的擬合優(yōu)化模型為：

(8)

(9)

考慮到模型中有10個(gè)變量進(jìn)行優(yōu)化時(shí)變量仍較多，故采用對(duì)多變量?jī)?yōu)化問(wèn)題表現(xiàn)良好的遺傳算法[8]進(jìn)行求解。遺傳算法將變量進(jìn)行編碼轉(zhuǎn)化為基因序列，通過(guò)模擬自然界繁衍、淘汰的規(guī)律對(duì)目標(biāo)函數(shù)進(jìn)行優(yōu)化，最終得到一個(gè)較優(yōu)的種群和最優(yōu)解。下面以第一組數(shù)據(jù)為例，對(duì)遺傳算法的求解過(guò)程進(jìn)行解釋。

求解時(shí)采用實(shí)數(shù)編碼，種群個(gè)數(shù)為200，迭代次數(shù)為500，利用MATLAB[9]遺傳算法工具包進(jìn)行求解。

1.4 求解結(jié)果

表1 標(biāo)定模型求解結(jié)果

x0=Δ0(-64)=-9.2773mm
y0=Δ0×95=6.2500mm
d=Δ0×2.8314=0.2765mm

式中，Δ0為一個(gè)小單元格的寬度。故CT旋轉(zhuǎn)中心在正方形托盤中的位置坐標(biāo)為(-9.2773,6.2500)mm，即距托盤中心點(diǎn)左偏9.2773mm，上偏6.2500mm，探測(cè)器單元之間的距離為0.2765mm。

CT系統(tǒng)的起始旋轉(zhuǎn)角度為29.7893°，基本以1°為間隔依次逆時(shí)針旋轉(zhuǎn)180°。

2 二維CT圖像重建模型

本模型基于傅里葉中心切片定理，用濾波反投影方法建立二維CT圖像重建模型。可以根據(jù)CT系統(tǒng)的X射線在180個(gè)方向上的衰減信息，重構(gòu)出測(cè)量物質(zhì)的吸收率和幾何形狀信息。

2.1 濾波反投影重建方法

在CT斷層成像技術(shù)中，最直接的方法是通過(guò)各個(gè)方向的掃描信息，直接通過(guò)累加法反投影計(jì)算出原圖像各個(gè)像素點(diǎn)的密度(衰減系數(shù))數(shù)值。但是由于是在旋轉(zhuǎn)的空間進(jìn)行重建，所以直接累加會(huì)引入星狀偽跡，即原圖像中密度為零的點(diǎn)，重建后密度不一定為零，圖像出現(xiàn)失真。為了克服這種現(xiàn)象，本文在投影以前先對(duì)投影數(shù)據(jù)進(jìn)行插值和濾波等修正，再把修正后的數(shù)據(jù)進(jìn)行反投影得到無(wú)偽跡的圖像[10-11]。

在利用濾波反投影對(duì)CT掃描投影數(shù)據(jù)進(jìn)行圖像重構(gòu)時(shí)，需要首先選擇濾波器的響應(yīng)函數(shù)。通過(guò)上面重建原理的分析，可以得到理論上的濾波頻率響應(yīng)函數(shù)為ω，但是實(shí)際系統(tǒng)存在一定帶寬，所以在設(shè)計(jì)濾波器時(shí)一般都有上限頻率。傳統(tǒng)的方法一般采用Ram-Lak(R-L)濾波器[12]，直接對(duì)前面的響應(yīng)函數(shù)設(shè)計(jì)一個(gè)帶限，如圖2(a)所示。但這種濾波器會(huì)在高頻產(chǎn)生吉布斯現(xiàn)象，使重構(gòu)的圖像邊緣產(chǎn)生振蕩。因此，這里采用一種更加平滑的方法，Shepp-Logan(S-L)濾波器[13]，響應(yīng)函數(shù)的圖像如圖2(b)所示，可以看作為ω加上了一個(gè)光滑約束。

圖2 兩種濾波器的頻率響應(yīng)

S-L濾波器的窗函數(shù)為：

(10)

因此，其頻率響應(yīng)函數(shù)為：

對(duì)應(yīng)的卷積函數(shù)為：

(12)

將180個(gè)方向上的投影信息(即X射線的吸收強(qiáng)度)通過(guò)濾波器，進(jìn)行濾波處理，將結(jié)果運(yùn)用反投影方法直接累加得到原圖像的吸收率分布，其反投影重建過(guò)程如圖3所示。

2.2 重建模型求解方法

為了得到探測(cè)器接收到的實(shí)際衰減信息，需要在濾波反投影重建前對(duì)接收信息的增益進(jìn)行還原。得到衰減信息的數(shù)據(jù)矩陣：

圖3 濾波反投影重建過(guò)程

(13)

式中，n=512,p=180，λ是CT系統(tǒng)對(duì)接收到的衰減信息進(jìn)行增益的比例，數(shù)據(jù)矩陣P中的元素pij表示在第j個(gè)旋轉(zhuǎn)角度θj下，探測(cè)器的第i個(gè)單元接收到的衰減信息，pij=p(si,θj)。這個(gè)衰減信息等效于該單元的探測(cè)直線上吸收率的線積分。

在得到了探測(cè)器的實(shí)際接收衰減信息后，利用濾波反投影重建方法對(duì)原始圖像的吸收率分布進(jìn)行重建。由于旋轉(zhuǎn)中心不一定在重建物體的幾何中心，所以在重建前需要在數(shù)據(jù)矩陣P每一列的首尾增加上100個(gè)單位的零填充值，以保證整個(gè)圖像能被完整還原。

(14)

最后，由于濾波反投影重建是一種數(shù)值方法，在原圖像上吸收率應(yīng)該為0的地方也會(huì)填充有微小數(shù)值，所以本文在處理最后結(jié)果時(shí)，將值小于0.1的點(diǎn)均置0，得到符合真實(shí)情況的吸收率分布。

2.3 求解結(jié)果

最后得到問(wèn)題二的吸收率分布情況如圖4所示。

同樣對(duì)問(wèn)題三進(jìn)行求解，得到的吸收率分布情況如圖5所示。

圖4 問(wèn)題二介質(zhì)的位置和幾何形狀

圖5 問(wèn)題三介質(zhì)的位置和幾何形狀

2.4 CT圖像重建模型的驗(yàn)證

最后，對(duì)本文中圖像重建的模型和算法進(jìn)行可行性驗(yàn)證。利用題目附件二的接收信息和已經(jīng)標(biāo)定的參數(shù)，對(duì)模板的位置、幾何形狀和吸收率進(jìn)行重建，并與實(shí)際結(jié)果進(jìn)行對(duì)比，說(shuō)明了該重建模型和算法的可行性。

圖6是利用該CT圖像重建算法重構(gòu)出的模板的位置和幾何形狀與原始圖像的對(duì)比。

圖6 CT圖像重建模型的驗(yàn)證

重建結(jié)果的位置和幾何形狀與原圖像的符合程度均非常好，只是重構(gòu)圖像在邊緣上比原始圖像略小一圈，這是因?yàn)閳D像邊緣的頻率太高，在通過(guò)濾波器時(shí)會(huì)被濾掉，但是這對(duì)分析模板的吸收率和確定位置和形狀并沒(méi)太大影響。此外，通過(guò)對(duì)比重建的模板各點(diǎn)的吸收率差異，可以看到重建吸收率與原始吸收率相比誤差極小，驗(yàn)證了該模型的準(zhǔn)確性。

3 結(jié)束語(yǔ)

對(duì)于CT系統(tǒng)參數(shù)的標(biāo)定，本文根據(jù)X射線衰減規(guī)律得到接受信息的理論計(jì)算式，基于最小二乘原理建立了理論與實(shí)際接受信息的損失函數(shù)，并以未知參數(shù)(中心點(diǎn)、探測(cè)單元間距、180個(gè)旋轉(zhuǎn)角)為決策變量，抓住題目信息隱含的約束條件，建立了單目標(biāo)優(yōu)化模型。

由于決策變量過(guò)多，采用單目標(biāo)優(yōu)化模型直接求解十分困難，甚至得不到正確解。因此按一定規(guī)則多次從180組數(shù)據(jù)中每次抽取6組數(shù)據(jù)進(jìn)行分組優(yōu)化，較好的利用了Radon變換正弦圖得到了旋轉(zhuǎn)角的大致初值范圍，減小了搜索范圍，增大了精度。抽取30次先行計(jì)算中心點(diǎn)和單元距離等數(shù)據(jù)，得到了精確統(tǒng)一的收斂解，然后利用該數(shù)據(jù)對(duì)旋轉(zhuǎn)角變量進(jìn)行逐組求解，求得了精確的180個(gè)CT系統(tǒng)旋轉(zhuǎn)角度值，較好地完成了系統(tǒng)的標(biāo)定。

對(duì)于問(wèn)題二和問(wèn)題三CT圖像的重建，采取了醫(yī)學(xué)上常用濾波反投影圖形重建方法，去除了直接反投影重建算法中會(huì)出現(xiàn)的星狀偽影，得到了精確的重建圖像以及未知介質(zhì)的具體信息，并利用問(wèn)題一中的已知模板驗(yàn)證了CT圖像重建模型的可行性。