基于圖像重建的CT系統(tǒng)參數(shù)標(biāo)定與成像模型

萬 犇，賈昀達(dá)，王宏曄，張 勇

(1.電子科技大學(xué) 數(shù)學(xué)科學(xué)學(xué)院，四川 成都 611731；2.電子科技大學(xué) 自動化工程學(xué)院，四川 成都 611731)

典型的二維CT系統(tǒng)是由探測器、待重建物體和X射線組成，平行入射的X射線垂直于探測器平面，每個探測器單元看成一個接收點，且等距排列。

X射線的發(fā)射器和探測器相對位置固定不變，整個發(fā)射-接收系統(tǒng)繞某固定的旋轉(zhuǎn)中心逆時針旋180次。對每一個X射線方向，在具有512個等距單元的探測器上測量經(jīng)位置固定不動的二維待檢測介質(zhì)吸收衰減后的射線能量，并經(jīng)過增益等處理后得到180組接受信息。

為了確定CT系統(tǒng)旋轉(zhuǎn)中心在正方形托盤中的位置、探測器單元之間的距離以及該CT系統(tǒng)使用的X射線的180個方向，在正方形托盤上放置固體介質(zhì)組成的標(biāo)定模板，模板的幾何信息可用一個256維的方陣來表示，其中每一點的數(shù)值反映了該點的吸收強(qiáng)度，本文中稱為“吸收率”。對應(yīng)于該模板的接收信息可用一個180維的向量來表示。利用這一模板及其接收信息，可對CT系統(tǒng)進(jìn)行參數(shù)標(biāo)定。

利用得到的標(biāo)定參數(shù)和某未知介質(zhì)的接收信息，可確定該未知介質(zhì)在正方形托盤中的位置、幾何形狀和吸收率等信息。

1 建模前準(zhǔn)備

根據(jù)CT系統(tǒng)的工作原理[1]，得到如圖1所示的CT系統(tǒng)示意圖。為了標(biāo)定CT系統(tǒng)的參數(shù)，這里采用的標(biāo)定模版由兩個均勻固體介質(zhì)組成，在二維平面上分別是一個橢圓和圓，具體的介質(zhì)形狀和參數(shù)如圖2所示。

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

圖2 模板示意圖(單位：mm)

在建立模型之前，先對已有數(shù)據(jù)進(jìn)行處理和觀察，以簡化相應(yīng)的物理機(jī)理，從而降低模型建立的難度。

1.1 CT系統(tǒng)旋轉(zhuǎn)過程的大致確定

根據(jù)標(biāo)定模版的吸收率方陣[2]，畫出標(biāo)定模版的形狀，如圖3所示，利用標(biāo)定模版得到180組探測器接收信息，三維圖如圖4所示。

圖3 依標(biāo)定模板吸收率方陣畫的標(biāo)定模版示意圖

圖4 180組探測器接收信息三維圖

從圖4得到，探測器接收信息的值有高有低，且經(jīng)過橢圓圓心附近的X射線對應(yīng)的探測器的值在同一組值中最高，說明探測器接收信息的值與X射線穿透物體的長度成正相關(guān)。

利用上述原理，可得到X射線方向經(jīng)由垂直于橢圓長軸到垂直于橢圓短軸，又由于CT系統(tǒng)是逆時針旋轉(zhuǎn)的，所以可以得到如圖5所示的CT系統(tǒng)旋轉(zhuǎn)過程。

圖5 CT系統(tǒng)旋轉(zhuǎn)過程

1.2 三個近似處理

從圖4容易看出，圓所對應(yīng)的探測器值的變化。根據(jù)圓所對應(yīng)的探測器的變化規(guī)律，得到CT系統(tǒng)旋轉(zhuǎn)了180°左右，每次旋轉(zhuǎn)在1°左右。由于偏轉(zhuǎn)角度只有1°左右的誤差以及探測器之間的間距較小，所以所有的探測器值中有最高值的那組可認(rèn)為是X射線的方向垂直于模板中橢圓的短軸，如圖6所示。

圖6 方向垂直橢圓短軸

圖7 方向垂直橢圓長軸

當(dāng)X射線垂直于橢圓長軸時，除去圓對應(yīng)的探測器，由于經(jīng)過橢圓圓心的弦長軸最長，原則上，此時探測器不為0的個數(shù)是最多的，所以橢圓對應(yīng)的探測器不為0的個數(shù)最多的時候大概是X射線方向垂直于橢圓長軸的時候，如圖7所示。

不考慮模板中圓對探測器值的貢獻(xiàn)，對于每一個方向下的探測器接收到的值都有一個最大值，根據(jù)平行線經(jīng)過橢圓原點的弦弦長最長的原理，可近似認(rèn)為每一組探測器接受數(shù)值中最大值對應(yīng)的X射線經(jīng)過橢圓的圓心，如圖8所示。

圖8 經(jīng)過橢圓圓心的X射線

由此建立CT系統(tǒng)參數(shù)標(biāo)定和成像模型前，做了如下3個近似處理[3]。

1)所有的探測器值中有最高值的那組為X射線的方向垂直于模板中橢圓的短軸情形下測得。

2)橢圓對應(yīng)的探測器不為0的個數(shù)最多的時候為X射線方向垂直于橢圓長軸。

3)每一組接收信息中，探測器的值最高的對應(yīng)的X射線經(jīng)過模板中橢圓的圓心。

1.3 消除圓對探測值的影響

根據(jù)圖4，有可能會因為圓和橢圓兩段長度的疊加，導(dǎo)致每一組最大值對應(yīng)的X射線并不穿過橢圓的圓心，為了能夠讓最大值對應(yīng)的X射線穿過橢圓圓心，必須消除圓對探測值的影響，所以可以采取濾波的方式對其進(jìn)行處理。

首先分析圓平均會影響多少個探測器的值，每個探測器影響的值大概是多少，然后對每一組數(shù)據(jù)進(jìn)行梯度檢測，如果相應(yīng)數(shù)據(jù)產(chǎn)生的梯度(或者二階差分值)滿足一定的條件，就減去圓對探測器所產(chǎn)生的影響。下文所有提到探測值都是指濾波之后的探測值。

2 CT系統(tǒng)參數(shù)標(biāo)定模型的建立

探測器單元之間的間距可利用某一段實際長度與對應(yīng)的探測器單元個數(shù)求解得到，CT旋轉(zhuǎn)中心和X射線180次方向可基于建立的直角坐標(biāo)系，經(jīng)機(jī)理分析得到的相應(yīng)方程求解得到。

2.1 確定探測器單元之間的間距

如果知道多個探測器單元的總長度以及對應(yīng)的探測器單元的個數(shù)，就可以確定探測器單元之間的間距。而對于實際長度，只知道模版以及正方形托盤的相應(yīng)參數(shù)。

首先找到X射線垂直于橢圓直徑透射時探測器的值。由之前的近似處理，橢圓對應(yīng)的探測器不為0的個數(shù)最多的時候為X射線方向垂直于橢圓長軸，所以橢圓對應(yīng)的探測器不為0的個數(shù)最多的時候就是X射線垂直于橢圓直徑透射的時候。

然后將橢圓的長軸長度作為實際長度，若要計算探測器單元之間的間距，就要知道此橢圓長軸對應(yīng)的探測器單元個數(shù)N，但是此時只能知道探測器值不為0的探測器單元個數(shù)Nl。Nl的求解只需要知道第一個不為0的探測器編號和最后一個不為0的探測器編號，對于第一個不為0的探測器編號由如下優(yōu)化模型確定：

min i

s.t.Rli≠0

其中，Rli表示第l次第i個探測器對應(yīng)的值。對于最后一個不為0的探測器編號由如下優(yōu)化模型確定：

max j

s.t.Rli≠0

假設(shè)上述問題求解結(jié)果為i0和j0，則Nl=j0-i0+1。

此時Nl個探測器單元之間的長度是小于橢圓長軸長度的，但是Nl個探測器單元對應(yīng)的探測器長度比它們之間的長度之和多出一個探測器單元間距d，而且探測器單元之間的間距較小，所以可直接將橢圓長軸對應(yīng)于Nl個探測器單元，由此可求解出探測器單元之間的間距：

(1)

式中，b為橢圓長半軸的長度。

2.2 確定CT系統(tǒng)的旋轉(zhuǎn)中心在正方形托盤的位置

若要確定旋轉(zhuǎn)中心，一定要有兩個方向，再建立直角坐標(biāo)系，通過中心旋轉(zhuǎn)的相關(guān)知識，設(shè)出未知量和參數(shù)，聯(lián)立方程，最后求解可得到CT系統(tǒng)的旋轉(zhuǎn)中心。

首先不妨選取X射線垂直于橢圓長軸和垂直于橢圓短軸時的兩種特殊情況。之前已經(jīng)得到X射線垂直于橢圓長軸透射時的情形，所以只需要得到X射線垂直于橢圓短軸透射時的情形。由探測器的探測值最大所對應(yīng)的那組數(shù)據(jù)即為X射線垂直于橢圓短軸透射時的接收信息，即：

max Rij

其中，Rij為X射線垂直橢圓長軸時第j個探測器的值。求解得到的i0對應(yīng)的Ri0為X射線垂直于橢圓短軸透射時的接收信息。

此時還可以得到兩組數(shù)據(jù)中各自的最大值Rlm和Rsm，對應(yīng)著的X射線分別與橢圓短軸和長軸重合，分別對應(yīng)的探測器編號為m與n。

然后以橢圓圓心為圓心，短軸為x軸方向，長軸為y軸方向建立直角坐標(biāo)系，如圖9所示。

圖9 求解旋轉(zhuǎn)中心示意圖

X射線垂直于橢圓長軸透射時，此時第m個探測器正對于x軸，縱坐標(biāo)為0，設(shè)橫坐標(biāo)為x0；X射線垂直于橢圓短軸透射時，此時第n個探測器正對于y軸，橫坐標(biāo)為0，設(shè)縱坐標(biāo)為y0。

將兩次旋轉(zhuǎn)中，對應(yīng)的探測器單元連接，得到一系列線段。根據(jù)中心旋轉(zhuǎn)的原理，這些線段的中垂線將相交于一點，該點就是旋轉(zhuǎn)中心。如圖9所示，為方便起見，將X射線垂直于橢圓長軸透射時，從上到下的探測器坐標(biāo)分別設(shè)為(x0,0)、(x0,y1)、(x0,y2)和(x0,y3)。將X射線垂直于橢圓短軸透射時，從左到右的探測器坐標(biāo)分別設(shè)為(x1,y0)、(x2,y0)、(x3,y0)和(0,y0)。

對于第m個探測器，兩點連線的斜率為y0/(x1-x0)，所以其中垂線的斜率為(x0-x1)/y0。中垂線還經(jīng)過兩點的中點((x0+x1)/2,y0/2)，所以可以得到對應(yīng)的中垂線方程：

(2)

同理，可以得到其他3條中垂線方程：

(3)

(4)

(5)

式中，x1=-k1d，x2=-k2d，x3=-k3d，y1=-k3d，y2=-k2d，y3=-k1d,0

將式(2)～式(5)進(jìn)行聯(lián)立，可得一個方程組：

(6)

求解方程組(6)，得到(x，y)、(x0，y0)，其中(x,y)即為CT系統(tǒng)旋轉(zhuǎn)中心在坐標(biāo)系中的坐標(biāo)，即可確定在正方形托盤中的位置。

2.3 確定X射線的180個位置

由于最大值與透射的長度有關(guān)，所以先應(yīng)用數(shù)值擬合得到Rim與對應(yīng)的透射長度lim的關(guān)系，再利用解析幾何知識得到投射長度lim與X射線的旋轉(zhuǎn)角度θi關(guān)系，從而求解出X射線的180個方向。

2.3.1 數(shù)值擬合求探測值和透射長度的關(guān)系

不妨假設(shè)探測值與透射長度之間存在線性關(guān)系，可利用數(shù)值擬合求探測值和透射長度之間的線性關(guān)系，先求出回歸函數(shù)的參數(shù)，再做出圖像，分析殘差值，如果殘差值較小，則探測值與透射長度之間存在線性關(guān)系，并且求得探測值與透射長度的解析表達(dá)式。

對著多組數(shù)據(jù)進(jìn)行擬合，不妨設(shè)擬合函數(shù)R=al+b，擬合函數(shù)中a、b使得殘差平方和最小，即：

利用最小二乘法[4-5]求解上述優(yōu)化問題，即對a、b求偏導(dǎo)，令其等于0，得到：

(7)

求解方程組(7)，可得到每組探測值中的最大值Rim與穿透長度lsi的線性關(guān)系：

Rim=alsi+b

(8)

2.3.2 基于解析幾何求解X射線的180個方向

本文為方便起見，X射線的方向θi定義為X射線方向與x軸正方向之間的夾角。首先通過解析幾何知識推導(dǎo)出X射線的方向θi與透射長度lsi的關(guān)系。如圖6所示，將

(9)

(10)

從而得到X射線方向的求解解析式：

(11)

式中，lsi為該組探測值中最大對應(yīng)的透射長度。

利用處理后的探測值和式(8)、式(11)，可以計算出X射線的180個方向。

2.4 CT系統(tǒng)參數(shù)標(biāo)定模型的求解結(jié)果

以正方形托盤左下角為原點建立直角坐標(biāo)系，通過MATLAB求解得到，CT系統(tǒng)的旋轉(zhuǎn)中心坐標(biāo)為(59.543 1,43.778 7)，位于正方形托盤中心右偏下33.101 0゜，距離11.391 9 mm處；探測器單元間的距離為0.276 8 mm；X射線的180個方向是繞旋轉(zhuǎn)中心從右偏下60.335 4゜,每次逆時針旋轉(zhuǎn)大約1゜到左偏上61.347 9゜。

3 CT系統(tǒng)成像模型的建立

首先基于傅立葉切片定理利用濾波反投射方法得到原來二維物體的圖像，利用灰度表示為矩陣。由于探測器單元的間隔大于正方形托盤上離散點間的距離，所以得到的矩陣維數(shù)是大于所求的矩陣維數(shù)，需對矩陣進(jìn)行壓縮。利用模版得到矩陣中值與吸收率之間的關(guān)系，從而求得二維物體的吸收率。

3.1 基于傅立葉切片定理的濾波反投射方法確定位置與形狀

傅立葉切片定理[6]說明了由一維的圖像信息可以恢復(fù)出二維的圖像信息，具體過程如算法1。

算法1：濾波反投射方法[7]

1)計算每個投影的一維傅立葉變換；

2)用濾波函數(shù)|ω|乘以每個傅立葉變換；

3)得到步驟2導(dǎo)致的每個濾波后的變換的一維傅立葉反變換；

4)對步驟3得到的一維反變換求和，得到二維函數(shù)。

此時，由于旋轉(zhuǎn)中心與正方形托盤中心并不重合，所以還需要對上述得到的圖像與矩陣進(jìn)行平移，然后就可以得到原來二維物體的位置和形狀。

3.2 矩陣的壓縮處理

由于探測器單元之間的間距為0.276 8 mm，而將正方形托盤離散化后，一共256×256個離散點，每個離散點之間的間距為0.39 mm，所以重建圖像中的離散點之間的間隔也為0.276 8 mm，若要使重建圖片與原圖相同，則需要362×362個離散點，所以經(jīng)濾波反投射方法得到的矩陣A為362維的，要轉(zhuǎn)化為256維的，需要對矩陣進(jìn)行降維。矩陣或者圖像的壓縮可直接使用MATLAB中的imresize函數(shù)[8-9]。

3.3 基于數(shù)值擬合確定吸收率

一個二維物體可以用一個吸收率矩陣X表示，該二維物體經(jīng)過CT系統(tǒng)可以得到180組探測數(shù)據(jù)，同樣也可以表示成一個矩陣R。矩陣R經(jīng)過濾波反投射方法可以變?yōu)?62維的矩陣A，經(jīng)過矩陣壓縮處理得到256維的矩陣Y，如圖10所示。

圖10 CT系統(tǒng)成像流程圖

此時可以通過上述方法得到矩陣Y，由矩陣Y推算原來二維物體的吸收率矩陣X，所以現(xiàn)在需要尋找矩陣Y與矩陣X之間的關(guān)系，也就是X中每一個元素對應(yīng)到Y(jié)是怎么變化的，只要求得矩陣X與矩陣Y之間的關(guān)系，即可確定原來物體的吸收率矩陣。

可以利用模板進(jìn)行求解矩陣X與矩陣Y之間的關(guān)系，由于矩陣X中的元素值都為0或1，所以最后得到的矩陣Y應(yīng)該也是兩個常數(shù)，由于一些“噪聲”的干擾，矩陣Y中的元素為在兩個常數(shù)之間波動的一些數(shù)。主要探究的是矩陣X中值為1的元素經(jīng)過上述一系列變換變?yōu)榱硕嗌?，所以需要對矩陣Y進(jìn)行預(yù)處理，將小于一定值的元素舍去，求剩下元素的期望，作為矩陣X中值為1的元素變換后的值。

矩陣X前乘上不同的系數(shù)，觀察變換后的值的變化情況，利用數(shù)值擬合求解出兩者之間的關(guān)系，從而可以通過矩陣Y得到吸收率矩陣X。具體過程如算法2。

算法2：基于數(shù)值擬合的確定吸收率算法

1)確定迭代倍數(shù)k，迭代次數(shù)i=1，終止條件n，以及已知吸收率0-1矩陣X。

2)如果i=n，停止迭代，轉(zhuǎn)步驟4)；否則，Xi=kXi-1，xi=ki，其中xi為Xi中最大元素。將X經(jīng)CT系統(tǒng)，濾波反投射和矩陣壓縮得到對應(yīng)的矩陣Yi。

3)找出矩陣Yi中最大元素max，舍去小于max/2的元素值，對剩下元素求平均得yi。轉(zhuǎn)步驟2)。

4)分別對xi和對應(yīng)值yi進(jìn)行數(shù)值擬合，從而確定矩陣X與矩陣Y之間的關(guān)系。

如果矩陣X與矩陣Y之間是線性關(guān)系，則存在常數(shù)K，使得X=KY，從而可根據(jù)探測值矩陣推算出吸收率矩陣。

3.4 CT系統(tǒng)成像模型的求解結(jié)果

對于濾波反投射方法，可直接利用MATLAB中的iradon[10-12]函數(shù)，由于根據(jù)CT系統(tǒng)參數(shù)標(biāo)定模型，求得CT系統(tǒng)的旋轉(zhuǎn)中心在橢圓圓心(正方形托盤中心)右偏下33.101 0，距離11.391 9mm的地方，所以對iradon函數(shù)得到的圖像進(jìn)行平移，然后再依據(jù)X射線的旋轉(zhuǎn)角度，可以得到對應(yīng)的二維物體的位置和形狀，根據(jù)文獻(xiàn)[1]提供的附件3和附件5數(shù)據(jù)，兩個未知介質(zhì)的形狀如圖11和圖12所示。

圖11 文獻(xiàn)[1]的附件3黑白形狀圖

圖12 文獻(xiàn)[1]的附件5黑白形狀圖

從圖中可以看出，文獻(xiàn)[1]中附件3的二維物體是一個“人腦”CT形狀，位于正方形托盤的中心；文獻(xiàn)[1]中附件5的二維圖片是某個“生物組織”CT形狀，位于正方形托盤的中心。

本文吸收率矩陣X選取標(biāo)定模板的吸收率矩陣，分別用不同的系數(shù)xi乘上矩陣X，再經(jīng)CT系統(tǒng)，濾波反投射和矩陣壓縮得到圖像矩陣Y，然后對矩陣Y中元素進(jìn)行篩選和求期望，得yi，利用線性函數(shù)數(shù)值擬合xi和yi，得如圖13所示圖像。

圖13 吸收率矩陣和圖像矩陣線性擬合圖像

從圖13中可以看出，xi和yi幾乎是呈線性關(guān)系的，也就是說吸收率矩陣X和圖像矩陣Y呈線性關(guān)系為：

X=1.944Y

(12)

從而通過上述操作可以將接收器接收到的180組一維接收信息的數(shù)據(jù)轉(zhuǎn)化成256維的吸收率矩陣。

4 結(jié)束語

本文由于近似處理，減小計算量的同時，也損失了精度；在求解CT系統(tǒng)的旋轉(zhuǎn)中心時，由于初值對方程組的求解影響特別大，所以穩(wěn)定性不高。但是本文所建立的模型對CT系統(tǒng)參數(shù)標(biāo)定和成像都具有一定的指導(dǎo)意義。