張 震， 宋衛(wèi)東， 張豐收

(河南科技大學 醫(yī)學技術與工程學院， 河南 洛陽 471000)

錐束X-CT系統(tǒng)相比傳統(tǒng)CT系統(tǒng)有許多優(yōu)點，隨其不斷發(fā)展，在醫(yī)學成像和工業(yè)無損檢測中占據(jù)著重要地位。其成像質(zhì)量受制于安裝誤差，安裝誤差如果不加以校正，會引起重建圖像出現(xiàn)環(huán)形偽影[1]，所以對錐束X-CT系統(tǒng)進行幾何參數(shù)標定是很有必要的。近年來各種標定方法層出不窮，其中Sun和侯穎等[2,3]設計了一個鉆有4個小孔的金屬板作為校準模板,并通過在一個投影角度下得到投影數(shù)據(jù)，利用解析法得出相關參數(shù)。此方法前提是需要將模板擺放到旋轉(zhuǎn)臺的中心位置，有一定的精度要求。劉曉鵬等[4]提出了利用細琴弦和帶豁口的鋼尺得到幾何參數(shù)。Meng[5]和Yang等[6]提出了無模體標定的方法，陳煉等[7]提出了利用一個簡單的校正模型得到少量的投影，計算出幾何參數(shù)的方法，此方法在計算之前運用了一定的假設。李曉晴等[8]提出了平行雙絲法來獲取幾何參數(shù)，此方法雖然模體制造簡單、計算速度較快，但是此模體求得的參數(shù)不全，而且在計算參數(shù)時默認探測器無傾斜角和扭轉(zhuǎn)角，實際安裝中難以實現(xiàn)。王玨等[9]在點模型求解方法上進行了改進，有效降低了在求解幾何參數(shù)時對標定模體安裝精度的要求。黃秋紅等[10]研究了在不考慮成像系統(tǒng)畸變的情況下，靜態(tài)系統(tǒng)射線源陣列的幾何參數(shù)標定方法。楊帥和周凌宏等[11,12]分別從圖像質(zhì)量以及成像系統(tǒng)方面研究幾何參數(shù)標定方法，楊帥等人提出的方法比較耗時，周凌宏等人提出的方法對標定環(huán)境因素要求較高。

1 錐束CT系統(tǒng)掃描結構及其關鍵幾何參數(shù)定義

錐束CT系統(tǒng)的掃描結構如圖1所示，左側(cè)O-XYZ坐標系為平板探測器，中間A-MNL坐標系為待檢測物體，右側(cè)S為射線源。將探測器中間行和中間列相交的焦點定義為中心，用右手法則建立O-XYZ坐標系，以轉(zhuǎn)軸為L軸建立A-MNL坐標系，L軸為旋轉(zhuǎn)中心軸。S(sx,sy,sz)定義為射線源焦點的坐標，并定義S在探測器上的投影坐標為P(px,O,pz)。DAO為探測器到旋轉(zhuǎn)中心的距離，DSA為旋轉(zhuǎn)中心到射線源的距離，DSO為射線源到探測器的距離。并且存在射線源主射線SO與旋轉(zhuǎn)軸垂直。理想狀態(tài)下，投影坐標P(px,O,pz)在O-XYZ坐標系的原點上。主射線SO是L軸和Z軸的公垂線，探測器不存在任何的繞X軸傾斜、繞Z軸扭轉(zhuǎn)以及在XOY平面上的旋轉(zhuǎn)。但是由于安裝誤差不可避免，制造誤差也會導致出現(xiàn)探測器傾斜、扭轉(zhuǎn)以及射線源焦點投影不在探測器中心等情況，這樣會造成重建結果出現(xiàn)偽影，直接影響重建結果的精確度。如圖2所示，定義α角為探測器繞Z軸旋轉(zhuǎn)的角度，β角為繞X軸旋轉(zhuǎn)的角度，γ角為探測器在XOZ面旋轉(zhuǎn)的角度，ΔX為探測器在X軸方向的誤差，ΔZ為探測器在Z軸上的誤差，需要求得的參數(shù)為DSA、DSO、DAO、α、β、γ、ΔX、ΔZ。

圖1 錐束CT系統(tǒng)結構圖及其相關幾何參數(shù)Cone beam CT system structure diagram and the related geometric parameters

圖2 探測器安裝誤差示意圖Detector installation error

2 四橢圓法基本原理

圖3 錐束CT系統(tǒng)標定體模三維模型1. 可調(diào)節(jié)小鋼珠，2. 有機玻璃管，3. 固定鋼珠，4. 固定底座 Cone beam CT system calibration phantom three-dimensional model1. Adjustable small steel ball, 2. plexiglass tube, 3. fixed steel ball, 4. fixed base

本文中提出的方法是基于圓軌道掃描為假設前提，其使用的標定模體如圖3所示，需要兩個有機玻璃管(2)，每個有機玻璃管內(nèi)包含兩個直徑為2 mm的小鋼珠(1)。由于2個有機玻璃管內(nèi)各有兩個小鋼珠，鋼珠的運動軌跡為4個橢圓，故稱為四橢圓法。在設計時，標定模體的2個有機玻璃管應盡量保證與旋轉(zhuǎn)平臺垂直，因為其垂直度會影響成像的質(zhì)量，如果誤差較大，會導致計算所得參數(shù)誤差增大。為保證垂直度要求，文中標定模體設有2個有機玻璃管的固定底座，底座上有2個和有機玻璃管外直徑相等的裝配孔，裝配孔采用基軸制加工，以保證玻璃管與固定底座的垂直度。其次應保證固定底座的平面度，平面度可以用機加工設備來保證。其中一個有機玻璃球管遠離旋轉(zhuǎn)中心，另外一個有機玻璃球管離旋轉(zhuǎn)中心距離較近，旋轉(zhuǎn)中心與兩有機玻璃球管三者構成不等邊三角形，且4個小鋼珠中有2個為固定，另外2個為可調(diào)節(jié)鋼珠，適用于任意錐束角不同的錐束CT系統(tǒng)的參數(shù)標定。四橢圓法獲取全參數(shù)的基本思想為：將制成的標定模體放在待測物體轉(zhuǎn)臺上，首先把標定模體上2個有機玻璃管旋轉(zhuǎn)到主射線的任意一側(cè)，固定鋼珠一端靠近射線源，調(diào)節(jié)可調(diào)節(jié)鋼珠位置，使4個鋼珠在探測器上的投影為2個點，并保證2個有機玻璃管與轉(zhuǎn)臺平面垂直，轉(zhuǎn)臺旋轉(zhuǎn)一周，獲取標定模體360°投影圖像，然后對投影圖像進行數(shù)據(jù)處理，計算出需要得到的各個標定參數(shù)。

2.1 探測器投影視圖參數(shù)預處理

(1)

直線AB的長度lAB:

(2)

圖4 理想探測器與實際探測器投影軌跡示意圖Schematic diagram of the projected path of the ideal detector and the actual detector

直線CD的長度lCD:

(3)

直線AC的長度lAC:

(4)

直線BD的長度lBD:

(5)

2.2 射線源焦點在探測器上投影橫坐標ΔX、射線源到探測器的距離DSO的計算方法

本方法所采用的標定體模為雙有機玻璃和4個小鋼珠所組成，由于2個有機玻璃管分別固定在固定底座相對于旋轉(zhuǎn)中心非對稱的不同位置，因此，在標定模體順時針和逆時針旋轉(zhuǎn)一周時，必定有且僅有2個位置使小鋼珠在探測器上運動軌跡相重合。如圖5所示，A0代表標定模體的初始位置，θ1和θ2分別表示在兩個軌跡投影位置重合時標定模體所旋轉(zhuǎn)的角度。S表示射線源的焦點，A1和A3表示在第一個位置軌跡重合的情況，D1表示在第一個位置重合時X射線源在探測器上形成的坐標，A2和A4表示在第二個位置上軌跡重合的情況，D2表示第二個位置重合時X射線源在探測器上形成的坐標位置。在探測器理想的情況下，A點、D點應與D1點重合，B點、C點應與D2點重合，但是由于安裝誤差，實際中并不相重合，為了減少計算誤差，D1、D2點的橫坐標我們分別取投影點A、D和B、C橫坐標的平均值。

圖5 px和Dso的計算模型示意圖Schematic diagram of the calculation model of px and Dso

過旋轉(zhuǎn)中心A分別做SD1和SD2的垂線，則幾何關系可得：

(6)

(7)

由于初始位置不同，φ角的表達式不同，但是計算思路相同，本文只寫出一種情況。

則射線源到平板探測器的距離：

(8)

2.3 探測器的扭轉(zhuǎn)角α、傾斜角β的計算方法

圖2中A、B、C、D代表兩次投影重合的4個點。由2.2節(jié)選取投影位置相重合處，得到側(cè)視圖(圖6)，由圖2空間幾何關系得知，直線AB、CD、AC、BD的長度以及傾斜度只與α和β角的大小有關，與其它參數(shù)無關。則可令γ=0°，ΔX=0，ΔZ=0。在此關系下推導α和β的表達式，如圖6所示，可得以下式：

(9)

(10)

由式(9)和式(10)相除可得到：

(11)

在式(11)中，∠MAS=π/2+∠A′MA-(1/2)∠S，∠AMS=π/2-∠A′MA，∠ADS=π/2-∠A′MA-(1/2)∠S。

其中，∠S可由標定模體中4個小鋼珠位置相對關系求出，AD為已知。

圖7為俯視圖，過M點做ME垂直于SO并交SO于E。在圖7中可得如下幾何關系：

(12)

(13)

(14)

(15)

由式(14)和式(15)相除可得到：

(16)

A″M=A″S-MS=OS/cosφ-MS

(17)

其中,SO已知，φ和MS已知，則可聯(lián)立式(11)、式(13)、式(16)、式(17)求得α角，同理可求得β角。

圖6 ADS面視圖View(seen from ADS)

圖7 α角計算示意圖Schematic diagram of α angle solution

2.4 探測器在ZOX平面內(nèi)旋轉(zhuǎn)γ角的計算方法

γ角為探測器在ZOY平面內(nèi)的旋轉(zhuǎn)角，在本方法所用的模體中，模體旋轉(zhuǎn)一周，得到運動模體運動軌跡參數(shù)。在探測器繞p點旋轉(zhuǎn)γ角度時，投影圖像如圖8所示。實線表示實際投影圖形，虛線表示理想投影圖形。J、K點分別為上下橢圓的中心為已知，則可由J、K點坐標求得γ角。

(18)

圖8 γ角計算示意圖Schematic diagram of γ angle solution

2.5 射線源到旋轉(zhuǎn)中心的距離DSA及旋轉(zhuǎn)中心到探測器的距離DAO計算方法

在旋轉(zhuǎn)過程中，標定體模4個小球所確定的平面存在總是兩次與錐束主射線垂直，小球在探測器上的投影總存在以下幾何關系，如圖9所示。

圖9 DSA計算模型示意圖Schematic diagram of the calculation model of DSA

直線A1A2表示在第一次4個小球確定的平面與折射線垂直時的情況，直線A4A3表示第二次4個小球確定的平面與折射線垂直時的情況。兩個小球之間的法向距離為d，其測量實現(xiàn)方法為已知鋼珠球的直徑為amm，有機玻璃管壁厚bmm。采用數(shù)碼游標卡尺多次測量2個玻璃管外壁之間的距離，記錄測量次數(shù)和每次測量值，最后求取平均值x，則d=x-a-2b。OA2=R，OA2=r，則由相似三角形定理可得：

(19)

DAO=DSO-DSA

(20)

2.6射線源焦點在探測器上投影縱坐標ΔZ的計算方法

在2.2中提出，在標定過程中存在兩處在不同有機玻璃管的2個小球投影重疊的情況，則以在投影重疊時為基準面建立平面直角坐標系時，可得到圖10中的幾何關系。

圖10 ΔZ計算模型示意圖Schematic diagram of the calculation model of ΔZ

Z為旋轉(zhuǎn)中心A在探測器的投影點，C點為過A點到直線SD1的垂足點，根據(jù)上圖幾何關系可知： △SAC≈△SD1Z。

(21)

由上可知需先求出SD1和AC的長度。

(22)

(23)

根據(jù)海倫-秦九韶公式可得：

l周長=AA3+A3A1+A1A

(24)

AA3、A3A1、A1A均可由標定模體位置關系得出：S△A3AA1=

(25)

(26)

由式(24)、式(25)可求得AC的長，由式(22)、式(26)可求得SD1的長，由式(21)可以求得D1Z的長。則：

(27)

3 標定結果實驗驗證

為了驗證本方法求得幾何參數(shù)的準確性，制作了標定模體，并在錐束CT系統(tǒng)上進行了相關實驗驗證。為了更精確擬合橢圓軌跡參數(shù)，小鋼珠的直徑不宜選擇太小，本文中小鋼珠采用直徑為2 mm的鋼球。實驗平臺主要器材及參數(shù)為：X射線源采用東芝E7242X型X射線球管，最大電壓為125 kV，角度為14°。平板探測器采用VARIAN PaxScan2520型數(shù)字影像探測器，有效面積為250 mm×250 mm，像元矩陣為1536×1920，單個像素尺寸為127 μm。旋轉(zhuǎn)中心臺由日本SANMOTION生產(chǎn)的R2AA08075FXP11W型交流伺服電動機和SYNTEC生產(chǎn)的CNC Controller數(shù)控系統(tǒng)構成。

3.1 具體實驗操作方案

1) 把標定模體垂直放置在旋轉(zhuǎn)臺上，旋轉(zhuǎn)一定角度，使兩個有機玻璃管位于主射線的同一側(cè)，固定位置鋼珠所在的有機玻璃管靠近X射線源一側(cè)，然后調(diào)整標定模體，使兩個不同有機玻璃管內(nèi)的小鋼珠投影重合。標定體模調(diào)整結束，開始標定。

2) 使旋轉(zhuǎn)臺順時針和逆時針旋轉(zhuǎn)一周，得到4個小球的運動軌跡。根據(jù)文中提到的計算方法擬合橢圓參數(shù)方程，并計算出4個相應的交點坐標。

3) 記錄旋轉(zhuǎn)臺旋轉(zhuǎn)到兩次投影重合時的角度，計算相應的ΔX和射線源到探測器之間的距離DSO。

4)根據(jù)文中計算α、β、γ的方法，計算出結果，并校準α、β、γ的角度。

5)記錄4個小鋼珠所在的平面與探測器平行時，標定體模上4個小鋼珠的投影橫坐標，并計算射線源到旋轉(zhuǎn)中心的距離DSA以及旋轉(zhuǎn)中心到探測器的距離DAO。

6)在小鋼珠投影重合時根據(jù)幾何關系，計算出最后一個幾何參數(shù)ΔZ。

圖11為搭建的實驗平臺。

圖11 實驗平臺實物圖Experimental platform

3.2 實驗結果

從表1中可知，由本文中提到的基于四橢圓錐束CT系統(tǒng)幾何參數(shù)的算法得出來的計算值與設計時的真實值相比較，求得的全參數(shù)誤差相對很小，滿足實際的應用需求。

為了驗證本方法的實際效果，利用河南科技大第一附屬醫(yī)院實際人體口腔數(shù)據(jù)進行重建，其重建算法為圓軌跡濾波反投影FDK重建算法，重建過程為：1)對得到的投影數(shù)據(jù)進行加權預處理，即對得到投影數(shù)據(jù)與濾波函數(shù)進行卷積操作；2)把預處理過的數(shù)據(jù)進行一維的斜坡濾波；3)對濾波后的數(shù)據(jù)做錐形束的加權反投影。圖12和圖13為校正前和校正后重建的圖像，從圖中可以看出校正前圖像環(huán)影較重，而且細節(jié)處模糊不清，輪廓處有偽影。校正后偽影基本消除，清晰度較高。為了更直觀地分析結果，計算出了兩個圖像的峰值信噪比。其計算流程為：1)把兩個重建的圖像分別導入MATLAB中，并添加相同的椒鹽噪聲；2)生成系統(tǒng)預定義的3×3濾波器對兩張?zhí)砑釉肼暤膱D像進行均值濾波；3)求得兩張圖像的均值誤差MSE；4)求得兩張圖像的峰值信噪比PSNR。最后求得參數(shù)未校正圖像的峰值信噪比為11.0487，參數(shù)校正后圖像的峰值信噪比為14.0413。顯然參數(shù)校正后得到的圖像比未校正得到的圖像質(zhì)量更好，可滿足正常醫(yī)用的要求。

表1 實際幾何參數(shù)與計算幾何參數(shù)比較

圖12 參數(shù)未校正重建的圖像Image reconstructed with uncorrected parameter

圖13 本文方法參數(shù)校正之后重建的圖像Image reconstructed after parameter correction in this paper

4 結論

基于小球的橢圓軌跡法以及平行雙絲法，本文設計了一種通用標定模體并提出基于本標定模體的幾何全參數(shù)解析算法，從成像敏感參數(shù)出發(fā)環(huán)環(huán)相扣，得到的幾何參數(shù)比較精確。改善了目前標定方法在解析計算標定過程中標定模體的單一性、標定參數(shù)不全、參數(shù)獲取過程中計算較為復雜等缺點。本文設計的這種可調(diào)節(jié)的通用標定模體，可以適用于任意不同錐角的幾何參數(shù)標定。在標定前只需要相應調(diào)整標定模體，標定過程中，通過簡單的幾何數(shù)學模型，即可獲取錐束X-CT系統(tǒng)中的幾何參數(shù)。實驗結果表明，該方法對標定錐束CT系統(tǒng)幾何參數(shù)是有效的，使用的標定模體通用性強、結構簡單，適用于各種CT參數(shù)標定。但是目前使用該方法存在一個前提：標定前需要調(diào)整模體，使兩個小鋼珠投影相重合，這一條件是非常難以精確滿足的。因此今后我們將進一步研究如何得到一種與標定體模初始狀態(tài)無關的幾何參數(shù)標定方法。