張永紅

(渭南師范學(xué)院數(shù)理學(xué)院，陜西渭南714099)

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

張永紅

(渭南師范學(xué)院數(shù)理學(xué)院，陜西渭南714099)

計(jì)算機(jī)斷層成像術(shù)(Computed Tomography，簡(jiǎn)稱CT)是用來獲取觀測(cè)目標(biāo)斷層圖像的一門技術(shù)。然而，在CT系統(tǒng)安裝時(shí)往往存在誤差，從而影響成像質(zhì)量，因此需要對(duì)安裝好的CT系統(tǒng)進(jìn)行參數(shù)標(biāo)定。在探測(cè)器單元等距的前提下，首先計(jì)算出探測(cè)器單元之間的距離，然后以正方形托盤中心為坐標(biāo)原點(diǎn)建立直角坐標(biāo)系，分別找到幾種特殊的旋轉(zhuǎn)方向，進(jìn)而找出旋轉(zhuǎn)中心坐標(biāo)(－9.27，6.27)，最后計(jì)算出每次的旋轉(zhuǎn)方向。數(shù)值實(shí)驗(yàn)表明，該模型簡(jiǎn)單易行，計(jì)算結(jié)果比較準(zhǔn)確。

計(jì)算機(jī)斷層成像術(shù);幾何校正;旋轉(zhuǎn)中心;X射線

1 問題重述

計(jì)算機(jī)斷層成像術(shù)CT(Computed Tomography)可以在不破壞樣品的情況下，利用樣品對(duì)射線能量的吸收特性對(duì)生物組織和工程材料的樣品進(jìn)行斷層成像，由此獲取樣品內(nèi)部的結(jié)構(gòu)信息。一種典型的二維CT系統(tǒng)如圖1所示(全國(guó)大學(xué)生數(shù)學(xué)建模競(jìng)賽2017年A題)，平行入射的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)過增益等處理后得到180組接收信息。

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

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

圖3 建立直角坐標(biāo)系(單位:mm)

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)行成像。請(qǐng)建立相應(yīng)的數(shù)學(xué)模型和算法，解決以下問題:

在正方形托盤上放置兩個(gè)均勻固體介質(zhì)組成的標(biāo)定模板，模板的幾何信息如圖2所示，相應(yīng)的數(shù)據(jù)文件見文獻(xiàn)［1］中的附件1，其中每一點(diǎn)的數(shù)值反映了該點(diǎn)的吸收強(qiáng)度，這里稱為“吸收率”。對(duì)應(yīng)于該模板的接收信息見文獻(xiàn)［1］中的附件2。請(qǐng)根據(jù)這一模板及其接收信息，確定CT系統(tǒng)旋轉(zhuǎn)中心在正方形托盤中的位置、探測(cè)器單元之間的距離以及該CT系統(tǒng)使用的X射線的180個(gè)方向。［1］

2 問題分析

2.1 確定探測(cè)器單元之間的距離

針對(duì)該問題，要求計(jì)算探測(cè)器單元之間的距離，該問題可通過查詢CT系統(tǒng)的工作原理，求出探測(cè)器吸收強(qiáng)度與探測(cè)器單元距離的對(duì)應(yīng)關(guān)系，從而求出探測(cè)器單元之間的距離，或者根據(jù)CT系統(tǒng)中，發(fā)射—接收系統(tǒng)中兩束射線可穿過固體介質(zhì)最大直徑對(duì)應(yīng)的單元格個(gè)數(shù)求解出探測(cè)器單元之間的距離。

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

針對(duì)該問題，需要確定CT系統(tǒng)旋轉(zhuǎn)中心在正方形托盤上的位置，第一步建立合適的平面直角坐標(biāo)系，由于在任意方向時(shí)，探測(cè)器所接收的橢圓和圓的最大吸收強(qiáng)度分別為通過橢圓中心和圓心的射線(如果沒有通過中心，可以選擇鄰近的射線)，利用Matlab軟件找出最大值及其所在的行和列，從而可以確定出幾種特殊的旋轉(zhuǎn)方向，進(jìn)而根據(jù)平面幾何相關(guān)知識(shí)求解出CT系統(tǒng)旋轉(zhuǎn)中心在正方形托盤上的位置。

2.3 確定CT系統(tǒng)使用的 X射線的180個(gè)方向

針對(duì)該問題，同2.2中所述，找出某些特殊點(diǎn)的坐標(biāo)及其旋轉(zhuǎn)角度，首先需要確定 X射線的初始入射方向，根據(jù)文獻(xiàn)［1］附件2中第一列數(shù)據(jù)，可以獲得X射線的初始入射時(shí)探測(cè)器的吸收強(qiáng)度，進(jìn)而查找出一些其他特殊的旋轉(zhuǎn)方向，從而最終確定出180個(gè)旋轉(zhuǎn)方向。

3 模型假設(shè)

(1)假設(shè)所有的平行光束都照到探測(cè)器上，探測(cè)器單元等距排列。［2－4］

(2)旋轉(zhuǎn)載物臺(tái)為完全水平放置。

(3)載物臺(tái)中軸線所確定的平面在平板探測(cè)器的某個(gè)垂面上。

(4)平板探測(cè)器表面完全平整光滑。

(5)附件中所給的數(shù)據(jù)真實(shí)可靠。

4 模型的建立和求解

4.1 模型準(zhǔn)備

建立坐標(biāo)系:以正方形托盤中心為坐標(biāo)原點(diǎn)O，橢圓的短軸為ox軸，橢圓的長(zhǎng)軸為oy軸建立平面直角坐標(biāo)系，如圖3所示。

4.2 探測(cè)器單元之間的距離模型

在圖3中，假設(shè)當(dāng)X射線與oy軸平行時(shí)，通過橢圓左右兩側(cè)的X射線對(duì)應(yīng)的探測(cè)器［5－6］分別為Ti、Tj，通過圓左右兩側(cè)的分別為Tk、Tl，則橢圓所對(duì)應(yīng)探測(cè)器個(gè)數(shù)為n，圓所對(duì)應(yīng)的探測(cè)器個(gè)數(shù)為m，圓的半徑為r，橢圓的長(zhǎng)短軸分別為a、b，探測(cè)器單元之間的距離為d，因?yàn)槠叫腥肷涞腦射線垂直于探測(cè)器平面，并且每個(gè)探測(cè)器單元看成一個(gè)接收點(diǎn)，且等距排列，但是由于X射線可能與橢圓或圓不是正好相切的，因此利用橢圓和圓同時(shí)計(jì)算，然后對(duì)其加權(quán)平均(權(quán)重取為α、β，α+β=1)，綜上所述可以建立出探測(cè)器單元之間的距離模型為:

其中:n= Ti－Tj，m= Tk－ Tl。

利用文獻(xiàn)［1］中的附件2中所給的數(shù)據(jù)進(jìn)行計(jì)算:通過Matlab在附件2的數(shù)據(jù)中分別求取最大值得到相關(guān)的數(shù)值為:Ti=277，Tj=169，Tk=430，Tl=402，r=4，a=40，b=15，α = β =，代入式(1)可以解出d=0.28，即探測(cè)器單元之間的距離為 0.28 mm。

4.3 旋轉(zhuǎn)中心模型

在模型準(zhǔn)備中的直角坐標(biāo)系里，假設(shè)CT系統(tǒng)旋轉(zhuǎn)中心在正方形托盤中的位置為(x0，y0)，假設(shè)旋轉(zhuǎn)中心在探測(cè)器單元中心對(duì)應(yīng)的X射線T1=256上，在附件2中找出兩種特殊的旋轉(zhuǎn)方向如圖4(c)和圖4(d)所示，假設(shè)圖4(c)和圖4(d)中X射線對(duì)應(yīng)的探測(cè)單元分別為Ti、Tj，則可以建立旋轉(zhuǎn)中心模型為:

其中:d為探測(cè)器單元之間的距離。

利用附件2中的數(shù)據(jù)進(jìn)行計(jì)算:Ti=216，Tj=223，d=0.28，則通過式(2)，可以計(jì)算出旋轉(zhuǎn)中心坐標(biāo)為( － 9.27，6.27)。

4.4 旋轉(zhuǎn)方向模型

附件2中的數(shù)據(jù)為512×180二維數(shù)字矩陣。共有180列:分別表示該CT發(fā)射—接收系統(tǒng)繞固定旋轉(zhuǎn)中心逆時(shí)針旋轉(zhuǎn)180次;共有512行:分別表示每次旋轉(zhuǎn)時(shí)512個(gè)探測(cè)單元接收的X射線吸收強(qiáng)度。例如:第一列數(shù)據(jù)表示CT發(fā)射—接收系統(tǒng)初始位置射入X射線所得的吸收強(qiáng)度，數(shù)值為0表示X射線未穿過任何介質(zhì)，直接到達(dá)探測(cè)器，數(shù)值不為0則表示經(jīng)過待檢測(cè)介質(zhì)吸收后衰減的射線吸收強(qiáng)度。通過Matlab中的最大值函數(shù)確定出吸收率最大時(shí)所對(duì)應(yīng)的旋轉(zhuǎn)次數(shù)，從而可以得出CT系統(tǒng)在此坐標(biāo)系中圍繞旋轉(zhuǎn)中心旋轉(zhuǎn)的幾種特殊的旋轉(zhuǎn)位置，如圖4(a)至圖4(e)。其中:圖4(a)給出了CT系統(tǒng)在坐標(biāo)系O－xy中的初始位置，圖4(b)給出了CT系統(tǒng)旋轉(zhuǎn)16次位置，圖4(c)給出了CT系統(tǒng)旋轉(zhuǎn)89次的位置，圖4(d)給出了CT系統(tǒng)旋轉(zhuǎn)151次的位置，圖4(e)給出了CT系統(tǒng)旋轉(zhuǎn)180次的位置。

圖4 幾種特殊的旋轉(zhuǎn)位置

由圖4(a)可以看出初始位置時(shí)X射線是由第四象限射到第二象限，由圖4(b)可以看出當(dāng)旋轉(zhuǎn)16次時(shí)橢圓和圓介質(zhì)開始重合，由圖4(c)可以看出X射線與ox軸平行，方向相反，由圖4(d)可以看出X射線與oy軸平行，方向相反，由圖4(e)可以看出旋轉(zhuǎn)180次后正好旋轉(zhuǎn)了180°，X射線是由第二象限射到第四象限。

由于在任意方向時(shí)，探測(cè)器所接收的橢圓和圓的最大吸收強(qiáng)度分別為通過橢圓中心和圓心的射線(如果沒有通過中心，可以選擇鄰近的射線)，如圖5所示。

圖5 通過中心的X射線

對(duì)附件2中的數(shù)據(jù)，利用Matlab中的最大值函數(shù)，查找出每一次旋轉(zhuǎn)過程中探測(cè)器所接收的橢圓和圓的最大吸收強(qiáng)度對(duì)應(yīng)的探測(cè)器單元，并分別記作Tei、Tci，l表示這兩個(gè)探測(cè)單元之間的距離，每次的旋轉(zhuǎn)角度記為θi，i=1，2，…，180(旋轉(zhuǎn)180次)，則可以建立旋轉(zhuǎn)方向模型:

這里，βi表示探測(cè)器所在直線與ox軸之間的夾角(如圖5所示)，初始時(shí)β0=0，l'表示橢圓中心和圓心之間的距離，d為探測(cè)器相鄰單元之間的距離。

利用附件2中的數(shù)據(jù)進(jìn)行計(jì)算，其中:d=0.28 mm，l'=45 mm，分別計(jì)算出一些旋轉(zhuǎn)角度和旋轉(zhuǎn)次數(shù)。

(1)當(dāng)i=1時(shí)，即初始位置時(shí)X射線與oy軸之間的夾角θ1，從附件2中第一列容易找出Te1=274，Tc1可得出初始位置時(shí)X射線與oy軸之間的夾角 θ1= β1=29.64°。

同理可以算出，當(dāng) i=2 時(shí)，θ2= β2=30.99°。

為了簡(jiǎn)便計(jì)算，假設(shè)由第1次到第89次旋轉(zhuǎn)方向均勻，則

(2)當(dāng)i=89時(shí)，顯然，X射線與oy軸之間的夾角為90°。

為了簡(jiǎn)便計(jì)算，假設(shè)由第90次到第151次旋轉(zhuǎn)方向均勻，則

(3)當(dāng)i=151時(shí)，X射線與oy軸之間的夾角顯然為180。

為了簡(jiǎn)便計(jì)算，假設(shè)由第152次到第180次旋轉(zhuǎn)方向均勻，則

(4)當(dāng)i=180時(shí)，X射線與oy軸之間的夾角顯然為180°+β1=209.64°。

綜上所述，可以給出附件2中的180個(gè)旋轉(zhuǎn)方向。

5 結(jié)語

利用Einstein旋轉(zhuǎn)圓盤幾何知識(shí)［7］給出了CT系統(tǒng)中探測(cè)器單元之間的距離計(jì)算方法、旋轉(zhuǎn)中心計(jì)算模型、旋轉(zhuǎn)方向模型。數(shù)值實(shí)驗(yàn)表明，該方法簡(jiǎn)單易行。

［1］中國(guó)工業(yè)與應(yīng)用數(shù)學(xué)學(xué)會(huì).2017年高教社杯全國(guó)大學(xué)生數(shù)學(xué)建模競(jìng)賽賽題［DB/OL］.(2017－09－14)［2017－09－15］.中國(guó)大學(xué)生在線網(wǎng)站，http://www.mcm.edu.cn/html_cn/node/460baf68ab0ed0e1e557a0c79b1c4648.html.

［2］李興東，楊民，李德紅，等.一種改進(jìn)的三維CT系統(tǒng)射線源焦點(diǎn)投影坐標(biāo)測(cè)量方法［J］.光學(xué)學(xué)報(bào)，2011，31(12):112－116.

［3］張聰哲，李曉葦，楊昆.扇束等角型CT與等距型CT的比較研究［J］.CT理論與應(yīng)用研究，2013，22(2):215－223.

［4］敖波，曾亞斌，吳偉，等.等角扇形束投影的CT重建算法研究［J］.失效分析與預(yù)防，2013，8(4):206－211.

［5］張順利，李衛(wèi)斌，唐高峰.濾波反投影圖像重建算法研究［J］.咸陽師范學(xué)院學(xué)報(bào)，2008，23(4):47－49.

［6］吳夢(mèng)秋，程正興.反投影重建算法的改進(jìn)算法［J］.西安工業(yè)學(xué)院學(xué)報(bào)，2004，24(1):76－81.

［7］韓雙明，白秀英.Einstein旋轉(zhuǎn)圓盤幾何性質(zhì)的研究［J］.渭南師范學(xué)院學(xué)報(bào)，2012，27(12):46－48.

A Geometry Calibration Method for Computed Tomography

ZHANG Yong－h(huán)ong
(School of Mathematics and Physics，Weinan Normal University，Weinan 714099，China)

CT，computed tomography，is a technology used to obtain the tomography image of the observation target.However，there are often errors in the installation of CT system，which may affect the imaging quality.It is necessary to calibrate the parameters of the installed CT system.On the premise of equidistance in detector unit，firstly，the distance between the detectors is calculated.Secondly，a rectangular coordinate system is set up with the square tray center as the origin of coordinates to find some special directions of rotation respectively and then find out the center of rotation.Finally，the direction of each rotation is calculated.Numerical experiments show that the model is with high measurement precision and can satisfy the demand of high precision reconstruction in practical application.

Computed Tomography;geometry calibration;rotation center;X－ray

TP391

A

1009－5128(2017)24－0038－05

2017－10－20

陜西省扶持學(xué)科基金資助項(xiàng)目:基于小波的圖像加密算法研究(9);渭南師范學(xué)院教育科學(xué)研究項(xiàng)目:以一個(gè)數(shù)學(xué)教研室為例研究中學(xué)數(shù)學(xué)教師的教育科研能力(2016JYKX018);渭南師范學(xué)院特色學(xué)科建設(shè)項(xiàng)目:數(shù)學(xué)方法在秦東經(jīng)濟(jì)社會(huì)發(fā)展中的應(yīng)用(14TSXK02)

張永紅(1976—)，男，陜西合陽人，渭南師范學(xué)院數(shù)理學(xué)院副教授，理學(xué)碩士，主要從事計(jì)算機(jī)圖形學(xué)及圖像處理研究。

【責(zé)任編輯 牛懷崗】