CN104282008A - 对图像进行纹理分割的方法和装置 - Google Patents

Abstract

本发明提供一种对图像进行纹理分割的方法，该图像具有多种纹理，该方法包括以下步骤：将所述图像转换为灰度图像；根据所述图像中最大的纹理单元的大小，将所述灰度图像划分为大小相同的多个区域，并提取与所述多个区域一一对应的多个梯度特征向量；对提取的所述多个梯度特征向量进行模糊聚类分析，将所述多个区域进行分类，从而将所述图像中具有相同纹理的部分分类到相同的类别中。本发明提供一种对图像进行纹理分割的装置。本发明的对图像进行纹理分割的方法和装置可以适应多方向、多灰度级，并且计算量小。

Description

对图像进行纹理分割的方法和装置

技术领域

本发明涉及对图像进行分割的方法和装置，尤其涉及对图像进行纹理分割的方法和装置。

背景技术

随着信息技术的发展，图像越来越多的被用来获取和分析各种信息。图像处理领域的图像分割也已经成为研究的重点。

目前图像分割的方法主要分为颜色分割和纹理分割两大类。由于纹理较颜色特征信息更为丰富，纹理分割的意义更为重要。纹理分割目前主要有基于灰度共生矩阵和基于小波分析的纹理分割方法。《基于灰度共生矩阵的IKONOS影像中竹林信息提取》主要是将图像分割为多个小块区域，提取每个小区域的灰度共生矩阵作为区域的纹理特征，再通过对特征进行分类达到纹理分割的目的。《基于小波变换的纹理图像分割》将图像进行小波分解，在多尺度以及垂直、水平、45度方向上提取纹理特征，再对特征进行分类，达到纹理分割的目的。

上述纹理分割方法，能够提取纹理特征并构造特征向量，但是对于以下几种情况，不能很好的适应。由于灰度对比度也包含了很多信息，而灰度共生矩阵随着灰度级别的增加，计算量呈几何级数增加，一般多用于8个灰度级的图像纹理特征提取。并且灰度共生矩阵多提取指定8个方向的特征，当纹理方向不一致时，该方法无法提取准确的纹理特征。小波变换方法同样在多个纹理计算方向的情况下，计算量和准确性变差。对于图像分割的子窗口的尺寸，也严重影响计算速度，窗口越大，计算量越大。当纹理单元大小不均匀时，以上方法也不能很好的判断分割出来。如此，现有的图像纹理分割方法存在无法适应多方向、多灰度级以及计算量大的问题。

发明内容

为了解决上述问题，本发明提供了一种图像分割方法和装置，可以解决无法适应多方向、多灰度级以及计算量大的问题。

一种对图像进行纹理分割的方法，所述图像具有多种纹理，包括以下步骤：

a)将所述图像转换为灰度图像；

b)根据所述图像中最大的纹理单元的大小，将所述灰度图像划分为大小相同的多个区域，并提取与所述多个区域一一对应的多个梯度特征向量；

c)对提取的所述多个梯度特征向量进行模糊聚类分析，将所述多个区域进行分类，从而将所述图像中具有相同纹理的部分分类到相同的类别中。

在步骤b）中，所述多个区域中的每个区域的大小大于等于所述最大的纹理单元的大小。

在步骤b）中，通过以下步骤提取所述多个梯度特征向量：

b1)以水平向右为梯度方向角的0度方向，计算所述多个区域中的每个区域中的4m×4m个像素中的每个像素的梯度值和第一梯度方向角，m是大于等于1的整数；

b2)对计算的所述多个区域中的任一个区域中的4m×4m个像素各自的梯度值和梯度方向角进行统计，得到梯度直方图，在所述梯度直方图中，多个直方柱与多个第一方向角一一对应，将所述多个直方柱中具有峰值g_max的直方柱所对应的第一方向角作为所述任一个区域的梯度主方向角；

b3)以所述任一个区域的中心为中心，以所述梯度主方向角为所述梯度方向角的0度方向，限定与所述任一个区域大小相同的新区域；

b4)根据步骤b1）的计算结果，取出位于所述新区域中的4m×4m个像素中的每个像素的梯度值和第一梯度方向角，其中取出的每个像素的第一梯度方向角变更为第二梯度方向角，所述第二梯度方向角是从所述第一梯度方向角顺时针旋转到所述梯度主方向角的角度；

b5)将所述新区域均匀划分为16个子区域，每个子区域中具有m×m个像素，以所述梯度主方向角为所述0度方向，将0-360度方向均匀划分为多个第二方向角，如果任一个子区域中的任一个像素的第二梯度方向角不等于任一个第二方向角，则将所述任一个像素的梯度值分别分解到所述多个第二方向角中与所述任一个像素的第二梯度方向角相近的两个第二方向角上，得到与所述两个第二方向角分别对应的两个分解梯度值，如果任一个子区域中的任一个像素的第二梯度方向角等于任一个第二方向角，则将所述任一个像素的梯度值作为与所述任一个第二方向角相对应的分解梯度值，从而对于每个子区域，获得多个分解梯度值，所述多个分解梯度值中的任一个分解梯度值与所述多个第二方向角中的一个第二方向角相对应；

b6）将所述多个分解梯度值中与任一个第二方向角对应的数个分解梯度值进行累加，获得与所述任一个第二方向角对应的累加梯度值，从而对于所述每个子区域，获得与所述多个第二方向角一一对应的多个累加梯度值，将所述多个累加梯度值分别对于所述峰值g_max进行归一化，得到多个归一化累加梯度值；

b7）将所述新区域中16个子区域各自的多个归一化累加梯度值进行排列，获得所述任一个区域的梯度特征向量；

其中，对所述每个区域进行步骤b2)-b7)，从而提取与所述多个区域一一对应的所述多个梯度特征向量。

在步骤b1)中，根据以下等式计算每个像素的梯度值和第一梯度方向角，

等式1：

$g(x,y)=\sqrt{(P(x+1,y)-P{(x-1,y))}^2+P(x,y+1)-P{(x,y-1))}^2}$

等式2：

θ（x,y)＝αtan2((P(x,y+1)-P(x,y-1))/(P(x+1,y)-P(x-1,y)))

其中，αtan2表示反正切函数，g(x,y)表示位于所述灰度图像中坐标为(x,y)的像素点的梯度值，θ(x,y)表示位于所述灰度图像中坐标为(x,y)的像素点的第一梯度方向角，P(x,y)表示位于所述灰度图像中坐标为(x,y)的像素点的灰度值。

在步骤b2)中，所述梯度直方图的角度范围是0-360度，将0-360度均匀划分成所述多个第一方向角。

在步骤b4)中，当所述新区域包含超出所述灰度图像的超出范围时，将所述任一个区域中距离所述超出范围最近的多个像素各自的梯度值和第一梯度方向角取出作为所述超出范围中的多个像素各自的梯度值和第一梯度方向角。

在步骤b5)中，根据以下等式得到与所述两个第二方向角分别对应的两个分解梯度值，

等式3：g_θ1=g×sin(θ₂-θ)/sin(θ₂-θ₁)

等式4：g_θ2=g×sin(θ-θ₁)/sin(θ₂-θ₁)

其中，g表示所述任一个像素的梯度值，θ表示所述任一个像素的第二梯度方向角，θ₁和θ₂表示与θ相近的两个第二方向角，g_θ1表示g被分解到第二方向角θ₁上的分解梯度值，g_θ2表示g被分解到第二方向角θ₂上的分解梯度值。

在步骤b6)中，对于所述每个子区域，根据以下等式得到所述多个归一化累加梯度值，

等式5： $g_{\mathrm{ej}}=\frac{\underset{i=1}{\overset{m\×m}{\mathrm{\Σ}}}{g}_{\mathrm{eij}}}{g_{\max }}$

其中，g_ej表示在第e个子区域中的所述多个第二方向角中的第j个第二方向角上的归一化累加梯度值，g_eij表示在第e个子区域中的第i个像素在第j个第二方向角上的分解梯度值，e，j是大于1的整数，1≤e≤16。

在步骤b7）中，将所述新区域划分为大小相同的4个次区域，每个次区域包含4个所述子区域，根据步骤b6）中获得的多个归一化累加梯度值，按照逐个子区域逆时针方向，分别得到对应于4个次区域的4个次梯度特征向量G’₁、G’₂、G’₃、G’₄，其中，

G’₁=｛g₁₁、g₁₂、…g_1a、g₂₁、g₂₂、…g_2a、g₃₁、g₃₂、…g_3a、g₄₁、g₄₂、…g_4a｝，

G’₂=｛g₅₁、g₅₂、…g_5a、g₆₁、g₆₂、…g_6a、g₇₁、g₇₂、…g_7a、g₈₁、g₈₂、…g_8a｝，

G’₃=｛g₉₁、g₉₂、…g_9a、g₁₀₁、g₁₀₂、…g_10a、g₁₁₁、g₁₁₂、…g_11a、g₁₂₁、g₁₂₂、…g_12a｝，

G’₄=｛g₁₃₁、g₁₃₂、…g_13a、g₁₄₁、g₁₄₂、…g_14a、g₁₅₁、g₁₅₂、…g_15a、g₁₆₁、g₁₆₂、…g_16a｝，，

其中，a表示所述多个第二方向角的个数，且a是大于等于1的整数，1≤j≤a，

将所述次梯度特征向量G’₁、G’₂、G’₃、G’₄依次排列组成所述任一个区域的所述梯度特征向量G_r=｛g₁₁、g₁₂…g_16a｝，G_r表示与第r个区域对应的梯度特征向量，1≤r≤n，n表示所述多个区域的个数，且r，n是整数。

在步骤c）中，通过以下步骤进行所述模糊聚类分析：

c1)根据所述图像中的纹理种类，设置所述模糊聚类分析的C个类别，从所述多个梯度特征向量中任意取出一个梯度特征向量作为某个类别的当前聚类中心，得到C个当前聚类中心；

c2)计算所述多个梯度特征向量分别到每个当前聚类中心的隶属度，对于所述每个当前聚类中心，得到与所述多个梯度特征向量一一对应的多个隶属度，将所述多个隶属度中大于等于预定阈值t的隶属度所对应的梯度特征向量所对应的区域归类到与所述每个当前聚类中心所对应的类别，

c3)根据对于所述每个当前聚类中心计算得到的所述多个隶属度，计算与所述C个类别一一对应的C个更新聚类中心，并将所述C个更新聚类中心作为所述C个当前聚类中心；

c4）当所述多个梯度特征向量被归类的类别发生变化的比率大于预定的阈值v或重复步骤c2)和c3)的次数小于次数阈值N时，重复进行步骤c2），否则进行步骤c4）；

c5）将与所述C个当前聚类中心的隶属度都小于预定阈值t的梯度特征向量对应的区域归类到未知类别。

进一步包括如下步骤：

c6）将与归类到所述未知类别的某个区域所对应的4个次梯度特征向量G’₁、G’₂、G’₃、G’₄分别组成4个新梯度特征向量G’’₁、G’’₂、G’’₃、G’’₄，

G’’₁=｛g₁₁、g₁₂…g_4a，g₁₁、g₁₂…g_4a，g₁₁、g₁₂…g_4a，g₁₁、g₁₂…g_4a｝，

G’’₂=｛g₅₁、g₅₂、…g_8a，g₅₁、g₅₂、…g_8a，g₅₁、g₅₂、…g_8a，g₅₁、g₅₂、…g_8a｝，

G’’₃=｛g₉₁、g₉₂、…g_12a，g₉₁、g₉₂、…g_12a，g₉₁、g₉₂、…g_12a，g₉₁、g₉₂、…g_12a｝，

G’’₄=｛g₁₃₁、g₁₃₂、…g_16a，g₁₃₁、g₁₃₂、…g_16a，g₁₃₁、g₁₃₂、…g_16a，g₁₃₁、g₁₃₂、…g_16a｝，

c7）计算任一个新梯度特征向量分别到步骤c3)中的所述C个当前聚类中心的多个隶属度，并将所述任一个新梯度特征向量所对应的次区域分类到与所述多个隶属度中的最大隶属度对应的类别，从而将4个次区域各自分类到相应的类别中。

步骤c2)中，根据以下等式计算所述隶属度，

等式6： ${\mathrm{\μ}}_t\left(G_r\right)=\frac{{(1/{\left|\right|G_r-p_t\left|\right|}^2)}^{1/(b-1)}}{{\mathrm{\Σ}}_{k=1}^C{(1/{\left|\right|G_r-p_k\left|\right|}^2)}^{1/(b-1)}},$

其中，μ_t(G_r)表示第r个区域的梯度特征向量到第t个类别的当前聚类中心的隶属度，G_r表示所述多个区域中的第r个区域的梯度特征向量，p_t表示第t个类别的当前聚类中心，1≤t≤C，t为大于等于1的整数，b是模糊因子，用于控制聚类分析结果的模糊程度。

步骤c3)中，根据以下等式计算所述C个更新聚类中心中的第t个类别的更新聚类中心p_t’，

等式7： ${p_t}^{,}=\frac{{\mathrm{\Σ}}_{r-1}^n{\left[{\mathrm{\μ}}_t\left(G_r\right)\right]}^b{G}_r}{{\mathrm{\Σ}}_{r-1}^n{\left[{\mathrm{\μ}}_t\left(G_r\right)\right]}^b}.$

在步骤b3）中，以所述任一个区域的中心为中心，将所述任一个区域顺时针旋转一角度，所述角度为所述梯度主方向角，旋转后的任一个区域所包含的区域为所述新区域。

本发明还提供一种对图像进行纹理分割的装置，该图像分割装置使用如上所述的对图像进行纹理分割的方法。

附图说明

图1是根据本发明的实施例的图像；

图2是根据本发明的实施例的对图像进行纹理分割的方法的流程图；

图3是根据本发明的实施例的区域r的示意图；

图4是根据本发明的实施例的梯度直方图；

图5是根据本发明的实施例的区域r在一部分灰度图像中的示意图；

图6是根据本发明的实施例的新区域B的示意图；

图7是根据本发明的实施例的新区域B的16个子区域；

图8是图2中的步骤s2的具体流程图；

图9是图2中的步骤s3的具体流程图；

图10是根据本发明的实施例的对图像进行分割的结果示意图；

图11是根据本发明的实施例的对图像进行进一步分割的结果示意图。

具体实施方式

下面结合实施例，对本发明进行详细说明。

本实施例中，采用IKNOS多光谱影像作为数据源，分辨率为1m，从一个地区中获取如如图1所示的图像，该图像具有512×512个像素，且该图像可以是彩色图。该图像具有多种纹理，本例中，例如具有该地区中的沙漠纹理和人工物纹理这两种纹理。

参考图2，本发明的对图像进行纹理分割的方法包括以下步骤：

在步骤S0中，获取该图像，在步骤S1中，将该图像转换为灰度图像。

接着，在步骤S2中，根据该图像中最大的纹理单元的大小，将灰度图像划分为大小相同的多个区域，并提取与多个区域一一对应的多个梯度特征向量。其中，多个区域中的每个区域的大小均大于等于该最大的纹理单元的大小。另外，一种纹理包括多个纹理单元，且多个纹理单元的大小可以一样或不一样。

在步骤S2中，具体进行如图8所示的步骤。

如图8所示，首先，在步骤S21中，以水平向右为梯度方向角的0度方向，计算多个区域中的每个区域中的4m×4m个像素中的每个像素的梯度值和第一梯度方向角，m是大于等于1的整数。较佳的是，m是大于等于4的整数。图3是根据本发明的实施例的一个区域r的示意图，其中m=4，即，区域r具有16×16个像素。这里，根据以下等式计算每个区域中的每个像素的梯度值和第一梯度方向角，

等式1：

$g(x,y)=\sqrt{(P(x+1,y)-P{(x-1,y))}^2+P(x,y+1)-P{(x,y-1))}^2}$

等式2：

θ（x,y)＝αtan2((P(x,y+1)-P(x,y-1))/(P(x+1,y)-P(x-1,y)))

其中，g(x,y)表示位于灰度图像中坐标为(x,y)的像素点的梯度值，θ(x,y)表示位于灰度图像中坐标为(x,y)的像素点的第一梯度方向角，P(x,y)表示位于灰度图像中坐标为(x,y)的像素点的灰度值。αtan2表示反正切函数，求得的角度值等于X轴与通过原点和给定坐标点(x,y)的直线之间的夹角，且其范围是0-360度。

这里，例如该图像具有512×512个像素，则相应的灰度图像具有512×512个像素，每个区域具有16×16个像素，那么本例中一共有32×32=1024个区域。对1024个区域中的每个区域中的每个像素进行计算，得到该灰度图像中的每个像素的梯度值和第一梯度方向角。

接着，在步骤S22中，对计算的多个区域中的任一个区域中的16×16个像素各自的梯度值和梯度方向角进行统计，得到梯度直方图，在梯度直方图中，多个直方柱与多个第一方向角一一对应，将多个直方柱中具有峰值g_max的直方柱所对应的第一方向角作为任一个区域的梯度主方向角。

本例中，以图3所示的区域r作为任一个区域进行说明。图3中的多个黑色小箭头代表计算的各个像素的梯度值和第一梯度方向角，对这些梯度值和第一梯度方向角进行统计，得到如图4所示的梯度直方图，其中，横轴表示第一方向角，纵轴表示在各个第一方向角上的梯度累加值。该梯度直方图的第一方向角的范围是0-360度，以10度为间隔，将360度均匀分成36个第一方向角，即，0度、10度、20度、…350度。其中，0度就是360度。

具体的，例如区域r中的某个像素的梯度值为g’，第一梯度方向角为θ’，如果θ’不等于上述任一个第一方向角，那么通过以下等式将该像素的梯度值g’分解到与θ’相近的两个第一方向角θ₁’、θ₂’上。

g_θ1’=g’×sin(θ₂’-θ’)/sin(θ₂’-θ₁’)

g_θ2’=g’×sin(θ’-θ₁’)/sin(θ₂’-θ₁’)

例如，某个像素的第一梯度方向角θ’为15度，那么就将g’分解到两个第一方向角10度和20度上，即，θ₁’=10度、θ₂’=20度，得到对应于10度的梯度值g₁₀和对应于20度的梯度值g₂₀。

如果θ’等于上述任一个第一方向角，例如θ’为10度，那么就将梯度值g’作为对应于10度的梯度值g₁₀，而不进行上述分解。

对区域r中的每个像素进行如上操作，对于每个第一方向角可以得到一个或多个梯度值，并将对应于每个第一方向角的一个或多个梯度值进行累加，得到与每个第一方向角对应的梯度累加值，作为图4中的多个直方柱。如图4所示，这些直方柱与36个第一方向角一一对应，另外，为了简化，图4中仅显示了与0度、10度、20度、30度、40度、50度、60度、70度一一对应的8个直方柱，且例如具有峰值g_max的直方柱对应于40度，那么就将40度作为区域r的梯度主方向角，如图5中箭头F所示的方向角为梯度主方向角。

在步骤S23中，以区域r的中心O为中心，以梯度主方向角40度为梯度方向角的0度方向，限定与区域r大小相同的新区域B，如图6所示。这里，以区域r的中心O为中心，将区域r顺时针旋转一角度，该角度为梯度主方向角40度，即，将区域r顺时针旋转40度，旋转后的区域r所包含的区域为新区域B。

在步骤S24中，取出位于新区域B中的16×16个像素中的每个像素的梯度值和第一梯度方向角，其中取出的每个像素的第一梯度方向角变更为第二梯度方向角，该第二梯度方向角是从第一梯度方向角顺时针旋转到梯度主方向角的角度。

例如，某个像素的第一梯度方向角θ’为15度，从15度顺时针旋转到梯度主方向角40度的角度是335度，即，某个像素的第二梯度方向角为335度。

当新区域B包含超出灰度图像的超出范围时，将区域r中距离超出范围最近的多个像素各自的梯度值和第一梯度方向角取出作为超出范围中的多个像素各自的梯度值和第一梯度方向角。这里，如果区域r位于灰度图像的边界附近，那么在区域r进行上述旋转后得到的新区域B有可能包含位于灰度图像之外的超出范围，而在超出范围中并没有任何像素。为此，本发明在超出范围中填充数个像素，并将区域r中距离某个填充像素最近的像素的梯度值和第一梯度方向角作为某个填充像素的梯度值和第一梯度方向角。这样可以确保新区域B中仍然有16×16个像素。

这里，取出的位于新区域B中的16×16个像素中的每个像素的梯度值就是前面计算出的每个像素的梯度值，但每个像素的第一梯度方向角变更为第二梯度方向角，如此，后续提取的区域r的梯度特征向量可以保证对纹理方向旋转具有不变性。

在步骤S25中，将新区域B均匀划分为16个子区域，即，子区域①、②、……如图7所示，每个子区域中具有4×4个像素，以梯度主方向角40度为0度方向，将0-360度方向均匀划分为多个第二方向角，如果任一个子区域中的任一个像素的第二梯度方向角不等于任一个第二方向角，则将所述任一个像素的梯度值分别分解到所述多个第二方向角中与所述任一个像素的第二梯度方向角相近的两个第二方向角上，得到与所述两个第二方向角分别对应的两个分解梯度值，如果任一个子区域中的任一个像素的第二梯度方向角等于任一个第二方向角，则将所述任一个像素的梯度值作为与所述任一个第二方向角相对应的分解梯度值，从而对于每个子区域，获得多个分解梯度值，所述多个分解梯度值中的任一个分解梯度值与所述多个第二方向角中的一个第二方向角相对应。

本例中，将0-360度方向均匀划分为8个第二方向角，即，0度、45度、90度、135度、180度、225度、270度、315度。此外，为便于说明，新区域B如图7中所示，但是该新区域B仍然以梯度主方向角40度为0度方向，且8个第二方向角是以40度为0度方向的8个角度。

当任一个像素的第二梯度方向角不等于上述任一个第二方向角时，根据以下等式得到某个像素在两个第二方向角上的两个分解梯度值，

等式3：g_θ1=g×sin(θ₂-θ)/sin(θ₂-θ₁)

等式4：g_θ2=g×sin(θ-θ₁)/sin(θ₂-θ₁)

其中，g表示任一个像素的梯度值，θ表示任一个像素的第二梯度方向角，θ₁和θ₂表示与θ相近的两个第二方向角，g_θ1表示g被分解到第二方向角θ₁上的分解梯度值，g_θ2表示g被分解到第二方向角θ₂上的分解梯度值。

以子区域①为例，有4×4=16个像素，例如，任一个像素的第二梯度方向角θ为335度，不等于8个第二方向角中的任一个，那么与335度相近的两个第二方向角分别为315度和0度，即θ₁为315度，θ₂为0度。根据上述等式3和等式4可以得到g₃₁₅和g₀，g₃₁₅是对应于第二方向角为315度的分解梯度值，g₀是对应于第二方向角为0度的分解梯度值。另外，若任一个像素的第二梯度方向角θ是8个第二方向角中的一个，例如，θ=0度，则不对这一个像素的梯度值进行分解，且将这一个像素的梯度值作为对应于第二方向角0度的分解梯度值。

对于子区域①，对其中的16个像素中的每个像素的梯度值进行上述分解，可以获得多个分解梯度值。

在步骤S26中，将多个分解梯度值中与任一个第二方向角对应的一个或数个分解梯度值进行累加，获得与任一个第二方向角对应的累加梯度值，从而对于每个子区域，获得与多个第二方向角一一对应的多个累加梯度值，将多个累加梯度值分别对于峰值g_max进行归一化，得到多个归一化累加梯度值。

对于每个子区域，根据以下等式得到多个归一化累加梯度值，

等式5： $g_{\mathrm{ej}}=\frac{\underset{i=1}{\overset{m\×m}{\mathrm{\Σ}}}{g}_{\mathrm{eij}}}{g_{\max }}$

其中，g_ej表示在第e个子区域中的多个第二方向角中的第j个第二方向角上的归一化累加梯度值，g_eij表示在第e个子区域中的第i个像素在第j个第二方向角上的分解梯度值，e，j是大于1的整数，1≤e≤16。本例中，如上所述m=4。

以子区域①为例，将分别对应于8个第二方向角，即，分别对应于0度、45度、90度、135度、180度、225度、270度、315度的一个或数个分解梯度值进行累加，例如，将对应于第二方向角为0度的一个或数个分解梯度值进行累加，得到子区域①中与0度对应的一个累加梯度值。如此，可以得到子区域①的8个累加梯度值。

其中，0度、45度、90度、135度、180度、225度、270度、315度这8个第二方向角分别是第1、2、…7、8个第二方向角。

对于16个子区域分别进行如上操作，可以得到16×8=128个累加梯度值。根据等式5，将128个累加梯度值分别除以峰值g_max，得到128个归一化累加梯度值，从而得到新区域B中的128个归一化累加梯度值。

在步骤S27中，将新区域B中的128个归一化累加梯度值进行排列，获得区域r的梯度特征向量。

其中，将新区域B划分为大小相同的4个次区域，每个次区域包含4个子区域，如图7所示，例如，第1个次区域包含标号①、②、③、④的4个子区域，第2个次区域包含标号⑤、⑥、⑦、⑧的4个子区域，第3个次区域包含标号⑨、⑩、的4个子区域，第4个次区域包含标号的4个子区域。

根据步骤S26中获得的128个归一化累加梯度值，按照逐个子区域逆时针方向，分别得到对应于4个次区域的4个次梯度特征向量G’₁、G’₂、G’₃、G’₄，其中，

G’₁=｛g₁₁、g₁₂、…g_1a、g₂₁、g₂₂、…g_2a、g₃₁、g₃₂、…g_3a、g₄₁、g₄₂、…g_4a｝，

G’₂=｛g₅₁、g₅₂、…g_5a、g₆₁、g₆₂、…g_6a、g₇₁、g₇₂、…g_7a、g₈₁、g₈₂、…g_8a｝，

G’₃=｛g₉₁、g₉₂、…g_9a、g₁₀₁、g₁₀₂、…g_10a、g₁₁₁、g₁₁₂、…g_11a、g₁₂₁、g₁₂₂、…g_12a｝，

G’₄=｛g₁₃₁、g₁₃₂、…g_13a、g₁₄₁、g₁₄₂、…g_14a、g₁₅₁、g₁₅₂、…g_15a、g₁₆₁、g₁₆₂、…g_16a｝，

其中，a表示多个第二方向角的个数，且a是大于等于1的整数，1≤j≤a，本例中a=8，

如图7所示，按照逐个子区域逆时针方向，子区域①、②、③、④构成的第1个次区域的梯度特征向量G’₁=｛g₁₁、g₁₂…g₁₈、g₂₁、g₂₂、…g₂₈、g₃₁、g₃₂、…g₃₈、g₄₁、g₄₂、…g₄₈｝，子区域⑤、⑥、⑦、⑧构成的第2个次区域的梯度特征向量G’₂=｛g₅₁、g₅₂…g₅₈、g₆₁、g₆₂、…g₆₈、g₇₁、g₇₂、…g₇₈、g₈₁、g₈₂、…g₈₈｝，子区域⑨、⑩、构成的第3个次区域的梯度特征向量G’₃=｛g₉₁、g₉₂…g₉₈、g₁₀₁、g₁₀₂、…g₁₀₈、g₁₁₁、g₁₁₂、…g₁₁₈、g₁₂₁、g₁₂₂、…g₁₂₈｝，子区域构成的第4个次区域的梯度特征向量G’₄=｛g₁₃₁、g₁₃₂…g₁₃₈、g₁₄₁、g₁₄₂、…g₁₄₈、g₁₅₁、g₁₅₂、…g₁₅₈、g₁₆₁、g₁₆₂、…g₁₆₈｝。

将4个次梯度特征向量G’₁、G’₂、G’₃、G’₄依次排列组成区域r的梯度特征向量G_r=｛g₁₁、g₁₂…g₁₆₈｝，G_r表示与第r个区域（即，区域r）对应的梯度特征向量，1≤r≤n，n表示多个区域的个数，且r，n是整数，本例中，n=1024。

对于1024个区域中的每个区域，进行如区域r一样的处理，即，对每个区域进行步骤S22-S27，从而提取1024个梯度特征向量，即提取与1024个区域一一对应的1024个梯度特征向量。

本发明中，区域的大小是以像素的个数为单位，两个区域的大小相同，即两个区域的像素个数相同。本例中，1024个区域中的每个区域的大小是16×16个像素。

如图2所示，在步骤S3中，对提取的多个梯度特征向量进行模糊聚类分析，将多个区域进行分类，从而将图像中具有相同纹理的部分分类到相同的类别中。

在步骤S3中，例如使用模糊C均值聚类法来进行模糊聚类分析，其具体步骤如图9所示。

参考图9，在步骤S31，根据图1所示的图像中的纹理种类，设置模糊聚类分析的C个类别，从多个梯度特征向量中任意取出一个梯度特征向量作为某个类别的当前聚类中心，得到C个当前聚类中心。

本例中，该图像具有两类纹理（两个纹理类别），沙漠纹理和人工物纹理，即，C=2。例如取出G₁作为第1个类别（即，对应于沙漠纹理的类别）的当前聚类中心，取出G₂作为第2个类别（即，对应于人工物纹理的类别）的当前聚类中心。如上所述，G₁表示第1个区域的梯度特征向量，G₂表示第2个区域的梯度特征向量。作为替换，也可以取G₁作为第2个类别的当前聚类中心，即2个当前聚类中心可以是相同的。

在步骤S32，计算多个梯度特征向量分别到每个当前聚类中心的隶属度，对于每个当前聚类中心，得到与多个梯度特征向量一一对应的多个隶属度，将多个隶属度中大于等于预定阈值t的隶属度所对应的梯度特征向量所对应的区域归类（分类）到与每个当前聚类中心所对应的类别。

这里，根据以下等式计算隶属度，

等式6： ${\mathrm{\μ}}_t\left(G_r\right)=\frac{{(1/{\left|\right|G_r-p_t\left|\right|}^2)}^{1/(b-1)}}{{\mathrm{\Σ}}_{k=1}^C{(1/{\left|\right|G_r-p_k\left|\right|}^2)}^{1/(b-1)}},$

其中，μ_t(G_r)表示第r个区域的梯度特征向量到第t个类别的当前聚类中心的隶属度，G_r表示多个区域中的第r个区域的梯度特征向量，p_t表示第t个类别的当前聚类中心，1≤t≤C，t为大于等于1的整数，b是模糊因子，用于控制聚类分析结果的模糊程度，b的取值范围通常与多个区域的数量n有关，较佳的，1.25≤b≤2.5。本例中，例如取模糊因子b=2。

本例中，一共是两类，则有两个当前聚类中心，即，p₁=G₁，p₂=G₂。预定阈值t例如为60%，且根据上述等式6计算的隶属度μ₁(G_r)≥60%，则将第r个区域归类到第1个类别。对于第1个类别，可以计算得到1024个隶属度，对于第2个类别，也可以计算得到1024个隶属度。

在步骤S33，根据对于每个当前聚类中心计算得到的多个隶属度，计算与C个类别一一对应的C个更新聚类中心，并将C个更新聚类中心作为C个当前聚类中心，

这里，根据以下等式计算2个更新聚类中心中的第t个类别的更新聚类中心p_t’，

等式7： ${p_t}^{,}=\frac{{\mathrm{\Σ}}_{r-1}^n{\left[{\mathrm{\μ}}_t\left(G_r\right)\right]}^b{G}_r}{{\mathrm{\Σ}}_{r-1}^n{\left[{\mathrm{\μ}}_t\left(G_r\right)\right]}^b}.$

在步骤S34中，当多个梯度特征向量被归类的类别发生变化的比率大于预定的阈值v或重复步骤S32、S33的次数小于次数阈值N时，重复进行步骤S32，否则进行步骤S35。

本例中，阈值v例如为10%，次数阈值N例如为10，在每次进行步骤S32时，多个梯度特征向量中的一些梯度特征向量被归类的类别会发生变化，例如从归类为第1个类别改变为归类为第2个类别。当归类的类别发生变化的一些梯度特征向量占全部多个梯度特征向量的变化比例大于10%或者重复步骤S32和S33的重复次数小于10时，回到步骤S32，否则进行步骤S35。当该变化比例小于等于10%时，各个区域的隶属度稳定，即各个区域被归类的类别稳定。

另外，设置次数阈值N是为了避免长时间重复步骤S32和S33。

在步骤S35，将与2个当前聚类中心的隶属度都小于预定阈值t的梯度特征向量所对应的区域归类到未知类别。

例如，第3个梯度特征向量到第1个类别的当前聚类中心p₁的隶属度小于60%，到第2个类别的当前聚类中心p₁的隶属度也小于60%，那么将第3个区域归类为未知类别。其中，该未知类别与上述2个类别不同。

这样，可以将1024个区域分别归类到上述2个类别以及未知类别中，也就是说，图像中与某个区域对应的部分就被归类到其中一个类别中，如此，可以将图像中具有相同纹理的部分划分到相同的类别中。划分到第1个类别的是图像中具有沙漠纹理的部分，划分到第2个类别的是图像中具有人造物纹理的部分。如图10所示，将图像中具有沙漠纹理的部分划分到第1个类别中，将图像中具有人造物纹理的部分划分到第2个类别中，而标有“？”的一些区域是被归类到未知类别的区域。

如此，本发明的对图像进行纹理分割的方法可以适应多方向、多灰度级，并将计算量小。

较佳的，本发明还包括步骤S36，在步骤S36中，将与归类到未知类别的例如第3个区域所对应的4个次梯度特征向量G’₁、G’₂、G’₃、G’₄分别组成4个新梯度特征向量G’’₁、G’’₂、G’’₃、G’’₄，

G’’₁=｛g₁₁、g₁₂…g₄₈，g₁₁、g₁₂…g₄₈，g₁₁、g₁₂…g₄₈，g₁₁、g₁₂…g₄₈｝，

G’’₂=｛g₅₁、g₅₂、…g₈₈，g₅₁、g₅₂、…g₈₈，g₅₁、g₅₂、…g₈₈，g₅₁、g₅₂、…g₈₈｝，

G’’₃=｛g₉₁、g₉₂、…g₁₂₈，g₉₁、g₉₂、…g₁₂₈，g₉₁、g₉₂、…g₁₂₈，g₉₁、g₉₂、…g₁₂₈｝，

G’’₄=｛g₁₃₁、g₁₃₂、…g₁₆₈，g₁₃₁、g₁₃₂、…g₁₆₈，g₁₃₁、g₁₃₂、…g₁₆₈，g₁₃₁、g₁₃₂、…g₁₆₈｝，

可以看出，本发明是第3个区域划分为4个次区域，每个次区域对应于一个次梯度特征向量，将一个次梯度特征向量扩展为原来的4倍，形成一个新梯度特征向量，共形成4个新梯度特征向量G’’₁、G’’₂、G’’₃、G’’₄，且4个新梯度特征向量G’’₁、G’’₂、G’’₃、G’’₄一一对应于4个次区域。

步骤S37中，计算任一个新梯度特征向量分别到步骤S32中的2个当前聚类中心的2个隶属度，并将任一个新梯度特征向量所对应的次区域分类到与2个隶属度中的最大隶属度所对应的类别，从而将4个次区域分类到相应的类别中。

根据上述等式6，计算G’’₁到第1个类别的当前聚类中心p₁的隶属度μ₁(G₃)以及G’’₁到第2个类别的当前聚类中心p₂的隶属度μ₂(G₃)，得到2个隶属度，例如两个隶属度中最大的是隶属度μ₁(G₃)，那么就将G’’₁所对应的次区域（如图7中的子区域①、②、③、④构成的第1个次区域）分类到第1个类别，即，图像中与该次区域对应的部分被归类为第1个类别（沙漠纹理）。对4个新梯度特征向量分别进行上述计算，从而将4个次区域分别归类到第1个类别（沙漠纹理）或第二个类别（人造物纹理）。如此，可以进一步将属于未知类别的区域中的4个次区域分别归类到第1个类别（沙漠纹理）或第二个类别（人造物纹理）。从而将图像中与未知类别的区域对应的部分进一步被归类到第1个类别（沙漠纹理）或第二个类别（人造物纹理）。

如图11所示，图10中的标有“？”的一些区域分别被进一步归类为第1个类别（沙漠纹理）或第二个类别（人造物纹理），使得图像被完全归类为沙漠纹理和人造物纹理两个类别的纹理。并且，与图10相比，图11中的沙漠纹理和人造物纹理的边界线更加清晰准确。也就是说，对于具有混合纹理的区域，可以进一精确细分，使得分割更加精确。

另外，本发明还提供一种对图像进行纹理分割的装置，其特征在于，该图像分割装置使用如上所述的对图像进行纹理分割的方法。

本发明是一种新的基于局部的梯度特征向量和聚类分析的纹理分割方法。其特点是通过对局部区域的纹理用一个多维梯度特征向量来表征其特征。实际地物纹理单元通常具有不规则的方向，本发明通过提取梯度主方向并旋转将其作为0度方向，可以使提取的梯度特征向量具有方向不变性。由于使用的是梯度特征向量，所以能够忽略纹理上的亮度差异。在纹理边界处的区域存在混合类型的情况下，可以将区域细分并重新对次区域生成新梯度特征向量，并进行分类，从而能够精确细分具有混合纹理的区域，使得分割更加精确。

虽然本发明的特定实施例已被描述，但这些实施例只通过实例的方式进行表述，并不意欲限制本发明的范围。实际上，本文描述的创新方法可以通过各种其他形式实施；此外，也可以进行对本文描述的方法和系统的各种省略、替代和改变而不背离本发明的精神。附后的权利要求及其等同内容的目的是涵盖落入本发明的范围和精神内的这样的各种形式或修改。

Claims (15)

1.一种对图像进行纹理分割的方法，所述图像具有多种纹理，其特征在于，包括以下步骤：

a)将所述图像转换为灰度图像；

b)根据所述图像中最大的纹理单元的大小，将所述灰度图像划分为大小相同的多个区域，并提取与所述多个区域一一对应的多个梯度特征向量；

c)对提取的所述多个梯度特征向量进行模糊聚类分析，将所述多个区域进行分类，从而将所述图像中具有相同纹理的部分分类到相同的类别中。

2.如权利要求1所述的图像分割方法，其特征在于，在步骤b）中，所述多个区域中的每个区域的大小大于等于所述最大的纹理单元的大小。

3.如权利要求2所述的图像分割方法，其特征在于，在步骤b）中，通过以下步骤提取所述多个梯度特征向量：

b1)以水平向右为梯度方向角的0度方向，计算所述多个区域中的每个区域中的4m×4m个像素中的每个像素的梯度值和第一梯度方向角，m是大于等于1的整数；

b2)对计算的所述多个区域中的任一个区域中的4m×4m个像素各自的梯度值和梯度方向角进行统计，得到梯度直方图，在所述梯度直方图中，多个直方柱与多个第一方向角一一对应，将所述多个直方柱中具有峰值g_max的直方柱所对应的第一方向角作为所述任一个区域的梯度主方向角；

b3)以所述任一个区域的中心为中心，以所述梯度主方向角为所述梯度方向角的0度方向，限定与所述任一个区域大小相同的新区域；

b4)根据步骤b1）的计算结果，取出位于所述新区域中的4m×4m个像素中的每个像素的梯度值和第一梯度方向角，其中取出的每个像素的第一梯度方向角变更为第二梯度方向角，所述第二梯度方向角是从所述第一梯度方向角顺时针旋转到所述梯度主方向角的角度；

b5)将所述新区域均匀划分为16个子区域，每个子区域中具有m×m个像素，以所述梯度主方向角为所述0度方向，将0-360度方向均匀划分为多个第二方向角，如果任一个子区域中的任一个像素的第二梯度方向角不等于任一个第二方向角，则将所述任一个像素的梯度值分别分解到所述多个第二方向角中与所述任一个像素的第二梯度方向角相近的两个第二方向角上，得到与所述两个第二方向角分别对应的两个分解梯度值，如果任一个子区域中的任一个像素的第二梯度方向角等于任一个第二方向角，则将所述任一个像素的梯度值作为与所述任一个第二方向角相对应的分解梯度值，从而对于每个子区域，获得多个分解梯度值，所述多个分解梯度值中的任一个分解梯度值与所述多个第二方向角中的一个第二方向角相对应；

b6）将所述多个分解梯度值中与任一个第二方向角对应的数个分解梯度值进行累加，获得与所述任一个第二方向角对应的累加梯度值，从而对于所述每个子区域，获得与所述多个第二方向角一一对应的多个累加梯度值，将所述多个累加梯度值分别对于所述峰值g_max进行归一化，得到多个归一化累加梯度值；

b7）将所述新区域中16个子区域各自的多个归一化累加梯度值进行排列，获得所述任一个区域的梯度特征向量；

其中，对所述每个区域进行步骤b2)-b7)，从而提取与所述多个区域一一对应的所述多个梯度特征向量。

4.如权利要求3所述的图像分割方法，其特征在于，在步骤b1)中，根据以下等式计算每个像素的梯度值和第一梯度方向角，

等式1：

$g(x,y)=\sqrt{(P(x+1,y)-P{(x-1,y))}^2+P(x,y+1)-P{(x,y-1))}^2}$

等式2：

θ（x,y)＝αtan2((P(x,y+1)-P(x,y-1))/(P(x+1,y)-P(x-1,y)))

其中，αtan2表示反正切函数，g(x,y)表示位于所述灰度图像中坐标为(x,y)的像素点的梯度值，θ(x,y)表示位于所述灰度图像中坐标为(x,y)的像素点的第一梯度方向角，P(x,y)表示位于所述灰度图像中坐标为(x,y)的像素点的灰度值。

5.如权利要求4所述的图像分割方法，其特征在于，在步骤b2)中，所述梯度直方图的角度范围是0-360度，将0-360度均匀划分成所述多个第一方向角。

6.如权利要求5所述的图像分割方法，其特征在于，在步骤b4)中，当所述新区域包含超出所述灰度图像的超出范围时，将所述任一个区域中距离所述超出范围最近的多个像素各自的梯度值和第一梯度方向角取出作为所述超出范围中的多个像素各自的梯度值和第一梯度方向角。

7.如权利要求6所述的图像分割方法，其特征在于，在步骤b5)中，根据以下等式得到与所述两个第二方向角分别对应的两个分解梯度值，

等式3：g_θ1=g×sin(θ₂-θ)/sin(θ₂-θ₁)

等式4：g_θ2=g×sin(θ-θ₁)/sin(θ₂-θ₁)

其中，g表示所述任一个像素的梯度值，θ表示所述任一个像素的第二梯度方向角，θ₁和θ₂表示与θ相近的两个第二方向角，g_θ1表示g被分解到第二方向角θ₁上的分解梯度值，g_θ2表示g被分解到第二方向角θ₂上的分解梯度值。

8.如权利要求7所述的图像分割方法，其特征在于，在步骤b6)中，对于所述每个子区域，根据以下等式得到所述多个归一化累加梯度值，

等式5： $g_{\mathrm{ej}}=\frac{\underset{i=1}{\overset{m\×m}{\mathrm{\Σ}}}{g}_{\mathrm{eij}}}{g_{\max }}$

其中，g_ej表示在第e个子区域中的所述多个第二方向角中的第j个第二方向角上的归一化累加梯度值，g_eij表示在第e个子区域中的第i个像素在第j个第二方向角上的分解梯度值，e，j是大于1的整数，1≤e≤16。

9.如权利要求8所述的图像分割方法，其特征在于，在步骤b7）中，将所述新区域划分为大小相同的4个次区域，每个次区域包含4个所述子区域，根据步骤b6）中获得的多个归一化累加梯度值，按照逐个子区域逆时针方向，分别得到对应于4个次区域的4个次梯度特征向量G’₁、G’₂、G’₃、G’₄，其中，

G’₁=｛g₁₁、g₁₂、…g_1a、g₂₁、g₂₂、…g_2a、g₃₁、g₃₂、…g_3a、g₄₁、g₄₂、…g_4a｝，

G’₂=｛g₅₁、g₅₂、…g_5a、g₆₁、g₆₂、…g_6a、g₇₁、g₇₂、…g_7a、g₈₁、g₈₂、…g_8a｝，

G’₃=｛g₉₁、g₉₂、…g_9a、g₁₀₁、g₁₀₂、…g_10a、g₁₁₁、g₁₁₂、…g_11a、g₁₂₁、g₁₂₂、…g_12a｝，

G’₄=｛g₁₃₁、g₁₃₂、…g_13a、g₁₄₁、g₁₄₂、…g_14a、g₁₅₁、g₁₅₂、…g_15a、g₁₆₁、g₁₆₂、…g_16a｝，，

其中，a表示所述多个第二方向角的个数，且a是大于等于1的整数，1≤j≤a，

将所述次梯度特征向量G’₁、G’₂、G’₃、G’₄依次排列组成所述任一个区域的所述梯度特征向量G_r=｛g₁₁、g₁₂…g_16a｝，G_r表示与第r个区域对应的梯度特征向量，1≤r≤n，n表示所述多个区域的个数，且r，n是整数。

10.如权利要求9所述的图像分割方法，其特征在于，在步骤c）中，通过以下步骤进行所述模糊聚类分析：

c1)根据所述图像中的纹理种类，设置所述模糊聚类分析的C个类别，从所述多个梯度特征向量中任意取出一个梯度特征向量作为某个类别的当前聚类中心，得到C个当前聚类中心；

c2)计算所述多个梯度特征向量分别到每个当前聚类中心的隶属度，对于所述每个当前聚类中心，得到与所述多个梯度特征向量一一对应的多个隶属度，将所述多个隶属度中大于等于预定阈值t的隶属度所对应的梯度特征向量所对应的区域归类到与所述每个当前聚类中心所对应的类别，

c3)根据对于所述每个当前聚类中心计算得到的所述多个隶属度，计算与所述C个类别一一对应的C个更新聚类中心，并将所述C个更新聚类中心作为所述C个当前聚类中心；

c4）当所述多个梯度特征向量被归类的类别发生变化的比率大于预定的阈值v或重复步骤c2)和c3)的次数小于次数阈值N时，重复进行步骤c2），否则进行步骤c4）；

c5）将与所述C个当前聚类中心的隶属度都小于预定阈值t的梯度特征向量对应的区域归类到未知类别。

11.如权利要求10所述的图像分割方法，其特征在于，进一步包括如下步骤：

c6）将与归类到所述未知类别的某个区域所对应的4个次梯度特征向量G’₁、G’₂、G’₃、G’₄分别组成4个新梯度特征向量G’’₁、G’’₂、G’’₃、G’’₄，

G’’₁=｛g₁₁、g₁₂…g_4a，g₁₁、g₁₂…g_4a，g₁₁、g₁₂…g_4a，g₁₁、g₁₂…g_4a｝，

G’’₂=｛g₅₁、g₅₂、…g_8a，g₅₁、g₅₂、…g_8a，g₅₁、g₅₂、…g_8a，g₅₁、g₅₂、…g_8a｝，

G’’₃=｛g₉₁、g₉₂、…g_12a，g₉₁、g₉₂、…g_12a，g₉₁、g₉₂、…g_12a，g₉₁、g₉₂、…g_12a｝，

G’’₄=｛g₁₃₁、g₁₃₂、…g_16a，g₁₃₁、g₁₃₂、…g_16a，g₁₃₁、g₁₃₂、…g_16a，g₁₃₁、g₁₃₂、…g_16a｝，

c7）计算任一个新梯度特征向量分别到步骤c3)中的所述C个当前聚类中心的多个隶属度，并将所述任一个新梯度特征向量所对应的次区域分类到与所述多个隶属度中的最大隶属度对应的类别，从而将4个次区域各自分类到相应的类别中。

12.如权利要求11所述的图像分割方法，其特征在于，步骤c2)中，根据以下等式计算所述隶属度，

等式6： ${\mathrm{\μ}}_t\left(G_r\right)=\frac{{(1/{\left|\right|G_r-p_t\left|\right|}^2)}^{1/(b-1)}}{{\mathrm{\Σ}}_{k=1}^C{(1/{\left|\right|G_r-p_k\left|\right|}^2)}^{1/(b-1)}},$

其中，μ_t(G_r)表示第r个区域的梯度特征向量到第t个类别的当前聚类中心的隶属度，G_r表示所述多个区域中的第r个区域的梯度特征向量，p_t表示第t个类别的当前聚类中心，1≤t≤C，t为大于等于1的整数，b是模糊因子，用于控制聚类分析结果的模糊程度。

13.如权利要求12所述的图像分割方法，其特征在于，步骤c3)中，根据以下等式计算所述C个更新聚类中心中的第t个类别的更新聚类中心p_t’，

等式7： ${p_t}^{,}=\frac{{\mathrm{\Σ}}_{r-1}^n{\left[{\mathrm{\μ}}_t\left(G_r\right)\right]}^b{G}_r}{{\mathrm{\Σ}}_{r-1}^n{\left[{\mathrm{\μ}}_t\left(G_r\right)\right]}^b}.$

14.如权利要求3所述的图像分割方法，其特征在于，在步骤b3）中，以所述任一个区域的中心为中心，将所述任一个区域顺时针旋转一角度，所述角度为所述梯度主方向角，旋转后的任一个区域所包含的区域为所述新区域。

15.一种对图像进行纹理分割的装置，其特征在于，所述图像分割装置使用如权利要求1-14中任一项所述的对图像进行纹理分割的方法。