王廣志 童 明 丁 輝

(清華大學(xué)醫(yī)學(xué)院，北京，100084)

?

基于幾何特征的X射線圖像中圓形標(biāo)志點(diǎn)自動(dòng)提取

王廣志 童 明 丁 輝

(清華大學(xué)醫(yī)學(xué)院，北京，100084)

從X射線圖像中快速準(zhǔn)確地提取圓形標(biāo)志點(diǎn)是基于C形臂的手術(shù)導(dǎo)航中進(jìn)行跟蹤的關(guān)鍵環(huán)節(jié)。由于在實(shí)際應(yīng)用中X射線圖像的標(biāo)記點(diǎn)半徑較小，并且有其他的背景噪聲和干擾物，使得一些傳統(tǒng)的圓形檢測方法往往魯棒性不夠高。針對這一問題，本文提出了一種基于幾何特征的圓形目標(biāo)檢測方法，利用這種方法能夠?qū)射線圖像中的圓形標(biāo)記點(diǎn)進(jìn)行自動(dòng)提取。實(shí)驗(yàn)測試表明該方法能夠有效地檢測不同半徑的圓，并且能夠通過調(diào)節(jié)參數(shù)檢測到不標(biāo)準(zhǔn)的圓形。實(shí)驗(yàn)結(jié)果顯示該方法在對C形臂X射線成像中的圓形標(biāo)志物的自動(dòng)提取具有較高的效率及魯棒性。

圓形檢測；曲線擬合；特征提??；標(biāo)志點(diǎn)提取

引 言

由于C形臂X射線機(jī)具有實(shí)時(shí)成像和使用便捷的特點(diǎn)，在心臟和骨科等介入式手術(shù)中被廣泛使用。近年來，基于C形臂的手術(shù)導(dǎo)航系統(tǒng)得到了更廣泛的使用[1-2]。在基于C形臂的手術(shù)導(dǎo)航系統(tǒng)中，為實(shí)現(xiàn)人體和器械的跟蹤，常常通過特殊形狀的標(biāo)記點(diǎn)來表征其位置，因此圖像中標(biāo)志點(diǎn)的自動(dòng)識(shí)別是實(shí)現(xiàn)跟蹤的關(guān)鍵。在實(shí)際系統(tǒng)設(shè)計(jì)中，一般采取不同直徑的金屬小球作為標(biāo)志物，通過檢測不同直徑小球的位置來進(jìn)行系統(tǒng)標(biāo)定和畸變矯正等。由于金屬小球的密度較大，X射線幾乎不能透過，因此在投影圖像中標(biāo)志物的像素灰度值接近0。同時(shí)因?yàn)闃?biāo)志物為球體，不同角度的射線穿過球體的厚度并不一致，所以標(biāo)志物在投影圖像中的影像灰度值并不一致，一般來講，靠近邊緣的灰度值較高，同時(shí)影像中噪聲和干擾目標(biāo)等經(jīng)常導(dǎo)致圖像中小球的邊緣不夠規(guī)則，給自動(dòng)監(jiān)測帶來困難。

圖像中的圓形目標(biāo)檢測可以借助圖像分割方法進(jìn)行，但圖像分割的計(jì)算效率通常無法滿足實(shí)時(shí)性的要求[3]。為快速檢測圖像中的圓形目標(biāo)，實(shí)用化的處理方法主要可以分為兩大類：一種是利用霍夫圓變換的方法[4-5]，另一種是基于圓形的幾何特征進(jìn)行檢測的方法[6-8]。第1類方法的基本思想是將邊緣圖像中的每一點(diǎn)都映射到參數(shù)空間中的一個(gè)區(qū)域并且將其貢獻(xiàn)累加，然后輸出參數(shù)空間的局部最大值作為檢測到的目標(biāo)。這種方法需要較大的存儲(chǔ)量和計(jì)算量，并且在圖像中目標(biāo)較多，圓形半徑較小時(shí)(例如小于10個(gè)像素)，正確的目標(biāo)圓輪廓在參數(shù)空間的累計(jì)量會(huì)較少，從而可能得不到正確穩(wěn)定的檢測結(jié)果。第2類方法往往基于邊緣像素點(diǎn)，利用圓形的幾何特征來檢測輪廓是否為圓形。例如，文獻(xiàn)[9]提出基于曲線擬合的兩步法檢測圓形，在圓形邊緣局部畸變的情況下，能夠正確地檢測到圓形且計(jì)算效率高?？紤]在手術(shù)導(dǎo)航中X射線圖像中目標(biāo)可能較多且其面積較小、像素不均勻，并且目標(biāo)之間無相互遮擋的特點(diǎn)，本文提出一種利用圓的幾何特征檢測輪廓是否為圓的新算法。該算法首先通過邊緣檢測算法得到圖像中的邊緣輪廓，然后計(jì)算輪廓的幾何特征，并利用幾何特征的值判斷輪廓是否為圓形，實(shí)驗(yàn)證明該方法具有較強(qiáng)的魯棒性，并且計(jì)算效率較高。

1 圓形目標(biāo)提取算法

1.1 圓的擬合與圓形特征

1.1.1 圓的擬合

給定一組點(diǎn){pi}={(xi,yi)}(i=1,…,n)，用一個(gè)圓擬合這組點(diǎn)可以通過求下式的最小值

(1)

通過將Q分別對x0,y0,r求偏導(dǎo)，并令其為0，可以求出其閉合解為[10]

(2)

(3)

1.1.2 圓度特征的計(jì)算

在數(shù)學(xué)上圓度定義為[10]

(4)

式中：I為輪廓的長度(pixel)；s為輪廓所包含的面積(pixel2)。由圓度的定義可以得出結(jié)論：當(dāng)輪廓為標(biāo)準(zhǔn)圓時(shí)，其圓度值為1，而輪廓為其他任意形狀，其圓度值均小于1，并且輪廓形狀越接近圓形，其圓度值越接近1。

基于本文所要解決的問題，定義一種變形的圓度特征為

(5)

式中：r為利用式(2)中圓擬合方法求出的半徑;S為輪廓內(nèi)所包含的像素點(diǎn)個(gè)數(shù)。

1.2 圓形目標(biāo)提取算法

圖1 算法流程圖Fig.1 Flow chart of the algorithm

文獻(xiàn)[9,11]分別利用圓度和曲線擬合的方法進(jìn)行圖像中圓形的檢測。利用圓度定位需要確定一個(gè)閾值，其中圓度大于閾值的輪廓被認(rèn)為是圓形。由于在實(shí)際應(yīng)用中閾值的確定受到圓半徑的影響，文獻(xiàn)[11]根據(jù)經(jīng)驗(yàn)來確定閾值，沒有提供一個(gè)自適應(yīng)選擇的方法；同時(shí)，利用圓度對凹形輪廓檢測也可能會(huì)出錯(cuò)，在實(shí)際應(yīng)用中將導(dǎo)致結(jié)果不夠魯棒。利用文獻(xiàn)[9]描述的方法，目標(biāo)輪廓為彎月形時(shí)，算法會(huì)將輪廓誤判為圓形。

考慮到C型臂X射線圖像中圓形目標(biāo)數(shù)量較多，其像素不均勻，并且目標(biāo)之間不相互遮擋的特點(diǎn)，本文設(shè)計(jì)了一種算法，綜合利用曲線擬合以及上述變形的圓度特征判別圖像中的輪廓是否為圓形，算法的流程如圖1所示。算法通過圖像預(yù)處理將原始圖像轉(zhuǎn)為二值圖像，然后通過查找輪廓的方式得到二值圖像中所有的輪廓，并檢測這些輪廓點(diǎn)是否構(gòu)成符合待檢測特征的圓形。預(yù)處理方法可以根據(jù)原圖像的特點(diǎn)選擇相應(yīng)的方法。本文針對C形臂X射線圖像的成像特點(diǎn)及目標(biāo)像素值比背景像素值小的特點(diǎn)，選擇如圖1右邊方框中所示的預(yù)處理流程。

1.2.1 圖像二值化

二值化操作將圖像從灰度圖變換為二值圖像，主要操作方法有：全局取閾值法和局部取閾值法。常用的全局自適應(yīng)閾值法有大津法[12]，此方法通過求類間方差最大得到一個(gè)最優(yōu)的全局閾值，適用于灰度直方圖形態(tài)為雙峰的圖像。局部取閾值法通過每個(gè)像素的鄰域特性確定其閾值。文獻(xiàn)[13]評估了40種不同的取閾值的方法，他們的測試集顯示Sauvola的局部自適應(yīng)二值化方法要比其他的二值化方法有更好的性能。其計(jì)算像素閾值的方法是

(6)

式中：m(x,y)，s(x,y)分別以像素(x,y)為中心，大小為ω×ω的窗中所覆蓋像素的平均值與方差；t(x，y)為像素(x,y)的閾值。然而由于該算法需要計(jì)算每個(gè)窗內(nèi)像素集合的均值和方差，計(jì)算量較大，在實(shí)時(shí)應(yīng)用中速度會(huì)受到限制。本文綜合考慮計(jì)算量和閾值化的性能，采取了一種簡化的局部自適應(yīng)閾值計(jì)算方法。考慮圖像中大小為ω×ω的窗，閾值可表示為

t(x,y)=m(x,y)-c

(7)

式中：c為常量，均值m(x,y)可以用積分圖的方法快速計(jì)算。窗的大小由目標(biāo)大小確定，在X射線圖像標(biāo)志點(diǎn)提取中，標(biāo)志點(diǎn)像的大小一般不超過10×10個(gè)像素，因此可以選擇窗口大小在15×15～50×50之間，對于面積大的目標(biāo)可以選擇更大的窗口。常量c與圖像中目標(biāo)和背景對比度有關(guān)，c取值越大，越多與背景對比度不高的像素點(diǎn)將被移除。在本文的應(yīng)用中，由于在手術(shù)導(dǎo)航的X射線圖像中小球目標(biāo)與背景對比度很高，為了去除其他噪聲點(diǎn)和物體的干擾，可以取c值為50左右。在提取標(biāo)志點(diǎn)的同時(shí)，抑制背景中物體的干擾。在二值化基礎(chǔ)上，可以提取成簇的亮像素的外輪廓。

1.2.2 圓形目標(biāo)檢測

通過二值化得到圖像中所有包含輪廓的待選像素坐標(biāo)后，對這些候選點(diǎn)依次檢測其輪廓是否形成圓形，檢測的具體方法如下：

(1) 利用1.1.1介紹的曲線擬合方法，求出參數(shù)x0,y0,r，定義擬合度為

(8)

式中：φ為輪廓上所有像素點(diǎn)偏離所擬合曲線的平均距離。容易計(jì)算得到φ10)個(gè)像素的正方形的擬合度約為0.33r2，可見正方形的擬合度隨邊長變大而成平方增大，ε2的值可以根據(jù)此設(shè)置，為了排除半徑較大的非圓形，取值為8～15時(shí)可以取得較好的結(jié)果，ε1的值是為了排除半徑較小的非圓形，取0.5～1時(shí)可以取得較好的結(jié)果。

圖2 彎月形輪廓示意圖Fig.2 Sketch of the curved shape

(3) 為進(jìn)一步提高算法的魯棒性，可利用關(guān)于目標(biāo)點(diǎn)的個(gè)數(shù)以及空間分布的先驗(yàn)知識(shí)，選出所提取的具有最好擬合度和圓度特征的目標(biāo)。

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

為了測試該方法對C形臂X射線圖像中圓形標(biāo)志物的識(shí)別準(zhǔn)確性，對圖3所示的從不同角度拍攝的X射線原始圖像，分別使用文獻(xiàn)[9]提出的檢測圓形輪廓的兩步法和本文的方法進(jìn)行圓的檢測及識(shí)別。原始圖片為有機(jī)玻璃模型在C形臂下的成像，其中有機(jī)玻璃模型上貼了不同半徑的球形鋼珠,該模型為了模擬心電射頻消融手術(shù)的動(dòng)態(tài)跟蹤而構(gòu)建，其中大、小圓筒分別模擬胸腔和心臟，里面插有消融導(dǎo)管頭端。測試所使用的圖像工作站參數(shù)為：DELL OPTIPLEX 380臺(tái)式計(jì)算機(jī)，Intel?CoreMT2 Due CPU，主頻2.93 GHz，內(nèi)存2 GB。使用Visual Studio 2008開發(fā)環(huán)境，利用OpenCV函數(shù)庫，使用C++語言實(shí)現(xiàn)這兩種算法。本文中沒有提供與霍夫圓變換的檢測結(jié)果對比，是因?yàn)槟繕?biāo)圓的半徑較小，利用霍夫變換方法無法成功檢測到圖3中的圓形標(biāo)志物。

圖4給出了對圖3中對應(yīng)兩幅圖二值化的處理結(jié)果。圖5為采用文獻(xiàn)[9]的方法對圖4中的圓目標(biāo)加以提取的結(jié)果，其中十字中心表示檢測到的圓心，小圓表示檢測到的圓形輪廓，可見圖4中的短線狀非圓目標(biāo)也被識(shí)別為圓形目標(biāo)。圖6為采用本文算法對圖4中圓目標(biāo)加以提取的結(jié)果(圖中的十字線和圓標(biāo)示與圖5相同)。從圖6可以看出二值化圖像中圓形的目標(biāo)被正確提取，而非圓的目標(biāo)被排除。從圖5，6算法的對比可以看出：利用本文算法能夠排除文獻(xiàn)[9]中算法對非圓目標(biāo)的誤判，比文獻(xiàn)[9]描述的算法具有更高的魯棒性。圖5中利用文獻(xiàn)[9]的方法誤將長條形(可以看作是半徑較大的彎月形目標(biāo))的目標(biāo)檢測為圓形目標(biāo)，而本文的方法能夠利用圓度的特征將此目標(biāo)排除。為驗(yàn)證實(shí)際應(yīng)用場景下的圓形標(biāo)志物提取效果，采用圖7(a)所示的圖像進(jìn)行測試，在人體胸腔表面粘貼了4個(gè)較大的球形標(biāo)志，表征人體的空間位置，可以看到在存在介入的探測導(dǎo)管圖像中，利用設(shè)置好的參數(shù)，能夠有效的排除場景中小的方形目標(biāo)，所設(shè)計(jì)的算法可以正確地提取感興趣的體表標(biāo)志點(diǎn)。

圖3 測試所用的X射線圖像Fig.3 X-ray images used to test the algorithms

圖4 預(yù)處理后二值化的圖像Fig.4 Binary images after thresholding

圖5 采用文獻(xiàn)[9]算法提取目標(biāo)的結(jié)果Fig.5 Extraction results using the algorithm of Ref. [9]

圖6 采用本文算法提取目標(biāo)的結(jié)果Fig.6 Extraction results using proposed algorithm

圖7 實(shí)際胸腔X射線圖像中圓形標(biāo)志點(diǎn)檢測結(jié)果Fig.7 Circular markers extraction of chest X-ray image

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

實(shí)驗(yàn)結(jié)果表明，本文采取的方法對X射線圖像中圓形標(biāo)記物自動(dòng)提取具有很好的魯棒性，解決了文獻(xiàn)[5，9]中檢測算法對于彎月形物體產(chǎn)生誤判的問題。從時(shí)間復(fù)雜度來看，假設(shè)圖像大小為n×n，二值化后得到的邊緣像素點(diǎn)數(shù)為m，那么采取霍夫圓變換的算法檢測圓的時(shí)間復(fù)雜度為O(mn2)[14]，而本文算法的復(fù)雜度為O(m)，因此本文的算法復(fù)雜度較低。為了測試算法的運(yùn)行效率，隨機(jī)選取6幅圖像，記錄算法所消耗的時(shí)間。實(shí)驗(yàn)結(jié)果表明，本文算法處理一幅圖像平均消耗時(shí)間小于70 ms，因此本文算法能夠?qū)崿F(xiàn)低幀率的實(shí)時(shí)在線檢測?？紤]到算法復(fù)雜度的降低，有可能采用GPU進(jìn)行加速，以達(dá)到更高的目標(biāo)提取速度，這對于在手術(shù)導(dǎo)航系統(tǒng)中實(shí)時(shí)應(yīng)用具有重要的意義。

在本文中對臨床實(shí)際獲得的胸腔X射線圖像序列，是對序列中每張圖像單獨(dú)進(jìn)行提取計(jì)算的，在提高本文算法速度后，有可能對序列運(yùn)動(dòng)圖像中的圓形目標(biāo)進(jìn)行跟蹤，同時(shí)，基于圖像中目標(biāo)運(yùn)動(dòng)的連續(xù)性，可以利用前一幀圖像中提取的目標(biāo)位置，定義下一幀圖像的搜索區(qū)域，從而排除對二值化圖像中所有候選點(diǎn)的重復(fù)檢測，進(jìn)一步提高檢測速度。通過修改本文算法的參數(shù)值，該算法也能很好地提取其他圖像中不同大小的圓形物體，而排除圖像中非圓的物體，圖8給出了兩個(gè)具有不同半徑圓形目標(biāo)，以及混有非圓目標(biāo)圖像的檢測結(jié)果，其中十字中心表示檢測到的圓心，圓表示檢測到的圓形輪廓。從圖8可以看出，本文算法正確地檢出了其中不同大小的圓形目標(biāo)。

圖8 本文算法對不同半徑圓形標(biāo)記和含有非圓標(biāo)記圖像的檢測結(jié)果Fig.8 Circular markers extraction of other test images using proposed algorithm

相對于Hough變換算法，本文所提出算法不能檢測兩個(gè)圓形相互遮擋的情景?？紤]在手術(shù)導(dǎo)航實(shí)際應(yīng)用中，所布置的圓形標(biāo)志點(diǎn)通常是用于系統(tǒng)標(biāo)定、跟蹤和畸變矯正等，因此其空間分布是設(shè)計(jì)好的，在系統(tǒng)設(shè)計(jì)中可以通過將標(biāo)志物置于特定的位置避免不同的標(biāo)志物投影重疊，并利用各個(gè)標(biāo)志點(diǎn)之間在投影圖像中的分布模式提高檢測的魯棒性。為了增強(qiáng)算法自身的魯棒性，也可以利用更復(fù)雜的策略，比如弧形檢測來對相互遮擋的圓形進(jìn)行檢測。

4 結(jié)束語

Hough變換方法用于檢測目標(biāo)圖像中的圓形具有一定局限性，其待檢測的圓形應(yīng)當(dāng)比較接近標(biāo)準(zhǔn)圓，且圓的半徑不能太小。在基于C形臂的手術(shù)導(dǎo)航系統(tǒng)中，圓形標(biāo)志物的投影半徑一般較小，Hough變換難以檢測到目標(biāo)圓形。本文提出的通過輪廓擬合和圓度判別的算法能夠有效地解決上述問題，更穩(wěn)定地檢測圖像中圓形的目標(biāo)，并且算法具有線性的時(shí)間復(fù)雜度，因此在X射線手術(shù)導(dǎo)航系統(tǒng)中具有實(shí)用價(jià)值。

[1] Pascal Fallavollita. Acquiring multiview C-arm images to assist cardiac ablation procedures[J]. EURASIP Journal on Image and Video Processing, 2010, 2010:1-10.

[2] Wang M, Ding H, Wang G. Target visibility enhancement for C-arm cone beam CT-fluoroscopy-guided hepatic needle placement: Implementation and accuracy evaluation [J]. Int J Comput Assist Radiol Surg, 2015(10): 263-273.

[3] 鄭馨，王勇，汪國有. EM 聚類和SVM自動(dòng)學(xué)習(xí)的白細(xì)胞圖像分割算法[J]. 數(shù)據(jù)采集與處理, 2013, 28(5):614-619.

Zheng Xin, Wang Yong， Wang Guoyou. White blood cell segmentation using expectation-maximization and automatic support vector mechine learning [J]. Journal of Data Acquisition and Processing, 2013, 28(5):614-619.

[4] Kim C, Ballard D H, Sklansky J. Finding circles by an array of accumulators[J]. Communications of the Association for Computing Machinery, 1975,18(2): 120-122.

[5] 雷芳.基于改進(jìn)GHT算法的圓形目標(biāo)物的檢測[J]. 電子設(shè)計(jì)工程, 2014, 22(18): 26-28.

Lei Fang. Circle objective detection based on improved GHF[J]. Electronic Design Engineering, 2014, 22(18)： 26-28.

[6] Dori Dov, Liu Wenyin. Stepwise recovery of arc segmentation in complex line environments[J]. International Journal on Document Analysis and Recognition, 1998, 1(1): 62-71.

[7] 王林波,王延杰,邸男,等. 基于幾何特征的圓形標(biāo)志點(diǎn)亞像素中心定位[J]. 液晶與顯示, 2014, 29(6): 1003-1008.

Wang Linbo, Wang Yanjie, Di Nan, et al. Subpixel location of circle target center based on geometric features[J]. Chinese Journal of Liquid Crystals and Displays, 2014, 29(6), 1003-1008.

[8] 吳一全,孟天亮,吳詩婳. 圖像閾值分割方法研究進(jìn)展20年(1994-2014)[J]. 數(shù)據(jù)采集與處理, 2015,30(1):1-23.

Wu Yiquan， Meng Tianliang, Wu Shihua. Research progress of image thresholding methods in recent 20 years(1994-2014)[J]. Journal of Data Acquisition and Processing, 2015,30(1):1-23.

[9] 于中, 徐文立, 陳峰. 邊緣圖像中圓形輪廓的兩步法檢測[J]. 計(jì)算機(jī)輔助設(shè)計(jì)與圖形學(xué)學(xué)報(bào)， 2005,17(2): 243-246.

Yu Zhong, Xu Wenli, Chen Feng. Two-step method for detection of circle profiles in edge images[J]. Journal of Computer Aided Design & Computer Graphics, 2015,17(2): 243-246.

[10]武仲科，焦海星，戴國忠．一種線段和圓弧的逼近方法及其在工程圖紙矢量化中的應(yīng)用[J]．計(jì)算機(jī)輔助設(shè)計(jì)與圖形學(xué)學(xué)報(bào)，1998，10(4)：328-332.

Wu Zhongke, Jiao Haixing, Dai Guozhong. An algorithm of approximating line segment and circular arcs and its application in vectorization of engineering darwings[J]. J CAD & CG, 1998, 10(4): 328-332.

[11]張之穩(wěn), 宋世軍. 利用圓度進(jìn)行圓定位的方法[J]. 山東大學(xué)學(xué)報(bào):工學(xué)版, 2006, 36(1): 19-22.

Zhang Zhiwen, Song Shijun. Circle locating method based on roundness[J]. Journal of Shandong University: Engieeering Science, 2006, 36(1): 19-22.

[12]Nobuyuki O. A threshold selection method from gray-level histograms[J]. IEEE Trans Sys Man Cyber， 1979, 9(1): 62-66.

[13]Sezgin M，Sankur B. Survey over image thresholding techniques and quantitative performance evaluation[J]. Journal of Electronic Imaging， 2004， 13(1)：146-165.

[14]Christopher H. Reduction of computational complexity of hough transforms using a convolution approach[C] ∥ 24th IVCNZ Donald Bailey. Wellington, New Zealand： IEEE Computer Society Press, 2009: 373-378.

Automatic Circular Marker Detection for X-ray Images Base on Geometric Parameter

Wang Guangzhi， Tong Ming， Ding Hui

(School of Medicine, Tsinghua University, Beijing, 100084，China)

Automatic detection of markers for fluoroscopic images is the key for fluoroscopy-based navigation system. Due to the small radius of markers and background noises, the robust and the precision of traditional circle detection algorithms are not satisfied. In order to solve this problem, a geometric-parameter-based algorithm, which can extract markers automatically from X-ray images, is proposed. The method can detect the circles of different radius in fluoroscopic images and non-standard circles by adjusting the parameter. The experimental results show that the proposed algorithm holds high efficiency and robustness.

detection of circles； conic fitting; feature extraction; extraction of markers

國家自然科學(xué)基金(81471759，81127003，51361130032)資助項(xiàng)目。

2015-08-19；

2015-09-15

TP391.41

A

王廣志(1959-)，男，博士，教授，研究方向：醫(yī)學(xué)圖像引導(dǎo)手術(shù)，E-mail: wgz-dea@tsinghua.edu.cn。

童明(1987-)，男，碩士研究生，研究方向：機(jī)器學(xué)習(xí)與醫(yī)學(xué)圖像處理。

丁輝(1967-)，女，碩士，高級(jí)工程師，研究方向：醫(yī)學(xué)圖像處理， E-mail: Dinghui@tsinghua.edu.cn。