高海東 楊 民 魏東波

孟凡勇

(北京航空航天大學(xué) 機(jī)械工程及自動(dòng)化學(xué)院，北京100191)

(中國(guó)科學(xué)院過(guò)程工程研究所，北京100190)

近些年來(lái)，隨著計(jì)算機(jī)技術(shù)的飛速發(fā)展和面陣探測(cè)器的出現(xiàn)，錐束X射線計(jì)算機(jī)斷層成像(XCT，X-ray Computed Tomography)技術(shù)日益成為現(xiàn)代無(wú)損檢測(cè)(NDT，Non-Destructive Testing)領(lǐng)域內(nèi)的研究熱點(diǎn)［1－3］.錐束 XCT技術(shù)因其掃描效率高的優(yōu)勢(shì)得到了快速發(fā)展，并成功的應(yīng)用到飛機(jī)渦輪葉片與火箭發(fā)動(dòng)機(jī)等領(lǐng)域的無(wú)損檢測(cè).錐束XCT系統(tǒng)利用Radon正逆變換對(duì)進(jìn)行重建成像［4］，這就要求重建坐標(biāo)系與掃描投影坐標(biāo)系保持一致，吻合的程度直接影響重建CT圖像的質(zhì)量［5］.但在實(shí)際系統(tǒng)中，探測(cè)器不可避免地存在安裝誤差，導(dǎo)致掃描投影坐標(biāo)系偏離重建坐標(biāo)系，影響著最終反投影地址的計(jì)算精度，造成重建圖像出現(xiàn)偽影，丟失圖像細(xì)節(jié)，影響檢測(cè)結(jié)果.因此，對(duì)錐束XCT系統(tǒng)重建幾何坐標(biāo)系進(jìn)行準(zhǔn)確標(biāo)定，是開(kāi)展工業(yè)CT成像研究一項(xiàng)基礎(chǔ)而又十分關(guān)鍵的工作.文獻(xiàn)［5］提出了XCT系統(tǒng)的6個(gè)幾何坐標(biāo)配準(zhǔn)參數(shù)，其中關(guān)于探測(cè)器的幾何配準(zhǔn)參數(shù)有:探測(cè)器繞中間行旋轉(zhuǎn)角η、探測(cè)器繞中間列旋轉(zhuǎn)角ξ及探測(cè)器繞中心射線的扭轉(zhuǎn)角α，并針對(duì)上述3個(gè)配準(zhǔn)參數(shù)進(jìn)行了數(shù)值模擬實(shí)驗(yàn)，發(fā)現(xiàn)探測(cè)器的扭轉(zhuǎn)角對(duì)重建圖像的影響最大，產(chǎn)生了嚴(yán)重的偽影.本文針對(duì)探測(cè)器扭轉(zhuǎn)角α，提出基于球投影橢圓軌跡的測(cè)量與校正方法，使重建坐標(biāo)系達(dá)到Radon變換對(duì)的要求，消除重建偽影，保證重建CT圖像的準(zhǔn)確性.

1 探測(cè)器扭轉(zhuǎn)角α測(cè)量

1.1 探測(cè)器位姿分析

FDK(Feldkamp)重建算法是目前錐束XCT成像技術(shù)的主流重建算法［6］.FDK重建算法是在圖1a中射線源和探測(cè)器構(gòu)造的重建坐標(biāo)系Oxyz中進(jìn)行的，探測(cè)器的投影坐標(biāo)系為xdOdzd，算法理論上要求兩個(gè)坐標(biāo)系的關(guān)系為:Odzd軸平行于Oz軸，Odxd軸平行于Ox軸.然而在實(shí)際的錐束XCT系統(tǒng)安裝中，不可避免地存在機(jī)械安裝誤差，導(dǎo)致Odzd軸不平行于Oz軸，Odxd軸不平行于Ox軸，相當(dāng)于坐標(biāo)系xdOdzd繞y軸旋轉(zhuǎn)了一定的角度，該角度為探測(cè)器的扭轉(zhuǎn)角α，如圖1b所示.扭轉(zhuǎn)角α的存在將給重建斷層圖像帶來(lái)偽影［7］，造成斷層圖像發(fā)生形變，甚至導(dǎo)致特征信息丟失，影響圖像的分辨力和細(xì)節(jié)的有效檢出.

圖1 投影與重建坐標(biāo)系示意圖

1.2 扭轉(zhuǎn)角α測(cè)量方法

本文測(cè)量探測(cè)器扭轉(zhuǎn)角α的方法原理如圖2所示.xOz為理想狀態(tài)下不存在安裝誤差的重建坐標(biāo)系，xdOdzd為探測(cè)器存在安裝誤差扭轉(zhuǎn)角α的探元坐標(biāo)系，標(biāo)志點(diǎn)P為非中心射束平面的空間一點(diǎn)，點(diǎn)P在轉(zhuǎn)臺(tái)的帶動(dòng)下繞轉(zhuǎn)軸旋轉(zhuǎn).點(diǎn)P在探測(cè)器上的投影軌跡為曲線C，由檢測(cè)系統(tǒng)的特性可知，曲線C呈現(xiàn)近似橢圓狀.MN為該橢圓長(zhǎng)軸所在的直線，直線l與探測(cè)器坐標(biāo)系的Odxd軸平行，直線MN與直線l的夾角即為探測(cè)器的扭轉(zhuǎn)角α.

圖2 探測(cè)器扭轉(zhuǎn)角α測(cè)量原理圖

令點(diǎn) P 第 i幀投影坐標(biāo)為 pi(xd－i，zd－i)，采集點(diǎn) P 的多幀投影得到點(diǎn)集 Q={pi(xd－i，zd－i)}，對(duì)點(diǎn)集Q進(jìn)行橢圓擬合求出探測(cè)器的扭轉(zhuǎn)角α.常用的橢圓擬合方法主要有3類(lèi):基于最小二乘的擬 合 法［8］、基 于 不 變 矩 的 方 法［9］以 及 基 于HOUGH變換的橢圓擬合方法［10］.基于最小二乘的橢圓擬合方法適用性強(qiáng)，是曲線擬合常用方法.本文擬選擇基于對(duì)稱(chēng)性的最小二乘法擬合橢圓方程，進(jìn)而求出長(zhǎng)軸與坐標(biāo)軸的夾角.橢圓存在2條正交對(duì)稱(chēng)軸，其邊界上至少存在2個(gè)點(diǎn)其外法線方向相反，這2個(gè)點(diǎn)互為對(duì)偶點(diǎn)，對(duì)偶點(diǎn)具有平移、旋轉(zhuǎn)、縮放不變性［11］.基于橢圓的上述特性，對(duì)擬合點(diǎn)的采樣進(jìn)行優(yōu)化，避免傳統(tǒng)算法的大量無(wú)效采樣，降低橢圓擬合的計(jì)算量.對(duì)點(diǎn)集Q，利用對(duì)偶點(diǎn)切線方向相同的性質(zhì)選取對(duì)偶點(diǎn)，增大采樣點(diǎn)落在同一個(gè)橢圓上的概率，得到采樣點(diǎn)集Q'，Q'滿足 Q'?Q.

橢圓的一般方程為

令向量 β=［x2，xy，y2，x，y］，則橢圓方程為

向量 θ=［a，b，c，d，e，f］決定了橢圓的形狀參數(shù)，通過(guò)最小二乘擬合來(lái)尋找適合的參數(shù)向量θ，使各點(diǎn)到候選橢圓邊緣的距離平方和最小.一般橢圓的5 個(gè)參數(shù)描述為:S={A，B，Xc，Yc，α}.其中A，B 表 示 半 軸 長(zhǎng)，(Xc，Yc)表 示 中 心 坐標(biāo)，α(－π/2＜α＜π/2)表示長(zhǎng)軸與坐標(biāo)軸正向的夾角.利用坐標(biāo)的平移及旋轉(zhuǎn)關(guān)系，得出 α與 θ=［a，b，c，d，e，f］的關(guān)系:

在上述思想的指導(dǎo)下，實(shí)驗(yàn)過(guò)程中可用鋼珠在探測(cè)器成像平面xdOdzd上的投影質(zhì)心代替標(biāo)志點(diǎn) P，第 i幀投影質(zhì)心(xd－i，zd－i)的計(jì)算方法為

式中，m，n為投影圖像尺寸;f(xd，zd)為坐標(biāo)(xd，zd)處的投影值.綜上所述，探測(cè)器位姿參數(shù)測(cè)算方法的流程為:采集多幀鋼珠投影Q，利用式(4)求取每個(gè)鋼珠投影的質(zhì)心坐標(biāo)(xd－i，zd－i)，根據(jù)鋼珠投影質(zhì)心坐標(biāo)采用基于對(duì)稱(chēng)性的最小二乘法進(jìn)行橢圓擬合，得到橢圓參數(shù)向量θ，利用式(3)求得探測(cè)器的扭轉(zhuǎn)角α.

2 校正方法

由于探測(cè)器的安裝誤差，導(dǎo)致探測(cè)器并非水平豎直，從而造成重建圖像偽影，在這種情況下需要進(jìn)行校正.如圖3所示為引入探測(cè)器扭轉(zhuǎn)角α錐束XCT重建坐標(biāo)系統(tǒng)示意圖.

錐束射線視為沿z軸方向不同傾斜角度的扇束射線堆積而成，中心平面sOt上的數(shù)據(jù)重構(gòu)屬于扇束FBP(Filtered Backprojection)精確重建，對(duì)于非中心平面的重建，通過(guò)對(duì)扇束FBP重建公式進(jìn)行修正，得到FDK重建公式.設(shè)待重建衰減系數(shù)分布函數(shù)為f(s，t，z)，視角 β下的投影值表示為 p(s，z，β).物體的衰減系數(shù)分布 f(s，t，z)可由錐束射線掃描FDK算法表示為

h(s)為斜坡卷積函數(shù);pe(s，z，β)為等效投影:

令(s1，z1)為經(jīng)過(guò)待重建點(diǎn)(s，t，z)的射線與探測(cè)器的交點(diǎn)，U為加權(quán)因子，由圖3所示幾何關(guān)系，引入探測(cè)器扭轉(zhuǎn)角α，則有

其中

當(dāng)探測(cè)器扭轉(zhuǎn)角α=0°時(shí)，式(8)～式(10)退化為原始方法計(jì)算的投影地址.

由于上述方法改善了投影地址的計(jì)算，使反投影過(guò)程更加準(zhǔn)確，從而消除了探測(cè)器扭轉(zhuǎn)角對(duì)重建圖像的影響.

3 實(shí)驗(yàn)結(jié)果

3.1 探測(cè)器扭轉(zhuǎn)角測(cè)量

實(shí)驗(yàn)條件:鋼珠直徑d=10 mm置于轉(zhuǎn)臺(tái)上，焦距D=1200 mm，轉(zhuǎn)臺(tái)距探測(cè)器200 mm，鋼珠距轉(zhuǎn)臺(tái)中心d=75 mm，高度h=57 mm，探測(cè)器尺寸為1920×1536 像素，探元尺寸:0.127 mm/像素，采集投影180幀.采用1.2節(jié)方法測(cè)量扭轉(zhuǎn)角，擬合得參數(shù)向量θ各元素:

由式(3)可得長(zhǎng)軸與 Odxd軸的夾角為:α=0.0227 rad，換算成角度為1.3006°，即扭轉(zhuǎn)角 α=1.3006°.圖4a所示為取對(duì)數(shù)后鋼珠疊加投影.圖4b所示為投影質(zhì)心軌跡，虛線為橢圓長(zhǎng)軸所在直線，實(shí)線與探測(cè)器的Odxd平行，放大部分示意橢圓長(zhǎng)軸與探測(cè)器坐標(biāo)軸的夾角α.

圖4 探測(cè)器位姿參數(shù)測(cè)算結(jié)果

3.2 重建校正

利用上述得到的探測(cè)器扭轉(zhuǎn)角，通過(guò)第2節(jié)改進(jìn)的FDK重建算法進(jìn)行校正.實(shí)驗(yàn)條件:射線源焦點(diǎn)尺寸0.2mm，管電壓150kV，管電流2mA，焦距D=1200 mm，轉(zhuǎn)臺(tái)距探測(cè)器200 mm，探測(cè)器探元尺寸0.127 mm，被檢對(duì)象為50 mm×50 mm×100 mm的Al，采集投影360幅.

圖5a所示為利用原始FDK重建算法重建第30層得到的斷層圖像，圖5b為在計(jì)入探測(cè)器扭轉(zhuǎn)角之后利用本文方法校正得到的重建斷層圖像.改進(jìn)算法對(duì)探測(cè)器的扭轉(zhuǎn)角進(jìn)行了有效校正，使得反投影地址得以精確計(jì)算，消除了由于探測(cè)器扭轉(zhuǎn)角產(chǎn)生的偽影.對(duì)比圖5a與圖5b可以發(fā)現(xiàn)，由探測(cè)器扭轉(zhuǎn)角引起的圖像邊緣模糊得到了很好的校正，斷層圖像的直角特征得到了保留.

圖5 校正前后重建斷層圖像對(duì)比

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

本文討論了錐束XCT掃描系統(tǒng)中投影坐標(biāo)系與重建坐標(biāo)系的配準(zhǔn)問(wèn)題，針對(duì)探測(cè)器的安裝誤差，提出了基于球投影的探測(cè)器扭轉(zhuǎn)角的測(cè)量方法，并在三維錐束FDK重建算法的基礎(chǔ)上，結(jié)合空間坐標(biāo)變換原理，推導(dǎo)了相應(yīng)的校正方法.校正結(jié)果表明了扭轉(zhuǎn)角測(cè)量結(jié)果及重建斷層圖像的準(zhǔn)確性.本文方法在實(shí)現(xiàn)過(guò)程中，為保證測(cè)量和校正精度，需要多幀標(biāo)志點(diǎn)的投影，另外探測(cè)器安裝誤差的扭轉(zhuǎn)角的旋轉(zhuǎn)中心也需要正確標(biāo)定，否則會(huì)導(dǎo)致校正效果不佳甚至重建圖像錯(cuò)誤.

References)

［1］Roberts J A，Drage N A，Davies J.Effective dose from cone beam CT examinations in dentistry［J］.British Journal of Radiology，2009，82(973):35－40

［2］Yan Guorui，Tian Jie，Zhu Shouping，et al.Fast cone-beam CT image reconstruction using GPU hardware［J］.Journal of X-Ray Science and Technology，2008，16:225－234

［3］Periago D R，Scarfe W C，Moshiri M，et al.Linear accuracy and reliability of cone beam CT derived 3-dimensional images constructed using an orthodontic volumetric rendering program［J］.Angle Orthodontist，2008，78(3):387－395

［4］Jiang Hsieh.計(jì)算機(jī)斷層成像技術(shù)原理、設(shè)計(jì)、偽像和進(jìn)展［M］.北京:科學(xué)出版社，2006:27－65

Jiang Hsieh.Computed tomography principle，design，artifacts and recent advance［M］.Beijing:Science Press，2006:27－65(in Chinese)

［5］侯穎.錐束XCT的定標(biāo)方法研究［D］.大連:大連理工大學(xué)信息與通信工程學(xué)院，2010:28－42

Hou Ying.Research on calibration method in cone-beam XCT［D］.Dalian:School of Information and Communication Engineering，Dalian University of Technology，2010:28－42(in Chinese)

［6］Miao Hui，Wang Tingting，Zhao Huijuan，et al.Improvement based on FDK reconstruction algorithm of conebeam CT［C］//Proc of SPIE.Belingham:SPIE，2010，7570:1－7

［7］戎軍艷.錐束微CT系統(tǒng)偽影校正與性能測(cè)試［D］.北京:中國(guó)科學(xué)院高能物理研究所，2009:62－64

Rong Junyan.Artifacts correction and performance tests of cone beam micro-CT system［D］.Beijing:Institute of High Energy Physics，Chinese Academy of Science，2009:62－64(in Chinese)

［8］鄒益民，汪渤.一種基于最小二乘的不完整橢圓擬合算法［J］.儀器儀表學(xué)報(bào)，2006，27(7):808－812

Zou Yimin，Wang Bo.Fragmental ellipse fitting based on least square algorithm［J］.Chinese Journal of Scientific Instrument，2006，27(7):808－812(in Chinese)

［9］Mehmet Ali Aktas.Measuering shape ellipticity［J］.Computer A-nalysis of Image and Patterns，2011，6854:170－177

［10］Hu Xinyu，Chen Zuobing，Zhang Daode，The application of randomized hough transform in ellipse image detection［J］.Advanced Material Research，2010，159:388－392

［11］呂洪赫，姚振杰，易衛(wèi)東.基于對(duì)稱(chēng)性的最小二乘擬合隨機(jī)橢圓檢測(cè)法［J］.電子測(cè)量技術(shù)，2011，34(5):37－41

Lü Honghe，Yao Zhenjie，Yi Weidong.Random least square fitting elliptic detection algorithm based on symmetry［J］.Electronic Measurement Technology，2011，34(5):37－41(in Chinese)