CN106228545A - 一种佛像类唐卡一致性破损区域的自动分割方法 - Google Patents

Abstract

本发明提供了一种佛像类唐卡一致性破损区域的自动分割方法，首先对唐卡图像的头光区域进行垂直投影，利用一维函数对称性检测方法得到图像的对称轴，采用基于对称轴的分块分割方法得到初始分割结果；接着获取破损区域所覆盖的图像；然后利用Gabor变换提取纹理特征，结合Lab空间颜色特征，构造多尺度多特征集合，采用KNN分类得到二次分割结果；最后利用形态学运算对破损区域进行进一步细化，通过将小的破损区域去除，最终获得一致性破损区域的模板。本发明方法可以自动分割出佛像类唐卡图像中出现的各类较大尺度的线状和块状脱落区域等，而且分割速度快、效率高且精确度高，实现了唐卡图像破损区域的快速自动分割。

Description

一种佛像类唐卡一致性破损区域的自动分割方法

技术领域

本发明属于图像处理技术领域，涉及一种佛像类唐卡一致性破损区域的自动分割方法。

背景技术

所谓佛像类唐卡一致性破损区域的自动分割，就是利用一定算法首先提取佛像类唐卡图像的对称轴，采用基于对称轴的分块分割算法获得初始分割结果；接着利用Gabor变换提取初始分割的纹理特征，然后结合与人眼视觉一致的Lab颜色空间特征信息，共同组成多尺度多特征集合，采用KNN分类器实现破损区域的再次分割；最后利用形态学运算、去除小面积区域后得到最终破损区域的技术。唐卡图像破损区域的自动分割是一个实际应用问题，是实施唐卡破损修复的前提，其分割质量的好坏对于修复结果有着重要的影响，因而在文化遗产的数字化保护等方面有着重要的实际意义和理论价值。

破损唐卡图像主要存在着颜色剥落、龟裂、折痕，画布下陷、扭曲形变，虫蛀、鼠咬等损害现象，为了研究的方便并对比常见的破损形式，将其归结为三种类型：任意形状划痕(包括水平和垂直划痕)、斑点和大面积缺损。本发明主要解决的是大面积缺损的这种形式，这实质上属于图像分割的问题。图像分割，就是把图像分成各具特性的区域并提取出感兴趣区域的技术。早在20世纪60年代，国内外学者就对其进行了大量的研究，例如经典的算法有：阈值分割法、边缘检测法以及区域生长技术。但是，对于如何利用领域知识则是采用了交互式的方式。当然，随着图像分割技术的蓬勃发展，也涌现出了一些新的方法和思路。例如：结合模糊集、形状特征、变分、水平集和图论等理论进行图像分割。

然而由于图像分割本身的不确定性和复杂性，还没有一种算法能够解决所有的分割问题。关于图像破损区域自动分割方法的专门讨论仍然较少，破损块一般都由用户选择，如直接确定待修复或待移走的对象区域、从待修复区域的种子点开始用区域生长等方法获得破损块等，但这些方法都只能解决简单的破损，对复杂破损问题难以处理。目前来看，还没有一种统一的方法可以解决各种类型破损的图像分割。

发明内容

本发明的目的在于解决佛像类唐卡图像中出现的一致性破损区域如何进行自动分割的问题。本发明方法可以自动分割出佛像类唐卡图像中出现的各类较大尺度的线状和块状脱落区域等，而且分割速度快、效率高、精确度高。

为实现上述目的，本发明采取如下技术方案：

一种佛像类唐卡一致性破损区域的自动分割方法，首先对唐卡图像的头光区域进行垂直投影，利用一维函数对称性检测方法得到图像的对称轴，采用基于对称轴的分块分割方法得到初始分割结果；接着获取破损区域所覆盖的图像；然后利用Gabor变换提取纹理特征，结合Lab空间颜色特征，构造多尺度多特征集合，采用KNN分类得到二次分割结果；最后利用形态学运算对破损区域进行进一步细化，通过将小的破损区域去除，最终获得一致性破损区域的模板。

初始分割希望尽可能地包含破损区域的全部信息，由于破损区域与周围像素颜色相差不多，如果在整幅图像中只采用一个灰度阈值，则势必会丢失部分破损区域的信息。分块分割法自左向右、自上而下的对图像进行矩形分块，然后在各个子块内再进行阈值化分割，能够 最大程度地保留破损信息。

此外唐卡图像在构图时存在明显的对称性，尤其是佛像类。为此，本发明创新性的以对称轴为基准采用分块采样的策略进行初始分割。步骤概括如下：首先选定头部区域，接着分割出头光区域，然后在头光区域中检测出对称轴，最后基于对称轴进行分块分割得到初始分割结果。具体步骤如下：

a.选定头部区域

经过数字化得到的待分割的唐卡原始图像是彩色图像，记为F(x，y)，采取手工选定或计算机选定的方式得到图像中佛像的头部区域t(x，y)；

b.检测对称轴

利用基于划分的聚类算法k-means对头部区域t(x，y)进行分割，分割得到佛像类唐卡图像的头光区域，通过垂直投影的方法得到头光区域的一维函数，利用一维函数对称性检测算法得到头光区域的对称轴，这个对称轴就是整幅图像的对称轴。

b.1图像的灰度化处理

在进行破损区域的自动分割时使用的是灰度图像，将上述步骤a得到的头部区域首先进行灰度处理。设头部区域t(x，y)为一幅包含有M×N个像素点的二维图像，采用平均值法，即对每个像素点R、G、B三分量的亮度求其平均值，然后将求得的这个平均值作为灰度图像的灰度值；

b.2提取头光区域

在佛像类唐卡中，大多以主尊为中心，四周密密麻麻对称地画出 大量繁杂的人物、建筑、山水和花卉等，而在主尊，如佛陀、菩萨、圣者的头顶和身后，都笼罩着一道辉煌的灵光圈，分别称为头光和背光。虽然不同年代的唐卡画像在处理头光、背光图案样式上略有不同，但是头光大多表现为颜色单一、稳定的环形区域。本发明采用基于划分的聚类算法k-means分割头光区域。具体步骤如下：

b.2.1分割头部区域

第一步：将步骤a得到的头部区域t(x，y)的所有像素点构成集合T＝{t_i|i＝1，2，...total},T中的样本点用三个属性(x_m，y_m，f_m)进行描述，其中x_m代表像素点的水平坐标的值，y_m代表像素点垂直坐标的值，f_m代表像素点(x_m，y_m)的灰度值，随机选择k个元素，作为k个簇的初始中心点A＝{μ_j|j＝1，2，...，k}；

第二步：分别计算剩下的元素到k个簇中心的欧式距离，其中元素到簇中心点的欧式距离可以通过下述公式计算：

式中，为头部区域像素点i的灰度值，表示本次迭代中簇j中心点的灰度值，s即为图像中像素点与簇中心像素点灰度值的差异；

将T中剩余元素分别划归到距离最近的簇中，得到k个不同的元素聚类，即：

若||t_fi-μ_fj||＜||t_fi-μ_fl||，则t_i∈μ_j，

其中j，l＝1，2，...，k，j≠l

第三步：计算每个聚类中所有点的灰度值平均值，并将这个平均值作为新的聚类中心，新的聚类中心表达式为：

其中N_j是类W_j中样本的数目，t_i为像素点i的灰度值；

反复执行步骤(2)、(3)，直到聚类中心点不再变化为止，最终得到头部区域的分割结果。

b.2.2获得头光区域

鉴于在头部区域中头光部分的颜色稳定，像素点数目众多，在得到的这k个结果类中，分别统计每个聚类的像素点数目，然后选择像素点数目最多的区域，这个包含像素点最多的颜色稳定的区域就是佛教类唐卡图像的头光区域。

b.3获取对称轴

b.3.1将头光区域所有像素点叠加投影得到S(x)

将头光区域所有像素点沿着垂直方向叠加投影，即：

其中，n表示图像的高度，t(x，y)为待处理的头部区域，mask(x，y)为头光区域的二值模板，若某一像素点位于头光区域，则该点的mask(x，y)值为1，若在头光区域外则mask(x，y)值为0；

b.3.2求S(x)的奇函数与偶函数

式(3)所得的垂直方向投影函数是一维函数，假设对称轴的坐标为x_s(w/2≤x_s≤W-w/2)，其中对称区域的宽度为w，W是式(3)所得一维函数的宽度；则S(u)是以x_s为坐标原点的函数，S(u)＝S(x_s-u)。

任意函数可以表示为奇函数与偶函数的和，计算S(u)的奇函数分量odd(u)和偶函数分量even(u)。

b.3.3对称轴l求取

由于唐卡图像具有对称性，偶函数所占比重大，因此取偶函数分量even(u)求取对称轴l。使用式(6)对偶函数分量even(u)进行归一化得：

式中，W为一维函数S(u)的宽度；

计算对称性测度，对称性测度symmetry可以通过下式计算：

式(7)中，当symmetry＝1时，说明函数even(u)完全对称；symmetry＝-1时，说明函数even(u)完全不对称。实验中，取0.6＜symmetry≤1时计算得到x_s的值，x_s的坐标值即为头光区域的对称轴的坐标值；

c.基于对称性的分块采样分割

c.1唐卡图像分块

将原始的待分割图像F(x，y)进行图像的灰度化处理，得到灰度化图像f(x，y)，灰度化处理过程同步骤b.1；唐卡图像一般为左右对称，故以对称轴l为轴，可以将灰度图像分为左半部分f_左和右半 部分f_右，即f＝f_左+f_右；

然后以对称轴l为基准从图像顶端开始分别向左右边界延伸，进行分块采样。步骤如下：

c.1.1选择左边区域图像块S_左

从f(x，y)的顶端开始，以对称轴为基准向左依次选取图像块S_左1、S_左2...S_左N，图像块的大小均为的像素块；

c.1.2选择右边区域图像块S_右

从f(x，y)的顶端开始，以对称轴l为基准向右依次选取图像块S_右1、S_右2...S_右N，图像块的大小均为的像素块；

c.1.3合并得到分块S

将关于对称轴l相对称的S_左i和S_右i进行合并可以得到大小为m×n的分块S_m×n＝[S_左iS_右i]，i＝1，2，...N；为了有效地处理，边界可采用补0来扩充。具体参见图2所示，图中浅色部分表示左边区域图像块S_左i，深色部分表示右边区域图像块S_右i。实验中，m和n都选取8。这里分块大小决定了分割的粒度，块大小较小时包含了较多的冗余信息，块大小太大则丢失较小的破损区域。

以对称轴为基准，从上到下，将以上步骤c.1.1至步骤c.1.3重复，直到将唐卡图像f(x，y)表示成M个大小相同的图像块S_m×n的集合。c.2对分块采样后的子图像块进行分割

经过步骤c.1，图像表示为M个大小相同的图像块S_m×n的集合，即f(x，y)＝{S_i|i＝1，2，...M}，然后在每个子图像块s_m×n中分别进行图像分割，本发明选用阈值法进行子图像块的分割。

最大类间方差法(OTSU)是一种自适应的阈值确定算法，它按图像的灰度特性，将图像分成背景和目标两个部分。如果背景和目标间的类间方差较大，则表明这两部分的差别越大。当部分目标或背景被错分时，就会导致这两部分的差别变小。因此，若某个分割阈值使得类间方差最大，则意味着该分割错分的概率最小，分割的效果最为理想。具体步骤如下：

c.2.1计算灰度级i在图像中出现的概率

设每个子图像块S_m×n构成的图像包含L个灰度级(0,1…,L-1)，灰度值为i的像素点数为n_i，图像总的像素点数为N＝n₀+n₁+...n_L-1。灰度级i在图像中出现的概率为：

c.2.2计算图像的目标C₁和背景C₂的概率

假定通过灰度级阈值t将每一个子图像块的像素点分成目标C₁和背景C₂两类，其中C₁表示所有灰度级在[0，1，...，t]范围内的像素点的集合，C₂表示所有灰度级在(t+1，...，L-1]范围内的像素点的集合。则C₁和C₂两个类的概率可表示为：

c.2.3计算图像的目标C₁和背景C₂两个类的均值和方差

两个类的均值分别为：

两个类的方差分别为：

c.2.4计算类内方差、类间方差和总的方差

类内方差类间方差和总的方差的计算如下：

其中，故

b.2.5采用最大类间方差得到分块分割阈值t

采用循环的方式，灰度级阈值t从0开始直到取得子图像块灰度级的最大值L-1，利用式(10)计算类间方差，当类间方差取得最大值时所对应的t就是分块分割的灰度级阈值；也可利用式(9)计算类内方差，当类内方差取得最小值时对应的t也即是分块分割的灰度级阈值；

b.2.6利用阈值t对子图像块S_m×n进行分割

对于子图像块S_m×n中的所有像素点，若其灰度值大于阈值t，则标记为前景即目标，并用1表示；若其灰度值小于或等于阈值t，则标记为背景，并用0表示。

d.获取覆盖区域

在完成对所有子图像块分割之后，将每个子图像块分割的结果合并得到破损唐卡图像的初始分割结果，这里记为h(x，y)，这时候破损区域的全部信息都已经包含在背景部分或者说物体部分中了，然而在一幅画作中，真实的破损形式众多，特征千差万别，欲获得破损区域的模板还需要对分块分割结果做进一步的处理。为此，根据步骤c得到的分割结果，获得基于对称性分块分割覆盖的灰度图像f_分割(x，y)和彩色图像f_分割(x，y)′。

f_分割(x，y)＝f(x，y)*h(x，y)

f_分割(x，y)′＝F(x，y)*h(x，y)

其中，如前所表示：F(x，y)为待分割的原始彩色图像，f(x，y)是其灰度化处理得到的灰度图像。

色彩美是唐卡绘画艺术的一大特色，因此，颜色特征信息应该作为识别破损区域的特征之一。线条美、构图美是唐卡绘画的另一大艺术特色因此，纹理特征信息也应该作为提高识别破损区域的特征之一。为此，本发明结合颜色、纹理特征构造多特征融合的特征集合，通过像素点分类的方法对一致性破损区域进行自动分割，得到二次分割结果。具体步骤如下：

e.图像特征提取

e.1提取图像Gabor纹理特征

利用Gabor滤波器提取图像的纹理特征，Gabor变换是一种有效的多尺度分析工具，具备分析图像局部细微变化的能力，可以很好的表征唐卡图像不同尺度下的纹理特征。图像Gabor特征可通过将输入图像与Gabor小波滤波器组卷积得到，即：

m(x，y)＝f_分割(x，y)*g(x，y) (11)

这里f_分割(x，y)是输入的步骤d得到的覆盖区域灰度图像，g(x，y)是Gabor滤波函数，定义为：

式中，λ是正弦函数波长，θ是Gabor核函数的方向，是相位偏移，σ为高斯函数的标准差，γ为空间的方面比，且γ的范围值在0.23～0.92之间，x,y是给定位置的图像坐标，x',y'是旋转之后的坐标。通过修改参数，便可以得到一组不同尺度和方向的纹理特征，并作为最终的图像纹理特征feature_texture。

e.2提取颜色特征

对于颜色特征的提取，考虑到像素点的颜色差距应该与人眼视觉感觉到的差别相一致，故首先将待分割破损唐卡彩色图像f_分割(x，y)′从RGB空间转换为CIE Lab颜色空间，然后再提取颜色特征进行处理。具体步骤如下：

e.2.1色彩空间转换

RGB空间转换为CIE Lab颜色空间具体转换过程为

则CIE Lab空间的亮度L分量，以及有关色彩的a，b分量分别表示为：

其中(X₀,Y₀,Z₀)是标准白色的刺激值。

e.2.2获得颜色特征

通过颜色特征转化得到颜色特征向量feature_color＝{L，a，b}。

e.3提取纹理和颜色特征

将步骤e.1获取的纹理特征feature_texture和步骤e.2获取的颜色特征feature_color最终构成融合纹理特征和颜色特征的特征集合feature＝{feature_texture,feature_color}。

f.多尺度多特征融合得到二次分割结果

f.1归一化特征

前面提取的纹理特征大小和取值范围都不一样，将它们称之为异质特征。异质特征需要归一化处理；本发明选取高斯归一化方法。

高斯归一化的方法概述如下：假设有M个样本，每个样本的特征被表示为：(u_Fi1,u_Fi2,u_Fi3,u_Fi4,u_Fi5,σ_Fi1,σ_Fi2,,σ_Fi3,σ_Fi4,σ_Fi5)，其中1≤i≤M；现在以其中的任一特征分量为例做归一化处理；设u_F1j、 u_F2j……u_FMj表示M个样本的第j维特征分量，u_Fj和σ_Fj分别表示第j维特征分量的期望和方差；高斯归一化公式是：

为了使更多u_Fij落到(-1,1)内，此处高斯归一化公式分母取3σ_Fj；

则本发明特征归一化按式(18)所示的高斯归一化公式进行：

式中，1≤i≤M，M为特征集合feature所包含的特征数量，u_Fij表示第i个特征的第j维分量，u_Fj和σ_Fj分别表示第j维分量的期望和方差；

它能够将特征u_Fij归一化到(0,1)区间，这样就使得各个特征分量所起的作用相同，不会因为某些取值较大的特征而淹没其它特征的作用。

f.2基于多尺度多特征融合的分类

采用基于划分的聚类算法k-means分割破损区域，且步骤f.2所采用的聚类算法k-means与步骤b.2.1相同；

破损区域的样本集合为步骤f.1得到的经过归一化后的特征集合feature＝(x₁，x₂，...，x_n)，且每个x_i都是一个特征向量，需要将这个特征集合分成K个不同的类S＝(s₁，s₂，...s_k)，虽然这里只需要判断哪些点是破损点，哪些点是非破损点，也就是说可以设置k为2，但是 由于唐卡图像的破损原因较多，存在诸如质地变脆、颜色脱落、霉斑、颜料龟裂等破损形式，因此不能简单的分为两类，故取k＝5。具体步骤如下：

第一步：随机选择k个元素，作为K个簇的初始中心点A＝{μ_i|i＝1，2，...，k}；

第二步：分别计算剩下的元素到K个簇中心的欧式距离，其中元素到簇中心点的欧式距离可以通过下述公式计算：

s＝（x_j-μ_i)²

其中，x_j和μ_i均为feature中的特征向量；

将这些元素分别划归到距离最近的簇中，即：

若||x_i-μ_j||＜||x_i-μ_l||，则x_i∈μ_j j，l＝1，2，...，K，j≠l

第三步：计算每个聚类中所有特征向量的质心，并将其作为新的聚类中心，新的聚类中心表达式为：

其中N_j是类s_j中样本的数目；

反复执行步骤(2)、(3)，直到聚类中心点不再变化为止，最终得到分类的结果。

g.形态学运算

g.1开运算

步骤f.2二次分割得到的结果存在大量小的破损区域、破损区域之间的粘连，采用先腐蚀后膨胀的过程即开运算进行处理，用来消除小物体、在纤细点处分离破损区域、平滑较大破损区域的边界，且同 时不明显改变其面积。

腐蚀操作为：

膨胀操作为：

其中B(x)为结构元素，结构元素在形态变换中的作用相当于信号处理中的“滤波窗口”，这里采用3×3的全1模板。E为二次分割得到的结果图像。用B(x)对E进行腐蚀的结果就是把结构元素B平移后使B包含于E的所有点构成的集合。用B(x)对E进行膨胀的结果就是把结构元素B平移后使B与E的交集非空的点构成的集合。

g.2闭运算

先膨胀后腐蚀的过程称为闭运算。用来填充破损区域内细小空洞、连接邻近破损区域、平滑其边界的同时并不明显改变其面积。且该步骤的膨胀和腐蚀过程同步骤g.1。

h.去除小面积区域

统计步骤g.2得到的不同破损区域中像素点的个数m_i，将m_i＜ε的区域去除，这里ε是阈值，可以通过反复实验的方式得到最佳值，本发明中，ε的取值为30～50之间。最后得到的结果就是唐卡图像破损区域模板。

本发明的有益效果在于：可以自动分割出佛像类唐卡图像中出现的各类较大尺度的线状和块状脱落区域等，而且分割速度快、效率高且精确度高，实现了唐卡图像破损区域的快速自动分割；而且该方法适用于各类型唐卡图像的破损问题，使得唐卡图像破损区域的分割获得了统一的处理方法；另外，本发明方法分割得到的破损区域界线明晰、分割结果准确。

附图说明

图1为本发明方法的流程图；

图2为本发明步骤c.1分块分割的示意图。

具体实施方式

一种佛像类唐卡一致性破损区域的自动分割方法，包括如下步骤：

a.选定头部区域

经过数字化得到的待分割的唐卡原始图像是彩色图像，记为F(x，y)，采取手工选定或计算机选定的方式得到图像中佛像的头部区域t(x，y)；其中，采用计算机选定的过程如下：

a.1唐卡图像归一化

通过对唐卡的数字图像实验统计，得出适合唐卡图像处理与识别的图像大小为：400×520，所以对破损唐卡图像都将大小归一化为400×520；

a.2头部区域的半自动选择

经过对唐卡的数字图像的实验统计，发现头部区域一般在归一化图像的中上部。故选择从点(120，80)到(240，200)范围内的像素点作为唐卡图像的头部区域t(x，y)。

b.检测对称轴

b.1图像的灰度化处理

将上述步骤a得到的头部区域首先进行灰度处理。设头部区域t(x，y)为一幅包含有M×N个像素点的二维图像，采用平均值法，即对每个像素点R、G、B三分量的亮度求其平均值，然后将求得的这个平 均值作为灰度图像的灰度值；

b.2提取头光区域

采用基于划分的聚类算法k-means分割头光区域。具体步骤如下：

b.2.1分割头部区域

第一步：将步骤a得到的头部区域t(x，y)的所有像素点构成集合T＝{t_i|i＝1，2，...total},T中的样本点用三个属性(x_m，y_m，f_m)进行描述，其中x_m代表像素点的水平坐标的值，y_m代表像素点垂直坐标的值，f_m代表像素点(x_m，y_m)的灰度值，随机选择k个元素，作为k个簇的初始中心点A＝{μ_j|j＝1，2，...，k}；试验中，k＝4；

第二步：分别计算剩下的元素到k个簇中心的欧式距离，其中元素到簇中心点的欧式距离可以通过下述公式计算：

式中，为头部区域像素点i的灰度值，表示本次迭代中簇j中心点的灰度值，s即为图像中像素点与簇中心像素点灰度值的差异；将T中剩余元素分别划归到距离最近的簇中，得到k个不同的元素聚类，即：

若||t_fi-μ_fj||＜||t_fi-μ_fl||，则t_i∈μ_j，

其中j，l＝1，2，...，k，j≠l

第三步：计算每个聚类中所有点的灰度值平均值，并将这个平均值作为新的聚类中心，新的聚类中心表达式为：

其中N_j是类W_j中样本的数目，t_i为像素点i的灰度值；

反复执行步骤(2)、(3)，直到聚类中心点不再变化为止，最终得到头部区域的分割结果。

b.2.2获得头光区域

鉴于在头部区域中头光部分的颜色稳定，像素点数目众多，在得到的这k个结果类中，分别统计每个聚类的像素点数目，然后选择像素点数目最多的区域，这个包含像素点最多的颜色稳定的区域就是佛教类唐卡图像的头光区域。

b.3获取对称轴

b.3.1将头光区域所有像素点叠加投影得到S(x)

将头光区域所有像素点沿着垂直方向叠加投影，即：

其中，n表示图像的高度，t(x，y)为待处理的头部区域，mask(x，y)为头光区域的二值模板，若某一像素点位于头光区域，则该点的mask(x，y)值为1，若在头光区域外则mask(x，y)值为0；

b.3.2求S(x)的奇函数与偶函数

式(3)所得的垂直方向投影函数是一维函数，假设对称轴的坐标为x_s(w/2≤x_s≤W-w/2)，其中对称区域的宽度为w，W是式(3)所得一维函数的宽度；则S(u)是以x_s为坐标原点的函数，S(u)＝S(x_s-u)。

任意函数可以表示为奇函数与偶函数的和，计算S(u)的奇函数分量odd(u)和偶函数分量even(u)。

b.3.3对称轴求取

使用式(6)对偶函数分量even(u)进行归一化得：

式中，W为一维函数S(u)的宽度；

计算对称性测度，对称性测度symmetry可以通过下式计算：

实验中，取0.6＜symmetry≤1时计算得到x_s的值，x_s的坐标值即为头光区域的对称轴的坐标值；

c.基于对称性的分块采样分割

c.1唐卡图像分块

将原始的待分割图像F(x，y)进行图像的灰度化处理，得到灰度化图像f(x，y)，灰度化处理过程同步骤b.1；唐卡图像一般为左右对称，故以对称轴为轴，可以将灰度图像分为左半部分f_左和右半部分f_右，即f＝f_左+f_右；

然后以对称轴为基准从图像顶端开始分别向左右边界延伸，进行分块采样。步骤如下：

c.1.1选择左边区域图像块S_左

从f(x，y)的顶端开始，以对称轴为基准向左依次选取图像块S_左1、S_左2...S_左N，图像块的大小均为的像素块；

c.1.2选择右边区域图像块S_右

从f(x，y)的顶端开始，以对称轴为基准向右依次选取图像块S_右1、S_右2...S_右N，图像块的大小均为的像素块；

c.1.3合并得到分块S

将关于对称轴相对称的S_左i和S_右i进行合并可以得到大小为m×n的分块S_m×n＝[S_左i S_右i]，i＝1，2，...N；为了有效地处理，边界可采用补0来扩充。具体参见图2所示，图中浅色部分表示左边区域图像块S_左i，深色部分表示右边区域图像块S_右i。实验中，m和n都选取8。

以对称轴为基准，从上到下，将以上步骤c.1.1至步骤c.1.3重复，直到将唐卡图像f(x，y)表示成M个大小相同的图像块S_m×n的集合。c.2对分块采样后的子图像块进行分割

经过步骤c.1，图像表示为M个大小相同的图像块S_m×n的集合，即f(x，y)＝{S_i|i＝1，2，...M}，然后在每个子图像块s_m×n中分别进行图像分割，本发明选用阈值法进行子图像块的分割。

具体步骤如下：

c.2.1计算灰度级i在图像中出现的概率

设每个子图像块S_m×n构成的图像包含L个灰度级(0,1…,L-1)，灰度值为i的像素点数为n_i，图像总的像素点数为N＝n₀+n₁+...n_L-1。灰度级i在图像中出现的概率为：

c.2.2计算图像的目标C₁和背景C₂的概率

假定通过灰度级阈值t将每一个子图像块的像素点分成目标C₁ 和背景C₂两类，其中C₁表示所有灰度级在[0，1，...，t]范围内的像素点的集合，C₂表示所有灰度级在(t+1，...，L-1]范围内的像素点的集合。则C₁和C₂两个类的概率可表示为：

c.2.3计算图像的目标C₁和背景C₂两个类的均值和方差

两个类的均值分别为：

两个类的方差分别为：

c.2.4计算类内方差、类间方差和总的方差

类内方差类间方差和总的方差的计算如下：

其中，故

b.2.5采用最大类间方差得到分块分割阈值t

采用循环的方式，灰度级阈值t从0开始直到取得子图像块灰度级的最大值L-1，利用式(10)计算类间方差，当类间方差取得最大值时所对应的t就是分块分割的灰度级阈值；也可利用式(9)计算类内方差，当类内方差取得最小值时对应的t也即是分块分割的灰度级阈值；

b.2.6利用阈值t对子图像块S_m×n进行分割

对于子图像块S_m×n中的所有像素点，若其灰度值大于阈值t，则标记为前景即目标，并用1表示；若其灰度值小于或等于阈值t，则标记为背景，并用0表示。

d.获取覆盖区域

在完成对所有子图像块分割之后，将每个子图像块分割的结果合并得到破损唐卡图像的初始分割结果，这里记为h(x，y)，这时候破损区域的全部信息都已经包含在背景部分或者说物体部分中了，然而在一幅画作中，真实的破损形式众多，特征千差万别，欲获得破损区域的模板还需要对分块分割结果做进一步的处理。为此，根据步骤c得到的分割结果，获得基于对称性分块分割覆盖的灰度图像f_分割(x，y)和彩色图像f_分割(x，y)′。

f_分割(x，y)＝f(x，y)*h(x，y)

f_分割(x，y)′＝F(x，y)*h(x，y)

其中，如前所表示：F(x，y)为待分割的原始彩色图像，f(x，y)是其灰度化处理得到的灰度图像。

e.图像特征提取

e.1提取图像Gabor纹理特征

利用Gabor滤波器提取图像的纹理特征，图像Gabor特征可通过将输入图像与Gabor小波滤波器组卷积得到，即：

m(x，y)＝f_分割(x，y)*g(x，y) (11)

这里f_分割(x，y)是输入的步骤d得到的覆盖区域灰度图像，g(x，y)是Gabor滤波函数，定义为：

式中，λ是正弦函数波长，θ是Gabor核函数的方向，是相位偏移，σ为高斯函数的标准差，γ为空间的方面比，且γ的范围在0.23～0.92之间，x,y是给定位置的图像坐标，x',y'是旋转之后的坐标。

这里选择6个方向，分别为(0，π/6，2π/6，3π/6，4π/6，5π/6)，波长λ分别为(1/0.1，1/0.2，1/0.3)，σ分别为(0.5,1,1.5,2)；

通过修改参数，便可以得到一组不同尺度和方向的纹理特征，并作为最终的图像纹理特征feature_texture。

e.2提取颜色特征

对于颜色特征的提取，考虑到像素点的颜色差距应该与人眼视觉感觉到的差别相一致，故首先将待分割破损唐卡彩色图像f_分割(x，y)′ 从RGB空间转换为CIE Lab颜色空间，然后再提取颜色特征进行处理。

具体步骤如下：

e.2.1色彩空间转换

RGB空间转换为CIE Lab颜色空间具体转换过程为

则CIE Lab空间的亮度L分量，以及有关色彩的a，b分量分别表示为：

其中(X₀,Y₀,Z₀)是标准白色的刺激值。

e.2.2获得颜色特征

通过颜色特征转化得到颜色特征向量feature_color＝{L，a，b}。

e.3提取纹理和颜色特征

将步骤e.1获取的纹理特征feature_texture和步骤e.2获取的颜色特征feature_color最终构成融合纹理特征和颜色特征的特征集合feature＝{feature_texture,feature_color}。

f.多尺度多特征融合得到二次分割结果

f.1归一化特征

本发明特征归一化按式(18)所示的高斯归一化公式进行：

式中，1≤i≤M，M为特征集合feature所包含的特征数量，u_Fij表示第i个特征的第j维分量，u_Fj和σ_Fj分别表示第j维分量的期望和方差；

f.2基于多尺度多特征融合的分类

采用基于划分的聚类算法k-means分割破损区域，且步骤f.2所采用的聚类算法k-means与步骤b.2.1相同；

破损区域的样本集合为步骤f.1得到的经过归一化后的特征集合feature＝(x₁，x₂，...，x_n)，且每个x_i都是一个特征向量，需要将这个特征集合分成K个不同的类S＝(s₁，s₂，...s_k)，取k＝5。具体步骤如下：

第一步：随机选择k个元素，作为K个簇的初始中心点A＝{μ_i|i＝1，2，...，k}；

第二步：分别计算剩下的元素到K个簇中心的欧式距离，其中元素到簇中心点的欧式距离可以通过下述公式计算：

s＝(x_j-μ_i)²

其中，x_j和μ_i均为feature中的特征向量；

将这些元素分别划归到距离最近的簇中，即：

若||x_i-μ_j||＜||x_i-μ_l||，则x_i∈μ_j j，l＝1，2，...，K，j≠l

第三步：计算每个聚类中所有特征向量的质心，并将其作为新的聚类中心，新的聚类中心表达式为：

其中N_j是类s_j中样本的数目；

反复执行步骤(2)、(3)，聚类中心点不再变化为止，最终得到分类的结果。

g.形态学运算

g.1开运算

步骤f.2二次分割得到的结果存在大量小的破损区域、破损区域之间的粘连，采用先腐蚀后膨胀的过程即开运算进行处理，用来消除小物体、在纤细点处分离破损区域、平滑较大破损区域的边界，且同时不明显改变其面积。

腐蚀操作为：

膨胀操作为：

其中B(x)为结构元素，结构元素在形态变换中的作用相当于信号处理中的“滤波窗口”，这里采用3×3的全1模板。E为二次分割得到的结果图像。用B(x)对E进行腐蚀的结果就是把结构元素B平移后使B包含于E的所有点构成的集合。用B(x)对E进行膨胀的结果就是把结构元素B平移后使B与E的交集非空的点构成的集合。

g.2闭运算

先膨胀后腐蚀的过程称为闭运算。用来填充破损区域内细小空洞、连接邻近破损区域、平滑其边界的同时并不明显改变其面积。且该步骤的膨胀和腐蚀过程同步骤g.1。

h.去除小面积区域

统计步骤g.2得到的不同破损区域中像素点的个数m_i，将m_i＜ε 的区域去除，这里ε是阈值，可以通过反复实验的方式得到最佳值，本发明中，ε的取值为30～50之间。最后得到的结果就是唐卡图像破损区域模板。

Claims (4)

1.一种佛像类唐卡一致性破损区域的自动分割方法，其特征在于，首先对唐卡图像的头光区域进行垂直投影，利用一维函数对称性检测方法得到图像的对称轴，采用基于对称轴的分块分割方法得到初始分割结果；接着获取破损区域所覆盖的图像；然后利用Gabor变换提取纹理特征，结合Lab空间颜色特征，构造多尺度多特征集合，采用KNN分类得到二次分割结果；最后利用形态学运算对破损区域进行进一步细化，通过将小的破损区域去除，最终获得一致性破损区域的模板，具体步骤如下：

a.选定头部区域

经过数字化得到的待分割的唐卡原始图像是彩色图像，记为F(x，y)，采取手工选定或计算机选定的方式得到图像中佛像的头部区域t(x，y)；

b.检测对称轴

利用基于划分的聚类算法k-means对头部区域t(x，y)进行分割，得到唐卡图像的头光区域，通过垂直投影的方法得到头光区域的一维函数，利用一维函数对称性检测算法得到头光区域的对称轴对称轴即为整幅图像的对称轴；

c.基于对称性的分块采样初次分割

c.1唐卡图像分块

将原始的待分割图像F(x，y)进行图像的灰度化处理，得到灰度化图像f(x，y)；唐卡图像一般为左右对称，故以对称轴为轴，可以将图像分为左半部分f_左和右半部分f_右，即f＝f_左+f_右；

然后以对称轴为基准从图像顶端开始分别向左右边界延伸，进行分块采样；步骤如下：

c.1.1选择左边区域图像块S_左

从f(x，y)的顶端开始，以对称轴为基准向左依次选取图像块S_左1、S_左2...S_左N，图像块的大小均为的像素块；

c.1.2选择右边区域图像块S_右

从f(x，y)的顶端开始，以对称轴为基准向右依次选取图像块S_右1、S_右2...S_右N，图像块的大小均为的像素块；

c.1.3合并得到分块S

将关于对称轴相对称的S_左i和S_右i进行合并可以得到大小为m×n的分块S_m×n＝[S_左iS_右i]，i＝1，2，…N；

以对称轴为基准，从上到下，将以上步骤c.1.1至步骤c.1.3重复，直到将唐卡图像f(x，y)表示成M个大小相同的图像块S_m×n的集合；c.2对分块采样后的子图像块进行分割

经过步骤c.1，图像f(x，y)可表示为M个大小相同的图像块S_m×n的集合，即f(x，y)＝{S_i|i＝1，2，…M}；然后采用阈值法对每个子图像块S_m×n进行图像分割，对于子图像块S_m×n中的所有像素点，若其灰度值大于阈值t，则标记为前景即目标，并用1表示；若其灰度值小于阈值t，则标记为背景，并用0表示；

d.获取覆盖区域

在完成对所有子图像块分割之后，将每个子图像块分割的结果合并得到破损唐卡图像的初始分割结果，这里记为h(x，y)，这时破损区域的全部信息就都已经包含在步骤c得到的背景部分中了，根据步骤c得到的分割结果，获得基于对称性分块分割覆盖的灰度图像f_分割(x，y)和彩色图像f_分割(x，y)′，其中：

f_分割(x，y)＝f(x，y)*h(x，y)

f_分割(x，y)′＝F(x，y)*h(x，y)

e.图像特征提取

e.1提取图像Gabor纹理特征

利用Gabor滤波器提取图像的纹理特征，将步骤d得到的灰度图像与Gabor小波滤波器组卷积得到：

m(x，y)＝f_分割(x，y)*g(x，y) (11)

其中，g(x，y)是Gabor滤波函数，定义为：

$\left\{\begin{array}{c}{x}^{\′}=x\cos \mathrm{\θ}-y\sin \mathrm{\θ}\\ y^{\′}=xcos\mathrm{\θ}+ysin\mathrm{\θ}\end{array}\right.---\left(13\right)$

式中，λ是正弦函数波长，θ是Gabor核函数的方向，是相位偏移，σ为高斯函数的标准差，γ为空间的方面比也即Gabor核高斯函数的纵横比，且γ的范围值在0.23～0.92之间，x,y是给定位置的图像坐标，x',y'是旋转之后的坐标；通过改变这些参数，可以得到一组不同波长、不同尺度和不同方向的纹理特征，并作为最终的图像纹理特征feature_texture；

e.2提取颜色特征

首先将待分割破损唐卡图像f_分割(x，y)′从RGB空间转换为CIE Lab颜色空间，然后再提取颜色特征进行处理：

e.2.1色彩空间转换

RGB空间转换为CIE Lab颜色空间的具体过程为

$\left(\begin{array}{c}X\\ Y\\ Z\end{array}\right)=\left(\begin{array}{ccc}0.607& 0.174& 0.200\\ 0.299& 0.587& 0.114\\ 0.000& 0.066& 0.116\end{array}\right)\left(\begin{array}{c}R\\ G\\ B\end{array}\right)---\left(14\right)$

则CIE Lab空间的亮度L分量，以及有关色彩的a，b分量分别表示为：

$L=116\×\sqrt{\frac{Y}{Y_0}}-16---\left(15\right)$

$a=500\×(\sqrt[3]{\frac{X}{X_0}}-\sqrt[3]{\frac{Y}{Y_0}})---\left(16\right)$

$b=200\×(\sqrt[3]{\frac{Y}{Y_0}}-\sqrt[3]{\frac{Z}{Z_0}})---\left(17\right)$

其中，(X₀,Y₀,Z₀)为标准白色的刺激值；

e.2.2获得颜色特征

通过颜色特征转化得到颜色特征向量feature_color＝{L，a，b}；

e.3提取纹理和颜色特征

将步骤e.1获取的纹理特征feature_texture和步骤e.2获取的颜色特征feature_color最终构成融合纹理特征和颜色特征的特征集合feature＝{feature_texture,feature_color}；

f.多尺度多特征融合得到二次分割结果

f.1归一化特征

对步骤e提取的融合纹理特征和颜色特征的特征集合采用高斯归一化方法进行归一化处理，

高斯归一化公式如式(18)所示：

$N\left(u_{Fij}\right)=\frac{\frac{u_{Fij}-u_{Fj}}{3{\mathrm{\σ}}_{Fj}}+1}{2}---\left(18\right)$

式中，1≤i≤M，M为特征集合feature所包含的特征数量，u_Fij表示第i个特征向量的第j维分量，u_Fj和σ_Fj分别表示第j维分量的期望和方差；

f.2基于多尺度多特征融合的分类

采用基于划分的聚类算法k-means再次分割破损区域，破损区域的样本集合为步骤f.1得到的经过归一化后的特征集合，以此得到图像的二次分割结果；

g.形态学运算

g.1开运算

二次分割得到的结果存在大量小的破损区域以及破损区域之间的粘连，采用先腐蚀后膨胀的过程即开运算对步骤f.2得到的分类结果进行处理，用以消除小的破损区域、在纤细点处分离破损区域、平滑较大破损区域的边界，并同时不明显改变其面积，其中：

腐蚀操作为：

膨胀操作为：

式中，B(x)为结构元素，结构元素在形态变换中的作用相当于信号处理中的“滤波窗口”，这里采用3×3的全1模板；E为步骤f.2得到的结果图像；

g.2闭运算

先膨胀后腐蚀的过程称为闭运算，用来填充破损区域内细小空洞、连接邻近破损区域、平滑破损区域边界，并同时不明显改变其面积；且步骤g.2中的膨胀和腐蚀操作同步骤g.1；

h.去除小面积破损区域

统计步骤g.2得到的不同破损区域中像素点的个数m_i，将m_i＜ε的分类区域去除，其中，ε的取值为30～50之间。

2.根据权利要求1所述的一种佛像类唐卡一致性破损区域的自动分割方法，其特征在于，步骤b的具体过程如下：

b.1图像的灰度化处理

将步骤a得到的头部区域首先进行灰度处理：设头部区域t(x，y)为一幅包含有M×N个像素点的二维图像，采用平均值法，即对每个像素点R、G、B三分量的亮度分别求其平均值，然后将求得的这个平均值作为灰度图像的灰度值；

b.2提取头光区域

b.2.1采用基于划分的聚类算法k-means分割头部区域

第一步：将步骤a得到的头部区域t(x，y)的所有像素点构成集合T＝{t_i|i＝1，2，...total},T中的样本点用三个属性(x_m，y_m，f_m)进行描述，其中x_m代表像素点的水平坐标的值，y_m代表像素点垂直坐标的值，f_m代表像素点(x_m，y_m)的灰度值，随机选择k个元素，则作为k个簇的初始中心点A＝{μ_j|j＝1，2，...，k}；

第二步：分别计算剩下的元素到k个簇中心的欧式距离，其中元素到簇中心点的欧式距离通过下述公式计算：

$s={(t_{f_i}-{\mathrm{\μ}}_{f_j})}^2---\left(1\right)$

式中，为头部区域像素点i的灰度值，表示本次迭代中簇j中心点的灰度值，s即为图像中像素点与簇中心像素点灰度值的差异；

将T中剩余元素分别划归到距离最近的簇中，得到k个不同的元素聚类，即：

若||t_fi-μ_fj||＜||t_fi-μ_fl||，则t_i∈μ_j，

其中j，l＝1，2，...，k，j≠l

第三步：计算每个聚类中所有点的灰度值平均值，并将这个平均值作为新的聚类中心，新的聚类中心表达式为：

${\mathrm{\μ}}_j=\frac{1}{N_j}{\mathrm{\Σ}}_{t_i\∈W_j}{t}_i,j=1,2,...k---\left(2\right)$

其中，N_j是聚类W_j中样本的数目，t_i为像素点i的灰度值；

反复执行步骤(2)、(3)，直到聚类中心点不再变化为止，最终得到头部区域的分割结果；

b.2.2获得头光区域

分别统计b.2.1中每个聚类的像素点数目，然后选择像素点数目最多的聚类区域，这个包含像素点最多的颜色稳定的区域即为唐卡图像的头光区域；

b.3获取对称轴

b.3.1将头光区域所有像素点叠加投影得到S(x)

将头光区域中的像素点沿着垂直方向叠加投影，即：

$S\left(x\right)={\mathrm{\Σ}}_{j=0}^{n}t(x,y)\×mask(x,y)---\left(3\right)$

其中，n表示图像的高度，t(x，y)为待处理的头部区域，mask(x，y)为头光区域的二值模板，若某一像素点位于头光区域，则该点的mask(x，y)值为1，若在头光区域外则mask(x，y)值为0；

b.3.2求S(x)的奇函数与偶函数

式(3)所得的垂直方向投影函数是一维函数，假设对称轴的坐标为x_s(w/2≤x_s≤W-w/2)，其中，w为唐卡图像对称区域的宽度，W是式(3)所得一维函数的宽度；则S(u)是以x_s为坐标原点的函数，且S(u)＝S(x_s-u)；

计算S(u)的奇函数分量odd(u)和偶函数分量even(u)：

$odd\left(u\right)=\frac{S(u+x_s)-S(x_s-u)}{2}---\left(4\right)$

$even\left(u\right)=\frac{S(u+x_s)+S(x_s-u)}{2}---\left(5\right)$

b.3.3对称轴求取

使用式(6)对偶函数分量even(u)进行归一化得：

${\mathrm{even}}_n\left(u\right)=even\left(u\right)-\frac{1}{w}{\∫}_{u=-w/2}^{w/2}even\left(u\right)du---\left(6\right)$

式中，W为一维函数S(u)的宽度；

计算对称性测度，对称性测度symmetry可以通过下式计算：

$symmetry=\frac{{\∫}_{u=-w/2}^{w/2}{\mathrm{even}}_n{\left(u\right)}^{2}du-{\∫}_{u=-w/2}^{w/2}odd{\left(u\right)}^{2}du}{{\∫}_{u=-w/2}^{w/2}{\mathrm{even}}_n{\left(u\right)}^{2}du+{\∫}_{u=-w/2}^{w/2}odd{\left(u\right)}^{2}du}---\left(7\right)$

式(7)中，symmetry的取值范围为0.6～1，并计算此条件下的x_s值，求得的x_s值即为头光区域的对称轴的坐标值。

3.根据权利要求1所述的一种佛像类唐卡一致性破损区域的自动分割方法，其特征在于，步骤c.2中，具体过程如下：

c.2.1计算灰度级i在图像中出现的概率

设每个子图像块S_m×n构成的图像包含L个灰度级(0,1…,L-1)，灰度值为i的像素点的数量为n_i，图像总的像素点数量为N＝n₀+n₁+…n_L-1，灰度级i在图像中出现的概率为：

$p_i=\frac{n_i}{N},p_i\≥0,{\mathrm{\Σ}}_{i=0}^{L-1}{p}_i=1---\left(8\right)$

c.2.2计算图像的目标C₁和背景C₂的概率

假定通过灰度级阈值t将每一个子图像块的像素点分成目标C₁和背景C₂两类，其中C₁表示所有灰度级在[0，1，...，t]范围内的像素点的集合，C₂表示所有灰度级在(t+1，...，L-1]范围内的像素点的集合，则C₁和C₂两个类的概率可表示为：

${\mathrm{\ω}}_1=\underset{i=0}{\overset{t}{\Σ}}{p}_i$

${\mathrm{\ω}}_2=\underset{i=t+1}{\overset{L-1}{\Σ}}{p}_i$

c.2.3计算图像的目标C₁和背景C₂两个类的均值和方差两个类的均值分别为：

${\mathrm{\μ}}_1=\underset{i=0}{\overset{t}{\Σ}}{\mathrm{ip}}_i/{\mathrm{\ω}}_1$

${\mathrm{\μ}}_2=\underset{i=t+1}{\overset{L-1}{\Σ}}{\mathrm{ip}}_i/{\mathrm{\ω}}_2$

两个类的方差分别为：

${\mathrm{\σ}}_1^2=\underset{i=0}{\overset{t}{\Σ}}{(i-{\mathrm{\μ}}_1)}^2{p}_i/{\mathrm{\ω}}_1$

${\mathrm{\σ}}_2^2=\underset{i=t+1}{\overset{L-1}{\Σ}}{(i-{\mathrm{\μ}}_2)}^2{p}_i/{\mathrm{\ω}}_2$

c.2.4计算类内方差、类间方差和总的方差

类内方差类间方差和总的方差的计算如下：

${\mathrm{\σ}}_W^2={\mathrm{\ω}}_1{\mathrm{\σ}}_1^2+{\mathrm{\ω}}_2{\mathrm{\σ}}_2^2---\left(9\right)$

${\mathrm{\σ}}_B^2={\mathrm{\ω}}_1{\mathrm{\ω}}_2{({\mathrm{\μ}}_1-{\mathrm{\μ}}_2)}^2---\left(10\right)$

${\mathrm{\σ}}_T^2=\underset{i=0}{\overset{L-1}{\Σ}}{(i-{\mathrm{\μ}}_T)}^2{p}_i$

其中，故

c.2.5采用最大类间方差或最小类内方差得到分块分割的灰度级阈值t

采用循环的方式，t从0开始直到取得子图像块灰度的最大值L-1，利用式(10)计算相应的类间方差，当类间方差取得最大值时所对应的t就是分块分割的灰度级阈值；或利用式(9)计算类内方差，当类内方差取得最小值时对应的t也就是分块分割的灰度级阈值；

c.2.6利用阈值t对子图像块S_m×n进行分割

对于子图像块S_m×n中的所有像素点，若其灰度值大于阈值t，则标记为前景即目标，并用1表示；若其灰度值小于或等于阈值t，则标记为背景，并用0表示。

4.根据权利要求2所述的一种佛像类唐卡一致性破损区域的自动分割方法，其特征在于，步骤f.2所采用的基于划分的聚类算法k-means与步骤b.2.1相同，具体步骤如下：

破损区域的样本集合为步骤f.1得到的经过归一化后的特征集合feature＝(x₁，x₂，...，x_n)，且每个x_i都是一个特征向量，需要将这个特征集合分成K个不同的类S＝(s₁，s₂，...s_k)：

第一步：随机选择k个元素，作为K个簇的初始中心点A＝{μ_i|i＝1，2，...，k}；

第二步：分别计算剩下的元素到K个簇中心的欧式距离，其中元素到簇中心点的欧式距离可以通过下述公式计算：

s＝(x_j-μ_i)²

其中，x_j和μ_i均为feature中的特征向量；

将这些元素分别划归到距离最近的簇中，即：

若||x_i-μ_j||＜||x_i-μ_l||，则x_i∈μ_j j，l＝1，２，...，K，j≠ｌ；

第三步：计算每个聚类中所有特征向量的质心，并将其作为新的聚类中心，新的聚类中心表达式为：

${\mathrm{\μ}}_j=\frac{1}{N_j}{\mathrm{\Σ}}_{x_i\∈s_j}{x}_i,j=1,2,...K$

其中N_j是类s_j中样本的数目；

反复执行步骤(2)、(3)，直到聚类中心点不再变化为止，最终得到分类的结果。