基于直觉模糊C均值聚类的图像分割方法与流程

本发明属于数字图像处理领域，具体涉及到一种图像分割方法，可用于图像识别和计算机视觉的预处理。

背景技术：

从上个世纪七八十年代起，众多学者持续关注图像分割，图像分割技术已经成为众多领域的基础技术，只要关乎提取图像中内容的问题就不能缺少图像分割技术，图像分割工作的质量和效果将直接或者间接的影响后续的图像工程。现有的图像分割方法种类繁多，可以归纳为基于阈值的图像分割方法、基于边缘的图像分割方法、基于区域的图像分割方法、基于聚类的图像分割方法等，现有图像分割技术以基于阈值的图像分割方法和基于聚类的图像分割方法为主，基于阈值的图像分割是一种采用统计学信息分割的方法，通常提取合适的阈值，把图像中感兴趣的目标部分分离出来。但是其只适用于灰度差异明显的图像，此外，噪声等干扰元素也对图像的精确分割工作造成困难。到目前为止依然没有形成一种统一的理论体系来应用于不同的领域，图像分割依旧是一个热门研究领域。

基于聚类的图像分割方法根据图像内像素点之间的相似性把图像划分为多个类别，其具有分割效果好，耗时短，方法运行稳定等特点，近年来成为应用普遍的方法之一。常用的聚类分割方法有：k均值、模糊c-均值聚类、分层聚类和谱聚类)等等，其中模糊c-均值聚类是应用最为广泛的方法之一。模糊c均值方法最早由邓恩提出，参见：邓恩.一种和迭代自组织分析算法相关的模糊图像处理方法及其应用于检测紧凑易分离聚类.控制论学报.1973,3(3):32-57.(j.c.dunn.afuzzyrelativeoftheisodataprocessanditsuseindetectingcompactwell-separatedclusters[j].journalofcybernetics,1973,3(3):32-57.)，该方法由于没有考虑图像中的任何空间信息，对图像中的噪声敏感，使得分割含噪图像得不到理想的分割效果。后来许多学者发表了结合图像空间邻域信息改进方法。艾哈迈德等人将空间邻域信息项加入模糊c均值的目标函数中，提出了fcm_s算法，该算法虽然提高了对噪声的鲁棒性，但是计算复杂度较高，为了降低fcm_s算法的计算复杂度，陈松灿和张道强通过均值滤波和中值滤波的方式将邻域信息引入到算法的目标函数中，提出了fcm_s1和fcm_s2算法，参见：陈松灿,张道强.一种稳定的基于核函数的结合空间信息模糊c均值图像分割算法.美国电子电器工程师协会系统控制处理汇刊.卷34,1907–1916,2004.(s.chenandd.zhang,“robustimagesegmentationusingfcmwithspatialconstraintsbasedonnewkernel-induceddistancemeasure,”ieeetrans.syst,man,cybern,vol.34,pp1907-1916,2004.)；这两种算法中fcm_s1算法对高斯噪声处理效果较好，但对椒盐噪声鲁棒性差，二fcm_s2算法虽对椒盐噪声处理效果较好，但却对高斯噪声鲁棒性差，因而这两种算法不能做到对多种类型噪声的鲁棒普适性。蔡维玲等人结合图像的空间信息和灰度信息构造了线性加权和图像，提出了快速生成fcm算法，参见：蔡维玲,陈松灿,张道强.一种用于图像分割的引入局部信息的快速鲁棒的模糊c均值聚类算法.模式识别.卷40,825-838,2007.(w.cai,s.chen,andd.zhang,“fastandrobustfuzzyc-meansclusteringalgorithmsincorporatinglocalinformationforimagesegmentation,”patternrecognit.,vol.40,no.3,pp.825-838,mar.2007.)，该方法在对高斯噪声有较强鲁棒性的同时对椒盐噪声也具有较好的鲁棒性；但是以上算法都没有考虑数据更多的模糊性，查尔热等人进一步发现了采用直觉模糊集理论可以考虑数据更多的模糊性，对数据的分类更加精确，提出了基于直觉模糊数据的模糊聚类方法，参见：一种新颖的直觉模糊c均值聚类算法及其在医学图像中的应用.应用软计算(t,chaira.“anovelintuitionisticfuzzycmeansclusteringalgorithmanditsapplicationtomedicalimages.”appl.softcomput.11(2):1711-1717,2011.)；由于基于直觉模糊数据的模糊聚类方法对噪声也比较敏感，韦尔马等人进一步将局部空间信息引入直觉模糊c均值算法中，参见：韦尔马,阿格拉沃尔,沙兰.针对脑部图像分割的结合局部空间信息的改进直觉模糊c均值算法.应用软计算.543-557,2016.(h.verma,r.k.agrawal,a.sharan,“animprovedintuitionisticfuzzyc-meansclusteringalgorithmincorporatinglocalinformationforbrainimagesegmentation,”appl.softcomput.,543–557,2016)。

虽然上述改进方法一定程度上优化了模糊聚类算法在抗噪性能上的表现，但依然存在对噪声鲁棒性不够强，对聚类中心的初始值很敏感，不能自适应分析图像聚类数目等不足。

技术实现要素：

本发明的目的在于针对上述现有技术的不足，提供一种基于直觉模糊c均值聚类的图像分割方法，以增强对噪声的鲁棒性，提高分割精确度，减小对聚类中心初始值的敏感性，自适应地确定图像聚类数目。

为实现上述目的，本发明的技术包括如下步骤：

(1)输入待分割图像，若待分割图像为彩色图像，先将其转换为灰度图像；

(2)初始化参数：设最大迭代次数t＝100，停止阈值ε＝10^-5，模糊加权指数m＝2，邻域窗半径ω＝3，初始迭代次数t＝1，聚类数目默认初始值y＝2；

(3)使用基于开关均值策略构造对噪声鲁棒的直觉模糊集其中，1≤i≤n，n表示像素总数目；

(4)针对输入图像的灰度值分布，使用最小二乘法进行曲线拟合，选择拟合曲线的峰值点所对应的灰度级作为初始化的聚类中心值，并设估计聚类数的范围为{2，…，ymax}，聚类中心初始值范围集合c为其中，cy表示第y个聚类中心的初始值，1≤y≤ymax，ymax表示聚类数目的最大值，2≤ymax≤n，t表示转置；

(5)迭代计算线性加权函数系数hir、直觉模糊距离和直觉模糊目标函数j：

(5a)利用像素的位置信息和灰度信息构造线性加权函数系数hir：

其中，(ai，bi)和(ar，br)分别表示中心像素i和邻域像素r的空间坐标值，xi和xr分别表示中心像素i和邻域像素r的灰度值，uki表示xi对于第k个聚类中心ck的隶属度，1≤k≤y，1≤r≤ω²；

(5b)计算直觉模糊集下像素xi到聚类中心ck的直觉模糊距离其中，和分别表示像素xi与聚类中心ck的直觉模糊集；

(5c)结合(5a)得到的线性加权函数系数hir与(5b)得到的直觉模糊距离计算直觉模糊目标函数j，并求解直觉模糊目标函数j中的隶属度函数uki和聚类中心ck的直觉模糊集

(6)判断迭代停止条件||c^(t+1)-c^(t)||＜ε或运算迭代次数t＞t是否成立：

如果不成立，则令t＝t+1，返回(5)；

如果成立，输出隶属度矩阵u并判断y＜ymax是否成立，若成立，则令y＝y+1，取聚类中心初始值范围集合c的前y个数值，令聚类中心初始值集合为令初始迭代次数t＝1，返回(5)，若y＜ymax不成立，执行(7)，其中，为0聚类中心数值集合，表示第q个聚类中心迭代第t次的直觉模糊集，表示第p个聚类中心的初始值，1≤q≤y，1≤p≤y，u＝{uki}；

(7)从(6)输入(ymax-1)个隶属度矩阵u，并使用分级距离指数评价指标gd评价这些隶属度矩阵；设gd数值最大的隶属度矩阵为umax，将umax对应的聚类数目yfin作为输入图像的最终聚类数目，其中，1≤yfin≤ymax；

(8)对(7)输出的隶属度矩阵umax做像素分类，并进行检错分类，输出分割图像：

(8a)根据最大隶属度原则和最终聚类数目yfin，对(7)输出的隶属度矩阵umax的像素进行分类，输出标签矩阵ul，其中，ul内仅含类别数值{1，2，…，yw，…，yfin}，1≤yw≤yfin；

(8b)设中心像素i的标签类别值为li，在中心像素i的5乘5邻域内，判断标签类别值li的像素数量nl是否满足nl≤25/(yfin+1)的条件：

如果满足，则判定中心像素i为错分像素，并将错分像素的标签类别值li修正为邻域内同种标签类别值数量最多的标签类别值，执行(8c)；

如果不满足，则判定中心像素i为正确分配像素，执行(8c)；

(8c)判断标签矩阵ul中的像素是否全部遍历，如果是，输出最终标签矩阵uf，执行(8d)，如果不是，返回(8b)继续遍历；

(8d)对最终标签矩阵uf中标签类别值为yw的像素赋予不同对应的灰度值xw，并将赋予不同灰度值的矩阵作为分割图像输出。

与现有技术相比，本发明具有以下有益的技术效果：

第一，本发明使用基于开关均值策略构造对噪声鲁棒的直觉模糊集，利用像素的位置信息和灰度信息构造线性加权函数系数并引入到直觉模糊目标函数中，解决了聚类算法对噪声敏感的问题。

第二，本发明使用最小二乘法对图像的灰度值分布进行曲线拟合，并选择拟合曲线所有的峰值点所对应的灰度级作为聚类中心初始值范围，根据分级距离指数评价指标评价隶属度矩阵，选择适合输入图像的聚类数目，实现图像聚类数目的自适应地确定，并克服了传统聚类算法容易陷入局部最优的缺点。

第三，本发明使用检错策略筛选错分像素，并对错分元素进行正确分配，提高了分割精确度，使得分割效果更为理想。

附图说明

图1为本发明的实现流程图；

图2为用本发明和现有方法对berkeley图像数据库中的#3096图像进行仿真分割的结果对比图；

图3为用本发明与现有方法对berkeley图像数据库中的#238011图像进行仿真分割的结果对比图。

具体实施方式

以下结合附图对本发明的实施例及效果作进一步描述：

参见图1，本发明的实现步骤如下：

步骤1：输入待分割图像。

输入待分割图像，若待分割图像为彩色图像，先将其转换为灰度图像。

步骤2：设置初始参数值。

设最大迭代次数t＝100，停止阈值ε＝10^-5，模糊加权指数m＝2，邻域窗半径ω＝3，初始迭代次数t＝1，聚类数目默认初始值y＝2。

步骤3：构造对噪声鲁棒的直觉模糊集

现有技术构造直觉模糊集的方法包括有ifcm算法和iifcm算法，其均使用yager算子来构造直觉模糊集。

本实例使用基于开关均值策略修正输入图像的像素灰度值，并通过yager算子来构造对噪声鲁棒的直觉模糊集，其中，开关均值策略是根据像素i邻域内的像素灰度值分布情况做灰度值修正，其实现如下：

3.1)首先判断在中心像素i的3乘3邻域内的像素灰度值是否不全为0或者255：

如果成立，则表示邻域内的像素灰度值除0或者255外还有别的灰度值，先将邻域内灰度值为0和255的像素取出，再取邻域内其余像素的灰度均值，并设为中心像素的修正数值，再执行3.2)；

如果不成立，则不修改中心像素的灰度值，直接执行3.2)；

3.2)判断输入图像中的像素是否全部遍历：如果是，输出灰度值修正后的图像，执行3.3)，如果不是，返回3.1)继续遍历；

3.3)将3.2)输出的修正图像的像素灰度值计算为直觉模糊集

3.3.1)设μ(xi)表示像素灰度值xi的隶属度，v(xi)和π(xi)分别是xi使用yager算子构造的非隶属度和犹豫度，μ(xi)、v(xi)和π(xi)的计算公式分别如下：

π(xi)＝1-μ(xi)-v(xi)，

其中，0≤μ(xi)≤1，0≤v(xi)≤1，0≤μ(xi)+v(xi)≤1，α＝0.8表示yager算子系数；

3.3.2)根据3.3.1)计算得到的μ(xi)、v(xi)和π(xi)构造直觉模糊集

步骤4：针对输入图像的灰度值分布，使用最小二乘法进行曲线拟合，选择拟合曲线的峰值点所对应的灰度级作为初始化的聚类中心值，估计最大聚类数目及聚类中心初始值范围。

现有选取聚类中心初始值的技术有直接选取图像灰度直方图的峰值点作为初始聚类中心。

本实例先对图像灰度值分布使用最小二乘法进行曲线拟合，然后在平滑的拟合曲线选择峰值点，其实现如下：

4.1)将灰度图像设为256个灰度级，灰度级范围为[0，255]，并绘制灰度图像的灰度直方图；

4.2)统计灰度直方图上256个灰度级所对应的像素数目，设为灰度值分布频率f(l)，并使用七阶最小二乘法对f(l)做曲线拟合，将灰度值数据点拟合为灰度值曲线，其中，0≤l≤255；

4.3)求f(l)的一阶导数f′(l)，并比较f′(l)与f(l)前后灰度级的一阶导数f′(l-1)和f′(l+1)；

4.4)判断灰度级l是否在l＝0或者l＝255的条件下：

若在l＝0或者l＝255的条件下且l＝0，则再判断f′(l)＜0是否成立，如果成立，将灰度级l＝0作为峰值点输出，再执行4.5)；

若在l＝0或者l＝255的条件下且l＝255，则再判断f′(l)＞0是否成立，如果成立，将灰度级l＝255作为峰值点输出，再执行4.5)；

若不是在l＝0或者l＝255的条件下，即表示在0＜l＜255的条件下，则判断f′(l)＝0且f′(l-1)＞0和f′(l+1)＜0是否成立，如果成立，将灰度级l作为峰值点输出，再执行4.5)；

4.5)判断256个灰度级是否全部遍历，如果是，输出峰值点，执行4.6)，如果不是，返回4.4)继续遍历；

4.6)对4.5)输出的峰值点按照灰度分布频率值从大到小的顺序排序，统计峰值点的数目为ymax，设聚类数范围为2至ymax个，并统计从大到小排序的峰值点对应的灰度级l，得到灰度级数值范围为将灰度级对应数值统计为聚类中心初始值，即设聚类中心初始值范围集合c为其中：

ly表示在峰值点灰度分布频率值按照从大到小顺序中，第y个峰值点灰度分布频率值所对应的灰度级，cy表示第y个聚类中心的初始值，1≤y≤ymax，ymax表示聚类数目的最大值，2≤ymax≤n，t表示转置。

步骤5：迭代计算线性加权函数系数hir、直觉模糊距离和直觉模糊目标函数j。

现有的fcm_s1算法和fcm_s2算法将图像的局部空间信息引入到其目标函数中，其中，fcm_s1算法考虑的是图像的均值滤波特征，fcm_s2算法考虑的是图像的中值滤波特征。

本实例先利用像素的位置信息和灰度信息构造线性加权函数系数hir，然后将包含图像局部空间信息的hir引入到直觉模糊目标函数j，其实现如下：

5.1)利用像素的位置信息和灰度信息构造线性加权函数系数hir：

其中，(ai,bi)和(ar,br)分别表示中心像素i和邻域像素r的空间坐标值，xi和xr分别表示中心像素i和邻域像素r的灰度值，uki表示xi对于第k个聚类中心ck的隶属度，1≤k≤y，1≤r≤ω²；

5.2)计算直觉模糊集下像素灰度值xi到聚类中心ck的直觉模糊距离

5.2.1)设μ(ck)表示聚类中心ck的隶属度，v(ck)和π(ck)分别表示聚类中心ck使用yager算子构造的非隶属度和犹豫度，μ(ck)、v(ck)和π(ck)的计算公式分别如下：

π(ck)＝1-μ(ck)-v(ck)，

其中，0≤μ(ck)≤1，0≤v(ck)≤1，0≤μ(ck)+v(ck)≤1；

5.2.2)根据(5.2.1)得到的结果计算直觉模糊距离

其中，μ(xi)、v(xi)和π(xi)分别表示像素灰度值xi的隶属度、非隶属度和犹豫度；

5.3)结合(5.1)得到的线性加权函数系数hir与(5.2)得到的直觉模糊距离计算直觉模糊目标函数j：

其中，uki是像素xi对聚类中心ck的隶属度函数，m表示模糊加权指数；

5.4)使用拉格朗日乘数法求解直觉模糊目标函数j，并得到隶属度函数uki和聚类中心ck的直觉模糊集

5.4.1)在直觉模糊目标函数j中引入拉格朗日乘子λi，构建具有拉格朗日乘子的目标函数l：

5.4.2)分别计算具有拉格朗日乘子的目标函数l对于隶属度函数uki、聚类中心ck的直觉模糊集和拉格朗日乘子λi的偏导数和

5.4.3)在偏导数的条件下，通过如下求得隶属度函数uki与聚类中心ck的直觉模糊集

其中，k＝1,2,…,y，i＝1,2,…,n，r＝1,2,…,ω²。

步骤6：判断是否满足迭代停止条件。

如果满足迭代停止条件||c^(t+1)-c^(t)||＜ε或运算迭代次数t＞t，则输出隶属度矩阵u，并执行7),其中，为聚类中心数值集合，表示第q个聚类中心迭代第t次的直觉模糊集；

如果不满足迭代停止条件，则令t＝t+1，返回5.1)。

步骤7：判断当前聚类数目y是否满足条件。

如果满足当前聚类数目y小于最大聚类数目ymax，即y＜ymax，则令y＝y+1，然后取聚类中心初始值范围集合c的前y个数值，输出到聚类中心初始值集合，并令聚类中心初始值集合c⁽¹⁾为令初始迭代次数t＝1，返回5.1)，其中，表示第p个聚类中心的初始值，1≤p≤y；

如果不满足y＜ymax，则执行(8)。

步骤8：使用分级距离指数评价指标gd评价隶属度矩阵。

现有对隶属度矩阵的评价指标有划分系数vpc和划分熵vpe，其均能反映隶属度矩阵的模糊程度，其中vpc数值越大，隶属度矩阵的模糊性越小，vpe数值越小，像素分类越准确。

本实例使用分级距离指数评价指标gd评价隶属度矩阵，即从7)中输入(ymax-1)个隶属度矩阵u，使用分级距离指数评价指标gd评价这些隶属度矩阵，输出gd数值最大的隶属度矩阵，并设gd数值最大的隶属度矩阵为umax，将umax对应的聚类数目yfin作为输入图像的最终聚类数目，其中：分级距离指数评价指标gd公式如下：

其中，ui1表示第i个像素最大的隶属度值，ui2表示第i个像素第二大的隶属度值，1≤yfin≤ymax。

步骤9：对隶属度矩阵umax进行像素分类。

根据最大隶属度原则和最终聚类数目yfin，对8)输出的隶属度矩阵umax的像素进行分类，输出标签矩阵ul，其中，最大隶属度原则是根据vk＝argk{max(uki)}将像素i划分到隶属度最大的类别vk，ul内仅含类别数值{1，2，…，yw，…，yfin}，1≤yw≤yfin。

步骤10：对标签矩阵ul进行检错分类。

10.1)设9)输出的标签矩阵ul的中心像素i的标签类别值为li，在中心像素i的5乘5邻域内，判断标签类别值li的像素数量nl是否满足nl≤25/(yfin+1)的条件：

如果满足，则判定中心像素i为错分像素，并将错分像素的标签类别值lj修正为邻域内同种标签类别值数量最多的标签类别值，执行10.2)；

如果不满足，则判定中心像素i为正确分配像素，执行10.2)；

10.2)判断标签矩阵ul中的像素是否全部遍历，如果是，输出最终标签矩阵uf，执行11)，如果不是，返回10.1)继续遍历。

步骤11：对最终标签矩阵uf赋值并输出图像。

对10.2)输出的最终标签矩阵uf中标签类别值为yw的像素赋予不同对应的灰度值xw，并将赋予不同灰度值的矩阵作为分割图像输出，最终标签矩阵uf内仅含类别数值{1，2，…，yw，…，yfin}，赋值公式如下：

xw＝[255/yw]，

其中，[·]表示取内部的整数。

以下结合仿真实验，对本发明的技术效果作进一步说明；

1.仿真条件：

仿真实验在计算机inter(r)core(tm)i5-3210m2.50ghzcpu，8g内存，matlabr2014a软件环境下进行。

2.仿真内容：

仿真1，选取berkeley图像数据库中的编号为3096的图像，用本发明与现有fcm方法、fcm_s1方法、fcm_s2方法、fgfcm方法、ifcm方法和iifcm方法分别对其进行分割，结果如图2所示，其中：

2(a)是3096图像的原图；

2(b)是3096图像的标准分割图；

2(c)是3096图像的高斯含噪图像，噪声强度为0.006；

2(d)是用现有fcm方法对3096图像的高斯含噪图像的分割结果；

2(e)是用现有fcm_s1方法对3096图像的高斯含噪图像的分割结果；

2(f)是用现有fcm_s2方法对3096图像的高斯含噪图像的分割结果；

2(g)是用现有fgfcm方法对3096图像的高斯含噪图像的分割结果；

2(h)是用现有ifcm方法对3096图像的高斯含噪图像的分割结果；

2(i)是用现有iifcm方法对3096图像的高斯含噪图像的分割结果；

2(j)是用本发明对3096图像的高斯含噪图像的分割结果；

2(k)是3096图像的椒盐含噪图像，噪声强度为0.01；

2(l)是用现有fcm方法对3096图像的椒盐含噪图像的分割结果；

2(m)是用现有fcm_s1方法对3096图像的椒盐含噪图像的分割结果；

2(n)是用现有fcm_s2方法对3096图像的椒盐含噪图像的分割结果；

2(o)是用现有fgfcm方法对3096图像的椒盐含噪图像的分割结果；

2(p)是用现有ifcm方法对3096图像的椒盐含噪图像的分割结果；

2(q)是用现有iifcm方法对3096图像的椒盐含噪图像的分割结果；

2(r)是用本发明对3096图像的椒盐含噪图像的分割结果；

从图2可以看出，本发明在抑制噪声的基础上，可以将背景和目标分离开，图像聚类数目可以自适应确定，对聚类中心的初始值不敏感，所以本发明对含噪图像的分割效果优于现有fcm方法、fcm_s1方法、fcm_s2方法、fgfcm方法、ifcm方法和iifcm方法。

仿真2，选取berkeley图像数据库中的编号为238011的图像，用本发明与现有fcm方法、fcm_s1方法、fcm_s2方法、fgfcm方法、ifcm方法和iifcm方法分别对其进行分割，结果如图3所示，其中：

3(a)是238011图像的原图；

3(b)是238011图像的标准分割图；

3(c)是238011图像的高斯含噪图像，噪声强度为0.006；

3(d)是用现有fcm方法对238011图像的高斯含噪图像的分割结果；

3(e)是用现有fcm_s1方法对238011图像的高斯含噪图像的分割结果；

3(f)是用现有fcm_s2方法对238011图像的高斯含噪图像的分割结果；

3(g)是用现有fgfcm方法对238011图像的高斯含噪图像的分割结果；

3(h)是用现有ifcm方法对238011图像的高斯含噪图像的分割结果；

3(i)是用现有iifcm方法对238011图像的高斯含噪图像的分割结果；

3(j)是用本发明对238011图像的高斯含噪图像的分割结果；

3(k)是238011图像的椒盐含噪图像，噪声强度为0.01；

3(l)是用现有fcm方法对238011图像的椒盐含噪图像的分割结果；

3(m)是用现有fcm_s1方法对238011图像的椒盐含噪图像的分割结果；

3(n)是用现有fcm_s2方法对238011图像的椒盐含噪图像的分割结果；

3(o)是用现有fgfcm方法对238011图像的椒盐含噪图像的分割结果；

3(p)是用现有ifcm方法对238011图像的椒盐含噪图像的分割结果；

3(q)是用现有iifcm方法对238011图像的椒盐含噪图像的分割结果；

3(r)是用本发明对238011图像的椒盐含噪图像的分割结果；

从图3可以看出，本发明在抑制噪声的基础上，可以将背景和目标分离开，图像聚类数目可以自适应确定，对聚类中心的初始值不敏感，其对含噪图像的分割效果优于现有fcm方法、fcm_s1方法、fcm_s2方法、fgfcm方法、ifcm方法和iifcm方法。