基于圖像輪廓曲線提取的CT系統(tǒng)參數(shù)標(biāo)定

， ， ，

(1.福建工程學(xué)院 數(shù)理學(xué)院， 福建 福州 350118； 2.福建省大數(shù)據(jù)挖掘與應(yīng)用技術(shù)重點(diǎn)實(shí)驗(yàn)室， 福建 福州 350118) 3.福建工程學(xué)院 數(shù)據(jù)分析研究中心， 福建 福州 350118；

計(jì)算機(jī)斷層成像技術(shù)(computed tomogra-phy，CT)可以在不破壞樣品的情況下，利用樣品對射線能量的吸收特性對生物組織和工程材料的樣品在多個(gè)不同角度下進(jìn)行斷層成像，由此獲取樣品內(nèi)部結(jié)構(gòu)信息，但CT系統(tǒng)安裝時(shí)存在誤差，影響成像質(zhì)量[1], 因此需要借助于已知結(jié)構(gòu)的樣品，對安裝好的CT系統(tǒng)進(jìn)行參數(shù)標(biāo)定。目前國內(nèi)對投影旋轉(zhuǎn)中心(center of rotation, COR)的高精度標(biāo)定已有許多研究[2-4]。但對探測單元之間距離

和旋轉(zhuǎn)角度標(biāo)定的文獻(xiàn)還較少。本文針對平行束CT系統(tǒng)，以圖像輪廓曲線作為出發(fā)點(diǎn)，運(yùn)用幾何知識和線性回歸、最小二乘法等方法，建立了一種有效的參數(shù)標(biāo)定模型，為CT參數(shù)的標(biāo)定提供了一種精度較高的新方法。

1 數(shù)據(jù)的獲取及假設(shè)

本文數(shù)據(jù)來源于2017年全國大學(xué)生數(shù)學(xué)建模競賽A題。①依據(jù)問題，提出如下假設(shè)：1)實(shí)驗(yàn)

① 2017年高教社杯全國大學(xué)生數(shù)學(xué)建模競賽題目.[2017-12-09].http://mcm.blyun.com/front/detailTopic.

所收集的數(shù)據(jù)能客觀反映實(shí)際情況；2)正方形托盤所在的平面平行于探測器平面；3)托盤的旋轉(zhuǎn)軸垂直于托盤所在平面；4)平行入射的X射線在旋轉(zhuǎn)過程中，始終垂直于探測器平面；5)所測量的物體的尺寸遠(yuǎn)大于探測器單元的距離；6)所測量的物體在旋轉(zhuǎn)過程中不會(huì)發(fā)生形變。

2 標(biāo)定模板參數(shù)和CT的接收信息

在正方形托盤上放置兩個(gè)均勻固體介質(zhì)組成的標(biāo)定模板，模板幾何信息如圖1，CT系統(tǒng)使用的X射線從180個(gè)方向照射模板得到的512×180的像素投影圖像如圖2。

圖1 模板幾何示意圖(單位：mm)Fig.1 The template’s geometric diagram (unit: mm)

圖2 模板的投影圖像Fig.2 Projection image of the template

3 CT系統(tǒng)參數(shù)的標(biāo)定

3.1 輪廓曲線的提取

3.1.1 提取圖像輪廓

首先對原始圖像進(jìn)行灰度歸一化，得到灰度值落在[0,1]內(nèi)的灰度圖像。然后利用圖像輪廓區(qū)域的灰度值變化幅度較大的特點(diǎn)，對圖像的灰度矩陣按列做二階差分處理。最后從二階差分后的灰度矩陣中提取出絕對值大于設(shè)定閾值的點(diǎn)，即可提取出圖像的輪廓。在觀察二階差分后圖像的數(shù)值分布和測試多種閾值后，本文將閾值設(shè)為0.01，提取出的輪廓如圖3?？梢灾庇^地看到圖像輪廓由4條曲線組成。

圖3 輪廓圖像Fig.3 Contour image

3.1.2 分離輪廓曲線

在3.1.1中得到輪廓圖像上，有的線粗，有的線細(xì)，這是由于一條輪廓線在一列上對應(yīng)點(diǎn)的數(shù)量不同造成的。當(dāng)一條輪廓線在一列上只對應(yīng)一個(gè)點(diǎn)時(shí)，可以直接以該點(diǎn)坐標(biāo)作為曲線在該列的幾何坐標(biāo)，而當(dāng)一條輪廓線在一列上對應(yīng)多個(gè)點(diǎn)時(shí)，還需進(jìn)一步利用多個(gè)點(diǎn)的信息求出該曲線在這列上的準(zhǔn)確幾何位置。對此本文采用聚類的方法，把同一列上距離小于一定閾值的輪廓點(diǎn)聚為一類，并以聚類的幾何中心作為新輪廓點(diǎn),使得每列只剩下4或3個(gè)點(diǎn)(當(dāng)發(fā)生曲線相交時(shí)，有兩條曲線會(huì)共用一個(gè)點(diǎn)，故一列上只會(huì)有3個(gè)點(diǎn))。在反復(fù)選擇不同的距離閾值進(jìn)行聚類測試后，本文選擇將每列上距離在3個(gè)像素單位長度以內(nèi)的點(diǎn)聚為一類。

最后分離曲線的關(guān)鍵在于判定每列上的3或4個(gè)點(diǎn)是屬于哪條輪廓線。在此使用編號標(biāo)記的方法，來判斷每列點(diǎn)的歸屬。算法利用了本文所用標(biāo)定模板輪廓曲線的兩條特殊性質(zhì)：1)每列最多只有一對曲線會(huì)發(fā)生相交；2)曲線相交后一定會(huì)互相穿過對方。故當(dāng)曲線發(fā)生相交時(shí)，只要交換兩條輪廓線的高度編號，其他情況下保持高度編號不變，即可完成對所有點(diǎn)的判定。

令四元組(elld,ellu,cird,ciru)分別表示橢圓下輪廓線、橢圓上輪廓線、小圓下輪廓線、小圓上輪廓線所對應(yīng)高度編號，p表示曲線交點(diǎn)的高度編號。則算法具體描述如下：

Step1將每列的輪廓點(diǎn)按從下到上的順序標(biāo)記上高度編號1,2,3,4。

Step2通過觀察圖像矩陣的第一列給四元組賦以對應(yīng)的初始值(本文賦值為(1,2,3,4))，并從第2列開始執(zhí)行Step3。

Step3判斷當(dāng)前列是否有4個(gè)輪廓點(diǎn)，若有，令p=5并轉(zhuǎn)入Step5,否則轉(zhuǎn)入Step4。

Step4判斷上一列是否有4個(gè)輪廓點(diǎn)，若有，先計(jì)算上一列中距離最接近的兩個(gè)點(diǎn)高度編號，再交換四元組中與這兩個(gè)點(diǎn)編號相同元素的位置，并將p更新為兩個(gè)高度編號中的最小值。執(zhí)行Step5。

Step5按順序輸出四元組中的元素，輸出時(shí)若元素x的數(shù)值大于p，則輸出x-1。否則輸出x。執(zhí)行Step6。

Step6判斷是否有下一列。若有一列，進(jìn)入下一列并再次執(zhí)行Step3。否則終止算法。

表1中給出使用該算法得到幾個(gè)比較關(guān)鍵列的編號。

表1 關(guān)鍵列的編號表Tab.1 Numbering table of the key columns

最終利用該算法得到的編號，對點(diǎn)進(jìn)行歸類。得到分離效果如圖4，部分坐標(biāo)數(shù)據(jù)見表2。

3.2 標(biāo)定探測器單元之間的距離

由圖1可知，橢圓與小圓的圓心距為45 mm。利用表2的數(shù)據(jù)，分別對elld,ellu和cird,ciru取算術(shù)平均值，即可得到每個(gè)角度下小圓投影的幾何中心Ocir和橢圓投影的幾何中心Oell的坐標(biāo)。再將兩個(gè)中心坐標(biāo)相減，得到兩個(gè)圓心投影的距離d。最后測得小圓中心與橢圓中心投影的最大距離dmax=163.25像素長度。根據(jù)幾何關(guān)系，最大距離dmax約等于圓心距d0，故探測器單元間的距離ds應(yīng)滿足式(1)

圖4 分離的輪廓線Fig.4 Separated contour curves

列號elldellucirdciru1192．0362．5404．5433．02189．0364．0401．0431．03188．0364．0400．5429．04185．5365．0398．0426．5?????178159．0320．080．5108．5179158．0322．082．0110．5180157．0324．584．0112．5

(1)

探測器單元間距離的標(biāo)定值為0.275 7 mm。

3.3 標(biāo)定系統(tǒng)的各個(gè)旋轉(zhuǎn)角度

3.3.1 初步求解旋轉(zhuǎn)角度

圖5 圓心連線投影Fig.5 Projection of the centroid connection line

如圖5，Ocir和Oell別對應(yīng)小圓中心與橢圓中心在投影面X′上的投影。X射線的垂線與圓心連線的夾角θ滿足式(2)，因此可用反余弦函數(shù)求解d對應(yīng)的θ。

(2)

表3 角度分布表Tab.3 Angle distribution table

(3)

3.3.2 旋轉(zhuǎn)角度的修正

在表3可看出，初步解出的角度在180°附近產(chǎn)生明顯抖動(dòng)，故需要對角度的異常值進(jìn)行處理。針對此問題本文利用CT系統(tǒng)掃描時(shí)幾乎是等距旋轉(zhuǎn)的特點(diǎn)，對計(jì)算得到的角度數(shù)據(jù)以列號為自變量，角度作為因變量，對一元線性方程y=kx+b做線性回歸?；貧w圖像如圖6，具體回歸參數(shù)見表4，修正后的部分角度見表5。

圖6 線性回歸效果圖Fig.6 Diagram of linear regression effect

參數(shù)估計(jì)值/(°)置信區(qū)間/(°)b28．7470[28．4631，29．0310]k0．9979[0．9952，1．0006]

注：表中顯著性水平α=0.05

表5 修正后的部分角度表Tab.5 The revised part of the angle table

3.4 標(biāo)定旋轉(zhuǎn)中心位置

如圖7，X射線的方向垂直投影面X′，橢圓的圓心為S，T點(diǎn)為旋轉(zhuǎn)中心所在位置。橢圓中心在投影面上的投影為Oell。r=|ST|表示橢圓圓心到旋轉(zhuǎn)中心的距離，r′代表r的投影長度。φ0為旋轉(zhuǎn)中心與橢圓中心的連線關(guān)于橢圓短軸的夾角。lc為旋轉(zhuǎn)中心到探測器左端點(diǎn)的水平距離。

圖7 橢圓的投影圖Fig.7 Projection diagram of an ellipse

則有橢圓中心的投影滿足式(5)

Oell=lc+r’=lc+rcos(φ0-θ)=

lc+Acosθ+Bsinθ

(5)

其中，A=rcosφ0；B=rsinφ0。

這是一個(gè)關(guān)于cosθ,sinθ的線性函數(shù)，故可對式(5)使用最小二乘法進(jìn)行擬合，擬合圖像如圖8，擬合具體結(jié)果見表6。

圖8 擬合效果圖Fig. 8 Fitting effect diagram

Tab.6Tableoftheellipse’sprojectiontrajectory’sfittingresults

參數(shù)數(shù)值lc259．5727A33．4776B-22．7829均方誤差(MSE)0．0495

若以橢圓中心為原點(diǎn)，橢圓圓心到小圓圓心的連線作為X軸的正方向建立直角坐標(biāo)系，則CT系統(tǒng)旋轉(zhuǎn)中心在正方形托盤中的坐標(biāo)(xc,yc)滿足：

xc=-ds·rcosφ0=-ds·A=-9.228 1 mm

(6)

yc=-ds·rsinφ0=-ds·B= 6.280 1 mm (7)

其中，ds為探測器單元之間的距離。

4 誤差分析

4.1 探測器單元間距離的標(biāo)定誤差

4.1.1 橢圓和小圓投影長度的標(biāo)定誤差

只考慮由測量單元之間間距造成的誤差，其他誤差[5]忽略不計(jì)。如圖9，在感應(yīng)器端，1表示在該位置上探測到了物體，0表示沒探測到物體。因?yàn)樘綔y單元之間是有間距的，若用1的數(shù)量來計(jì)算投影長度，則測得的投影長度l與真實(shí)投影長度lr之間在左右端點(diǎn)上會(huì)存在誤差。設(shè)物體投影右端點(diǎn)誤差為D1，左端點(diǎn)誤差為D2則有l(wèi)r=l+D1+D2。由圖9可得，左、右端點(diǎn)誤差均在一個(gè)像素長度以內(nèi)，即D1、D2∈(0,1)。故在最壞情況下lr=l+2像素長度。

圖9 平行束投影圖Fig.9 Diagram of parallel beam projections

4.1.2 圓心的投影距離的誤差

因?yàn)樽笥叶它c(diǎn)的誤差，真實(shí)投影中心Cpr與測得的投影中心Cp滿足公式(8)，所以在最壞情況，計(jì)算得到的投影中心誤差也僅有0.5像素長度。

(8)

若用D11D21表示橢圓投影的左右端點(diǎn)誤差，D12D22表示小圓投影的左右端點(diǎn)誤差。則兩圓圓心投影的真實(shí)距離dr與測得的距離d滿足式(9)

(9)

由D1,D2∈(0,1)可得

(10)

4.1.3 測量單元間距的標(biāo)定誤差

在投影過程中，除了長度測量存在誤差外，不可能也恰好出現(xiàn)180°或者0°的投影角度，所以真實(shí)圓心距d0與測得兩圓心投影的最大距離dmax應(yīng)滿足如下不等式：

d0·cosφmin-1≤dmax≤d0+1

(11)

其中，φmin為投影方向與圓心連線的絕對值最小的夾角，而CT系統(tǒng)角度測量間距在1°左右，所以有

φmin≤0.5°

(12)

此時(shí)最大相對誤差可近似表示為：

(13)

4.2 旋轉(zhuǎn)角度的誤差

在初步計(jì)算旋轉(zhuǎn)角度時(shí)，先計(jì)算d與dmax的比值，再計(jì)算比值的反余弦函數(shù)值得到角度。其中d的測量誤差是影響角度精度的主要因素。由式(10)、(11)可知，d最大可能有1像素長度的誤差，所以角度誤差Dφ滿足式(14)。

?

(14)

根據(jù)式(14)可用MATLAB做圖，得到結(jié)果如圖10。當(dāng)θ=180°或θ=0°時(shí)，角度誤差上界達(dá)到最大值，此時(shí)Dφ≤ 6.345 0°，所以圖6會(huì)在180°附近發(fā)生明顯地抖動(dòng)。而從總體上來看，80%的角度下Dφ<1°。

圖10 角度誤差的上界與角度的關(guān)系圖 Fig.10 Diagram of the relation between the maxi-mum of the angle error and the angle

4.3 旋轉(zhuǎn)中心的誤差

旋轉(zhuǎn)中心的計(jì)算過程較為復(fù)雜，誤差受到太多因素的影響，難以直接估計(jì)。故只給出實(shí)際測得的精確值和標(biāo)定值之間的誤差，以此來大致衡量旋轉(zhuǎn)中心精度。

表7 旋轉(zhuǎn)中心誤差對照表Tab.7 Comparison table of rotation center errors

5 結(jié)論

通過提取CT掃描圖像的輪廓，分離出標(biāo)定模板上的橢圓和小圓輪廓曲線。然后利用得到的輪廓曲線，計(jì)算出小圓和橢圓的圓心坐標(biāo)。結(jié)合幾何關(guān)系對得到的圓心數(shù)據(jù)進(jìn)行擬合，得到了探測單元間的間距，180個(gè)投影角度和旋轉(zhuǎn)中心的標(biāo)定值。最后對標(biāo)定的參數(shù)誤差進(jìn)行定量分析，得出標(biāo)定值的誤差均在可接受的較小范圍內(nèi)，說明本文采取的標(biāo)定方法具有一定的魯棒性和可行性，可在CT參數(shù)的標(biāo)定方面應(yīng)用推廣。

[1] 王召巴.旋轉(zhuǎn)中心偏移對CT重建圖像質(zhì)的影響[J].兵工學(xué)報(bào),2001(3):323-326.

[2] 孟凡勇.基于投影原始數(shù)據(jù)的CT旋轉(zhuǎn)中心精確確定方法[C].中國體視學(xué)學(xué)會(huì).第十三屆中國體視學(xué)與圖像分析學(xué)術(shù)會(huì)議論文集.北京：中國體視學(xué)學(xué)會(huì),2013:6.

[3] 李保磊,傅健,黃巧珍，等.一種基于正弦圖的工業(yè)CT系統(tǒng)轉(zhuǎn)臺旋轉(zhuǎn)中心自動(dòng)確定方法[J].航空學(xué)報(bào),2009,30(7):1341-1345.

[4] 趙雨晴,叢長虹,胡春紅，等.基于正弦圖的投影旋轉(zhuǎn)中心校正方法[J].圖學(xué)學(xué)報(bào),2017,38(4):596-602.

[5] 鄒慧.工業(yè)CT系統(tǒng)誤差源—角度及尺寸誤差的校準(zhǔn)方法研究[D].杭州:中國計(jì)量學(xué)院,2014.