改进分水岭与K-means++结合的人体脊柱图像分割方法

郎成洪，华云松，张嘉棋，钟雪莲，王麒翔

（上海理工大学光电信息与计算机工程学院，上海 200093）

0 引言

人体脊柱图像的准确分割对于脊柱微创手术十分重要，图像分割效果会直接影响其三维重建图像效果以及微创手术前病灶信息的获取。由于不同成像设备和环境条件的影响，使得CT图像常常包含噪声。为提高图像分割效果，人们提出了许多图像分割方法，主要分为基于阈值、边缘、区域、聚类分析、能量泛函和人工神经网络的分割方法。

基于阈值的分割方法是根据一个或几个阈值对图像中的灰度值进行分类，常用的有自适应阈值法、最大熵法、模糊阈值法、类间阈值法等。此类方法计算简单、高效、速度快，但对噪声比较敏感，且对灰度差异不明显的目标分割效果不佳。如李长春等使用基于多阈值分割与B样条插值的方法去除CT图像中的金属伪影。

基于边缘的分割方法通过检测不同区域边缘来分割不同目标，边缘检测效果好且速度快，但无法得到好的区域结构，且抗噪性差。常见的边缘检测算子有Robert算子、Sobel算子、Prewitt算子、Canny算子等。如李庆忠等提出基于改进Canny算子的图像弱边缘检测算法，该算法通过改进梯度倒数加权滤波的方法检测平滑块中细小的弱边缘。

基于聚类分析的图像分割方法根据相似性准则，将图像像素在特征空间进行分类聚集并进行分割，然后将分割结果映射回原图像空间。此类方法主要包括K均值和模糊C均值聚类算法，但其主要缺点是灰度值不均匀、对噪声比较敏感，且运算开销大。如姜文斌等提出基于自适应权重法的K-means模型对遥感图像进行分割的方法，该方法利用二维高斯函数对遥感图像进行滤波平滑处理，然后根据收敛度和自适应值找到最优权重作为初始聚类中心，以提高图像分割精度。

基于能量泛函的分割方法是通过连续曲线表示目标边缘，然后定义一个能量泛函使其自变量包括边缘曲线，从而进行图像分割，最具代表性的是Osher提出的水平集模型。如黄国鹏等利用偏移场校正和局部区域信息拟合来构建能量函数，并通过加性算子分裂策略提高图像分割精度和效率。

基于人工神经网络的分割方法是通过训练多层感知机得到线性决策函数，然后通过决策函数对图像像素进行分类，从而分割出目标的过程。其优点是有效解决了图像灰度不均匀和噪声问题，但需要大量数据集，且速度很慢。如唐思源等提出一种基于改进粒子群优化算法的BP神经网络医学图像分割方法，提高了图像分割精度。

基于区域的分割方法主要采用分裂与合并两种方法，把具有相似性的像素连通，构成最终的分割区域。其优点是区域特征较好，且能有效克服图像分割空间小的缺点，但容易出现过分割现象。其中分水岭变换分割是最强大、应用最广泛的基于区域的图像分割方法之一。针对其过分割的缺点，人们提出一些改进的分水岭分割方法。如周小明采用最大熵算法对梯度图像进行校正，然后进行图像分割，取得了良好的图像分割效果。

随着机器学习算法的快速发展，越来越多机器学习算法被应用于解决分水岭算法的过分割问题。机器学习算法是一种监督算法，需要进行大量训练和相应标记才能对测试图像进行准确分割。但在大多数情况下，训练数据和测试数据的数量不足，无法保证训练图像的质量和相似度，故不需要标记数据的无监督方法在数据短缺情况下表现得更好。无监督的聚类方法不需要训练及标记数据，可自动将相似数据点分组到同一聚类中。模糊Cmeans聚类和K-means聚类算法是最著名的无监督聚类算法，已应用于各种领域的数据聚类。本文采用一种比K-means聚类算法更优秀的聚类算法——K-means++聚类算法对人体脊柱图像进行去噪，并修改区域最小值。相较于K-means，该算法分割效率更高。为了得到更好的分割效果，采用在传统分水岭算法基础上进行改进的标记符控制的分水岭算法对K-means++聚类算法处理后的初始聚类图像进行分割。

1 K-means++聚类算法

在原始人体脊柱CT图像中往往存在灰度相似的区域，这些区域会产生许多不必要的局部极小值，如图1所示。

Fig.1 CT image of human spine图1 人体脊柱CT图像

如果先用K-means对图像进行初始聚类处理，可减少这种不必要的局部极小值，从而改善分水岭算法的过分割现象。K-means++算法为解决K-means算法中聚类中心选取过于随意，且结果对初始聚类中心过于敏感的问题，对初始聚类中心选取方式进行改进，具体步骤如下：

（1）从数据集中随机选取一个样本作为初始聚类中心

c

。（2）计算数据集中每个样本

x

与已选取聚类中心的距离，用

D

(

x

)表示。

（3）计算每个样本被选为下一个聚类中心的概率，与当前已有聚类中心距离较远的样本有更大概率被选为下一个聚类中心。具体概率计算公式如下：

其中，

x

表示样本点，

D

(

x

)表示样本点到聚类中心的距离。（4）重复步骤（2）、（3），直到选择出

K

个聚类中心。（5）计算每个样本分别到

K

个聚类中心的距离，将各样本划分到与其距离最近的聚类中心所属类别中。

（6）针对每个类别计算其所有样本特征均值，将该均值作为该类别新的聚类中心。其均值计算公式如下：

其中，

i

表示该类别中的样本点数，

x

表示样本点。

（7）重复步骤（5）、（6），直到聚类中心位置不再发生变化。

原始人体脊柱图像的像素经过K-means++处理后分成了

K

类，由于在医学图像中，脊柱区域是聚类中心值最大且面积最大的连通区域，故可将聚类中心值最大且面积最大的连通分量

S

(

i

=1，2，…，

k

)作为脊柱图像区域。

其中，argmax表示括号内的最大值，即为脊柱区域面积值。

2 形态学处理与标记符控制的分水岭算法

分水岭变换是一种基于形态学的图像分割工具，其基本思想是把图像看作测地学上的拓扑地貌，把每个像素点灰度值表示为该点的海拔高度，然后模拟泛洪过程。假设在图像的局部极小值面刺穿一个小孔，并让水以均匀的速率从小孔中涌出，随着水深逐步增加，为防止不同集水盆地汇聚到一起，在两个集水盆地间修建大坝，即形成分水岭。

2.1 形态学技术去除噪声

利用K-means++聚类算法初始分割后的人体脊柱图像虽然过分割率下降了很多，但是仍然存在过分割现象，且脊柱区域连通分量会参杂部分背景像素，故对Kmeans++聚类算法初始分割后的脊柱图像进行基于重建的开操作和闭操作，在保留图像重要区域轮廓的同时去除细节和噪声。

2.1.1 开操作与闭操作

形态学处理有两种基本运算，分别为腐蚀和膨胀，开操作和闭操作就是结合腐蚀和膨胀两种基本形态学运算而形成的。开操作是先腐蚀再膨胀，闭操作是先膨胀再腐蚀。

形态学闭操作如下：

其中，

A

表示图像，

B

表示结构元。

形态学开操作如下：

其中，

A

表示图像，

B

表示结构元。

膨胀可填充图像中的小孔洞和图像边缘处的一些小凹陷，对图像的外部进行滤波；腐蚀可去除图像中一些小的“突刺”，对图像内部进行滤波，并缩小图像。开操作和闭操作是一种去噪处理，开操作可滤掉小的“突刺”，闭操作可填充孔洞缺口。

2.1.2 基于重建的开操作与闭操作

采用基于重建的开操作和闭操作清理使用Kmeans++聚类算法处理后的人体脊柱图像，先对图像采用基于重建的开操作，在开操作结果上利用基于重建的闭操作进行滤波重建。在简化初始分割图像的同时，准确定位轮廓分水线，从根源上消除过分割现象。在分割前需要选一个结构元素，用来对图像进行形态学重建。利用该结构元素处理图像后，图像中的灰度跃变急剧增加，可消除梯度对边缘方向的依赖性。

基于重建的开操作：

基于重建的闭操作：

2.2 标记符控制的分水岭算法

在标记符控制的分水岭算法中，首先将使用Kmeans++聚类算法初始分割后的人体脊柱图像转换成灰度图像，应用形态学算子对图像进行去噪处理，确定前景标记并与原始图像叠加，然后确定背景标记，获得原始图像的梯度图像并对其进行修改，最后应用分水岭变换生成分割图像。基于标记符控制的分水岭算法流程如图2所示。

Fig.2 Flow of watershed algorithm based on control marker图2 基于控制标记符的分水岭算法流程

具体步骤如下：

（1）由于脊柱骨骼区域在人体脊柱CT图中一般是灰度值较高的连通像素块，将形态学处理后的图像记为

f

。（2）计算

f

的局部区域最大值，并将局部区域最大值作为前景标记，然后把前景标记的灰度图像转换成二值图像，记为

f

。

（3）将前景图像叠加到原始图像上，某些对象中的前景标记会一直覆盖到对象边缘。

（4）通过对前景图像执行先闭操作后侵蚀的形态学去噪处理，清理标记斑点的边界，使得前景标记向内收缩。

（5）使用阈值分割方法对图像

f

进行分割，并将分割结果作为背景标记图像，记为

f

，以识别背景中的暗像素。（6）通过对背景标记图

f

进行距离变换，使背景标记不会非常靠近想要分割的物体边界，然后使用经过距离变换后的分水岭算法产生图像的分水岭脊线。对图像中所有像素(

i

，

j

)求距离变换公式如下：

（7）使用梯度幅值作为分割函数，并使用Sobel算子计算图像梯度。梯度值总是在物体边缘处高，而在物体内部低。

其中，

I

、

I

分别表示图像在垂直和水平方向的梯度图像，∇I表示图像梯度。

（8）梯度图像被修改，使得在前景和背景标记处仅存在区域最小值。经过修改后的梯度图像表示为：

其中，

IMMIN

(*)表示形态学极小值标定操作，∇I表示图像梯度，

f

为前景标记的二值化图。（9）将修改的梯度图像∇I′作为分水岭算法的输入，计算得到理想的物体分割结果，其过程可表示为：

其中，

WST

(*)表示为分水岭变换操作，∇I表示图像梯度。

3 实验数据及结果分析

3.1 实验数据

本文实验数据采用合作医院的双源64排CT机对患者进行扫描获取的人体脊柱CT图像，并通过DICOM Converter软件将人体脊柱CT图像的DICOM格式数据转换成JPG格式的原始实验图像，二维图像分辨率大小为512×512。本文实验数据集中共有186张不同部位的人体脊柱图像，部分数据集如图3所示。

Fig.3 Experimental data set of partial human spine images图3 部分人体脊柱图像实验数据集

3.2 实验结果分析

本文实验在PyCharm Community Edition 2020.3.5x64和Matlab R2018a软件平台上进行，首先直接应用分水岭变换分割方法对人体脊柱图像进行分割，并查看结果。图4显示了直接对人体脊柱原始图像进行分水岭分割的结果。如图所示，直接应用分水岭变换方法的分割结果显示出很高的过分割率，并导致错误的结果。为提高分割效果，采用K-means++聚类方法进行初始聚类。选用不同聚类数

K

对原始图像进行聚类，不同

K

值的K-means++聚类效果如图5所示。

Fig.4 Effect drawing of direct watershed图4 直接分水岭分隔效果

Fig.5 Effect of K-means++clustering with different K values图5 不同K值的K-means++聚类效果

由先验知识可知，在医学图像中，骨骼区域一般是完整且灰度值较高的连通像素块。通过分析以上不同

K

值的初始聚类图像可知，当

K=

3时，能得到原始图像的最佳聚类效果，在分割出脊柱骨的同时，尽量减少局部极小值点。对使用

K=3

的K-means++聚类算法处理后的初始聚类图像进行分水岭变换，可发现其错误的局部极小值明显较少，过分割率明显降低，如图6所示。

为进一步减少分水岭变换的过分割率，对K-means++算法聚类后的初始分割图像（见图5（b））应用5×5的圆形结构体进行腐蚀，然后对腐蚀图像进行形态学重建，再使用相同的结构体对重建后的脊柱图像进行膨胀处理，最后对膨胀后的图像进行形态学重建处理，使得K-means++处理后的人体脊柱图像变得平滑，如图7所示。

Fig.6 Effect of direct watershed transformation after K-means++processing图6 K-means++处理后的直接分水岭变换效果

Fig.7 Result after morphological denoising图7 形态学去噪处理后结果

人体脊柱医学图像中灰度值较大的区域一般为脊柱区域，故提取经过形态学处理后的结果图像（见图6）的区域最大值，得到前景标记。标记的前景图如图8所示。

Fig.8 Tagged foreground图8 标记的前景图

上述经过形态学处理后的结果图像（见图7）中比较亮的是人体脊柱骨区域，比较暗的是背景，故采用阈值分割法对形态学处理后的结果图像进行分割，将分割结果标记为背景图像，如图9所示。

Fig.9 Background image of tag图9 标记的背景图

为了不让背景标记太接近目标物体，对背景标记图（图9）进行距离变换，以达到细化目标物体的目的。由于梯度总是在边缘处高，而在物体内部低，利用Sobel算子计算原始图像梯度。修改梯度图像，使得梯度图像的区域最小值只出现在标记的前景和背景处，如图10所示。

Fig.10 Gradient graph calculated by Sobel operator图10 Sobel算子计算后的梯度图

最后将修改后的梯度图像（见图10（b））利用分水岭变换算法进行分割，分割结果如图11所示。

Fig.11 Effect drawing of watershed transformation for modified gradient image图11 修改后的梯度图像进行分水岭变换效果

3.3 比较实验分析

对实验结果采用

Dice

、

Jaccard

和

Precision

3种图像分割评价指标进行分析，3种指标取值范围均在0～1之间。该值越大，表示分割效果越好，具体定义如下：

其中，

SEG

表示算法对人体脊柱图像自动分割的结果图，

GT

表示采用MITK软件手动分割的结果图。MTIK软件手动分割结果如图12所示。

Fig.12 MITK manual segmentation result图12 MITK手动分割结果

对原始人体脊柱CT图像分别采用传统分水岭算法与K-means++结合的分割算法（K-means++and Watershed，KM++-W）、基于控制标记符分水岭与K-means结合的分割算法（K-means and Marker-Controlled Watershed，KMMCW）、基于神经网络的U-NET算法及本文提出的控制标记符分水岭算法与K-means++结合的分割算法（Kmeans++and Marker-Controlled Watershed，KM++-MCW）进行分割实验，然后对分割结果图进行二值化处理。4种算法分割效果如图13所示。

Fig.13 Segmentation rendering of four algorithms图13 4种算法分割效果

将分割结果分别代入式（12）-式（14）进行计算，得到4种分割算法的

Dice、Jaccard

和

Precision

评价指标，并对其进行比较，如表1所示。

Table 1 Three evaluation coefficient values of four segmentation algorithms表1 4种分割算法的3种评价系数值

由表1可以看出，KM++-MCW算法的3种图像分割评价指标数值均比另外3种算法高，从而证明了本文提出的人体脊柱图像分割算法具有较高的准确度及良好的图像分割效果。

对100幅不同部位的人体脊柱CT切片图像分别采用3种算法进行实验，并计算其

Dice

、

Jaccard

和

Precision

指标，指标数值分布如图14所示。从图14可以看出，4种算法的

Dice

系数均值分别为0.072 1（KM++-W）、0.726（KM-MCW）、0.740（U-NET）、0.892（KM++-MCW），

Jaccard

系数均值分别为0.052 5（KM++-W）、0.593（KM-MCW）、0.696（U-NET）、0.823（KM++-MCW），

Precision

系数均值分别为0.062 8（KM++-W）、0.736（KM-MCW）、0.751（U-NET）、0.854（KM++-MCW）。从4种算法的各评价系数均值可以看出，本文算法的分割效果较另外3种算法有明显提高。其中，U-NET是深度学习算法，需要大量数据进行训练并打上标签，在数据量较少且无标签的情况下，无监督学习的图像分割算法更具有优势。本文提出的人体脊柱图像分割算法相比U-NET算法，3种评价指标平均值分别高出了约15.2%、12.7%、10.3%，对人体脊柱图像的分割具有良好的精确度，能有效分割出人体脊柱图像中的脊柱骨骼区域。

Fig.14 Index value distribution图14 指标数值分布

4 结语

本文提出一种结合K-means++聚类算法与改进分水岭算法的人体脊柱图像分割方法，如果直接对人体脊柱原始图像应用分水岭变换方法进行分割会导致很高的过分割率。为防止这种情况，首先利用K-means++聚类算法对原始图像进行初始分割，然后应用形态学技术进行去噪处理，标记出前景图像和背景图像，计算原始图像的梯度图像并对其进行修改，使得前景和背景标记处仅存在区域最小值，最后对修改后的梯度图像应用分水岭变换得到分割结果。实验结果证明，该方法能很好地降低分水岭变换的过分割率，并能有效将原始图像中的人体脊柱骨骼区域分割出来。但是本文只针对人体脊柱图像进行了实验，对于如何将该算法运用于人体其它组织和器官的分割尚未进行研究，这也是未来研究的重点方向。