耿瑞旭，吳迪遠(yuǎn)，楊廣緒，程光輝

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

CT系統(tǒng)可利用樣品對射線能量的吸收特性對生物組織和工程材料進(jìn)行斷層成像[1]。本文對于一個(gè)CT系統(tǒng),借助于已知結(jié)構(gòu)的樣品標(biāo)定CT系統(tǒng)的參數(shù),并據(jù)此對未知結(jié)構(gòu)樣品進(jìn)行成像。

本文在建立了非線性規(guī)劃模型,Radon反變換復(fù)原模型和切線區(qū)域驗(yàn)證模型后對CT系統(tǒng)進(jìn)行標(biāo)定和成像；并利用具體的模板信息和數(shù)據(jù)吸收信息對CT系統(tǒng)的參數(shù)進(jìn)行了具體的標(biāo)定,確定了CT系統(tǒng)旋轉(zhuǎn)中心在正方形托盤中的位置、探測器單元之間的距離以及該CT系統(tǒng)使用的X射線的180個(gè)方向。在標(biāo)定完成后,利用上述系統(tǒng)對未知介質(zhì)的接收信息,確定該未知介質(zhì)在托盤中的位置、幾何形狀和吸收率等信息。

1 分析與模型假設(shè)

系統(tǒng)三個(gè)參數(shù)的求解可以分為兩個(gè)部分進(jìn)行:首先求解探測器單元之間的距離及系統(tǒng)使用的X射線的180個(gè)方向；其次利用得到的方向變化規(guī)律和距離數(shù)值求解旋轉(zhuǎn)中心。

為了利用接收器信息對輪廓規(guī)則的外形介質(zhì)進(jìn)行復(fù)原,可將信息復(fù)原分為三部分。

1)未知介質(zhì)的位置,每次旋轉(zhuǎn)時(shí)的兩條切線都會(huì)得到一個(gè)可行域,對180次可行域求交集,可以得到未知介質(zhì)的位置。

2)未知介質(zhì)的形狀,對接收信息進(jìn)行反Radon變換,即可得到未知介質(zhì)的形狀,并可驗(yàn)證得到的介質(zhì)形狀未發(fā)生畸變,僅發(fā)生了縮放。之后,根據(jù)比例關(guān)系,得到在正方形板上實(shí)際的介質(zhì)形狀,并將其所在位置與第一步中得到的位置結(jié)果進(jìn)行對比,從而驗(yàn)證形狀和位置的正確性。

3)根據(jù)Radon逆變換求出介質(zhì)每一點(diǎn)的衰減系數(shù),再根據(jù)標(biāo)定模板,求出衰減系數(shù)與吸收度的線性關(guān)系,從而進(jìn)行求解,即可得到未知介質(zhì)的位置、形狀和吸收度。

本文假定:

1)重建CT圖像時(shí)旋轉(zhuǎn)中心始終不發(fā)生偏移；

2)根據(jù)CT成像原理[2-3],Radon反變換得到的衰減系數(shù)與吸收強(qiáng)度成比例關(guān)系；

3)標(biāo)定模型的兩部分分別為橢圓和圓。

2 參數(shù)標(biāo)定模型的分析與建立

2.1 確定單元距離與射線方向的非線性優(yōu)化模型

參數(shù)標(biāo)定過程示意圖如圖1所示。

在此坐標(biāo)系中,添加了第i次旋轉(zhuǎn)時(shí)平行入射的X射線,以及橢圓吸收介質(zhì)E,圓形吸收介質(zhì)C,記射線l為X射線中與坐標(biāo)原點(diǎn)距離向量最大的射線(以向上為正方向)；β(i)為X射線方向與坐標(biāo)系x軸正方向的夾角；s(i)表示第i次旋轉(zhuǎn)時(shí)坐標(biāo)原點(diǎn)到512個(gè)接收單元對應(yīng)的射線的最大距離；d表示接收器相鄰接收單元的距離。對于圖1中第i次旋轉(zhuǎn)時(shí)的示意圖,設(shè)其中X射線的方程為:

圖1 標(biāo)定模型探測器與射線簇示意圖

首先求解X射線在橢圓中的弦長。橢圓E的方程為:

由于橢圓長軸a和短軸b的取值分別為15和40,故可得到第i次旋轉(zhuǎn)時(shí)從左到右的第j條X射線穿過的橢圓形材料長度為:

式中,

之后,求解X射線在圓中的弦長。圓C的方程為:

將X射線的方程代入式(4),可得:

式中,

最后,構(gòu)造目標(biāo)函數(shù)。CT系統(tǒng)標(biāo)定時(shí)的介質(zhì)對X射線的吸收是均勻的,故探測器每一個(gè)接收單元接收到的信息應(yīng)與該單元接收到的信號所穿過的介質(zhì)長度成正比,即:

式中:w0(i,j)為實(shí)際中給出的第i次旋轉(zhuǎn)時(shí)第j個(gè)探測器的接收信息,是已知的；K為比例系數(shù),L(i,j)為計(jì)算得到的第i次旋轉(zhuǎn)時(shí)第j個(gè)探測器穿過的介質(zhì)長度,即L(i,j)＝L1(i,j)＋L2(i,j)。為了消除接收到的信息與穿過的介質(zhì)長度之間的比例關(guān)系K,本文使用了比值法,即:

式中,等號左邊為實(shí)際值,等號右邊為由D(i,j)和β(i)兩個(gè)決策變量構(gòu)成的計(jì)算值。根據(jù)式(7),利用最小二乘原理可構(gòu)造出非線性規(guī)劃模型的優(yōu)化函數(shù)。即:

式中:i為旋轉(zhuǎn)的次數(shù),對每一個(gè)i,均可求出一組最優(yōu)的β(i)、s(i)和d；β(i)為第i次旋轉(zhuǎn)時(shí)射線與x軸正方向的夾角；s(i)為第i次旋轉(zhuǎn)時(shí)坐標(biāo)原點(diǎn)到512個(gè)接收單元的最大距離；d為探測器上相鄰兩個(gè)等距單元的距離,應(yīng)為與i無關(guān)的常量。

綜上所述,本文給出非線性優(yōu)化目標(biāo)的優(yōu)化函數(shù)如式(8)所示。

式中,

決策變量為β(i)、s(i)和d。

考慮到i有180個(gè)取值,故最后應(yīng)得到180組優(yōu)化結(jié)果,每組結(jié)果包含3個(gè)決策變量的最優(yōu)值。

2.2 確定旋轉(zhuǎn)中心位置的非線性規(guī)劃模型

在圖1所示的坐標(biāo)系(坐標(biāo)系1)中,設(shè)旋轉(zhuǎn)中心點(diǎn)坐標(biāo)為P(xo,yo),則在題目和附件給出的默認(rèn)坐標(biāo)系中,其坐標(biāo)為P1(xo＋50,yo＋50)。計(jì)算時(shí)統(tǒng)一在坐標(biāo)系1中進(jìn)行計(jì)算,即以P(xo,yo)為坐標(biāo)。

根據(jù)2.1節(jié)的求解結(jié)果和點(diǎn)到直線的距離公式,可得到第i次旋轉(zhuǎn)時(shí),旋轉(zhuǎn)中心點(diǎn)P(xo,yo)到第一條射線的距離為:

式中:k(i)為第i次旋轉(zhuǎn)時(shí)X射線的斜率,即k(i)＝tan[β(i)]；c1為第一條射線的截距； (xo,yo)為旋轉(zhuǎn)中心點(diǎn)P的坐標(biāo)。

由此,即可求得每次旋轉(zhuǎn)時(shí)旋轉(zhuǎn)中心點(diǎn)P(xo,yo)到第一條射線的距離r(i)。而優(yōu)化目標(biāo)便是使這180個(gè)距離值盡可能保持恒定,即方差盡可能小。因此,根據(jù)方差的定義,可以得到非線性規(guī)劃模型的目標(biāo)函數(shù)為:

對于約束條件,主要是對(xo,yo)的坐標(biāo)約束,即旋轉(zhuǎn)中心應(yīng)在正方形托盤中,故有如下約束:

綜上所述,該非線性規(guī)劃模型為:

其中,旋轉(zhuǎn)中心點(diǎn)坐標(biāo)(xo,yo)為決策變量。優(yōu)化目標(biāo)中,r(i)的計(jì)算公式為:

3 成像模型的分析與建立

3.1 基于切線交集的介質(zhì)輪廓與未知復(fù)原模型

首先,設(shè)第i次旋轉(zhuǎn)的X射線簇的斜率為ki,與待復(fù)原物質(zhì)相切的兩條切線的截距從大到小分別為ci,1和ci,2,則對于第i次旋轉(zhuǎn),其可行域的約束為:

因此,對于全部180次旋轉(zhuǎn),可得如下可行域的交集表示形式,這種通過求切線所圍區(qū)域的交集求出幾何圖形位置與輪廓的方法也可看作是一種線性規(guī)劃,即:

做出上述180組不等式的幾何圖形,即為180組切線圍成的圖形,也就得到了未知介質(zhì)圖形的輪廓與位置。

3.2 基于Radon逆變換的未知介質(zhì)幾何形狀復(fù)原模型

Radon及其逆變換是對物體的平行束投影及其復(fù)原的數(shù)學(xué)描述。一個(gè)N維函數(shù)的Radon變換是(N-1)維超平面上的積分值,相應(yīng)地,Radon逆變換就是指根據(jù)這些積分值反求得到原始被積函數(shù)[4]。在本題中,僅從二維對物體進(jìn)行投影,故N＝2,即函數(shù)f(x,y)為二維空間上的函數(shù)[5],其沿著任意投影射線防線的線積分就為其Radon變換。Radon變換表示為[6-7]:

Radon逆變換的公式為[8]:

用lQ,θ表示投影射線[3],其方程為:

式中,θ代表投影射線的法線和x軸正向之間的夾角,Q代表投影射線與坐標(biāo)原點(diǎn)之間的距離。

由此,通過對接受信息進(jìn)行Radon逆變換復(fù)原原始圖像是可行的。在MATLAB中,Radon逆變換的函數(shù)為iradon。

4 模型的求解

4.1 單元距離與射線方向的求解

求解該非線性約束模型的步驟如下:

1)求出由初始輸入值構(gòu)成的初始輸入x0；

2)將x0作為初始迭代點(diǎn),由fminseach函數(shù)求出極值x；

3)令x0＝x,重復(fù)步驟2)；

4)迭代180次后求解完成。

其中,fminsearch是MATLAB中所帶的函數(shù)[9-10],但使用之前需首先求出一個(gè)合適的初始迭代點(diǎn)[11]。初值估計(jì)主要對第一次旋轉(zhuǎn)時(shí)的角度、探測器相鄰單元之間的距離和第一次旋轉(zhuǎn)時(shí)第一條射線距坐標(biāo)原點(diǎn)的距離。

首先,利用標(biāo)定給出的第一次旋轉(zhuǎn)時(shí)的數(shù)據(jù)。對該列數(shù)據(jù)進(jìn)行繪圖,如圖2(a)所示；對數(shù)據(jù)表格進(jìn)行二值繪圖,如圖2(b)所示。

圖2 初始時(shí)刻接收的信息分布圖和接收信息的可視化

通過對數(shù)據(jù)進(jìn)行分析,圓形介質(zhì)的高度為28,即圓形介質(zhì)約8 mm的直徑占用了接收器大約28個(gè)單元,進(jìn)而得到d的取值。最終,考慮到旋轉(zhuǎn)為逆時(shí)針方向及擬合后的效果,將β(1)、d和s(1)作為初始值,代入MATLAB中fminsearch函數(shù)進(jìn)行求解,可得:

4.2 旋轉(zhuǎn)中心點(diǎn)的求解

對有約束非線性規(guī)劃模型,使用MATLAB優(yōu)化工具箱中的fmincon函數(shù)進(jìn)行求解[1]。該模型的求解復(fù)雜度較低,故不需要對初始迭代值進(jìn)行較準(zhǔn)確的估計(jì)。根據(jù)系統(tǒng)的旋轉(zhuǎn)中心常常與其幾何中心接近的生活經(jīng)驗(yàn),在實(shí)際求解時(shí),設(shè)置坐標(biāo)系1中的原點(diǎn)(0,0)作為初始迭代點(diǎn)。模型的求解步驟如下:

1)將初始值(0,0)作為初始迭代點(diǎn)；

2)由MATLAB優(yōu)化工具箱中的fmincon函數(shù)求出極值(x,y)；

3)將極值(x,y)作為新的迭代點(diǎn),重復(fù)步驟2)；

4)求解成功后退出程序。

5 計(jì)算結(jié)果

5.1 參數(shù)標(biāo)定計(jì)算結(jié)果

對探測器單元間距離和射線方向的求解結(jié)果如表1所示。

表1 CT系統(tǒng)參數(shù)標(biāo)定結(jié)果

按照4.2節(jié)中所述的方法進(jìn)行求解,可得到旋轉(zhuǎn)中心P在圖1所示坐標(biāo)系中的位置坐標(biāo)是P(-9.272 8,6.288 6)。如圖3所示,給出了180個(gè)方向中的四個(gè)方向,即第一次旋轉(zhuǎn)時(shí)的方向(與x軸正方向夾角),第60次旋轉(zhuǎn)的方向(與x軸正方向夾角),第120次旋轉(zhuǎn)的方向(與x軸正方向夾角),第180次旋轉(zhuǎn)的方向(與x軸正方向夾角)。

圖3 四次旋轉(zhuǎn)過程示意圖

5.2 未知介質(zhì)復(fù)原結(jié)果

對于已有的反投影數(shù)據(jù),可以通過切線可行域交集模型求得待測介質(zhì)的位置,并通過反Radon變換得到未知介質(zhì)的真實(shí)形狀。

圖4 圖像復(fù)原結(jié)果

如圖4所示,圖(a)為Radon逆變換得到的直接結(jié)果；圖(b)為進(jìn)行按比例縮放和位置驗(yàn)證后的最終結(jié)果。圖中心部分的大橢圓即為物體外輪廓,其中下部有兩個(gè)孔洞,上部有兩個(gè)橢圓形物體。由圖4可知,通過Radon逆變換得到的復(fù)原結(jié)果的位置與通過取切線交集區(qū)域得到的位置幾乎完全吻合。

6 結(jié)束語

本文針對CT系統(tǒng)參數(shù)標(biāo)定的問題建立了無約束非線性規(guī)劃模型,并用直接搜索法對模型進(jìn)行了求解,從而得到了較為精確的CT系統(tǒng)參數(shù)估計(jì)模型。針對CT系統(tǒng)的圖像重建問題,分別利用Radon反變換和切線可行域取交集法得到圖像的形狀和外部輪廓。