一种新颖的MR脊柱图像自动分割算法研究

叶伟，陶晶，陈小宇，林千早，魏玲，张小兵，陈万洪

1.徐州医科大学附属淮安医院 影像科，江苏 淮安 223002；2.淮安市第四人民医院 影像科，江苏 淮安 223002

引言

脊柱是人体中最复杂的承载结构，成人大约长70 cm，分为5个椎体段[1]。椎间盘突出、椎管狭窄和退行性椎间盘是常见的脊柱疾病，MRI是一种行之有效的影像检查手段，但鉴别病变部位和定量分析往往需要借助图像分割技术。由于MR图像自身分辨率和对比度较差，易受噪声、伪影和局部体积效应等影响，使得MR脊柱图像自动分割成为一项具有挑战性的任务。

常用的分割方法有全局阈值法、边缘检测法、区域增长方法等[2-4]。灰度阈值由灰度直方图决定，由于MR图像的灰度分布具有随机性，往往导致全局阈值法无法确定最佳分割阈值。此外，阈值分割法忽略图像像素位置信息，加剧了图像空间不确定性。边缘检测法可用于识别感兴趣区域的轮廓边界，但对边界的连续性和光滑度要求较高，且需要复杂的后处理操作。区域增长法是阈值法的延伸[5]，将满足一定灰度相似性测度准则的像素通过邻域扩展连通成区域，但对初始种子点选取要求较高，且对噪声敏感。

模糊C均值聚类区别于K均值、C均值等硬性分类算法，是指通过优化目标函数得到每个样本点对所有类中心的隶属度，从而决定样本点的类属以达到自动对样本数据进行分类的目的[6]。本文采用基于模糊C均值聚类的分割算法快速准确的提取MR脊柱部分，辅助医生对脊柱退化状态分类和畸形的诊断。

1 方法

首先，运用各向异性扩散滤波对原始图像进行预处理，再采用模糊C均值聚类算法进行图像分割，最后利用形态学方法进行后处理（图1）。

图1 MR脊柱图像分割流程

1.1 图像预处理

各向异性滤波是一种基于偏微分方程的滤波技术，建立于热量的各向异性扩散理论[7]。本研究采用各向异性扩散滤波平滑原始MR脊柱图像，在图像平坦区域选择大尺度平滑，边缘区域选择小尺度评卷，能在抑制噪声的同时保持图像边缘信息。各向异性扩散滤波将图像看作热量场，每个像素看作热流，根据当前像素和领域像素关系，判断是否要向周围扩散[8]。如果某个邻域像素和当前像素差别较大，表示此邻域像素可能处于边界，则当前像素停止向该方向继续扩散，确保图像边缘完整清晰。主要步骤如下：

首先，设置初始值和边界值，然后采用公式(1)进行迭代运算得到预处理图像[9]。其中I表示源图像，t表示迭代次数。

N、S、E、W代表4个散度，在4个方向上对当前像素求偏导数，由公式(2)给出[10]。

其中N、S、E、W 4个方向的导热系数由公式(3)给出[11]，其中导热相关系数k取值越大，平滑效果越好，边缘保留越差。

根据预实验结果，综合考虑计算量和平滑效果，迭代次数取20，导热系数k取15，λ取为0.15。

1.2 核密度估计的模糊C均值聚类算法优化

模糊C均值聚类算法属于无监督学习方法，将普通C均值算法对于数据的硬性划分改进为柔性的模糊划分[12]。主要思想是使得被划分到同一簇的对象之间相似度最大，而不同簇之间的相似度最小，已广泛应用于模式识别和图像处理研究等领域[13]。模糊C均值聚类算法是通过使得已知数量的聚类目标函数达到最小化来获得数据集分类，目标函数由公式(4)给出[14]。

其中，m是一个加权指数，0≤m≤∞，通常m取2.0是合理的[15]；uij是群集的隶属度，0≤uij≤1，且xi是第i个数据点；vj是模糊组j的聚类中心；‖xi-vj‖2是表示数据点与聚类中心间的距离；N表示像素数，而M表示簇中心数，根据预实验结果取值为3。

模糊聚类是公式(4)所示的目标函数的迭代优化，通过更新隶属度函数uij和簇中心vj进行。使目标函数达到最小的必要条件由公式(5)给出[16]：

首先选取(0,1)间的随机数初始化uij，再利用公式(6)计算M个聚类中心vj，和目标函数，如果前后两次目标函数值的改变量小于某个阈值ε，则迭代停止，否则采用公式(5)更新uij，重复上述操作。ε是迭代停止的临界值，本文 ε取 0.01。

模糊C均值聚类算法虽然简洁高效，但初始的M个初始聚类中心vj很难确定，选择不当，容易陷入局部最优解，本文采用核密度估计来确定模糊C均值聚类所需的M个初始聚类中心，核密度估计是一种非参数估计的建模方法，可直接从图像的像素连续变化值中估计出概率密度函数，此时图像像素概率密度估计曲线的峰值点即是初始聚类中心[17]。本文选取高斯函数作为核函数，等概率的选出t个格点的像素灰度值作为观察值。则高斯核密度估计函数关系式如公式(7)所示[18]：

用公式(8)在图像中等概率的选出t个格点，窗宽h由公式(9)确定[19]。

其中Ixy为第x行，第y列的像素点灰度值，yj是第j个格点灰度值，max(Ixy)是图像像素最大值，min(Ixy)是图像像素灰度最小值，h为窗宽，STD是图像像素的标准差，IQR为图像像素的四分位差，N是图像的像素数量。

1.3 图像后处理

本研究采用一系列形态学操作从聚类输出图像中提取椎体。主要步骤如下：

第一步，先进行2次图像膨胀，采用3×3的结构元素扫描图像中的每一个元素，并与其覆盖的二值图像进行“与”运算，填充椎体内细小空洞；再进行2次图像腐蚀，此时结构元素与二值图像做“或”运算，能够在不明显改变物体面积的同时平滑边界，去除孤立点。

第二步，区域滤波。调用区域滤波函数roifilt2(H,I,BW)对二值图像BW中灰度值为1（白色）区域进行滤波处理，忽略灰度值为0（黑色）的区域，滤波结果返回到输出矩阵J中，其中I为目标图像数据矩阵。

第三步，采用区域测度提取椎体。根据椎体和周围组织形态分析，椎体纵横比为1.5～2.0，将此作为区域测度有助于剔除脊柱周围的韧带和肌肉。

1.4 图像质量评价指标

使用 Dice系数（Dice Similarity Coefficient，DC）和HD距离（Hausdorff Distance，HD）评价图像分割结果，将放射科医师对图像的手动分割结果作为参考标准。

Dice系数用于度量两个集合的一致性[20]，由公式(10)给出。其中X，Y表示两个集合，分子表示两个集合的相交操作后的点集合，分母表示两个集合总和，Dice系数越接近于1，图像分割结果越精确。

HD距离用于评价两个二进制点集之间的相似性[21]，见公式(11)。其中，h(X,Y)称为从集合X～Y的单向HD，见公式(12)。HD是测量点的相似位置，能够抵抗点位置的微小波动，准确测量点集之间非相似性。HD越接近于0，图像分割越精确。

2 结果与讨论

为了定性、定量评估本文提出算法的优越性，强健性和实用性，选取5幅MR脊柱T1加权矢状位平面图像作为实验对象，图像矩阵大小均为512×512，具有512个灰阶，所用算法均在MATLAB 2013a编程环境下仿真实现。

2.1 不同算法脊柱分割结果

应基于模糊C均值分割算法、Otsu阈值分割算法和K均值分割算法的MR脊柱提取过程，见图2～4。首先对源图像进行预处理，接着采用不同算法进行分割，然后采用腐蚀、滤波和基于纵横比的形态学处理。结果表明，模糊C均值聚类分割方法能完整的保留脊柱信息，而基于Otsu阈值算法和K均值聚类方法均只分割出部分椎体，脊柱信息严重不完整，这是由于Otsu法难以确定最佳阈值，导致分割不足；K均值算法分聚类中心和数目对分割结果影响很大，容易出现分割不足和过分割。

2.2 临床实例分析

为了验证本文分割算法的临床应用价值，选取4例轻度腰痛患者的矢状位MR脊柱T1加权图像，其中包括2名女性和2名男性，年龄为45～60岁。4名患者的脊椎分割结果，见图5，由模糊C均值聚类分割后进行标记着色与原始图像叠加所得。可见模糊C均值聚类分割算法能够适用于多种不同的脊柱类型，为临床医生诊断提供便利。

图2 基于模糊C均值聚类分割结果

图3 基于Otsu阈值分割结果

图4 基于K-均值聚类分割结果

图5 4例腰痛患者脊柱分割结果

此外，分别采用Otsu阈值算法、K均值聚类算法和模糊C均值聚类算法对4例患者的脊柱图像进行分割，并以放射科主任医师手动分割结果作为验证参考标准，分别计算不同分割算法所得的DC值和HD值，结果见图6。DC值大小排序依次为模糊C均值聚类算法、K均值分割算法和Otsu阈值法，其中基于本文提出算法的DC值分别为0.835，0.900，0.852和0.830（接近于1），HD值在模糊C均值分割算法中最小，分别为0.040，0.040，0.036和0.040，其次是K均值聚类分割算法，Otsu阈值分割算法的HD值最大。表明基于模糊C均值分割算法所得的分割结果与参考标准相似性最强，分割结果最稳定，其次是K均值聚类算法，Otsu阈值分割效果最差。

图6 不同分割算法的Dice系数（a）和HD距离（b）比较

3 结论

本文提出了一种基于改进的模糊C均值聚类的无监督自动分割算法，并应用于脊柱MR冠状位平面图像的分割。其中各向异性滤波函数能平滑原始图像，高斯核密度估计可以为模糊C均值聚类算法选出最佳的聚类的中心和聚类数目，模糊C均值算法得到二值化初始分割图像，一系列的形态学操作成功提取出椎体感兴趣区域，对于退行性腰椎病变等腰椎疾病的定量分析具有重要指导意义。定量和定性评估均表明模糊C均值聚类算法能快速准确提取出脊椎区域，明显优于K均值聚类算法和Otsu阈值法，具有很强的鲁棒性，是一种可行的MR脊柱图像分割算法。

[参考文献]

[1] 刘松涛,殷福亮.基于图割的图像分割方法及其新进展[J].自动化学报,2012,38(6):911-922.

[2] 王江涛,练煜,石红岩,等.结合模糊C均值聚类和边缘检测算法的彩色图像分割[J].兰州文理学院学报(自然科学版),2015,29(2):61-65.

[3] 刘俊,马燕,陈坤.一种改进的基于区域生长的彩色图像分割方法[J].计算机应用与软件,2012(11):288-291.

[4] 曾文权,何拥军,崔晓坤.基于各向异性滤波和空间FCM的MRI图像分割方法[J].计算机应用研究,2014,31(1):316-320.

[5] 付丽娟,姚宇,付忠良.中值滤波与各向异性扩散相结合的医学图像滤波方法[J].计算机应用,2014,34(1):145-148.

[6] Kannan SR,Ramathilagam S,Sathya A,et al.Effective fuzzy c-means based kernel function in segmenting medical images[J].Comput Biol Med,2010,40(6):572-579.

[7] Li BN,Chui CK,Chang S,et al.Integrating spatial fuzzy clustering with level set methods for automated medical image segmentation[J].Comput Biol Med,2011,41(1):1-10.

[8] Gong M,Liang Y,Shi J,et al.Fuzzy C-means clustering with local information and kernel metric for image segmentation[A].IEEE Transactions on Image Processing A Publication of the IEEE Signal Processing Society[C].2013,22(2):573-584.

[9] 杨素华,陈琼,罗艳芬.基于Graph-Cuts的脑部MRI图像脑组织提取方法[J].中国生物医学工程学报,2014,33(5):525-531.

[10] 黄会营,庄淑君.改进的正交匹配追踪图像快速重建算法[J].光学技术,2014,40(6):515-519.

[11] 谢永华,徐赵飞,范文晓.基于高斯尺度空间粗糙度描述子的花粉图像分类识别[J].计算机应用,2015,35(7):2039-2042.

[12] 陈金位,吴冰.二维直方图重建和降维的Otsu阈值分割算法[J].图学学报,2015,36(4):570-575.

[13] 徐超,黄风华,毛政元.一种改进的二维Otsu阈值分割算法[J].电子技术应用,2016,42(12):108-111.

[14] Ananthi VP,Balasubramaniam P,Lim CP.Segmentation of gray scale image based on intuitionistic fuzzy sets constructed from several membership functions[J].Pattern Recogn,2014,47(12):3870-3880.

[15] 周啸虎,高伟,张子齐.基于质子密度和弛豫时间的大脑MR图像分割新算法[J].中国医疗设备,2016,31(10):25-28.

[16] 宋国权,李金锋.基于聚类算法的脑部MR图像分割[J].中国医疗设备,2017,32(1):26-29.

[17] Gómez D,Yáñez J,Guada C,et al.Fuzzy image segmentation based upon hierarchical clustering[J].Knowl-Based Syst,2015,87(C):26-37.

[18] Dubey YK,Mushrif MM.FCM clustering algorithms for segmentation of brain MR images[J].Adv Fuzzy Syst,2016(3):1-14.

[19] 朱然,李积英.几种优化FCM算法聚类中心的方法对比及仿真[J].计算机技术与发展,2015(5):17-20.

[20] Š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.

[21] 李俊林,符红光.改进的基于核密度估计的数据分类算法[J].控制与决策,2010,25(4):507-514.