陳志祥 須穎 石錦洋

摘要：工業(yè)計(jì)算機(jī)斷層成像技術(shù)（computer tomography，CT）系統(tǒng)幾何參數(shù)不準(zhǔn)確會(huì)影響重建圖像的質(zhì)量和尺寸測(cè)量的精度。為校正工業(yè) CT 幾何位姿失配參數(shù)，提出基于標(biāo)準(zhǔn)球板模體的幾何校正方法。首先，利用球板在兩個(gè)不同位置的投影圖求解出射線源到探測(cè)器的距離以及射線源到運(yùn)動(dòng)平臺(tái)的距離。再在360°范圍內(nèi)間隔一定角度均勻采集球板投影信息，擬合球心投影軌跡，解析求解探測(cè)器偏角η、φ和探測(cè)器中心點(diǎn)坐標(biāo)。最后根據(jù)重建出的體素模型的球心距測(cè)量結(jié)果，利用牛頓迭代法，以工業(yè) CT 尺寸測(cè)量最大允許誤差為約束條件迭代求解探測(cè)器偏角θ。實(shí)驗(yàn)結(jié)果顯示：該文幾何校正方法可以求解出所有工業(yè) CT 幾何位姿失配參數(shù)。參數(shù)校正后，提升圖像質(zhì)量，消除重建圖像偽影，且工業(yè) CT 球心距測(cè)量最大誤差由44.8μm 降低到16.3μm，表明該幾何校正方法可以滿足實(shí)際工業(yè) CT 系統(tǒng)校正的要求。

關(guān)鍵詞：工業(yè) CT;幾何校正;校正模體;圖像重建

中圖分類號(hào)： TP391.4文獻(xiàn)標(biāo)志碼： A文章編號(hào)：1674–5124（2021）12–0023–06

Geometric parameter calibration of industrial CT system

CHEN Zhixiang，XU Ying，SHI Jinyang

（School of Mechanical and Electrical Engineering， Guangdong University of Technology， Guangzhou 510006， China）

Abstract： Incorrect geometric parameters of the industrial CT （computer tomography） system will affect the qualityof thereconstructedimageandtheaccuracyof dimensionalmeasurement. Inordertocalibrate geometric mismatching parameters of industrial CT， a geometric calibration method based on the standard spherical plate phantom was proposed. Firstly， the distance from the ray source to the detector and the distance from the ray source to the rotating platform were calculated by using the projection of the spherical plate at two different positions. Then， the projection information of the spherical plate was uniformly collected at a certain angle interval within 360° range， fitted the sphere center projection trajectory， and the detector deflection angle η ， φ and the detector center point coordinates were calculated analytically. Finally， based on the spherical- center distance measurement results of the voxel model， the Newton iteration method was used to solved the detectordeflectionangle θ bytakingtheMaximumPermissibleErrorof theindustrialCTdimensional measurementastheconstraintcondition. Theexperimentalresultsshowedthatthegeometriccorrection method could solved all the geometric position mismatch parameters of industrial CT. After the parameters were calibrated， the image quality was improved， and the reconstructed image artifacts were eliminated. The maximum error of industrial CTspherical distance measurement was reduced from 44.8μm to 16.3μm，indicated that the geometric calibration method could met the requirements of actual industrial CT system calibration.

Keywords： industrial computer tomography; geometry calibration; calibration phantom; image reconstruction

0引言

工業(yè)計(jì)算機(jī)斷層成像技術(shù)（ computer tomography，CT）是一門集光學(xué)、機(jī)械、控制和計(jì)算機(jī)科學(xué)于一體的數(shù)字化 X 射線成像技術(shù)，其利用物體對(duì) X 射線吸收率的不同來(lái)獲取物體的多角度投影圖，并根據(jù)這些投影圖重建出物體的三維圖像[1]。目前這項(xiàng)技術(shù)正在向尺寸測(cè)量領(lǐng)域發(fā)展，但因?yàn)楣I(yè) CT 通常采用 FDK 算法對(duì)投影數(shù)據(jù)進(jìn)行重建，而 FDK（ Feldkamp-Davis-Kress）算法要求工業(yè) CT系統(tǒng)幾何位姿滿足理想條件[2]，實(shí)際情況下，工業(yè) CT 系統(tǒng)難于滿足理想條件，系統(tǒng)幾何參數(shù)會(huì)存在偏差。在圖像重建過(guò)程中，如果不進(jìn)行校正，會(huì)影響工業(yè) CT 尺寸測(cè)量的精度。

針對(duì)這個(gè)問(wèn)題，Noo[3]利用其設(shè)計(jì)的兩點(diǎn)模體，通過(guò)獲取模體多角度的投影數(shù)據(jù)，擬合橢圓運(yùn)動(dòng)軌跡，根據(jù)橢圓的相關(guān)參數(shù)求解系統(tǒng)幾何參數(shù)，在忽略探測(cè)器偏角θ對(duì)其他參數(shù)的影響的情況下，精確求解其余參數(shù)，但實(shí)際工業(yè) CT 系統(tǒng)中偏角θ是一定存在的，且也會(huì)影響圖像的重建質(zhì)量。Jacobson[4]設(shè)計(jì)了一種金屬絲校準(zhǔn)模型，利用兩條不共面的直線在不同角度的投影平面內(nèi)的交點(diǎn)位置的不同，根據(jù)這些交點(diǎn)的軌跡來(lái)計(jì)算出系統(tǒng)幾何參數(shù)，但該方法需要準(zhǔn)確知道校正模體中每一根金屬絲準(zhǔn)確的空間位置，在實(shí)際過(guò)程中，很難精確定位空間中金屬絲的位置，且對(duì)空間中金屬絲線位置測(cè)量的難度明顯高于球板球心位置測(cè)量的難度。Zhao[5]利用含有九個(gè)標(biāo)準(zhǔn)球的正方形模體，采集模體單個(gè)角度的投影數(shù)據(jù)，根據(jù)標(biāo)準(zhǔn)球的投影位置關(guān)系迭代求出系統(tǒng)幾何參數(shù)，雖能通過(guò)迭代求解出系統(tǒng)所有的失配參數(shù)，但該方法僅根據(jù)球板一個(gè)位置的投影關(guān)系求解，其求解精度難以保證。目前，工業(yè) CT 系統(tǒng)參數(shù)校正的方法[6-14]都是利用特定校正模體，采集模體不同角度的投影數(shù)據(jù)，根據(jù)投影數(shù)據(jù)建立起與工業(yè)CT 系統(tǒng)幾何位姿參數(shù)的映射關(guān)系，再求解出相應(yīng)的參數(shù)。

本文采用高精度球板模體，采集球板模體在不同角度，不同位置下的二維投影圖，求出不同情況下球板各球體球心在探測(cè)器上的投影坐標(biāo)，擬合運(yùn)動(dòng)軌跡，再根據(jù)投影關(guān)系，利用解析法精確求解除θ外的失配參數(shù)。之后，用帶參數(shù)的 FDK 重建算法軟件重建出的球板三維圖像，將其導(dǎo)入可視化軟件 VG Studio MAX 3.0中并提取邊緣，測(cè)量球體球心距，分析工業(yè) CT 測(cè)量誤差，以工業(yè) CT 尺寸測(cè)量最大允許誤差為約束條件，利用牛頓迭代法求解系統(tǒng)最后一個(gè)失配參數(shù)。相較于其他校正方法，本方法可求解出工業(yè) CT 系統(tǒng)所有的失配參數(shù)，且用工業(yè) CT 尺寸測(cè)量誤差來(lái)對(duì)比校正前后的效果，更為直觀。

1工業(yè) CT 系統(tǒng)和幾何參數(shù)

工業(yè) CT 系統(tǒng)主要由射線源、旋轉(zhuǎn)平臺(tái)和探測(cè)器構(gòu)成，根據(jù) FDK 算法的要求，理想的工業(yè) CT 系統(tǒng)幾何結(jié)構(gòu)應(yīng)滿足：射線源焦點(diǎn)和探測(cè)器中心點(diǎn)的連線經(jīng)過(guò)旋轉(zhuǎn)軸，且垂直于探測(cè)器平面;旋轉(zhuǎn)軸與探測(cè)器中心列平行。工業(yè) CT 系統(tǒng)理想位姿關(guān)系如圖1所示。其中f 為射線源到旋轉(zhuǎn)軸的距離;d 為射線源到探測(cè)器中心的距離;（u0 ， v0）為射線源焦點(diǎn)在探測(cè)器上的垂直投影點(diǎn)。

為更好地研究工業(yè) CT 系統(tǒng)的幾何位姿關(guān)系，建立三維坐標(biāo)系O? XYZ，如圖1所示。X 軸沿射線方向;Y 軸沿水平方向垂直于射線方向且與探測(cè)器中心行平行;Z 軸由右手坐標(biāo)系原則確定。同樣在探測(cè)器所在平面上建立二維坐標(biāo)系O?uv，u 軸沿探測(cè)器中心行方向;v 軸沿探測(cè)器中心列方向。工業(yè) CT 系統(tǒng)在安裝調(diào)試過(guò)程中，系統(tǒng)位姿關(guān)系一定會(huì)存在偏差，工業(yè) CT 的校正即校正這些參數(shù)。經(jīng)國(guó)內(nèi)外眾多研究人員對(duì)工業(yè) CT 校正的研究，需要校正的參數(shù)可概括為：探測(cè)器繞 X、Y、Z 軸的轉(zhuǎn)動(dòng)角度η、θ、φ;射線源焦點(diǎn)在探測(cè)器上的垂直投影點(diǎn)（u0 ， v0）和在 u、v 方向上的偏移Δu、Δv;源到旋轉(zhuǎn)中心的距離以及源到探測(cè)器中心的距離。其中和探測(cè)器有關(guān)的5個(gè)參數(shù)如圖2所示。

2參數(shù)校正

2.1校正模體

本文使用的校正模體是球板標(biāo)準(zhǔn)器，如圖3所示，由12個(gè)高精度的紅寶石球體和碳纖維平板組成，部分球體之間的球心距已知。

2.2計(jì)算d、f

在保持管電壓、管電流和曝光時(shí)間等掃描參數(shù)不變的情況下，僅改變探測(cè)器到運(yùn)動(dòng)平臺(tái)的距離，采集兩組不同位置時(shí)的球板二維投影圖。擬合球板各球體在二維投影圖上的圓心位置，再計(jì)算圓心距 L，L 為像素?cái)?shù)量與像素大小的乘積，即

式中：L——被測(cè)量，mm;

n——像素?cái)?shù)量;

M——放大比;

D——像元尺寸，mm。

根據(jù)式（1），利用上述采集的兩組投影圖數(shù)據(jù)，利用球板各球體間球心距，即可建立多個(gè)關(guān)于 d 和 f 的二元一次方程，在精確知道運(yùn)動(dòng)平臺(tái)沿射線源方向上的運(yùn)動(dòng)距離的情況下，可求解出多個(gè) d、f值，最后取 d、f 的平均值作為最終校正結(jié)果。

2.3計(jì)算η、u0、v0

采集 N（N 為偶數(shù)）張1號(hào)球和12號(hào)球每隔一定角度的投影圖，計(jì)算出1號(hào)球和12號(hào)球在不同角度投影圖下的球心位置（ui ， vi）（i =1， 2，3， ··· ， N）。并根據(jù)這些球心坐標(biāo)擬合運(yùn)動(dòng)軌跡，可以得到兩個(gè)橢圓（橢圓1和橢圓2），擬合出的橢圓可用橢圓方程式表示如下：

式中：（ ， ）——橢圓中心點(diǎn)坐標(biāo);

a、b、c——橢圓方程參數(shù)。

由式（3）分別計(jì)算出1號(hào)球和12號(hào)球球心連線在每隔180°投影圖下的中心點(diǎn)位置（1 ， 1）、（2 ， 2），如圖4所示。

其中j = i+ N/2;i =1， 2， ··· ， N/2。

探測(cè)器面內(nèi)扭轉(zhuǎn)角η即可由式（4）計(jì)算得出。探測(cè)器平面內(nèi)射線源焦點(diǎn)的垂直投影點(diǎn)（u0 ， v0）可由式（5）、（6）計(jì)算出。

2.4計(jì)算φ

Sun[15]等利用四球模體校正系統(tǒng)幾何參數(shù)，提出平行于轉(zhuǎn)軸的兩球投影圓心距之間滿足以下關(guān)系：

式中：——平行于轉(zhuǎn)軸的兩球投影圓心距，

——實(shí)際圓心距，mm。

根據(jù)這一原理，將 1 號(hào)球到 10 號(hào)球的球心距和 3 號(hào)球到 12 球的球心距 SD3-12 代入式 （7），可計(jì)算出。

2.5 計(jì)算

根據(jù)上述計(jì)算出的探測(cè)器失配參數(shù)對(duì)在理想坐標(biāo)系下的 1 號(hào)球球心坐標(biāo)進(jìn)行坐標(biāo)變換得，則：

其中：

經(jīng)過(guò)坐標(biāo)變換后得到：

射線源焦點(diǎn)坐標(biāo)S =（0， 0， 0）經(jīng)坐標(biāo)變換后，射線源焦點(diǎn)在理想坐標(biāo)系下的坐標(biāo)應(yīng)為S＼=（xs ，ys ，zs）;將 S＼與p 相連交探測(cè)器于點(diǎn)p = v ＼），p ＼即為1號(hào)球球心經(jīng) X 射線投射后在探測(cè)器上的投影點(diǎn)。其中：

同理可計(jì)算出其他球體球心經(jīng) X 射線投射后在探測(cè)器上的投影點(diǎn)坐標(biāo)。

此時(shí)使用帶參數(shù)的 FDK 重建算法軟件，利用上面所計(jì)算出的6個(gè)系統(tǒng)失配參數(shù)，重建出此時(shí)球板模體的三維圖像，再利用 VG Studio MAX 3.0圖像可視化軟件測(cè)量各球體間的球心距，找出此時(shí)球心距誤差最大的兩球體。以1號(hào)球和12號(hào)球?yàn)槔?，此時(shí)兩球球心在探測(cè)器上的投影點(diǎn)間距離為：

式（12）中，除θ外的其他參數(shù)均已知，故S D1?12= f （θ）是一個(gè)關(guān)于θ的函數(shù)，解析求解過(guò)于復(fù)雜，可用牛頓迭代法求解，以工業(yè) CT 尺寸測(cè)量最大允許誤差為迭代約束條件，迭代求解出θ。

3實(shí)驗(yàn)驗(yàn)證

3.1實(shí)驗(yàn)設(shè)備及參數(shù)

為驗(yàn)證上述方法的可行性，在實(shí)際工業(yè) CT 系統(tǒng)中進(jìn)行了實(shí)驗(yàn)，實(shí)驗(yàn)所使用的設(shè)備為天津三英精密儀器有限公司生產(chǎn)的 nanoVoxel-4000，如圖5所示。射線源為開(kāi)管發(fā)射式，最大管電壓為220 kV，最大管電流為1.0 mA;探測(cè)器成像面積為427 mm ×427 mm，像素矩陣為3072×3072;像素細(xì)節(jié)分辨力為500 nm，空間分辨率為2?m。

3.2校正過(guò)程

首先驗(yàn)證 d 和f，在 d 固定的情況下，改變f其他參數(shù)如表1所示，采集兩組投影圖。表1 初始參數(shù)

通過(guò)計(jì)算可得，射線源到探測(cè)器的距離誤差為–5.2029 mm，射線源到旋轉(zhuǎn)軸的距離誤差為–4.3907 mm。

在調(diào)整工業(yè) CT 系統(tǒng) d 和f后，再次掃描校正模體，保持 d=1500 mm，f=500 mm，在其他參數(shù)不變的情況下360°范圍內(nèi)采集球板模體的投影圖，然后根據(jù)本文方法求解圓心，擬合橢圓軌跡如圖6所示，橢圓參數(shù)計(jì)算出工業(yè) CT 其余幾何參數(shù)，最后利用迭代法經(jīng)4次迭代后計(jì)算出工業(yè) CT 需校正的最后一個(gè)參數(shù)。

至此，所有參數(shù)校正完畢，最終校正結(jié)果如表2所示。

3.3校正前后效果對(duì)比

實(shí)際重建過(guò)程中，利用帶參數(shù)的 FDK 重建算法軟件重建球板模體，校正前后的重建圖像對(duì)比如圖7所示，可以看出，校正后球板模體的重建圖像質(zhì)量有了明顯提升，消除了圖像偽影，提升了圖像清晰度。

利用 VG Studio MAX 3.0圖像可視化軟件分別對(duì)校正前后重建的球板模體三維圖像進(jìn)行尺寸測(cè)量，兩次測(cè)量結(jié)果如圖8所示。可知校正后工業(yè) CT 尺寸測(cè)量的最大誤差由44.8?m 降低到16.3?m，效果明顯，本文工業(yè) CT 系統(tǒng)參數(shù)校正的方法降低了工業(yè) CT 尺寸測(cè)量的誤差。

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

本文基于高精度球板模體，提出了一種工業(yè) CT 幾何參數(shù)校正的方法，利用球板在不同角度下的投影圖像，獲取球心投影運(yùn)動(dòng)軌跡，擬合橢圓，再根據(jù)橢圓參數(shù)和工業(yè) CT 尺寸測(cè)量的最大誤差，結(jié)合解析法與迭代法，可求解出工業(yè) CT 系統(tǒng)所有幾何失配參數(shù)。實(shí)驗(yàn)表明，校正后重建圖像的質(zhì)量有了明顯的提升，減少了圖像偽影，校正后工業(yè) CT 尺寸測(cè)量最大誤差由44.8?m 降低到16.3?m，能滿足實(shí)際工業(yè) CT 系統(tǒng)參數(shù)的校正要求。但是，本文幾何校正方法需要高精度的校正模板，需要準(zhǔn)確知道球板模體各球心距，也受限于運(yùn)動(dòng)平臺(tái)回轉(zhuǎn)精度和沿射線源方向的運(yùn)動(dòng)精度，否則會(huì)影響校正參數(shù)的準(zhǔn)確性。

參考文獻(xiàn)

[1]張朝宗.工業(yè) CT 技術(shù)和原理[M].北京：科學(xué)出版社， 2009.[2] FELDKAMP L A， DAVIS L C， KRESS J W. Practical cone-beam algorithm[J]. Journal of the Optical Society of America A， 1984， 1（6）：612-619.

[3] NOO F， CLACKDOYLE R， MENNESSIER C， et al. Analyticmethodbasedonidentificationofellipseparametersfor scanner calibration in cone-beam tomography[J]. Physics in Medicine and Biology， 2000， 45（11）：3489-3508.

[4] JACOBSON M W， KETCHA M D， CAPOSTAGNO S， et al.A line fiducial method for geometric calibration of cone-beam CTsystemswithdiversescantrajectories[J]. Physicsin Medicine and Biology， 2018， 63（2）：025030.

[5] ZHAO J T， HU X D， ZOU J，etal. Geometric parametersestimationandcalibrationincone-beammicro-CT[J]. Sensors， 2015， 15（9）：22811-22825.

[6]周凌宏， 李翰威， 徐圓， 等.錐束 CT 圓軌道掃描的幾何校正[J].光學(xué)精密工程， 2014， 22（10）：2847-2854.

[7] SMEKALLV，KACHELRIEM，STEPINAE，etal.Geometricmisalignmentandcalibrationincone-beam tomography[J]. Medical Physics， 2004， 31（12）：3242.

[8] CHO Y， MOSELEY D J，SIEWERDSEN J H， et al. Accuratetechniqueforcompletegeometriccalibrationof cone-beam computedtomographysystems[J]. MedicalPhysics， 2005，32（4）：968-983.

[9] XIANG Q， WANG J， CAI Y F. A geometric calibrationmethodforconebeamCTsystem[C]//EighthInternational ConferenceonDigitalImageProcessing （ICDIP 2016）. International Society for Optics and Photonics， 2016.

[10] PATEL V， CHITYALA R N， HOFFMANN K R， et al. Self- calibrationofacone-beammicro-CTsystem[J]. MedicalPhysics， 2008， 36（1）：48-58.

[11]張才鑫， 黃魁東， 陳華.基于標(biāo)定量具的快速錐束 CT 幾何校準(zhǔn)方法[J].儀器儀表學(xué)報(bào)， 2019， 40（1）：175-182.

[12] ZHOU K， HUANG Y X， MENG X X， et al. A new method for cone-beamcomputedtomographygeometricparameters estimation[J]. JournalofComputerAssistedTomography， 2016， 40（4）：639-648.

[13] LI G， LUO S H， YOU C Y， et al. A novel calibration method incorporatingnon-linearoptimizationandball-bearing markers for cone-beam CT with a parameterized trajectory[J]. Medical Physics， 2019， 46（1）：152-164.

[14] LI H C， HONG W， LIU Y， et al. A novel method of micro- tomographygeometricanglecalibrationwithrandom phantom[J]. Journal of X-ray science and technology， 2017，25（4）：641-652.

[15] SUN Y， HOU Y， ZHAO F Y， et al. A calibration method for misalignedscannercogeometryinne-beamcomputed tomography[J]. NDT and E International， 2006， 39（6）：499-513.

（編輯：莫婕）