基于高斯混合模型和核密度估计的全身骨骼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