基于高斯混合模型和核密度估計的全身骨骼SPECT圖像分割算法研究

徐磊，孟慶樂，楊瑞，田書暢，蔣紅兵

南京醫科大學附屬南京醫院（南京市第一醫院） a.醫療設備處；b.核醫學科，江蘇 南京 210006

基于高斯混合模型和核密度估計的全身骨骼SPECT圖像分割算法研究

徐磊a，孟慶樂b，楊瑞b，田書暢a，蔣紅兵a

南京醫科大學附屬南京醫院（南京市第一醫院） a.醫療設備處；b.核醫學科，江蘇 南京 210006

目的提出一種基于高斯混合模型的骨掃描圖像分割算法，可自動識別全身骨骼SPECT圖像中的病變區域。方法首先對二維全身骨骼SPECT圖像進行銳化、平滑、灰度變換等預處理；然后采用核密度估計方法擬合出預處理圖像的像素概率密度函數曲線，并根據曲線的峰值點確定期望最大值（EM）算法的初始值；再應用高斯混合模型對圖像進行分割；最后使用模板匹配算法排除誤識別的區域。結果應用本研究提出的圖像分割算法所得到的圖像清晰度和對比度優于其他圖像分割算法，且本研究提出的圖像分割算法的 相似性測度明顯高于其他圖像分割算法，平均耗時最短。結論基于高斯混合模型和核密度估計的全身骨骼SPECT圖像分割算法是一種高效、實用的骨掃描圖像分割算法。

全身骨骼顯像；高斯混合模型；核密度估計；EM算法

全身骨骼顯像是臨床上一種重要的檢查項目，因靈敏度高、能夠全身成像、成本低等優點廣泛應用于病灶檢測和腫瘤骨轉移識別等方面[1-3]。正常人的全身骨骼顯像放射性分布對稱，圖像均勻清晰；而患者的腫瘤骨轉移或骨骼外傷部位則會出現放射性核素濃聚現象，在圖像上表現為團點狀的亮點[4]。目前全身骨骼顯像診斷主要依靠醫生的個人經驗，主觀因素影響較大，容易出現誤診、漏診的情況。因此開發針對全身骨骼顯像的計算機輔助診斷系統具有重要意義，該系統可運用圖像分割算法提取出骨掃描圖像上的感興趣區域，便于對病灶進行客觀的定性、定量分析。

圖像分割是利用圖像特征集合的相似性準則對圖像像素進行聚類分組，將圖像劃分成不相交區域的方法[5]，而骨掃描圖像分割結果的優劣直接影響后續病灶的識別與定量分析。目前針對骨掃描圖像分割的研究工作主要在國外開展：Sadik等[6]采用的基于標準圖譜的圖像分割算法具有較高的系統敏感性和特異性；Huang等[7]提出的基于模糊集的圖像閾值化的分割算法主要用于圖像分區；Sajn等[8]提出的基于先驗知識的特征點檢測算法的敏感性和特異性易受病人差異性影響。相比之下，國內相關研究較少，朱春媚等[9]嘗試采用BP神經網絡算法進行圖像分割，但該算法存在收斂速度較慢、網絡層次較多等缺點。由此可知，目前針對全身骨骼顯像的圖像分割尚沒有標準的分割算法。

本研究聯合使用基于期望最大值（Expectation Maximization，EM）算法的高斯混合模型（Gaussian Mixture Model，GMM）和核密度估計算法對全身骨骼顯像的圖像進行分割。其中高斯混合模型算法簡單、高效，可以處理大像素的圖像，且對團點狀區域的圖像聚類效果顯著；EM算法能估算出高斯混合模型中的參數值；核密度估計可以根據圖像灰度的概率密度曲線上的峰值點確定EM所需的初始值。

1 資料與方法

1.1 資料

本研究選取南京市第一醫院核醫學科2015年1月～6月間收治的200例患者的全身骨骼顯像的圖像作為實驗樣本，其中無腫瘤骨轉移和發生腫瘤骨轉移的患者圖像各100例。每例患者均在注射99Tcm之后2～3小時采集前身和后身兩張SPECT骨掃描圖像，圖像采集矩陣大小為256×1024，灰度級為256，所用算法均在Matlab2013a開發環境下仿真實現。

1.2 方法

1.2.1 圖像分割流程

全身骨骼顯像的圖像信噪比低、分辨率差、差異性大，采用單純高斯混合模型對其進行圖像分割并不能達到理想的分割效果。本研究將圖像的處理流程分為預處理、高斯混合模型進行圖像分割和后處理3部分（圖1）。其中，預處理流程主要負責銳化圖像、去除噪聲、擴大灰度動態范圍；高斯混合模型則負責將圖像劃分成若干區域，并提取出感興趣區域；后處理流程主要通過對稱程序和模板匹配法去除誤識別的區域。

圖1 全身骨骼顯像的圖像處理流程

1.2.2 圖像預處理

受成像原理限制，骨掃描圖像的灰度范圍很窄并伴隨很高的噪聲。為了盡可能保留圖像的骨骼信息，突出細節信息，需要通過圖像預處理來改善圖像質量。本研究采取的預處理策略如下：首先采用拉普拉斯變換突出圖像的細節部分；然后采用梯度法突出其邊緣，經均值濾波器平滑后的梯度圖像掩蔽拉普拉斯圖像；最后通過灰度變換擴展圖像的灰度動態范圍（圖2）。

圖2 骨掃描圖像的增強處理流程

1.2.3 高斯混合模型

高斯混合模型是一種經典的無監督聚類模型，能對大型數據進行高效聚類，廣泛用于圖像分割、模式識別、數據挖掘等領域[10]。高斯混合模型的基本思想[11]如下：用高斯概率密度函數精確量化事物，將一個事物分解為若干基于高斯概率密度函數的模型。高斯混合模型應用于圖像分割的原理可以理解為：假設圖像被分割成K個區域，每個區域的像素服從均值為μ、方差為σ2的正態分布，則整個圖像特征分布可以用高斯混合模型來描述。

第K類的單高斯分布如式(1)所示，由此可知，整個圖像的混合高斯分布可由式(2)表示。

其中，Ixy是圖像第x行、第y列像素的灰度值；πk是Ixy屬于第K類的概率權重，πk滿足式(3)；fK(Ixy；θK)是Ixy在第K類時的概率密度函數，θ=（μ1，μ2，...，μK；σ21，σ2

2，...，σ2

K）；所有參數形成的未知矢量Φ=(π1，π2，...，πK；θ1，θ2，...，θK)。

1.2.4 基于EM算法求解高斯混合模型參數

在求解高斯混合模型參數時，直接采用最大似然估計高斯混合模型的參數比較困難，本研究選用的EM算法是一種快捷、簡便的最大似然估計方法，可以優化聚類過程。

首先，構造Ixy的對數似然估計函數，表達形式如式(4)所示。

1.2.5 基于核密度估計確定EM算法初值

EM算法簡單、穩定，但對初值點比較敏感，隨機選取初值點容易導致其陷入局部最優解[13]。本研究選用核密度估計確定EM算法初始值，該方法可直接從圖像的像素連續變化值中擬合出概率密度函數曲線[14]，并根據曲線峰值點確定初始值。一般選用高斯函數作為核函數，需要等概率選出 個格點的像素灰度值作為觀察值。高斯核密度估計函數關系式如式(9)[15]所示。

在圖像中等概率選出m個格點可由式(10)完成，此處m取256，窗寬h由式(11)確定[16]。

其中，yj是第j個格點的像素灰度值，Max(Ixy)、Min(Ixy)分別是圖像像素的灰度最大值和最小值，h為窗寬，STD是圖像像素的標準差，IQR為圖像像素的四分位差。

1.2.6 圖像后處理

在嘴巴、鼻子、膀胱等部位常出現核素聚集現象，特征與病變區域相似，利用模板匹配法和對稱程序可以檢測出這些誤識別的區域，然后對這些區域的灰度值進行翻轉，使之變為正常區域。

1.2.7 圖像分割結果評價

采用主觀視覺對圖像預處理的效果進行定性評價，選用Tanimoto相似性測度和平均所耗時間對高斯混合模型的圖像分割效果進行定量評估，0≤Tanimoto≤1，Tanimoto相似性測度越接近于1，平均耗時越短，說明圖像分割結果越精確，算法收斂越快[17]。Tanimoto相似性測度可用式(12)表示。

其中，XK是分割后第k類區域的像素數量，YK是源圖像中第k類目標區域的像素數量，

圖3 骨掃描圖像預處理結果

2 結果

2.1 骨掃描圖像的分割步驟

綜上所述，本研究確定的骨掃描圖像的分割步驟如下所示。

（1）輸入骨掃描圖像X=Ixy，進行圖像預處理。

（2）使用式(10)等概率選出256個格點，使用式(11)計算窗寬h，然后使用式(9)擬合出骨掃描圖像像素的核密度估計曲線。

（3）根據核密度估計曲線的峰值點確定EM算法所需的初始值，峰值點hpi即是初始值μi，σ即是所需的σi，σ的計算公式如式(13)所示。

（4）利用EM算法迭代計算最大似然估計參數。設置初始值Фold，利用公式(5)～(8)更新參數值，如果logL(Ф(new)) -log L(Ф(old))＜ε，迭代停止，否則繼續使用公式(5)～(8)迭代。ε是臨界值，本研究選取ε=0.01。

（5）輸出圖像分割結果，并對其進行模板匹配處理，得到最終骨掃描圖像分割結果。

2.2 圖像預處理結果

圖像預處理包括銳化、平滑、灰度變換等，具體結果，見圖3。從視覺效果看，拉普拉斯突出了圖像的細節部分，梯度變換提取出了圖像的感興趣區域，均值濾波器提高了圖像信噪比，灰度變換提高了圖像對比度。

2.3 基于核密度估計確定EM算法初值結果

本研究選用高斯核密度估計和直方圖來確定EM算法所需的初始值見圖4。相比于圖像直方圖，高斯核密度估計曲線具有更明顯的峰值點，圖中帶箭頭的峰值點可選為初始值。為了兼顧計算量以及分割精度，EM初始值選取規則如下[14]：當峰值點個數＜10時，所有峰值點均作為初始值；當峰值點＞10時，將峰值點由大到小排列，取前10個峰值點作為初始值。由此可得初始值目K，選為10，具體值為μ1=1.6，μ2=12.1，μ3=17.3，μ4=25.1，μ5=30.1，μ6=38.1，μ7=45.9，μ8=51.5，μ9=56.2，μ10=64.0。

圖4 骨掃描圖像的直方圖和高斯核密度估計曲線

2.4 不同算法的圖像分割效果比較

本研究選取不同算法對骨掃描圖像進行了分割，結果見圖5。由圖可知，采用本研究提出的算法分割的圖像清晰度、對比度明顯高于其他算法。

對骨掃描圖像的不同分割算法Tanimoto的定量評價結果見表1，可以看出，神經網絡算法的 相似性測度最小，且平均耗時最長；K均值聚類算法Tanimoto可以在一定程度上提升分割精度；基于高斯混合模型算法的相似性測度相對較大，平均耗時最短；其他三種算法雖然分割精度較高，但是都存在耗時較長的缺陷。

表1 骨掃描圖像的不同分割算法的定量評價結果

3 討論

本研究提出了基于高斯混合模型的全身骨骼顯像的圖像分割算法。整個圖像分割流程簡潔明了，其中圖像預處理突出了圖像的骨骼部分和感興趣區域，提高了圖像的對比度和清晰度；高斯核密度估計根據峰值點確定了EM算法所需要的初始值，避免了隨機選取初始值的盲目性；EM算法更準確地估計了高斯混合模型的參數，達到了快速分割圖像的效果；圖像后處理排除了誤識別區域，使分割結果更為準確。

綜上所述，本研究提出的全身骨骼顯像的圖像分割算法能快速、準確、穩定地提取出感興趣區域，識別出微小可疑病變，便于醫生診斷；其Tanimoto相似性測度和平均耗時都最優，說明該算法是一種可行的、較優的骨掃描圖像分割算法。在以后的工作中，本研究將聯合使用Matlab和C++設計并開發出一套骨掃描圖像分割輔助診斷系統。

圖5 不同算法的骨掃描圖像分割結果

[1] 王慶文,左麗,朱乙聲,等.單光子發射計算機斷層顯像全身骨掃描在脊柱關節病診斷中的價值[J].中華風濕病學雜志,2014, 18(12):815-818.

[2] 魯珊珊,王德杭,李永軍,等.全身MRI與核素骨掃描對骨轉移瘤診斷價值的對照研究[J].中華放射學雜志,2011,45(5):459462.

[3] ?ajn L,Kononenko I,Mil?inski M.Computerized segmentation and diagnostics of whole-body bone scintigrams[J].Comput Med Imag Grap,2007,31(7):531-541.

[4] 張鵬,王國申,孔維慧,等.骨肉瘤99Tcm-MDP全身骨掃描的特點及定量分析[J].中國誤診學雜志,2012,12(4):859-860.

[5] 李杰,陳國棟.基于改進區域生長算法的肝臟管道圖像分割方法[J].中國醫療設備,2014,29(10):19-23.

[6] Sadik M,Hamadeh I,Nordblom P,et al.Computer-assisted interpretation of planar whole-body bone scans[J].J Nucl Med,2008,49(12):1958-1965.

[7] Huang JY,Kao PF,Chen YS.A set of image processing algorithms for computer-aided diagnosis in nuclear medicine whole body bone scan images[J].IEEE T Nucl Sci,2007, 54(3):514-522.

[8] ?ajn L,Kukar M,Kononenko I,et al.Computerized segmentation of whole-body bone scintigrams and its use in automated diagnostics[J].Comput Meth Prog Bio,2005,80(1):47-55.

[9] 朱春媚,田聯房,陳萍,等.基于BP神經網絡的全身骨SPECT圖像分割[J].生物醫學工程學雜志,2007,24(5):1050-1053.

[10] 陳允杰,張建偉,韋志輝,等.基于高斯混合模型的活動輪廓模型腦MRI分割[J].計算機研究與發展,2015,44(9):1595-1603.

[11] 朱峰,羅立民,宋余慶,等.基于自適應空間鄰域信息高斯混合模型的圖像分割[J].計算機研究與發展,2015,48(11):2000-2007.

[12] Portela NM,Cavalcanti GDC,Ren TI.Semi-supervised clustering for MR brain image segmentation[J].Expert Syst Appl,2014,41(4):1492-1497.

[13] 孫殿柱,劉華東,史陽,等.基于核密度估計的散亂點云邊界特征提取[J].農業機械學報,2013,44(12):275-279,268.

[14] 宋余慶,謝從華,朱玉全,等.基于近似密度函數的醫學圖像聚類分析研究[J].計算機研究與發展,2015,43(11):1947-1952.

[15] Chen TB,Chen JC,Lu HS.Segmentation of 3D microPET images of the rat brain via the hybrid gaussian mixture method with kernel density estimation[J].J X-ray Sci Technol,2012,20(3):339-349.

[16] 李毅,王遠弟.基于核密度估計的圖像平滑的最優停止[J].上海大學學報(自然科學版),2011,17(1):103-110.

[17] 施明登,姚江河,白鐵成.Gabor 傳遞函數相似性測度圖像質量評價算法[J].科技通報,2014,30(10):133-135.

A Research on the Segmentation Algorithm for the Whole Body SPECT Image via the Gaussian Mixture Model with Kernel Density Estimation

XU Leia, MENG Qing-leb, YANG Ruib, TIAN Shu-changa, JIANG Hong-binga
a.Department of Medical Equipment; b.Department of Nuclear Medicine, Nanjing First Hospital, Nanjing Medical University, Nanjing Jiangsu 210006, China

ObjectiveTo propose a novel segmentation algorithm for the whole-body bone scan image based on the Gaussian mixture model (GMM) which is used for the automatic recognition of the lesion area.MethodsFirst, we sharpened and smoothed the 2D SPECT whole-body scan image for preprocessing. Second, Gaussian kernel density estimation was adapted to obtain the initial value of the expectation-maximization (EM) algorithm by fitting the curve of probability density function. Then we segmented the image using the GMM algorithm. Finally, the template match method was used to eliminate the wrong recognized areas.ResultsFrom subjective evaluation, the presented segmentation method can provide clearer and more detailed activity structures and improve the image quality. Quantitatively experimental results indicate that the GMM algorithm can generate a higher degree of Tanimoto similarity than other methods, and has a less running time.ConclusionKernel density estimation can effectively prevent the blindness of the initial value selection in the EM algorithm. Thus the lesion areas will be segmented accurately by combination of the GMM and EM method. Therefore, the proposed method is a feasible algorithm for the whole-body bone scan image segmentation.

TP391.41

A

10.3969/j.issn.1674-1633.2016.02.011

1674-1633(2016)02-0048-04

2015-09-21

南京市醫學科技發展資金“青年工程”人才培養專項經費資助項目（QRX11033）。

蔣紅兵，研究員級高級工程師，主要從事生物醫學工程研究。

通訊作者郵箱： cmdjhb@126.com

Abstract:: whole-body bone scan; Gaussian mixture model; kernel density estimation; expectationmaximization algorithm