灰度共生矩阵（GLCM）原理与Python实现：从纹理分析到特征提取

1. 从“看山是山”到“看山是纹理”：为什么需要灰度共生矩阵

我们每天都在处理图像，无论是手机拍照、刷短视频，还是做设计、搞科研。很多时候，我们评价一张图片，会说“这张图很清晰”、“那张图很模糊”，或者“这个材质看起来很有质感”。这里的“质感”，在图像处理领域，有一个更专业的词，叫做纹理。

纹理是什么？它不是颜色，也不是亮度，而是像素之间的一种空间排列规律。想象一下，一块光滑的大理石板和一块粗糙的砂纸，即使它们的平均亮度相同，我们也能一眼分辨出来，靠的就是纹理。再比如，一片森林的卫星图像和一片沙漠的，颜色可能相近，但森林的像素点分布杂乱、有规律，沙漠则相对均匀，这种差异也是纹理。

那么，计算机如何“看见”并量化这种纹理呢？这就是灰度共生矩阵要解决的问题。GLCM，全称Gray-Level Co-occurrence Matrix，是上世纪70年代由Haralick等人提出的一种经典纹理分析方法。它的核心思想非常直观：不孤立地看单个像素的灰度值，而是看一对像素（共生对）的灰度值组合，在特定空间关系下出现的频率。

这就像我们看一幅点彩画。如果只看一个色点，你只能知道它是红色还是蓝色。但如果你统计“在它右边1厘米处出现另一个红色点”这种情况在整个画布上出现了多少次，你就能量化出这幅画在水平方向上的“红色连续性”或“红色聚集程度”。GLCM做的就是这个统计工作，只不过把“色点”换成了“灰度级”，把“1厘米”换成了特定的距离和方向。

为什么这个方法历经半个世纪依然经典？因为它将抽象的“纹理感觉”转化为了一个具体的、可计算的矩阵。从这个矩阵中，我们可以进一步提取出诸如“对比度”、“能量”、“熵”、“同质性”等十几个量化特征。这些特征，就是让计算机理解图像纹理的“语言”。

在遥感领域，GLCM特征用于区分林地、农田、水体；在医学影像中，它帮助鉴别肿瘤组织的纹理异质性；在工业检测里，它用来发现产品表面的划痕、瑕疵；甚至在艺术分析中，也能辅助鉴定画作的笔触风格。可以说，只要你的问题涉及到“模式”、“结构”、“粗糙度”，GLCM很可能就是一个有力的工具。

很多人初次接触GLCM会觉得公式复杂、概念抽象，从而望而却步。其实，它的内核非常朴素。接下来，我将抛开复杂的数学外衣，用最直白的思路和可运行的Python代码，带你从零构建GLCM，并亲手计算出那些强大的纹理特征。你会发现，它不过是一个“数数”的游戏，但数出来的结果，却能揭示图像的深层秘密。

2. 拆解GLCM：一个“数对子”的统计游戏

理解GLCM，最关键的一步是忘掉矩阵，先理解“共生”这个概念。我们通过一个极度简化的例子来建立直觉。

假设我们有一张4x4的微型图像，其灰度值已经被量化为只有0, 1, 2三个级别（0最黑，2最白）。图像数据如下：

1, 1, 0, 0 1, 1, 0, 0 2, 2, 1, 1 2, 2, 1, 1

现在，我们定义一种空间关系：一个像素和它正右侧的像素。用参数表示就是：距离d=1，方向θ=0°（水平向右）。我们要为这种关系建立一个GLCM。

怎么建？就是遍历图像中的每一个像素（作为“参考像素”），找到它右边距离为1的那个像素（作为“邻居像素”），记录下这一对像素的灰度值(i, j)，其中i是参考像素的灰度值，j是邻居像素的灰度值。然后，在整个图像上统计每一对(i, j)出现的次数。

让我们手动数一下：

• 从左上角(0,0)的像素1开始，它右边的像素是1，所以我们记录到一对(1, 1)。
• 接着是(0,1)的像素1，右边是0，记录(1, 0)。
• (0,2)的像素0，右边是0，记录(0, 0)。
• 第一行最后一个像素(0,3)的0右边没有像素了，跳过。
• 第二行同理，会得到(1,1),(1,0),(0,0)。
• 第三行：(2,2)的2右边是2，记录(2,2)；2右边是1，记录(2,1)；1右边是1，记录(1,1)。
• 第四行同理，得到(2,2),(2,1),(1,1)。

现在，我们来统计所有可能的(i, j)对出现的次数。因为灰度级有0,1,2，所以我们的矩阵是3x3的。统计结果如下：

• (0,0)出现了 2次（来自第一、二行的第三个位置）
• (0,1)出现了 0次
• (0,2)出现了 0次
• (1,0)出现了 2次（来自第一、二行的第二个位置）
• (1,1)出现了 4次（来自第一、二行的第一个位置，第三、四行的第三个位置）
• (1,2)出现了 0次
• (2,0)出现了 0次
• (2,1)出现了 2次（来自第三、四行的第二个位置）
• (2,2)出现了 2次（来自第三、四行的第一个位置）

把这个统计结果填入矩阵，其中行索引i代表参考像素灰度值，列索引j代表邻居像素灰度值，就得到了我们的GLCM：

i\j	0	1	2
0	2	0	0
1	2	4	0
2	0	2	2

这个矩阵就是GLCM。它的大小由图像的最大灰度级决定。矩阵中的每个元素P(i, j)，表示在定义的方位关系下，灰度值为i和j的像素对出现的次数。

注意：在实际应用中，我们通常会将这个“次数矩阵”转换为“概率矩阵”，即每个元素除以所有元素之和，使得矩阵所有元素之和为1。这消除了图像大小对统计量的影响，便于不同图像间的比较。上面的矩阵总和是12，所以概率矩阵中P(1,1)=4/12≈0.333。

空间关系的多样性：我们刚才只用了(d=1, θ=0°)这一种关系。实际上，d和θ是可以变化的。常见的θ有0°, 45°, 90°, 135°四个方向。对于每个方向，我们都可以计算出一个GLCM。因此，描述一个图像的纹理，我们往往会得到一组GLCM（例如4个方向各一个）。

矩阵的对称性：在上面的统计中，我们只考虑了“参考像素→邻居像素”这一个顺序。但有时我们也会考虑“邻居像素→参考像素”，即(j, i)对。更常见的做法是，在统计时同时考虑这两种顺序（即无方向性），这样得到的GLCM是对称矩阵。具体实现时，可以在统计完(i, j)后，也累加一次(j, i)到矩阵中（或者直接使用P(i, j) + P(j, i)）。很多库的默认设置就是生成对称的GLCM。

理解了这个“数对子”的过程，GLCM就不再神秘。它本质上是一个条件概率的统计表，反映了在某种空间构型下，某种灰度值组合出现的可能性。图像纹理越粗糙、对比越强烈，GLCM的数值就越容易集中在主对角线两侧；图像纹理越均匀、平滑，GLCM的数值就越容易集中在主对角线上。

3. 从矩阵到特征：Haralick特征的物理意义与计算

得到了GLCM概率矩阵P(i, j)之后，我们面对的是一个Ng x Ng的矩阵（Ng是灰度级数）。直接使用这个矩阵作为特征维度太高，且不直观。因此，Haralick等人定义了14个从GLCM中推导出的标量特征，用来量化纹理的不同属性。这里我们重点讲解最常用、物理意义最明确的5个。

在计算之前，有一个至关重要的预处理步骤：灰度级压缩。如果原图是8位灰度图，灰度级有256级，那么GLCM就是256x256的矩阵，不仅计算量大，而且矩阵会非常稀疏（很多位置是0），统计不稳定。通常，我们会通过线性或非线性变换，将灰度级压缩到更少的级别，例如8级、16级或32级。这是GLCM应用中非常关键且容易被忽略的一步，直接影响到特征的稳定性和区分度。

假设我们已经得到了一个归一化的（元素之和为1）、灰度级为Ng=8的GLCM概率矩阵P。下面我们来计算特征，我会给出数学公式和对应的Python实现思路。

3.1 对比度

公式：Contrast = Σ Σ |i-j|² * P(i, j)物理意义：衡量图像的清晰度和纹理的沟壑深浅。对比度值大，说明纹理反差大，图像局部变化剧烈，视觉效果清晰、棱角分明。想象一下锐利的岩石边缘与平滑的水面，岩石的对比度特征值会远大于水面。计算解读：|i-j|是像素对的灰度差。这个公式对灰度差大的像素对赋予了更大的权重（平方）。因此，如果图像中有很多亮度差异很大的像素紧挨着（比如黑白棋盘格），那么对比度值就会非常高。Python计算核心：

contrast = 0.0 for i in range(Ng): for j in range(Ng): contrast += ((i - j) ** 2) * P[i, j]

3.2 能量

公式：Energy = Σ Σ P(i, j)²物理意义：也称为“角二阶矩”。它反映了图像纹理的均匀程度和粗糙度。能量值高，表明GLCM中的元素分布非常集中，可能大量集中于主对角线附近，意味着图像灰度分布均匀、纹理粗糙、变化缓慢。一块纯色布料的能量值会接近1。计算解读：由于P(i, j)是概率值，平方之后会更小。只有当某些P(i, j)值本身很大时，其平方和才会大。所以能量是“主导元素”强度的度量。Python计算核心：

energy = 0.0 for i in range(Ng): for j in range(Ng): energy += P[i, j] ** 2

3.3 熵

公式：Entropy = - Σ Σ P(i, j) * log(P(i, j))(约定0*log(0)=0)物理意义：衡量图像纹理的随机性或混乱程度。熵值越大，表示GLCM中元素分布越均匀，图像纹理越复杂、越无序。一片杂乱无章的灌木丛的熵值会很高。它与能量特征通常呈负相关。计算解读：信息论中的概念。概率分布越平均，信息熵越大。如果GLCM中只有一个元素为1，其余为0（极度均匀的图像），熵为0；如果所有元素都相等（极度混乱的图像），熵最大。Python计算核心：

entropy = 0.0 for i in range(Ng): for j in range(Ng): if P[i, j] > 0: # 避免log(0) entropy -= P[i, j] * np.log(P[i, j])

3.4 同质性

公式：Homogeneity = Σ Σ P(i, j) / (1 + |i-j|)物理意义：也称为“逆差矩”。衡量图像纹理的局部均匀性。同质性值高，说明纹理局部变化小，GLCM中的元素主要集中在对角线附近（|i-j|小）。光滑的表面同质性高。计算解读：分母中的(1+|i-j|)使得灰度差越大的像素对，其贡献权重越小。因此，这个特征对主对角线附近的元素非常敏感。Python计算核心：

homogeneity = 0.0 for i in range(Ng): for j in range(Ng): homogeneity += P[i, j] / (1 + abs(i - j))

3.5 相关性

公式：Correlation = Σ Σ [ (i - μ_i)(j - μ_j) * P(i, j) ] / (σ_i * σ_j)其中，μ_i = Σ i * P(i, j),μ_j = Σ j * P(i, j),σ_i² = Σ (i - μ_i)² * P(i, j),σ_j² = Σ (j - μ_j)² * P(i, j)物理意义：衡量GLCM中行元素和列元素之间的线性依赖关系。可以理解为图像在指定方向上，像素灰度值的线性关联程度。如果图像在该方向上有明显的线性结构（如条纹），相关性会较高。计算解读：这是概率论中标准相关系数公式在二维联合概率分布上的应用。计算稍复杂，需要先求出行和列的边际分布及它们的均值和标准差。Python计算核心：

# 计算边际概率和均值 px = np.sum(P, axis=1) # 行和， P(i) py = np.sum(P, axis=0) # 列和， P(j) ux = np.sum(np.arange(Ng) * px) uy = np.sum(np.arange(Ng) * py) # 计算标准差 sigmax = np.sqrt(np.sum((np.arange(Ng) - ux) ** 2 * px)) sigmay = np.sqrt(np.sum((np.arange(Ng) - uy) ** 2 * py)) # 计算相关性 correlation = 0.0 if sigmax * sigmay > 0: # 避免除零 for i in range(Ng): for j in range(Ng): correlation += ((i - ux) * (j - uy) * P[i, j]) / (sigmax * sigmay) else: correlation = 1 # 如果标准差为0，说明灰度恒定，定义为完全相关

这五个特征是最常被使用的，它们从不同侧面刻画了纹理。在实际项目中，我们通常会计算多个方向（如0°, 45°, 90°, 135°）的GLCM，然后对每个特征，取各个方向上的均值或最大值作为最终的特征值，以得到旋转不变的纹理描述。

4. 手把手实现：从图片到特征向量的完整Python代码

理论说得再多，不如亲手跑一遍代码来得实在。这一部分，我们将抛开任何高级图像处理库中可能封装好的GLCM函数，从最底层开始，用NumPy和PIL（或OpenCV）实现整个流程。我会详细解释每一步的意图和注意事项。

4.1 环境准备与图像读取

首先，确保你的环境安装了必要的库。我们使用Pillow读取图片，NumPy进行矩阵运算。

pip install Pillow numpy

然后，我们写一个函数来读取图像并转换为灰度图。这是所有计算的基础。

import numpy as np from PIL import Image def load_and_grayscale(image_path): """ 加载图像并转换为灰度图。 参数: image_path: 图像文件路径。 返回: gray_img: 二维NumPy数组，代表灰度图像，值范围[0, 255]。 """ img = Image.open(image_path).convert('L') # 'L'模式表示8位灰度图 gray_img = np.array(img) return gray_img

4.2 核心：GLCM计算函数

这是整个教程最核心的部分。我们将实现一个函数，根据指定的距离、方向和灰度级，计算GLCM。

def compute_glcm(image, d=1, theta=0, gray_levels=16, symmetric=True, normalized=True): """ 计算灰度共生矩阵。 参数: image: 二维NumPy数组，输入灰度图像。 d: 整数，像素对之间的距离（像素单位）。 theta: 浮点数，方向角度（度）。0（右），45（右上），90（上），135（左上）。 gray_levels: 整数，将灰度压缩到的级别数（通常为8, 16, 32）。 symmetric: 布尔值，是否生成对称GLCM。 normalized: 布尔值，是否将GLCM归一化为概率矩阵。 返回: glcm: 二维NumPy数组，形状为(gray_levels, gray_levels)的GLCM。 """ # 1. 灰度级压缩 # 将[0, 255]的灰度值线性映射到[0, gray_levels-1] max_val = image.max() # 避免除零，同时进行线性缩放和取整 compressed = ((image.astype(np.float32) / max_val) * (gray_levels - 1)).astype(np.int32) # 确保值在有效范围内 compressed = np.clip(compressed, 0, gray_levels - 1) # 2. 根据角度计算偏移量(dx, dy) # 注意图像坐标系：原点在左上角，y轴向下，x轴向右。 # 因此，角度需要对应转换。 theta_rad = np.deg2rad(theta) dx = int(round(d * np.cos(theta_rad))) # x方向偏移 dy = int(round(d * np.sin(theta_rad))) # y方向偏移 # 由于y轴向下，我们通常用 -dy，但这里根据常见定义，我们保持dy为正表示向下。 # 在图像处理中，0度通常指向右，90度指向下（正y方向）。 # 为了符合习惯（0度右，90度下），我们直接使用计算出的dx, dy。 # 3. 初始化GLCM矩阵 glcm = np.zeros((gray_levels, gray_levels), dtype=np.float64) # 4. 获取图像尺寸，并计算有效的像素对区域 rows, cols = compressed.shape # 参考像素的有效范围：要保证邻居像素仍在图像内 start_row = max(0, -dy) if dy < 0 else 0 end_row = rows - dy if dy > 0 else rows start_col = max(0, -dx) if dx < 0 else 0 end_col = cols - dx if dx > 0 else cols # 5. 遍历有效区域，统计共生对 for i in range(start_row, end_row): for j in range(start_col, end_col): ref_val = compressed[i, j] neighbor_val = compressed[i + dy, j + dx] glcm[ref_val, neighbor_val] += 1 # 6. 如果要求对称，则加上转置 if symmetric: glcm = glcm + glcm.T # 7. 如果要求归一化，则转换为概率 if normalized: glcm_sum = glcm.sum() if glcm_sum > 0: glcm = glcm / glcm_sum else: # 如果矩阵全零（理论上不会发生），则返回零矩阵 pass return glcm

关键点解析与避坑指南：

1. 灰度级压缩：(image / max_val) * (gray_levels - 1)这个线性映射是常用方法。但这里有个坑：如果图像对比度很低，max_val可能很小，导致压缩后动态范围依然很窄。更稳健的做法是使用np.percentile或直接设定固定的缩放区间（如2%和98%分位数之间的值映射到[0, gray_levels-1]），或者使用直方图均衡化后再压缩。这里为了简单，采用了线性映射。

2. 坐标偏移计算：角度到(dx, dy)的转换需要小心图像坐标系（y轴向下）。我们使用标准数学坐标系（x向右，y向上）计算，然后应用到图像上。对于theta=90°，(dx=0, dy=1)，表示参考像素和它正下方的像素组成一对，这符合“90度方向”的常见定义。

3. 遍历边界处理：start_row, end_row等变量的计算是为了确保(i+dy, j+dx)这个索引不会超出图像边界。这是实现中容易出错导致索引越界的地方。

4. 对称化处理：glcm = glcm + glcm.T这一步将矩阵与其转置相加，使得P(i, j) = P(j, i)。这相当于在统计时同时考虑了(i, j)和(j, i)这对顺序，使得矩阵对称，得到的纹理特征具有方向无关性（对于某些特征）。这是很多应用中的默认做法。

4.3 特征提取函数

基于第3章讲解的公式，我们实现特征提取函数。这里我们计算对比度、能量、熵、同质性四个最常用的特征。

def glcm_features(glcm): """ 从GLCM中提取Haralick纹理特征。 参数: glcm: 归一化的灰度共生矩阵。 返回: features: 字典，包含计算出的特征值。 """ # 检查glcm是否全零（理论上归一化后不会，但安全起见） if np.sum(glcm) == 0: return {'contrast': 0, 'energy': 0, 'entropy': 0, 'homogeneity': 0} Ng = glcm.shape[0] # 灰度级数 i, j = np.ogrid[:Ng, :Ng] # 创建网格索引，用于向量化计算 # 注意：i和j是二维数组，i代表行索引（参考像素灰度），j代表列索引（邻居像素灰度） # 1. 对比度 contrast = np.sum((i - j) ** 2 * glcm) # 2. 能量 / 角二阶矩 energy = np.sum(glcm ** 2) # 3. 熵 # 避免log(0)，只对非零元素计算 entropy_elements = glcm * np.log(glcm + (glcm == 0)) # glcm==0的位置，log参数为1，结果为0 entropy = -np.sum(entropy_elements) # 4. 同质性 / 逆差矩 homogeneity = np.sum(glcm / (1 + np.abs(i - j))) features = { 'contrast': contrast, 'energy': energy, 'entropy': entropy, 'homogeneity': homogeneity } return features

代码优化技巧：

• 使用np.ogrid创建网格索引i和j，然后利用NumPy的广播机制进行向量化运算，这比用双重for循环快几个数量级。
• 熵的计算中，glcm + (glcm == 0)是一个小技巧。(glcm == 0)得到一个布尔矩阵，在数值运算中True为1，False为0。这样，当glcm为0时，log的参数变为1，log(1)=0，整个项为0，避免了log(0)的警告或错误。

4.4 整合与实战：分析一张纹理图片

现在，我们把所有函数整合起来，用一张真实的图片做测试。我们使用一张布纹图片。

def analyze_texture(image_path, distances=[1], angles=[0, 45, 90, 135], gray_levels=16): """ 分析图像纹理的主函数。 参数: image_path: 图片路径。 distances: 距离列表。 angles: 角度列表（度）。 gray_levels: 灰度级数。 返回: all_features: 列表，包含每个(d, angle)组合下的特征字典。 avg_features: 字典，所有方向和距离下各特征的平均值。 """ # 1. 加载图像 img_gray = load_and_grayscale(image_path) print(f"图像加载成功，尺寸: {img_gray.shape}， 灰度范围: [{img_gray.min()}, {img_gray.max()}]") all_features = [] feature_sums = {'contrast': 0.0, 'energy': 0.0, 'entropy': 0.0, 'homogeneity': 0.0} count = 0 # 2. 遍历所有距离和角度组合 for d in distances: for angle in angles: # 计算GLCM glcm = compute_glcm(img_gray, d=d, theta=angle, gray_levels=gray_levels) # 提取特征 features = glcm_features(glcm) features['distance'] = d features['angle'] = angle all_features.append(features) # 累加特征值用于求平均 for key in feature_sums.keys(): feature_sums[key] += features[key] count += 1 # 3. 计算平均特征 avg_features = {key: feature_sums[key] / count for key in feature_sums} return all_features, avg_features # 使用示例 if __name__ == "__main__": # 替换为你的图片路径 image_path = "fabric_texture.jpg" # 示例：一张布纹图片 try: detailed_features, avg_result = analyze_texture(image_path) print("\n--- 详细特征（每个方向） ---") for feat in detailed_features: print(f"距离 {feat['distance']}, 角度 {feat['angle']}°: " f"对比度={feat['contrast']:.4f}, " f"能量={feat['energy']:.4f}, " f"熵={feat['entropy']:.4f}, " f"同质性={feat['homogeneity']:.4f}") print("\n--- 平均纹理特征 ---") for key, value in avg_result.items(): print(f"{key}: {value:.4f}") except FileNotFoundError: print(f"错误：找不到图片文件 '{image_path}'，请检查路径。") except Exception as e: print(f"处理过程中发生错误: {e}")

运行这段代码，你会得到四个方向上的纹理特征值，以及它们的平均值。对于纹理明显的图片（如木材、织物、云层），你会看到不同方向的特征值可能有差异（各向异性），而平均特征则给出了一个整体的纹理描述。

5. 参数调优与实战经验：让GLCM真正为你所用

代码跑通了，只是第一步。要让GLCM特征在实际项目（如分类、分割、检索）中发挥最大效用，参数的选择和预处理至关重要。这部分是我在实际项目中踩过坑后总结的经验。

5.1 关键参数深度解析

1. 灰度级数：这是最重要的参数之一。

   • 值太小（如4）：会丢失大量灰度细节，导致纹理信息过度简化，不同纹理可能被映射到相同的GLCM，区分度下降。
   • 值太大（如64或256）：GLCM矩阵会变得非常大且稀疏，计算量激增，并且由于统计样本相对不足（图像像素数有限），矩阵中会出现很多零值，使得特征估计不稳定，噪声敏感。
   • 经验值：对于大多数自然图像，16级是一个很好的起点，它在计算复杂度和信息保留之间取得了平衡。对于非常平滑或纹理简单的图像，可以尝试8级；对于高分辨率、纹理极其丰富的图像（如显微图像），可以考虑32级。务必在你的数据集上进行交叉验证，选择表现最好的级别。

2. 距离：

   • d=1：考察的是最相邻像素之间的关系，描述的是微观纹理、细粒度结构。
   • d增大：考察的是更大尺度上的像素关系，描述的是宏观纹理、粗粒度结构。例如，对于周期性纹理（如格子布），当d等于纹理周期时，GLCM会表现出特定的模式。
   • 如何选择：没有固定答案。通常的做法是尝试一组距离，如[1, 2, 3, 4]，然后计算每个距离下的特征，将所有特征拼接成一个长特征向量。这叫做多尺度纹理分析，能同时捕捉不同尺度的纹理信息。

3. 方向：

   • 默认计算[0°, 45°, 90°, 135°]四个方向，然后对每个特征取平均值。这是为了获得旋转不变性，即无论图像如何旋转，纹理特征值大致不变。
   • 在某些应用中，纹理的方向性本身就是重要信息（如检测木材纹理走向）。这时，可以保留各个方向的特征，或者取最大值/最小值来表征最强的方向性响应。

4. 图像区域：

   • GLCM统计的是全局信息。如果图像包含多种纹理（如一幅风景画中有天空、树木、地面），全局GLCM会混合所有信息，失去区分力。
   • 标准做法是使用滑动窗口：将图像划分为许多重叠或非重叠的小块（如64x64像素），对每个小块计算GLCM特征。这样，每个小块得到一个特征向量，整张图就变成了一个特征图或一组特征向量，可用于像素级或区域级的纹理分析。

5.2 预处理：比想象中更重要

1. 灰度标准化/均衡化：直接使用原始灰度值计算GLCM，其值受光照影响极大。同一纹理，在亮光和暗光下，GLCM会完全不同。因此，必须先进行光照归一化。常用方法有：

   • 直方图均衡化：增强对比度，使灰度分布更均匀，对光照变化有一定鲁棒性。
   • 自适应直方图均衡化：效果更好，能增强局部对比度。
   • 更简单的方法：对每个图像块，进行零均值单位方差标准化，即(patch - mean(patch)) / std(patch)。这能有效消除整体亮度和对比度的影响。

2. 滤波去噪：图像噪声（尤其是椒盐噪声、高斯噪声）会严重污染GLCM，因为噪声点会创造出原本不存在的灰度值组合。在计算GLCM前，进行适度的平滑滤波（如高斯滤波、中值滤波）是很有必要的。但滤波强度不宜过大，以免模糊掉真正的纹理边缘。

5.3 一个完整的、鲁棒的纹理分析流程示例

结合以上经验，一个更健壮的纹理特征提取流程如下：

def robust_texture_feature_extraction(image_path, window_size=64, step=32, gray_levels=16, distances=[1, 3]): """ 鲁棒的纹理特征提取流程：滑动窗口 + 预处理 + 多尺度。 参数: image_path: 图片路径。 window_size: 滑动窗口大小。 step: 滑动步长。 gray_levels: 灰度级。 distances: 多尺度距离列表。 返回: feature_matrix: 二维数组，每行是一个窗口的特征向量。 window_locations: 每个窗口的左上角坐标列表。 """ from scipy.ndimage import gaussian_filter from skimage.exposure import equalize_hist # 可选，需要scikit-image img_gray = load_and_grayscale(image_path) rows, cols = img_gray.shape # 1. 全局预处理：轻度高斯滤波去噪 img_filtered = gaussian_filter(img_gray.astype(np.float32), sigma=0.5) features_list = [] locations = [] angles = [0, 45, 90, 135] # 2. 滑动窗口遍历 for y in range(0, rows - window_size + 1, step): for x in range(0, cols - window_size + 1, step): window = img_filtered[y:y+window_size, x:x+window_size] # 3. 窗口内预处理：局部对比度归一化 (Local Contrast Normalization) window_mean = np.mean(window) window_std = np.std(window) if window_std > 1e-6: # 避免除零 window_normalized = (window - window_mean) / window_std else: window_normalized = window * 0.0 # 4. 将归一化后的窗口值线性映射到[0, gray_levels-1] # 由于经过了零均值单位方差，值域大致在[-3,3]，我们将其缩放到[0, gray_levels-1] win_min, win_max = window_normalized.min(), window_normalized.max() if (win_max - win_min) > 1e-6: window_scaled = ((window_normalized - win_min) / (win_max - win_min) * (gray_levels - 1)).astype(np.int32) else: window_scaled = np.zeros_like(window_normalized, dtype=np.int32) window_scaled = np.clip(window_scaled, 0, gray_levels - 1) window_feature_vector = [] # 5. 多尺度、多方向GLCM特征计算 for d in distances: for angle in angles: glcm = compute_glcm(window_scaled, d=d, theta=angle, gray_levels=gray_levels) feats = glcm_features(glcm) # 将四个特征值按顺序加入向量 window_feature_vector.extend([feats['contrast'], feats['energy'], feats['entropy'], feats['homogeneity']]) # 特征向量长度 = len(distances) * len(angles) * 4 features_list.append(window_feature_vector) locations.append((y, x)) return np.array(features_list), locations

这个流程提取出的feature_matrix，每一行代表图像中一个局部区域的纹理特征。这个矩阵可以直接用于训练分类器（如SVM、随机森林）进行纹理分类，或者通过聚类来分析图像中不同的纹理区域。

5.4 常见问题与排错

• 特征值全是0或NaN：检查你的图像数据。可能是图像加载失败（全黑或全白），或者灰度压缩后所有像素值相同，导致GLCM除主对角线外全为0。添加数据有效性检查。
• 计算速度慢：对于大图像或小窗口、多尺度多方向的情况，循环计算GLCM会非常慢。解决方案：
  1. 使用scikit-image库的greycomatrix和greycoprops函数，它们是高度优化的C扩展，速度极快。
  2. 如果必须自己实现，确保使用NumPy向量化操作，避免Python层级的循环。
  3. 考虑减少灰度级数、距离和方向的数量。
• 特征区分度不高：可能是参数不适合你的数据。尝试调整灰度级、距离，并务必进行局部归一化。也可以考虑使用其他纹理特征描述子（如LBP、Gabor滤波器）作为补充或替代。

GLCM是一个强大的工具，但它不是万能的。它对于具有明显周期性或结构性纹理的图像效果很好，但对于完全随机或语义性很强的“纹理”（比如人脸），效果可能有限。理解其原理和局限，合理地选择参数和预处理方法，才能让它在你具体的项目里真正发光发热。