X線相干散射圖像環(huán)圓心準(zhǔn)確定位的內(nèi)相關(guān)算法

劉太輝，劉景鑫，李敬筠，宋建中

1.中國科學(xué)院 長春光學(xué)精密機(jī)械與物理研究所，吉林 長春 130013；2.北華大學(xué) 計(jì)算機(jī)學(xué)院，吉林 吉林132031；3.吉林大學(xué)中日聯(lián)誼醫(yī)院信息中心，吉林 長春 130033

X線相干散射圖像環(huán)圓心準(zhǔn)確定位的內(nèi)相關(guān)算法

劉太輝1，2，劉景鑫3，李敬筠3，宋建中1

1.中國科學(xué)院 長春光學(xué)精密機(jī)械與物理研究所，吉林 長春 130013；2.北華大學(xué) 計(jì)算機(jī)學(xué)院，吉林 吉林132031；3.吉林大學(xué)中日聯(lián)誼醫(yī)院信息中心，吉林 長春 130033

X射線相干散射成像分析中，能量的分布為一系列同心圓環(huán)。準(zhǔn)確地確定圓心是精確提取能量譜線的關(guān)鍵。本文采用K均值聚類結(jié)合形態(tài)學(xué)濾波確定內(nèi)環(huán)區(qū)域，使用形心法近似確定圓心，對圓心使用標(biāo)準(zhǔn)差內(nèi)相關(guān)運(yùn)算逼近標(biāo)準(zhǔn)圓心。實(shí)驗(yàn)結(jié)果證明，該算法實(shí)用、有效。

X線；相干散射；能量分布；K均值聚類；形態(tài)學(xué)；內(nèi)相關(guān)

專欄——醫(yī)學(xué)成像新技術(shù)的研究與應(yīng)用

編者按：醫(yī)學(xué)影像學(xué)科是醫(yī)院發(fā)展及醫(yī)院收入的重要支柱，與臨床醫(yī)學(xué)工程學(xué)科有著極其密切的聯(lián)系。醫(yī)學(xué)影像設(shè)備的良好應(yīng)用需要臨床工程師與醫(yī)療、技術(shù)和護(hù)理等多方面人員的密切合作。加強(qiáng)醫(yī)學(xué)影像的質(zhì)量控制，提高輻射防護(hù)手段和水平，開發(fā)現(xiàn)有醫(yī)學(xué)影像設(shè)備的新功能，開展新設(shè)備、新技術(shù)的應(yīng)用和推廣，進(jìn)而提高全院醫(yī)療設(shè)備的應(yīng)用水平、利用率特別是設(shè)備功能利用率，是我們臨床醫(yī)學(xué)工程人員的重要工作內(nèi)容。

欄目主編：劉景鑫

吉林大學(xué)中日聯(lián)誼醫(yī)院信息中心主任，吉林省醫(yī)學(xué)影像工程技術(shù)研究中心主任，吉林省醫(yī)學(xué)影像質(zhì)量控制中心副主任，中華醫(yī)學(xué)會醫(yī)學(xué)工程學(xué)分會常委，《中國醫(yī)療設(shè)備》常務(wù)編委。教授，碩士生導(dǎo)師，主持國家自然科學(xué)基金及省市等多項(xiàng)科研課題。

郵箱： liujingxin@126.com

醫(yī)學(xué)X射線相干散射技術(shù)的應(yīng)用研究是近些年在國際上興起的一種新的醫(yī)用病理和物質(zhì)成分分析技術(shù)，可以應(yīng)用于病理分析、疾病的早期診斷和輔助治療等方面。目前，國際上在醫(yī)學(xué)相干散射研究上進(jìn)行了多種形式的研究工作。Puumalainen 等人利用求X射線的相干散射與Compton散射的比率的方式，分析在樣品中骨骼礦物質(zhì)成份的平均電荷密度起伏，用來確定骨骼礦物質(zhì)密度[1]。骨礦含量（BMC）的測定采用的是單純的測量相干散射截面與空間角之間的關(guān)系[2]。Westmore等人也進(jìn)行了相干散射成像及物質(zhì)成份之間關(guān)系的研究，通過對不同的單一物質(zhì)進(jìn)行成像，可以獲得不同空間分布的對稱圓環(huán)圖像[3]。

以臨床影像設(shè)備成像的方式進(jìn)行相干散射研究，是臨床研究的較好途徑。能量分布曲線（能譜）可以準(zhǔn)確地反映散射線的空間特征，因此，準(zhǔn)確地提取圖像的能譜曲線，是相干散射成像技術(shù)研究的一個關(guān)鍵環(huán)節(jié)。相干散射圖像具有噪聲影響嚴(yán)重、圖像分布彌散、沒有確定的邊界、內(nèi)部信息互相混疊不易提取的問題。本文綜合多種圖像處理技術(shù)，提出自相關(guān)優(yōu)化方法，準(zhǔn)確確定圖像環(huán)的圓心，為提取圖像能量譜提供支持。

1 X射線相干散射圖像特點(diǎn)

如圖1所示，在樣品圖中，目標(biāo)是黑色彌散狀的圓環(huán)分布，其沒有確定的邊界。內(nèi)部的不規(guī)則白色區(qū)是射線阻擋器的偽影，形狀不規(guī)則，與圓環(huán)不共心。在外標(biāo)定環(huán)（非圖像內(nèi)容）外，還有一些彌散的暗區(qū)，分布不具有對稱性。同時，通過直方圖可見，灰度的分布不具有明顯的分類特征，通過直方圖進(jìn)行分割不可行。因此，目標(biāo)圓環(huán)的圓心只能夠通過內(nèi)標(biāo)定環(huán)以內(nèi)的區(qū)域進(jìn)行標(biāo)定。

圖 1 散射圖像及其直方圖

2 近似圓心算法

2.1 K均值聚類結(jié)合形態(tài)學(xué)濾波的圖像二值化

K均值聚類的基本含義，對于給定的一組數(shù)據(jù)（x1,x2,x3……xn）,如果每一個數(shù)據(jù)都是一個d維的實(shí)向量，則K均值聚類就是把n個數(shù)據(jù)分割為k組，（k

這里ui是集合Si的均值。這里以圖像的灰度值作為分類的依據(jù)，即

在此，K均值聚類的目標(biāo)是有效地劃分目標(biāo)區(qū)和背景區(qū)。對于散射圖像而言，由于噪聲的強(qiáng)度與有效值的強(qiáng)度非常接近，而且其分布的灰度范圍很窄，所以聚類數(shù)K不宜選擇過大， 3≤K≤10。重置聚類集內(nèi)像素點(diǎn)的灰度值，實(shí)現(xiàn)二值化：

經(jīng)過聚類分割后的圖像，背景中散落著噪聲，圓環(huán)的邊界模糊不清，目標(biāo)圓環(huán)的內(nèi)部也散落著白色的背景點(diǎn)。仔細(xì)分析背景和目標(biāo)圓環(huán)的圖像特征后，背景中的噪聲點(diǎn)的分布雜散，沒有區(qū)域性，而圓環(huán)內(nèi)的背景散落點(diǎn)稀少，不影響圓環(huán)的區(qū)域特征。這種分布特點(diǎn)恰好可以使用形態(tài)學(xué)濾波方法進(jìn)行進(jìn)一步的平滑處理。

選用3×3濾波器[6]

首先進(jìn)行一次膨脹運(yùn)算，加強(qiáng)圓環(huán)的區(qū)域性，然后進(jìn)行連續(xù)多次的腐蝕操作，剪除背景中的噪聲點(diǎn)。

2.2 正交能量投影及形心計(jì)算

二值化的圖像，需要準(zhǔn)確判別環(huán)內(nèi)區(qū)域的分布和大小，采用一般的梯度邊緣提取算法繁瑣而且結(jié)果不準(zhǔn)確。這里提出基于正交能量投影的算法。

首先，沿著圖像的x軸方向和y軸方向，分別作圖像的橫向和縱向的投影。投影的結(jié)果為該方向上單一行（列）內(nèi)像素點(diǎn)能量的累積。由于圖像本身是一個正交的平面，因此投影后的能量累積在特定坐標(biāo)方向上滿足圖像分辨率的要求，而且任意的兩個x、y方向上獨(dú)立累積能量所在的坐標(biāo)位置恰好唯一標(biāo)定一個圖像空間的像素點(diǎn)。

設(shè)圖像為P={p1,p2,……pmxn}(m,n)∈N，m、n分別為坐標(biāo)x、y軸方向的圖像分辨率。每一個數(shù)據(jù)pi是二維空間中的一個標(biāo)量p(m,n)。Em={e1,e2,……en}n∈N, Em為沿著m方向上能量投影集合，同理，En={e1,e2,……em}m∈N 為n方向上的能量投影集合。其中：

投影后的曲線中有明顯的波谷存在，波谷是該投影方向上圓環(huán)的能量導(dǎo)致的，也就是在兩個波谷之間的區(qū)域是圓環(huán)的內(nèi)部區(qū)域。求得兩條譜線上波谷之間區(qū)域的中點(diǎn)位置（xc，yc），其對應(yīng)的圖像像素必在圓環(huán)內(nèi)部區(qū)域上。以這個點(diǎn)為起始點(diǎn)進(jìn)行區(qū)域增長處理，可以獲得圓環(huán)內(nèi)部區(qū)域的準(zhǔn)確分布，進(jìn)而對這個區(qū)域使用形心算法，提取其形心。這個形心就是圓環(huán)的近似圓心。

圖 2 求得的形心

3 圓心準(zhǔn)確定位算法

相干散射圖像在理論上成像圓環(huán)完全共心，因此能否準(zhǔn)確地確定圓心的位置，是進(jìn)一步分析的關(guān)鍵。

圖像中，根據(jù)前面的分析可以知道，除在過圓心的水平軸(x)和垂直軸(y)上的像素點(diǎn)可以被明確地確定為某一個圓環(huán)上的點(diǎn)以外，其他的像素點(diǎn)都是通過歐幾里德距離計(jì)算后，近似歸類為某一個圓環(huán)。由于在成像探測器對空間中x射線的分布進(jìn)行成像的量化過程中，成像的結(jié)果受空間中的噪聲、成像探測器量化誤差和采樣算法的計(jì)算誤差的綜合影響，只有過圓心的坐標(biāo)軸上的點(diǎn)相對準(zhǔn)確地反映圖像對x射線能量分布。

假定當(dāng)前的圓心為理想圓心，那么當(dāng)前圓心對應(yīng)的圓環(huán)上的像素分布一定具有強(qiáng)相關(guān)性，這種相關(guān)是環(huán)內(nèi)像素間的相關(guān)，稱作內(nèi)相關(guān)性，即像素點(diǎn)的灰度分布在一個相對比較均勻的范圍內(nèi)。而當(dāng)圓心離開理想圓心位置的時候，由于空間射線分布的高斯特性，新的圓環(huán)上的像素點(diǎn)之間的相關(guān)性將減低。偏心越嚴(yán)重減低的程度越大。如圖3所示，假設(shè)理想圓心的位置為xc、yc。圓心向下偏移（xc+3、yc-2）時，所有當(dāng)前圓心對應(yīng)的圓環(huán)上的點(diǎn)全部改變，相關(guān)性會大大降低。

上面的形心方法獲得的圓心是一個偏離理想圓心的結(jié)果，通過上面相關(guān)性分析，可以通過相關(guān)性計(jì)算使圓心逐漸逼近理想圓心。

圖 3 圖像中的圓形分布示意圖

3.1 相關(guān)性測度

樣品點(diǎn)是可以準(zhǔn)確代表能量分布的像素點(diǎn)。通過前面的分析，在圖像中圓環(huán)的選取只能夠以像素單位進(jìn)行。無疑以理想圓心為原點(diǎn)的笛卡爾坐標(biāo)系是分析的標(biāo)準(zhǔn)坐標(biāo)系，那么x軸和y軸上的像素點(diǎn)就是圓環(huán)的標(biāo)稱點(diǎn)，這些點(diǎn)也最適合作為樣品點(diǎn)。其他的圓環(huán)上的點(diǎn)都具有不同的近似性，不適合作為樣品點(diǎn)使用 。

同一個圓環(huán)上的樣品點(diǎn)的灰度值分布，理想的情況下，應(yīng)該在一個標(biāo)準(zhǔn)值（平均值）附近擺動，相互間灰度差值很小。這種特性恰好符合標(biāo)準(zhǔn)偏差的理論特征。

那么樣品點(diǎn)的標(biāo)準(zhǔn)偏差就可以作為同心圓環(huán)上樣品點(diǎn)之間灰度分布均勻性的測度。

計(jì)算過程中，樣品點(diǎn)選擇為每一個圓環(huán)上以圓心為原點(diǎn)的笛卡爾坐標(biāo)系中，x軸上圓心對稱的兩個點(diǎn)，y軸上圓心對稱的兩個點(diǎn)，總計(jì)4個像素點(diǎn)。樣品圓環(huán)選定為進(jìn)入有效能量分布區(qū)中距離圓心30～60像素距離的一組圓環(huán)。

每一個樣品圓環(huán)上的樣品點(diǎn)之間（圓心為理論圓心）具有強(qiáng)相關(guān)性，它們的分布標(biāo)準(zhǔn)偏差應(yīng)該最小（分布在均值附近的很小的范圍內(nèi)）。由于單一的圓環(huán)的測度不能夠規(guī)避掉奇異點(diǎn)的影響，因此可以選擇多個圓環(huán)作為樣品圓環(huán)，這里選取為40個樣品環(huán)。在每一個樣品環(huán)都計(jì)算出標(biāo)準(zhǔn)偏差后，對組標(biāo)準(zhǔn)偏差求均值，作為當(dāng)前的圓心偏離理想圓心的測度。如果計(jì)算過程中，出現(xiàn)2個以上的圓心測度值相等的情況，則以這些樣品點(diǎn)作為新的樣品圓心，增加或減小樣品環(huán)的數(shù)量，重新計(jì)算每個樣品圓心的測度，直至沒有相同取值存在為止。

3.2 標(biāo)準(zhǔn)偏差

設(shè)當(dāng)前圓心對應(yīng)的Gi={g1,g2,……,gn}n∈N是第i個圓環(huán)上的樣品點(diǎn)的集合，那么：樣品點(diǎn)的均值為

此組樣品點(diǎn)的標(biāo)準(zhǔn)差為：

那么樣品圓環(huán)組的標(biāo)準(zhǔn)偏差集為{δ1,δ2,δ3,……,δi} i ∈N，則當(dāng)前圓心的標(biāo)準(zhǔn)差測度為：

3.3 計(jì)算方法

標(biāo)準(zhǔn)的圓心必然在圖像的內(nèi)環(huán)以內(nèi)，根據(jù)文獻(xiàn)資料[7，8]X射線前向小角度相干散射的成像角度在5～10°之間。因此，以成像距離（源到探測器）40cm計(jì)算，探測器尺寸（43cm×35cm）高x寬，成像板分辨率（2688×2208）像素。每一個像素點(diǎn)的空間角為0.0219°，那么有效的區(qū)域在200個像素點(diǎn)以外。極限情況下，形心算法算得的圓心位于距理想圓心100個像素點(diǎn)的位置上的時候，逼近算法可以準(zhǔn)確地搜索到理想圓心，這時候搜索范圍應(yīng)該為100個像素點(diǎn)。而近似圓心如果距離理想圓心的距離超過100個像素點(diǎn)，那么搜索有可能進(jìn)入能量分布區(qū)，相關(guān)性測度的計(jì)算將無法控制，這個圓心結(jié)果應(yīng)歸為壞圓心，需要重新計(jì)算。

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

如圖4所示，為典型的四個樣品的相關(guān)計(jì)算結(jié)果。在對15個樣品圖片進(jìn)行算法的實(shí)驗(yàn)計(jì)算后，可以得到基本的結(jié)論：① 形心算法獲得的形心不同程度地偏離理想圓心。② 相關(guān)強(qiáng)度隨著估計(jì)圓心偏離理想圓心而減低。③ 理想圓心具有全局域中最小的標(biāo)準(zhǔn)偏差，具有最大的相關(guān)強(qiáng)度。這說明上述的相關(guān)算法收斂而且有效。

圖 4 內(nèi)相關(guān)強(qiáng)度曲線及求精結(jié)果

由于沒有準(zhǔn)確的邊界，也沒有評價標(biāo)準(zhǔn)的支持，這里只是簡單地以目視可以區(qū)分的方式顯示圓心準(zhǔn)確定位的結(jié)果，進(jìn)一步的理論支持將在以后的科研工作中，通過穩(wěn)定的算法分析進(jìn)行驗(yàn)證。

圖 5 算法仿真比較圖像

如圖5所示，圖像中b為原始圖像，a為經(jīng)過分析后的放大圖像。為了對比方便，在圖像中，有兩對十字交叉線，前面的一對展示的是精確定位后的圓心，對應(yīng)白色的圓形區(qū)域，該圓形是通過圓形算法和圖像區(qū)域增長算法結(jié)合計(jì)算的仿真結(jié)果。后面被遮擋的十字交叉線是形心算法的估計(jì)圓心，對應(yīng)的圓形區(qū)域?yàn)榛疑浯蟛糠直徽趽踝?。可以通過兩個區(qū)域的疊加看出兩個圓心之間的差異，也可以通過十字交叉線展示兩種算法的結(jié)果的區(qū)別。通過目視可以明顯區(qū)別出精確計(jì)算獲得的圓心相比較情況下，更加接近于理想圓心（如圖6）。

圖 6 仿真結(jié)果三組對比圖

5 結(jié)論

通過算法結(jié)果的分析，本論文中的內(nèi)相關(guān)精確圓心計(jì)算方法是有效的，不同的仿真結(jié)果也證明了精確定位后的圓心相對于近似估計(jì)圓心更加接近于理想圓心。這對后期進(jìn)行相干散射圖像中提取和分析相干散射X射線的能量譜線來說是非常關(guān)鍵的。進(jìn)一步的工作應(yīng)該對噪聲和奇異點(diǎn)對圓環(huán)分布的影響進(jìn)行分析總結(jié)，并完善內(nèi)相關(guān)算法，得到穩(wěn)定、優(yōu)化的求精算法。

[1] P.Puumalainen, A.Uimarihuhta, E. Alhava,et al.A new photon scattering method for bone mineral density measurements[J]. Radiology,1976,120:723-724.

[2] M.Newton,D.W.L.Hukins,G.Harding.Bone composition measured by x-ray scattering[J].Phys.Med.Biol,1992,37: 1339-1347.

[3] M.S.Westmore,A.Fenster,I.A.Cunningham.Angulardependent coherent scatter measured with a diagnostic x-ray image intensifier-based imaging system[J].Med.Phys,1996, 23:723-733.

[4] Inaba M.,Katoh N.,Imai H.Applications of weighted Voronoi diagrams and randomiza-tion to variance-based k-clustering[C].Proceedings of 10th ACM Symposium on Computational Geometry,1994:332-339.

[5] Kanungo T,Mount D.M., Netanyahu N.S.,et al.An efficient k-means clustering algorithm:Analysis and implementation[J]. IEEE Trans. Pattern Analysis and Machine Intelligence,2002,24: 881-892.

[6] Pierre Soille.Morphological Image Analysis: Principles and Applications[M].2 ed. Santa Clara: Springer-Verlag Telos,2003.

[7] MacKay, David. Information Theory, Inference and Learning Algorithms[M]. Cambridge:Cambridge University Press,2003.

[8] Liu taihui, Sun junxi, Liu jingxin.Medical X-ray Coherent Scatter Images Processing based on the MRF Model[J]. Proc of SPIE,2006,6026:1-9.

Inner Correlation Algorithm on Determining Ring Center of X-Ray Coherent Scatter Image

LIU Tai-hui1,2,LIU Jing-xin3, LI Jingjun3, SONG Jian-zhong1
1.Changchun Institute of Optics, Fine Mechanics and Physics,Chinese Academy of Sciences, Changchun Jilin 130013, China;2.Beihua University, Computer College,Jilin Jilin 132031, China;3. Information Center,China-Japan Union Hospital of Jilin University, Changchun Jilin 130033, China

R814

B

10.3969/j.issn.1674-1633.2010.07.001

1674-1633(2010)07-0001-04

2010-06-21

國家自然科學(xué)基金(60372062)。

本文作者：劉太輝，副教授。

劉景鑫，教授,碩士生導(dǎo)師。

作者郵箱： luke_ciom@yahoo.com.cn

Abstract:The energy distribution of X-ray coherent scatter images is a series of concentric rings. It is the key step to determine the center of rings accurately. In this paper, the inside region of image was detected by K-mean and mathematical morphology methods. The center was determined with the centroid algorithm. Then the center was recalculated with the standard error inner correlation algorithm to approach the theory circle center. The experiment result shows that the inner correlation algorithm is effective and practical.

Key words:X-ray; coherent scatter; energy distribution; K-Mean cluster; mathematical morphology; inner correlation