改进的分水岭图像分割算法

孙惠杰，邓廷权，李艳超

(1.哈尔滨工程大学计算机科学与技术学院，黑龙江哈尔滨150001;2.哈尔滨师范大学计算机科学与信息工程学院，黑龙江哈尔滨150025;3.哈尔滨工程大学 理学院，黑龙江 哈尔滨150001;4.哈尔滨师范大学 黑龙江省智能教育与信息工程重点实验室，黑龙江哈尔滨150001)

图像分割是把图像分成各具特性的区域并从中 提取感兴趣目标的技术和过程，它是图像处理的重要环节。分割结果的好坏直接影响后续的图像分析、理解及识别。图像分割的方法很多，其中阈值分割、区域生长、分水岭分割算法、聚类算法、神经网络分割算法等都是一些经典方法。

分水岭算法起源于Digabel等对简单二值图像的分析，是一种基于二值数学形态学的分割算法［1］。为了增加模型的适用范围，Beucher等［2］建立了灰度图像分割的分水岭理论和算法。由于图像中噪声点和过多局部极小点的存在，使得分水岭分割算法存在严重的过分割现象。因此降低过分割现象是当今研究热点之一。目前，分水岭算法的改进算法［3-4］很多，大体可以分为3类:1)对原始图像在空间域进行预处理(区域平滑等)以消除过多的区域极值，进而减少连通分割类的个数［5］;2)对原始空间域图像进行特定的有意义的变换，根据其性质在变换后的图像上运用分水岭算法进行图像的连通分割，如在标记图像上进行分水岭分割［6］，将小波变换与分水岭算法相结合对分水岭算法进行改进［7］等;3)在原始图像进行分水岭变换并利用合适的区域合并方法对邻接的同质连通区域进行融合以消除过分割现象。

本文旨在结合区域生长原理和图像的视觉效果，基于优化方法对分水岭算法进行改进，提出一种新型的图像分割算法。

1 分水岭算法与区域生长

1.1 分水岭算法

分水岭算法［8］是基于模拟浸水过程实现的。将图像看成3D地形表示，即2D的地基(对应图像空间)加上第3维的高度(对应图像灰度值)。图像中存在一些局部极小点。假设在每个极小点打一个孔，水将从这些孔中慢慢浸入表面，从最低的极小值点开始，水逐渐淹没图像的积水盆。当来自2个不同极小值点区域的水面不断升高要汇集到一起时，在此筑起一道堤坝。在流溢过程的最后，每个极小值点被相应积水盆的堤坝所包围，整个堤坝集合构成分水岭。不同的积水盆代表图像的不同分区，实现图像分割。

实验表明在梯度图像上进行分水岭分割比在原始图像上做分割得到的结果更加准确，因此本文利用如下多尺度形态学梯度算子求取梯度图像［9］进行分水岭变换:

式中:GRAD为形态学梯度，Bi结构元素可以是任何形状，只要满足关系B0⊆B1⊆…⊆Bn即可。

1.2 区域生长

区域生长是一种经典的图像分割方法，其基本思想是将具有相似性质的像素集合起来构成区域。具体来说就是先对每个需要分割的区域找一个种子像素作为生长的起始点(种子点)，将种子点邻域中与种子像素有相同或相似性质的像素(根据某种事先确定的生长或相似准则来判定)合并到种子像素所在区域。将该区域中合并的像素当作新的种子点继续进行上面的过程，直到没有满足条件的像素被包括进来为止，这样一个区域就生长成了。这样的过程一般用四叉树算法可很快实现。

2 基于区域生长的分水岭分割算法

经典分水岭图像分割算法会在每2个相邻的极小点区域之间设定一个分水线，导致过分割现象的出现，影响图像分割的质量和视觉效果。

2.1 基于区域生长的分水岭算法

本文将区域生长算法和分水岭算法融合实现图像分割。借助于地形学解释，当某2个相邻的极小点区域的相邻灰度值小于某一给定阈值T时，不在这2个积水盆之间筑堤坝，而把它们合并成一个区域。

设灰度图像f的最小灰度值为hmin，最大灰度值为hmax，对于某一确定的阈值T，详细的融合图像分割算法可描述如下。

水流最先到达图像的灰度最小点hmin，将这样的点记为Mi(i=1，2，…，n)，n为像素值为hmin的点的个数。以Mi为种子点进行区域生长，用N8(Mi)表示点Mi的8邻域，经过区域生长之后，点Mi所对应的区域为

CBT(Mi)={p|p∈N8(Mi)，|f(p)-f(Mi)|≤T}式中:CB为积水盆;将CBT(Mi)中每个点p又作为种子点进行区域生长，生长准则和上述方法一样，直到没有满足生长条件的点为止。假设共对N个种子点进行了区域生长，hmin对应的积水盆为

从而得到hmin对应的积水盆。

将当前灰度值增加1，记为h，则h=hmin+1，重复以上步骤，直到h=hmax。得到所有积水盆，将积水盆的边缘作为分水线。

在该融合算法中，区域生长阈值T的选取决定了融合算法的效果。一般情况下，结合图像的具体特点和人的视觉特性，通过经验和试验获得。这种处理方法带有很大的主观性和盲目性。

为了获得区域生长算法中的最优阈值T，本文结合图像分割常用评价指标确定一个目标函数，通过寻求该目标函数的最优值确定其阈值。

2.2 目标函数的确定

目前评价图像分割结果好坏的准则很多［10-12］，其中区域间对比度，区域内部均匀性等是经典且常用的图像分割质量评判指标。

设Ri、Rj表示运用阈值T进行区域生长分水岭分割之后的2个相邻连通部分，Ai表示区域Ri的大小(面积)，μ 表示区域Ri的平均灰度，ki表示与Ri相邻的区域个数，定义归一化系数:

它们分别表示分割图像连通部分的最大方差和相邻连通部分间的灰度平均值的最大差。记

式中:Ti表示Ri与其相邻区域的平均对比度，而Mi表示区域Ri内部的均匀性。

为说明阈值对分割质量的影响，令区域生长阈值在［1，50］逐渐变化，对图1(a)进行区域生长分水岭算法分割，利用式(3)、(4)分别计算分割结果图像的区域内部均匀性和区域间对比度，以及它们各自的平均值，所得结果如图1(b)、(c)所示。从图中可以看出，区域间的对比度均值随着分割阈值的增大呈现增大的趋势，而区域内部均匀性则随着分割阈值的增大呈现减小的趋势。

经过实验和分析发现，当生长阈值增大时，分割出的连通部分面积将越来越大，其区域内部的均匀性将越来越小，而区域间的对比度越来越大。

好的图像分割效果意味着连通区域内部的一致性较好(均匀性较大)，相邻区域间的对比度也较大。如何限制合适的阈值进行区域增长以实现上述目的变得十分关键。最优阈值T的选择正是基于一致性和对比度这对矛盾评价指标的平衡和折中。本文运用熵来刻画这些评价指标的平衡，将分割图像区域内灰度值的一致性和区域间灰度对比度的香农熵作为评判图像分割质量优劣的目标函数，以此确定区域生长的最佳阈值。

熵是信息论中一个非常基本并有重要应用的概念，它刻画一个对象所蕴含的平均信息量。当分割区域内灰度值的一致性和区域间灰度对比度大时，每一个分割区域都是有意义的，蕴含的信息最多，从而可以用香农熵来量化分割结果所蕴含的信息量，进而评价图像分割结果的好坏。对于图像分割而言，其蕴含的信息量越多，分割效果越好。香农熵定义如下:

当pi为0或1时，熵值取最小值0;pi为1/e时熵值取得最大值，此时研究对象含有最多信息。香农熵在［1，1/e］之间为单增区间。因此为了方便且不引起歧义，将Ti和Mi规范化到得到［1，1/e］之间，得到Mi’=Mi/e，Ti’=Ti/e。设K表示分割图像连通区域的个数，则分割图像连通区域内部均匀性的熵和对比度的熵为

为了更准确地度量分割效果，确保分割图像连通区域不致过大，也不致过小，更加符合人的视觉特性，将分割区域的内部均匀性和区域间对比度的极端度量(最大值的熵和最小值的熵)也作为评判标准的一部分，形成如下目标函数H(T):

使H(T)达到最大的T所对应的分割结果蕴含的信息最多，分割结果是最有意义的。T即为区域生长的最优阈值。

2.3 基于粒子群的分水岭图像分割算法

基于熵最大的目标函数(6)不具有解析解，通常需要优化算法而非穷举法在有效的时间范围内可求得最优值。本文采用粒子群优化算法确定式(6)最优解，从而得到图像分割算法的最优阈值参数T。

粒子群优化算法［13-14］是一种基于群智能的优化算法，其把优化问题的解抽象成粒子运动的最优状态。1999年Clerc［15］提出带有收缩因子的粒子群优化算法，使得粒子群算法具有更好的收敛性。本文正是采用这种优化算法实现式(6)的最优解。该算法需要设定初始种群规模M，初始种群的速度:

初始位置:

式中:Vi(0)和Xi(0)分别表示第i个粒子的初始速度和初始位置。粒子的速度向量的更新公式为

结合图像分割特点，其粒子的位置向量的每个分量都应该是整数，其更新公式采用:

式中:［*］表示最接近* 的整数，rand表示［0，1］之间的随机值，Vi(t)为粒子i在第t次迭代的更新速度，要求在［-Vmax，Vmax］，Vmax为事先给定的最大更新速度。参数k为收缩因子，C1、C2为学习因子，C1控制粒子向自身学习的能力，C2控制粒子向群体学习的能力。依实验结果，取k=0.729，C1=C2=2.05。粒子群的初始状态X(0)和初始速度V(0)可随机选取，但要求X(0)在解空间内，V(0)在［-Vmax，Vmax］内。Pi(t)表示粒子i在第t次迭代时经过的最佳位置，由

确定，G(t)表示所有粒子在寻找最优解的过程中第t次迭代时所经过的最佳位置，其中:

基于粒子群分水岭图像分割算法的步骤如下:

1)粒子群初始化。确定阈值T:A＜T＜B，选取初始种群规模M，设定最大迭代次数tmax。设定粒子运动的最大速度为Vmax，在解空间中选择M个数作为粒子群的初始位置X(0)，并在［-Vmax，Vmax］内随机选择初始速度V(0)。

2)对于初始种群X(0)中的每个分量，依据2.1节所描述的方法求得其分割图像，利用式(1)、(2)计算Xi(0)对应的CT和DT。取Cmax=max{CT}，Dmax=max{DT}。以C，D为归一化系数，根据式(6)分别计算Xi(0)对应的目标函数H(Xi(0))。

3)根据式(9)、(11)，得到每个粒子的初始最优位置Pi(0)和所有粒子的最优位置G(0)。

4)令t=1。

5)依据式(7)、(8)更新速度和位置V(t)，X(t)。

6)依据2.1节所描述的方法求得Vi(t)对应的分割图像，利用式(1)、(2)计算Vi(t)对应的CT和DT。取Cmax=max{CT，C}，Dmax=max{DT，D}。以C、D为归一化系数，根据式(6)分别计算Xi(t)对应的目标函数H(Xi(t))与Pi(t-1)的目标函数H(Pi(t-1))。

7)根据式(10)，得到每个粒子在第t次迭代时所经过的最佳位置Pi(t)。根据式(12)，得到所有粒子在第t次迭代时所经过的最佳位置G(t)。

8)若t≤tmax，t=t+1，返回 7);若t＞tmax，则将G(t)作为区域生长的最优阈值，对应的图像为最优分割结果，每个区域的边界为分水线。算法结束。

经该算法后，图像的过分割现象会得到很大程度的缓解，但是分割区域仍然可能较多，主要原因是图像存在噪声，算法会将噪声点及其邻域像素聚类形成一个面积很小的区域。因此，依据实际情况，结合人的视觉感知特性对分割结果进行后处理。

设定一个临界值L，对于面积小于临界值L的区域Ri，计算Ri与其相邻区域间的灰度差异度:

式中:dif为灰度差异度，将Ri合并到差异度最小区域中，得最终分割结果。

设图像共有N个像素点，在基于区域生长分水岭算法中，因每个像素点都要和其周围的8个像素点比较大小，令k=8，则这一部分的时间复杂度为O(k×N)。当粒子群的规模为M，迭代次数为tmax时，本文算法总的时间复杂度为:O(M×tmax×k×N)。

3 实验结果与分析

为验证本文图像分割算法的有效性，在MATLAB 7.1平台上实现了该算法，并对实验结果进行分析。

3.1 参数的设置对实验结果的影响

本文算法中涉及的参数包括:种群规模M、粒子的最大运行速度Vmax、生长参数阈值T的取值范围和迭代次数tmax以及临界值L。

M越大代表种群包含粒子越多，寻优能力越强，但计算量越大［16］。一般选取M为一个较小的整数值即可达到寻优目的。由于区域生长阈值参数T(灰度值，整数)的取值范围有限，T很小时会使得区域生长的效果很不明显，较大时有可能会使图像产生欠分割现象，因此T不能太小，也不能太大。根据图像分割的特点，如果图像灰度变化比较剧烈，区域生长阈值参数T的取值范围［A，B］应稍大一些，而对于灰度值变化比较平缓的图像，T的取值范围应稍小一些。Vmax控制粒子的运动速度，Vmax太大，粒子可能会飞过最优区域，太小可能导致粒子无法跳出局部最优。Vmax的值也与T的值有很大关系，由于T的值有限，而Vmax控制粒子从一个值向另一个值转变的速度，因此Vmax也取的比较小。迭代次数tmax取的越大，计算时间越长，效果越好。初始种群V(0)是进行寻优的初始值，在解空间确定的情况下，标准粒子群算法规定可以随机确定初始种群。如果选择得当，会使算法具有更好的收敛性。在种群规模M和T的取值范围［A，B］确定的情况下，为了尽量避免粒子群算法陷入局部最优，使初始粒子的值均匀分布在T的取值范围内，依据:

确定位置向量中第j个粒子的初始值。

大量实验表明，对于很多图像，当T=50时，已经产生欠分割现象，如图2(e)～(h)所示。另一方面，当T＜10时，大多数图像存在过分割现象，如图2(i)～(l)所示(其为图2(a)～(d)在生长阈值为10时的分割结果)。

图2 阈值分别为50和10时，lena、pepper、house、brain图的分割结果Fig.2 Segmented results for images at different thresholds

因此一般可取10 ＜T＜50，在［11，49］求取最优值，当然也可以视不同的情况做一些调整。

当10 ＜T＜50 时，在式(7)中取k=0.729，C1=C2=2.05。令M=3，Vmax=4，tmax=40，利用式(13)计算初始值，基于本文算法，得到图2中的4幅图像的收敛情况如图3所示。从图中可以看到，对于house和brain图像，利用本文算法可以很快的收敛到最优值，但是pepper图像的收敛速度比较慢，lena图并没有收敛到最优值，其主要原因为house与brain图像的目标相对单一，且图像中不存在纹理区域，而pepper图像目标相对比较复杂，lena图像中存在纹理区域。

在M=8，Vmax=2，tmax=40的情况下，图2中的4幅图像的收敛情况如图4所示。从图3、4可以看出，对于比较复杂的图像(比如说pepper图像，lena图像)，应该将种群规模取的大一些，以免算法陷入局部最优。当种群规模比较大时，应相应的将最大速度取的小一些，以免粒子跳过全局最优值，而对于比较简单的图像(如house、brain图像)，种群规模取的较大，反而会增加运行时间，因此应将这一类图像的种群规模取的较小一些。一般M在［3，8］取值即可。同时根据M调整Vmax的大小。迭代次数tmax也不用过大，一般迭代5～10次即可收敛到全局最优值。

基于如上所述分割算法及确定的参数值，得到图5中图像(a)～(d)的最优分割结果，对应的最优生长阈值分别为35、27、29、26。对于临界值L的选择，主要根据分割区域的分布及大小并兼顾实际需要设定L的值，将无意义的小区域合并到与之相邻的区域中。图5(e)～(h)为合并后的图像分割结果，对应的临界值分别为 200、200、50、100。

图3 优化变量敏感性分析Fig.3 The sensitivity analyses of the optimization variables

图 4 lena、pepper、house、brain 图的收敛结果Fig.4 Convergence for lena，pepper，house and brain images

图5 合并之后的最终分割结果Fig.5 Final results of segmengtation

3.2 与现有算法的比较

利用改进分水岭算法对图2(a)进行分割。图6(a)为利用文献［17］所提出的改进分水岭算法对lena图像进行分割得到的结果，该算法将lena图像分成了64个区域，较传统分水岭和梯度分水岭算法有了很大程度的改进。图6(b)为利用本文算法对lena图像进行分割的分割结果，将lena图像分成了23个区域，其分割结果较传统分水岭算法、梯度分水岭算法以及文献［17］中提到的算法都有很大程度的改进，而且分割结果与人的直观视觉相一致。图6(d)为利用本文算法对6(c)的分割结果，(e)为在原图上叠加分割线的结果，分割效果不错，6(f)为利用文献［18］中提到的算法进行分割的结果，可以看出，本文算法的分割结果更加符合人的直观视觉。

图6 不同分割方法的比较Fig.6 Comparison of results of different segmentation methods

4 结束语

本文改进了现有的图像分割算法，提出一种区域生长与分水岭分割算法相结合的灰度图像自动分割算法。1)将灰度值接近的极小区域合并成为一个区域，进而改善分水岭算法的分割效果。2)引入生长阈值，使得分割区域内部具有较好的均匀性和对比度。3)结合香农熵提出一种图像分割结果评判准则，依此寻找最优分割结果。

实验证明，新算法有效地解决了过分割现象，是一种有效实用的图像分割和目标提取方法。

［1］DIGABEL H，LANTUEJOUL C.Iterative algorithms［C］//Proc 2nd European Symp Quantitative Analysis of Microstructures in Material Science，Biology and Medicine.Sturrgart，West Germany:Riederer Verlag，1978:85-99.

［2］VINCENT L，SOILLE P.Watersheds in digital spaces:an efficient algorithm based on immersion simulation ［J］.IEEE Transactions on Pattern Analysis and Machine Intelligence，1991，13(6):583-598.

［3］NG H P，ONG S H，FOONG K W C，et al.Masseter segmentation using an improved watershed algorithm with unsupervised classification［J］.Computers in Biology and Medicine，2008，38:171-184.

［4］RAO A R，SRINIVAS V V.Regionalization of watersheds by fuzzy cluster analysis［J］.Journal of Hydrology，2006，318:57-79.

［5］ZHANG M，ZHANG L，CHEN H D.A neutrosophic approach to image segmentation based on watershed method［J］.Signal Processing，2012，90:1510-1517.

［6］FLORES F C，LOTUFO R A.Watershed from propagated markers:An interactive method to morphological object segmentation in image sequences［J］.Image and Vision Computing，2010，28:1491-1514.

［7］JUNG C R.Combining wavelets and watersheds for robust multiscale image segmentation［J］.Image and Vision Computing，2007，25:24-33.

［8］王小鹏.形态学图像分析原理与应用［M］.2版.北京:清华大学出版社，2008:58.WANG Xiaopeng.Morphological image analysis principles and applications［M］.2nd ed.Beijing:Tsinghua University Press，2008:58.

［9］WANG D.A multiscale gradient algorithm for image segmentation using watershed［J］.Pattern Recognition，1997，30(12):2043-2052.

［10］章毓晋.图像分割评价技术分类和比较［J］.中国图像图形学报，1996，1(2):151-158.ZHANG Yujin.A classification and comparison of evaluation techniques for image segmentation［J］.China Journal of Image and Graphics，1996，1(2):151-158.

［11］POLAK M，ZHANG H，PI M H.An evaluation metric for image segmentation of multiple objects［J］.Image and Vision Computing，2009，27:1223-1227.

［12］CORDENESA R，LUIS-GARCIA R，BACH-CUADRAB M.A multidimensional segmentation evaluation for medical image data［J］.Computer Methods and Programs in Biomedicine，2009，96(2):108-124.

［13］王科俊，郭庆昌.基于粒子群优化算法和改进的Snake模型的图像分割算法［J］.智能系统学报，2007，2(1):53-58.WANG Kejun，GUO Qingchang.Image segmentation algorithm based on the PSO and improved Snake model［J］.CAAI Transactions on Intelligent Systems，2007，2(1):53-58.

［14］EBERHART R C，KENNEDY J.A new optimizer using particle swarm theory［J］.Sixth International Symposium on Micro Machine and Human Science，1995(3):39-43.

［15］CLERC M.The swarm and queen:towards a deterministic and adaptive particle swarm optimization［J］.IEEE International Congress on Evolutionary Computation，1999(3):1951-1957.

［16］李艳超.基于分割区域的图像压缩方法研究［D］.哈尔滨:哈尔滨工程大学，2011.LI Yangchao.Research on image compression based on image segmentation［D］.Harbin:Harbin Engineering University，2011.

［17］LUIS P.Fuzzy relations applied to minimize over segmentation in watershed algorithms［J］.Pattern Recognition Letters，2005，26:819-828.

［18］李苏祺，张广军.基于邻接表的分水岭变换快速区域合并算法［J］.北京航空航天大学学报，2008，34(11):1327-1348.LI Suqi，ZHANG Guangjun.Fast region merging algorithm for watershed transform based on adjacency list［J］.Journal of Beijing University of Aeronautics and Astronautics，2008，34(11):1327-1348.