基于灰度共生矩阵的纹理方向精细估算方法

摘要

本发明公开了一种基于灰度共生矩阵的纹理方向精细估算方法，包括：S1：选取图像；将相对位置的方向θ和距离r进行等间隔划分，计算这些相对位置对应的GLCM；S2：计算参数Z的等角度和等距离间隔的分布；S3：计算参数Z的积分Z′；S4：搜索Z′的极小值，确定纹理方向等步骤。它解决了现有技术中纹理方向估算方法不够灵活，抗噪声能力差等技术问题。本发明具有抗噪声性能好，可估算多个纹理方向，且精确度高的优点。

权利要求

1.基于灰度共生矩阵的纹理方向精细估算方法，包括如下步骤：

S1：选取图像，将相对位置的方向和距离r分别进行等间隔划分，方向以角度θ表示，角度划分范围为-90°到90°，计算这些等角度和等距离间隔的相对位置对应的GLCM；

这些GLCM可通过下述三种方法中的任一一种可计算图像关于任意相对位置的GLCM计算方法得到：方法一：

将图像以-θ角度进行旋转，再由公式(1)从旋转后的图像中计算得到相对位置(r,θ)对应的图像的GLCM：其中 表示向下取整；G(m,n；r,θ)表示旋转前的图像的GLCM中第m行第n列的元素值；

G′(m,n；r′,θ′)表示旋转后的图像或包含了旋转后的图像的矩形区域的GLCM中第m行第n列的元素值，它可以由公式(2)计算得到；

其中“card”表示集合中的元素数目；表示矩形区域中的像素集合；f(x,y)表示位置坐标为(x,y)处的灰度值，这里位置以矩形坐标系的坐标表示，该坐标系以旋转前的图像的水平方向和垂直方向作为坐标系的两轴；

归一化因子Q′(r′,θ′)为：

它表示矩形区域中能与相对位置(r′,θ′)重合上的像素对的数目；

在计算过程中，需要去除包含了旋转后的图像的矩形区域在图像外的四个角；这可以通过如下标记来实现：假设图像中的灰度范围为1到N，令四个角区域的灰度值为N+1，然后通过公式(2)和(3)计算整个矩形区域的GLCM，那么旋转后的图像的GLCM即该GLCM的第1行到第N行第1列到第N列的矩阵块；在所述纹理方向精细估算方法中，还可以将GLCM的参数“相对位置两端的间距r”划分为像素间距的整数倍，那么公式(1)可简化为：G(m,n；r,θ)＝G′(m,n；r,0°)公式(4)

方法二：在待计算的GLCM对应的相对位置(r,θ)的附近找到一些相对位置(r1,θ1),…,(rK,θK)，这些相对位置与图像中的像素是可以重合的；如此，这些相对位置(r1,θ1),…,(rK,θK)对应的GLCM的矩阵元素，即G(m,n；r1,θ1),…,G(m,n；rK,θK)，便可以通过公式(5)和(6)计算得到，而待计算的GLCM的矩阵元素G(m,n；r,θ)则可由这些GLCM的矩阵元素插值得到；

其中G(m,n；rk,θk)表示相对位置(rk,θk)对应的图像的GLCM中第m行第n列的元素值；

“card”表示集合中的元素数目；表示图像的像素集合；f(x,y)表示位置坐标为(x,y)处的灰度值，这里位置以矩形坐标系的坐标表示，该坐标系以图的水平方向和垂直方向作为坐标系的两轴；

归一化因子Q(rk,θk)为：

方法三：相对位置(r,θ)对应的图像的GLCM中第m行第n列的矩阵元素可由下述公式计算得到：其中，G(m,n；r,θ)表示相对位置(r,θ)对应的GLCM中第m行第n列的元素值；“card”表示集合中的元素数目；表示图像的像素集合；表示图像所在区域；f(i,j)表示位置坐标为(i,j)处的像素灰度值，f(i+rcosθ,j+rsinθ)由(i+rcosθ,j+rsinθ)邻近像素的灰度值插值和取整得到；这里位置以矩形坐标系的坐标表示，该坐标系以图的水平方向和垂直方向作为坐标系的两轴；

这里归一化因子Q(r,θ)为

S2：假设图像的灰度范围为1到N；通过公式(9)或(10)计算这些等角度和等距离间隔的相对位置对应的GLCM的Z(r,θ)，得到该图像的Z(r,θ)的等角度和等距离间隔的分布：其中w(m,n)是关于GLCM的矩阵元素在矩阵中的位置和矩阵对角线之间距离的递增函数；显然，较大的矩阵元素越集中于GLCM的对角线附近，那么Z(r,θ)的值就越小；

如果采用(m-n)2作为w(m,n)，那么公式(9)就变为：

此即对比度；

S3：计算参数Z的积分Z′：通过公式(11)来计算Z′(θ)；

Z′(θ)＝∫Z(r,θ)dr公式(11)

S4：搜索Z′(θ)的极小值，确定纹理方向：

Z′(θ)的极小值所对应的方向即图像中的纹理方向。

2.根据权利要求1所述的基于灰度共生矩阵的纹理方向精细估算方法，其特征在于：在步骤S1和S2之间，还包括如下步骤：将通过步骤S1计算得到的GLCM，加上该矩阵的转置，再将矩阵和的所有矩阵元素都除以2，得到对称的GLCM，将该对称的GLCM用于步骤S2中。

3.根据权利要求1所述的基于灰度共生矩阵的纹理方向精细估算方法，其特征在于：在步骤S3中，通过公式(13)来计算Z′(θ)；

其中，L表示距离r划分的数目，rl表示第l个距离。

说明书

基于灰度共生矩阵的纹理方向精细估算方法

技术领域

[0001] 本发明涉及一种纹理方向的估算方法，尤其涉及一种基于灰度共生矩阵(GLCM，Gray-level co-occurrence matrix)的纹理方向精细估算方法。

背景技术

[0002] 纹理能从许多可视表面中被人的视觉自然地识别出来，并给人一些特别的感觉，如方向感、周期感和粗糙感等。作为可视表面的内在属性，纹理在众多领域都是非常重要的研究课题。许多围绕纹理的研究工作都是在探索如何更好地提取纹理的特征来对其进行客观描述，GLCM便是其中一种发展起来的统计分析工具。

[0003] 传统的GLCM本质上是满足特定相对位置关系的像素对的灰度值联合概率分布。该矩阵可以通过对图像中的像素对进行统计得到。从该矩阵可以导出图像的二阶统计参数，这些参数分别表征了图像纹理的不同特性。有研究表明纹理和纹理之间的区别极大地依赖于纹理的二阶统计上的差异，因此由GLCM导出的二阶统计参数对于图像纹理的识别、分类等而言都是非常重要的。

[0004] 国外研究显示纹理的方向性在人的纹理感知中扮演了极其重要的角色，是人用来进行纹理识别和分类的一个重要特征。纹理方向估算在遥感图像分析中具有重要的应用价值，因为遥感图像中的纹理方向往往暗示了地貌或一些物理现象的相关信息。例如从合成孔径雷达遥感图像中获取海面风向信息的方法便是基于图像纹理的方向估算。海面风场引起的条纹以及海气边界层涡旋会在合成孔径雷达遥感图像中产生具有方向性的纹理，该纹理的方向通常与风向一致，因此便可以通过估算纹理方向来获得海面风向。除了海面风场之外，其他一些海洋现象也会在合成孔径雷达或光学遥感图像中产生具有方向性的纹理，如海浪、内波及沙波等。除遥感领域外，纹理方向估算在许多其他领域也有着广泛而重要的应用价值，如地震图像解译，医学成像分析以及工业中机器加工表面的评估等。

[0005] 基于傅立叶变换的方法是纹理方向估算方面一种较为典型的方法，但该方法通常需要一块包含许多纹理样本的矩形区域数据来精确地计算纹理方向估算所需的图像傅立叶谱。基于灰度梯度的方法是另一种比较典型的方法，但是该方法需要对图像进行差分运算，而差分运算是容易受噪声影响的；另外，在图像灰度的脊部和谷部区域，由于灰度梯度的幅度非常小而难以准确计算，因而需要增加额外的操作。

[0006] 从GLCM中提取的部分参数对于纹理方向是敏感的。但是传统GLCM图像纹理分析受图像像素位置的限制，一般只考虑特定方向上的GLCM(通常为图像的水平方向、垂直方向和对角线方向)，因而只能非常粗糙地在这些方向上估计纹理的走势，无法精确地估算纹理的方向。因此传统GLCM及其导出的参数通常不应用于纹理方向的精细估算。

发明内容

[0007] 本发明针对现有技术中，纹理方向估算方法不够灵活，抗噪声能力差等技术问题，提供一种基于灰度共生矩阵的纹理方向精细估算方法，该方法抗噪声性能好，可估算多个纹理方向，且精确度高。

[0008] 为此，本发明采用如下技术方案：

[0009] 一种基于灰度共生矩阵的纹理方向精细估算方法，包括如下步骤：

[0010] S1：选取图像，将相对位置的方向和距离(相对位置两端的间距)r分别进行等间隔划分，方向以角度θ表示，角度划分范围为-90°到90°，计算这些等方向(角度)和等距离间隔的相对位置对应的GLCM；

[0011] 这些GLCM可通过下述三种方法中的任一一种可计算图像关于任意相对位置的GLCM计算方法得到：

[0012] 方法一：

[0013] 将图像以-θ角度进行旋转，再由公式(1)从旋转后的图像中计算得到相对位置(r,θ)对应的图像(旋转前的图像)的GLCM：

[0014]

[0015] 其中 表示向下取整；G(m,n；r,θ)表示旋转前的图像的GLCM中第m行第n列的元素值；G′(m,n；r′,θ′)表示旋转后的图像或包含了旋转后的图像的矩形区域的GLCM中第m行第n列的元素值(从后文可以看到：旋转后的图像的GLCM即包含了旋转后的图像的矩形区域的GLCM的第1行到第N行第1列到第N列的矩阵块)，它可以由公式(2)计算得到；

[0016]

[0017] 其中“card”表示集合中的元素数目；P′表示矩形区域中的像素集合；f(x,y)表示位置坐标为(x,y)处的灰度值，这里位置以矩形坐标系的坐标表示，该坐标系以旋转前的图像的水平方向和垂直方向作为坐标系的两轴；

[0018] 归一化因子Q′(r′,θ′)为：

[0019]

[0020] 它表示矩形区域中能与相对位置(r′,θ′)重合上的像素对的数目。

[0021] 本发明中，提及的“card”均表示集合中的元素数目。

[0022] 在计算过程中，需要去除包含了旋转后的图像的矩形区域在图像外的四个角；这可以通过如下简单的标记来实现：不失一般性地，假设图像中的灰度范围为1到N，令四个角区域的灰度值为N+1，然后通过公式(2)和(3)计算整个矩形区域的GLCM，那么旋转后的图像的GLCM即该GLCM的第1行到第N行第1列到第N列的矩阵块；在本发明给出的精细估算纹理方向的方法中，为简便起见，可以将GLCM的参数“相对位置两端的间距r”划分为像素间距的整数倍，那么公式(1)可简化为：

[0023] G(m,n；r,θ)＝G′(m,n；r,0°)  公式(4)

[0024] 方法二：在待计算的GLCM对应的相对位置(r,θ)的附近找到一些相对位置(r1，θ1)，…，(rK，θK)，这些相对位置与图像中的像素是可以重合的；如此，这些相对位置(r1，θ1)，…，(rK，θK)对应的GLCM的矩阵元素，即G(m，n；r1，θ1)，…，G(m，n；rK，θK)，便可以通过公式(5)和(6)计算得到，而待计算的GLCM的矩阵元素G(m,n；r,θ)则可由这些GLCM的矩阵元素插值得到；

[0025]

[0026] 其中G(m，n；rk，θk)表示相对位置(rk，θk)对应的图像的GLCM中第m行第n列的元素值； 表示图像的像素集合；f(x,y)表示位置坐标为(x,y)处的灰度值，这里位置以矩形坐标系的坐标表示，该坐标系以图的水平方向和垂直方向作为坐标系的两轴；

[0027] 归一化因子Q(rk，θk)为：

[0028]

[0029] 方法三：相对位置(r,θ)对应的图像的GLCM中第m行第n列的矩阵元素可由下述公式计算得到：

[0030]

[0031] 其中， 表示图像的像素集合； 表示图像所在区域；f(i,j)表示位置坐标为(i,j)处的像素灰度值，f(i+rcosθ,j+rsinθ)由(i+rcosθ,j+rsinθ)邻近像素的灰度值插值和取整得到；这里位置以矩形坐标系的坐标表示，该坐标系以图的水平方向和垂直方向作为坐标系的两轴；

[0032] 这里归一化因子Q(r,θ)为

[0033]

[0034] S2：不失一般性，假设图像的灰度范围为1到N；通过公式(9)或(10)计算这些等角度和等距离间隔的相对位置对应的GLCM的Z(r,θ)，得到该图像的Z(r,θ)的等角度和等距离间隔的分布：

[0035]

[0036] 其中w(m,n)是关于GLCM的矩阵元素在矩阵中的位置和矩阵对角线之间距离的递增函数；显然，较大的矩阵元素越集中于GLCM的对角线附近，那么Z(r,θ)的值就越小；

[0037] 如果采用(m-n)2作为w(m,n)，那么公式(9)就变为：

[0038]

[0039] 此即对比度；

[0040] S3：计算参数Z的积分(求和)Z′：通过公式(11)来计算Z′(θ)；

[0041] Z′(θ)＝∫Z(r,θ)dr  公式(11)

[0042] S4：搜索Z′(θ)较为明显的极小值，确定纹理方向：

[0043]

[0044] Z′(θ)的极小值所对应的方向即图像中的纹理方向。

[0045] 作为一种替代方案，在步骤S1和S2之间，还包括如下步骤：将通过步骤S1计算得到的GLCM，加上该矩阵的转置，再将矩阵和的所有矩阵元素都除以2，得到对称的GLCM，将该对称的GLCM用于步骤S2中。

[0046] 作为一种具体的实施方案：在步骤S3中，可通过公式(13)来计算Z′(θ)；

[0047]

[0048] 其中，L表示距离r划分的数目，rl表示第l个距离。

[0049] 本发明的有益效果在于：如前文所述，基于傅立叶变换的方法和基于灰度梯度的方法是目前较为典型的两种纹理方向精细估算方法，两者分别从图像的傅立叶谱和灰度梯度出发。基于傅立叶变换的方法需要包含许多纹理样本的矩形区域数据，才能准确地提取出估算纹理方向所需的图像傅立叶谱，这使得该方法在估算非矩形区域或小区域的纹理的方向时不太灵活。基于灰度梯度的方法需要对图像进行差分运算，而差分运算对噪声是比较敏感的；另外，在图像灰度的脊部和谷部区域，由于灰度梯度的幅度非常小变得难以准确计算，因而需要增加额外的操作来估算纹理方向。而本发明从图像的GLCM出发来精细估算纹理方向，在原理上与现有技术中的两种典型的纹理方向精细估算方法是完全不同的，因而不存在上述问题，具有非常良好的噪声抗性，可以估算多个纹理方向，且精确度高，是精细估算纹理方向的一种全新的方法。

附图说明

[0050] 图1为本发明步骤S1中方法一的示意图；

[0051] 图2为本发明步骤S1中方法三的示意图；

[0052] 图3-1至图3-6为几种纹理图像和采用本发明的纹理方向精细估算方法估算的纹理方向；

[0053] 图4-1至图4-6给出了分别采用本发明步骤S1中的方法一至方法三得到的图3-1在取纹理方向和纹理方向的垂直方向作为相对位置的方向时的GLCM的元素值分布；

[0054] 图5-1至图5-3给出了由公式(10)计算得到的图3-1的Z(r，θ)的分布，图中纵轴表示r(单位：像素间距)，横轴表示角度θ；

[0055] 图6给出了图3-1对应的Z′(θ)曲线，图中，纵轴表示Z′，横轴表示角度θ；

[0056] 图7为本发明的流程示意图；

[0057] 图8-1至图8-6为采用本发明的方法对引入了加性高斯噪声后的图像的纹理方向的估算结果；

[0058] 图9-1至图9-6为采用本发明的方法对引入了椒盐噪声后的图像的纹理方向的估算结果；

[0059] 图10-1至图10-6为采用本发明的方法对引入了斑点噪声(乘性噪声)后的图像的纹理方向的估算结果。

具体实施方式

[0060] 为了使本技术领域的人员更好地理解本发明方案，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整的描述。

[0061] 实施例1：

[0062] 本发明的基于灰度共生矩阵的纹理方向精细估算方法，包括如下步骤：

[0063] S1：选取图像，将相对位置的方向和距离(相对位置两端的间距)r分别进行等间隔划分，其中，方向以角度θ表示，角度划分范围为-90°到90°，计算这些等角度和等距离间隔的相对位置对应的GLCM；

[0064] 本发明的纹理方向精细估算方法需要计算这些在角度和距离上等间隔的不同相对位置对应的GLCM；其中“角度θ”为相对位置中的方向与图像水平方向之间的夹角，“距离r”为相对位置两端的间距；为此，本发明首先提出三种可用于计算任意相对角度和距离的GLCM的方案。

[0065] 这里以非对称的GLCM为例，而对称的GLCM仅仅只是“非对称的GLCM与它的转置的和，再将矩阵和的所有矩阵元素都除以2”，因此类似。图像可以看作为某种场在像素点位置上的采样值的集合，那么一个非对称的GLCM本质上可以看作是满足特定相对位置的“位置对”上的场值(图像中用灰度值标定)的联合概率分布。在传统的GLCM分析中，仅用到了那些相对位置可以和图像的像素重合的GLCM。

[0066] 通过本发明提出的下述三种方法中的任一一种可计算图像关于任意相对位置的GLCM的计算方法来计算S1中等角度和等距离间隔的这些相对位置对应的GLCM：

[0067] 方法一，如图1所示：

[0068] 将图像以-θ角度进行旋转，再由公式(1)从旋转后的图像中计算得到相对位置(r,θ)对应的图像(旋转前的图像)的GLCM：

[0069]

[0070] 其中 表示向下取整；G(m,n；r,θ)表示旋转前的图像的GLCM中第m行第n列的元素值，G′(m,n；r′,θ′)表示旋转后的图像或包含了旋转后的图像的矩形区域(如图1所示)的GLCM中第m行第n列的元素值(从后文可以看到：旋转后的图像的GLCM即包含了旋转后的图像的矩形区域的GLCM的第1行到第N行第1列到第N列的矩阵块)，它可以由公式(2)计算得到；

[0071]

[0072] 其中 表示矩形区域中的像素集合；f(x,y)表示位置坐标为(x,y)处的灰度值，这里位置以矩形坐标系的坐标表示，该坐标系以旋转前的图像的水平方向和垂直方向作为坐标系的两轴；

[0073] 归一化因子Q′(r′,θ′)为：

[0074]

[0075] 它表示矩形区域中能与相对位置(r′,θ′)重合上的像素对的数目；

[0076] 在计算过程中，需要去除包含了旋转后的图像的矩形区域在图像外的四个角；这可以通过如下简单的标记来实现：不失一般性地，假设图像中的灰度范围为1到N，令四个角区域的灰度值为N+1，然后通过公式(2)和(3)计算整个矩形区域的GLCM，那么旋转后的图像的GLCM即该GLCM的第1行到第N行第1列到第N列的矩阵块；在本发明给出的精细估算纹理方向的方法中，为简便起见，可以将GLCM的参数“相对位置两端的间距r”划分为像素间距的整数倍，那么公式(1)可简化为：

[0077] G(m,n；r,θ)＝G′(m,n；r,0°)   公式(4)

[0078] 方法二：在待计算的GLCM对应的相对位置(r,θ)的附近找到一些相对位置(r1，θ1)，…，(rK，θK)，这些相对位置与图像中的像素是可以重合的；如此，这些相对位置(r1，θ1)，…，(rK，θK)对应的GLCM的矩阵元素，即G(m，n；r1，θ1)，…，G(m，n；rK，θK)，便可以通过公式(5)和(6)计算得到，而待计算的GLCM的矩阵元素G(m,n；r,θ)则可由这些GLCM的矩阵元素插值得到；

[0079]

[0080] 其中G(m，n；rk，θk)表示相对位置(rk，θk)对应的图像的GLCM中第m行第n列的元素值； 表示图像的像素集合；f(x,y)表示位置坐标为(x,y)处的灰度值，这里位置以矩形坐标系的坐标表示，该坐标系以图的水平方向和垂直方向作为坐标系的两轴；

[0081] 归一化因子Q(rk，θk)为：

[0082]

[0083] 方法三，如图2所示：

[0084] 相对位置(r,θ)对应的图像的GLCM中第m行第n列的矩阵元素可由下述公式计算得到：

[0085]

[0086] 其中， 表示图像的像素集合； 表示图像所在区域；f(i,j)表示位置坐标为(i,j)处的像素灰度值，f(i+rcosθ,j+rsinθ)由(i+rcosθ,j+rsinθ)邻近像素的灰度值插值和取整得到；这里位置以矩形坐标系的坐标表示，该坐标系以图的水平方向和垂直方向作为坐标系的两轴；

[0087] 这里归一化因子Q(r,θ)为

[0088]

[0089] 为清楚起见，图2示出了方法三中待统计的“位置对”的几何说明，图中：圆点表示图像的像素位置，矩形区域ABCD即图像所在区域 方点表示按照相对位置关系(r,θ)对像素位置(圆点)进行平移以后的位置；矩形区域AEFG中的圆点即公式(7)和公式(8)中待统计的满足相对位置(r,θ)的“位置对”的起始端位置；矩形区域A'E'F'G'中的方点即公式(7)和公式(8)中待统计的满足相对位置(r,θ)的“位置对”的末尾端位置。

[0090] 类似地，也可以令相对位置中的末尾端与图像像素重合而起始端上的值由像素上的值插值和取整得到，或者两端上的值都由像素上的值插值和取整得到。这两种方式本质上和方法二的方式没有区别，因此不再赘述。

[0091] S2：通过公式(9)或(10)计算这些等角度和等距离间隔的GLCM的Z(r,θ)，得到该图像的Z(r,θ)的等角度和距离间隔的相对位置对应的分布：

[0092] 灰度沿着纹理方向变化是最小的。这可以反映到图像的GLCM上。当相对位置的方向与纹理方向一致时，该相对位置对应的GLCM中值较大的元素集中在该矩阵的对角线附近。极端地，对于一个理想的条纹图像(即在纹理方向上没有灰度值变化的图像)，其与纹理方向一致的那些相对位置对应的GLCM中非零元素仅存在于矩阵的对角线上。

[0093] 作为一个例子，图3-1显示了一幅典型的纹理图像；而图4-1至图4-6给出了分别采用前面三种GLCM计算方法得到的图3-1在取纹理方向和纹理方向的垂直方向作为相对位置的方向时的GLCM的元素值分布，其中相对位置中的距离r取10倍像素间距。

[0094] 图4-1由GLCM计算方法一得到的沿纹理方向的GLCM；图4-2由GLCM计算方法二得到的沿纹理方向的GLCM，图4-3由GLCM计算方法三得到的沿纹理方向的GLCM，图4-4由GLCM计算方法一得到的与纹理方向垂直的方向的GLCM，图4-5由GLCM计算方法二得到的与纹理方向垂直的方向的GLCM，图4-6由GLCM计算方法三得到的与纹理方向垂直的方向的GLCM。

[0095] GLCM的这一特性可以被利用于图像纹理方向的精细估算。首先这一特性需要被定量地参数化。为此，定义如下参数来衡量GLCM中元素关于对角线的集中程度；不失一般性，假设图像的灰度范围为1到N。

[0096]

[0097] 其中w(m,n)是关于GLCM的矩阵元素在矩阵中的位置和矩阵对角线之间距离的递增函数；显然，较大的矩阵元素越集中于GLCM的对角线附近，那么Z(r,θ)的值就越小；如果采用(m-n)2作为w(m,n)，那么公式(9)就变为：

[0098]

[0099] 此即对比度；

[0100] S3：计算参数Z的积分Z′：通过公式(11)来计算Z′(θ)；

[0101] Z′(θ)＝∫Z(r,θ)dr  公式(11)

[0102] S4：搜索Z′(θ)较为明显的极小值，确定纹理方向：

[0103]

[0104] Z′(θ)的极小值所对应的方向即图像中的纹理方向。对于Z′(θ)有小幅波动起伏干扰的情况，可以先对Z′(θ)进行平滑去除小幅波动起伏干扰后，再搜索其明显的极小值确定纹理方向。步骤S2-S4的详细说明如下：

[0105] 在采用前面给出的计算GLCM的三种方法计算了GLCM后，Z(r,θ)的等角度(方向)和等距离间隔的分布图就可以通过公式(9)或(10)得到。图5-1至图5-3分别给出了由步骤S1中三种GLCM计算方法得到的等方向(角度)和等距离间隔的相对位置对应的GLCM通过公式(10)计算得到的图3-1的Z(r,θ)的分布。图中可以看到：Z(r,θ)的值在纹理方向(即与图3-1水平轴成近90°的方向)上总体较小，因此沿着纹理方向的GLCM相比其他方向的GLCM而言，值较大的矩阵元素总体上更加集中于矩阵对角线附近。图3-1中的纹理在其他方向(非纹理方向)上有近似的周期性。因此当两个位置之间的距离接近于近似的周期时，这两个位置上的灰度值差异是很小的。这导致了Z(r,θ)关于θ的曲线(保持r为恒值)在非纹理方向上也有极小值出现。因此，不能仅通过搜索Z(r,θ)的极小值来估算纹理方向。为了解决这一问题，本发明采用了Z(r,θ)关于r的积分(以Z′(θ)表示)来替代Z(r,θ)。

[0106] Z′(θ)＝∫Z(r,θ)dr  公式(11)

[0107] 实际应用中，可使用Z′(θ)的数值积分(离散求和)的形式，例如

[0108]

[0109] 其中，L表示距离r划分的数目，rl表示第l个距离。

[0110] 采用Z(r,θ)的积分(求和)而不采用Z(r,θ)的另一个好处在于：Z(r,θ)的积分(求和)可以看作是Z(r,θ)在距离区间内的综合结果，因而随机波动更小。

[0111] Z′(θ)较为明显的极小值对应的方向可以从Z′(θ)中搜索出来作为对图像纹理方向的精细估算结果。

[0112]

[0113] 图6给出了图3-1对应的Z′(θ)曲线，其中： 表示在纹理方向估算过程中步骤S1采用了方法一的GLCM计算方案， 表示在纹理方向估算过程中步骤S1采用了方
法二的GLCM计算方案， 表示在纹理方向估算过程中步骤S1采用了方法三的GLCM
计算方案。图中可以看到图3-1仅含一个主要的纹理方向。由此估算的纹理方向在图3-1中作了标记，可以看到估算的方向与图像的纹理方向非常一致。

[0114] 图7示出了本发明的基于灰度共生矩阵的纹理方向精细估算方法的流程。

[0115] 作为测试，对图3-2至图3-6中其他纹理的方向也用本发明的纹理方向精细估算方法进行了提取，并在图3-2至图3-6中作了标注。同样地，在图中： 表示在纹理方向估算过程中步骤S1采用了方法一的GLCM计算方案， 表示在纹理方向估算过程中步
骤S1采用了方法二的GLCM计算方案， 表示在纹理方向估算过程中步骤S1采用了方
法三的GLCM计算方案。可以看到通过本发明的估算方法能够准确地提取出纹理的方向。

[0116] 为了验证本发明的纹理方向精细估算方法对各类噪声的良好抗性，在图3-1至图3-6中分别引入了大的加性高斯噪声、椒盐噪声及斑点噪声(乘性噪声)。图8-1至图8-6为采用本发明的方法对引入了加性高斯噪声(信噪比为-10dB)后的图像的纹理方向的估算结果。图9-1至图9-6为采用本发明的方法对引入了椒盐噪声(被椒盐噪声污染的像素占图像的50％)后的图像的纹理方向的估算结果。图10-1至图10-6为采用本发明的方法对引入了斑点噪声(乘性噪声)后的图像的纹理方向的估算结果。斑点噪声模型如下

[0117] f′(i,j)＝(1+α)·f(i,j)  公式(14)

[0118] 其中f(i,j)和f′(i,j)分别表示原图像和引入了斑点噪声后的图像在像素(i,j)处的灰度值。α为均值为0，方差为2的均匀分布的随机噪声。可以看到图8-1至图8-6、图9-1至图9-6以及图10-1至图10-6中的纹理方向估算结果和图3-1至图3-6的纹理方向估算结果基本一致，所引入的大的加性高斯噪声、椒盐噪声及斑点噪声并没有在本方法估算图中纹理方向时造成明显的影响。

[0119] 同样地，在图8-1至图8-6、图9-1至图9-6以及图10-1至图10-6中： 表示在纹理方向估算过程中步骤S1采用了方法一的GLCM计算方案， 表示在纹理方向估算
过程中步骤S1采用了方法二的GLCM计算方案， 表示在纹理方向估算过程中步骤S1
采用了方法三的GLCM计算方案。

[0120] 实施例2：

[0121] 上文所提及的GLCM为不对称的GLCM。不对称的GLCM加上该矩阵的转置，再将矩阵和的所有矩阵元素都除以2，得到对称的GLCM；

[0122] 本实施例与实施例1的区别之处在于：在步骤S1和S2之间，还包括如下步骤：将通过步骤S1计算得到的GLCM，加上该矩阵的转置，再将矩阵和的所有矩阵元素都除以2，得到对称的GLCM，将该对称的GLCM用于步骤S2中替代不对称的GLCM。

[0123] 显然，所描述的实施例仅仅是本发明的一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都应当属于本发明保护的范围。