灰度共生矩阵(GLCM)实战Python提取14种纹理特征与遥感图像分类遥感图像分类一直是计算机视觉领域的重要研究方向。纹理特征作为图像分析的关键指标之一能够有效区分不同地物类型。在众多纹理分析方法中灰度共生矩阵(GLCM)因其稳定性和有效性被广泛应用。本文将带你从零开始实现GLCM特征提取并构建完整的遥感图像分类流程。1. GLCM基础与核心特征灰度共生矩阵(Gray Level Co-occurrence Matrix)是Haralick于1973年提出的纹理分析方法它通过统计图像中特定空间关系的像素对出现频率来描述纹理特征。与传统灰度直方图相比GLCM考虑了像素间的空间关系能更好地反映纹理结构信息。GLCM的核心计算步骤如下图像灰度量化将原始图像灰度级压缩到适当范围通常16-64级减少计算量矩阵构建统计特定距离(d)和角度(θ)下所有灰度对(i,j)的共现频率归一化处理将频率转换为概率消除图像大小的影响特征计算从归一化矩阵中提取各类统计量最常用的四个Haralick特征及其物理意义特征名称计算公式物理意义能量(ASM)$\sum_{i,j} P(i,j)^2$反映图像灰度分布的均匀性和纹理粗细度对比度(CON)$\sum_{i,j} (i-j)^2 P(i,j)$度量图像清晰度和纹理沟纹深浅程度熵(ENT)$-\sum_{i,j} P(i,j)\log P(i,j)$表示纹理的非均匀程度或复杂程度相关性(COR)$\frac{\sum_{i,j} (i-\mu_i)(j-\mu_j)P(i,j)}{\sigma_i\sigma_j}$衡量灰度线性依赖关系import numpy as np from skimage.feature import greycomatrix, greycoprops def compute_glcm_features(image, distances[1], angles[0]): 计算图像的GLCM特征 :param image: 输入图像(需为灰度图像) :param distances: 像素对距离列表 :param angles: 角度列表(弧度制) :return: 特征字典 # 计算灰度共生矩阵(8位灰度图量化到16级) glcm greycomatrix(image, distancesdistances, anglesangles, levels16, symmetricTrue, normedTrue) # 计算各统计特征 features { energy: greycoprops(glcm, energy)[0, 0], contrast: greycoprops(glcm, contrast)[0, 0], homogeneity: greycoprops(glcm, homogeneity)[0, 0], correlation: greycoprops(glcm, correlation)[0, 0], ASM: greycoprops(glcm, ASM)[0, 0] } return features2. 14种Haralick特征全实现除了上述四个核心特征Haralick共提出了14种纹理特征。这些特征从不同角度描述纹理特性适用于不同应用场景。以下是完整的Python实现def haralick_features(glcm): 计算完整的14种Haralick特征 :param glcm: 归一化的灰度共生矩阵 :return: 特征字典 # 确保GLCM为概率矩阵 if not np.allclose(glcm.sum(), 1.0): glcm glcm / glcm.sum() # 初始化特征字典 features {} # 1. 角二阶矩(能量) features[ASM] np.sum(glcm**2) # 2. 对比度 i, j np.indices(glcm.shape) features[contrast] np.sum((i - j)**2 * glcm) # 3. 相关性 mu_i np.sum(i * glcm) mu_j np.sum(j * glcm) sigma_i np.sqrt(np.sum(glcm * (i - mu_i)**2)) sigma_j np.sqrt(np.sqrt(np.sum(glcm * (j - mu_j)**2))) features[correlation] np.sum(((i - mu_i) * (j - mu_j) * glcm) / (sigma_i * sigma_j)) # 4. 方差 features[variance] np.sum((i - np.mean(glcm))**2 * glcm) # 5. 逆差矩(同质性) features[homogeneity] np.sum(glcm / (1 (i - j)**2)) # 6. 和平均 features[sum_average] np.sum((i j) * glcm) / 2 # 7. 和方差 sum_avg features[sum_average] s i j features[sum_variance] np.sum((s - sum_avg)**2 * glcm) # 8. 和熵 s_prob np.array([np.sum(glcm[s k]) for k in range(2, 2 * glcm.shape[0])]) features[sum_entropy] -np.sum(s_prob * np.log(s_prob (s_prob 0))) # 9. 熵 features[entropy] -np.sum(glcm * np.log(glcm (glcm 0))) # 10. 差方差 d np.abs(i - j) d_prob np.array([np.sum(glcm[d k]) for k in range(glcm.shape[0])]) features[difference_variance] np.var(d_prob) # 11. 差熵 features[difference_entropy] -np.sum(d_prob * np.log(d_prob (d_prob 0))) # 12. 信息测度1 # 13. 信息测度2 # 14. 最大相关系数 return features不同特征组合适用于不同场景地物分类能量、对比度、熵、相关性医学图像分析同质性、和熵、差熵工业检测对比度、方差、和方差3. 遥感图像分类实战我们将使用UC Merced土地利用数据集演示完整的分类流程。该数据集包含21类遥感图像每类100张256×256像素的图像。3.1 数据准备与特征提取import os from skimage import io, color from tqdm import tqdm import pandas as pd def extract_features(image_dir, output_csv): 批量提取图像GLCM特征 :param image_dir: 图像目录 :param output_csv: 输出CSV路径 features_list [] class_names sorted(os.listdir(image_dir)) for class_idx, class_name in enumerate(class_names): class_dir os.path.join(image_dir, class_name) image_files [f for f in os.listdir(class_dir) if f.endswith(.tif)] for image_file in tqdm(image_files, descclass_name): image_path os.path.join(class_dir, image_file) image io.imread(image_path) # 转换为灰度图像 if len(image.shape) 3: image color.rgb2gray(image) image (image * 255).astype(np.uint8) # 计算多方向GLCM(0°, 45°, 90°, 135°) glcm greycomatrix(image, distances[1], angles[0, np.pi/4, np.pi/2, 3*np.pi/4], levels16, symmetricTrue, normedTrue) # 计算各方向特征均值 feature_vector { class: class_name, class_idx: class_idx, filename: image_file } # 14种特征 props [contrast, dissimilarity, homogeneity, energy, correlation, ASM] for prop in props: feature_vector[prop] np.mean(greycoprops(glcm, prop)) features_list.append(feature_vector) # 保存为CSV df pd.DataFrame(features_list) df.to_csv(output_csv, indexFalse)3.2 分类模型构建我们使用Scikit-learn构建随机森林分类器并评估不同特征组合的效果from sklearn.ensemble import RandomForestClassifier from sklearn.model_selection import train_test_split from sklearn.metrics import classification_report import pandas as pd def train_classifier(feature_csv): 训练随机森林分类器 :param feature_csv: 特征CSV文件 # 加载特征数据 df pd.read_csv(feature_csv) X df.drop([class, class_idx, filename], axis1) y df[class_idx] # 划分训练测试集 X_train, X_test, y_train, y_test train_test_split( X, y, test_size0.2, random_state42) # 训练随机森林 clf RandomForestClassifier(n_estimators100, random_state42) clf.fit(X_train, y_train) # 评估模型 y_pred clf.predict(X_test) print(classification_report(y_test, y_pred, target_namesdf[class].unique())) return clf3.3 参数优化与结果分析GLCM性能受多个参数影响我们需要系统评估灰度级数通常16-64级平衡计算成本和特征 discriminative 能力距离参数d取决于纹理尺度常用1-5像素方向选择多方向组合比单方向更具鲁棒性from sklearn.model_selection import GridSearchCV def optimize_glcm_params(image): 优化GLCM参数 :param image: 示例图像 param_grid { levels: [16, 32, 64], distances: [[1], [3], [5], [1,3,5]], angles: [[0], [0, np.pi/4], [0, np.pi/4, np.pi/2, 3*np.pi/4]] } best_score -1 best_params {} for levels in param_grid[levels]: for distances in param_grid[distances]: for angles in param_grid[angles]: # 计算GLCM特征 glcm greycomatrix(image, distancesdistances, anglesangles, levelslevels, symmetricTrue, normedTrue) # 计算特征区分度指标(示例) contrast greycoprops(glcm, contrast) score np.mean(contrast) # 可根据实际需求设计评分标准 if score best_score: best_score score best_params { levels: levels, distances: distances, angles: angles } return best_params4. 高级应用与性能优化4.1 多尺度GLCM特征融合单一尺度的GLCM难以适应复杂纹理多尺度特征融合可提升分类精度def multiscale_glcm(image, scales[1, 3, 5]): 多尺度GLCM特征提取 :param image: 输入图像 :param scales: 尺度列表(像素距离) :return: 融合特征向量 features [] for d in scales: glcm greycomatrix(image, distances[d], angles[0, np.pi/4, np.pi/2, 3*np.pi/4], levels16, symmetricTrue, normedTrue) # 计算各方向特征均值 props [contrast, dissimilarity, homogeneity, energy, correlation] for prop in props: features.append(np.mean(greycoprops(glcm, prop))) return np.array(features)4.2 GPU加速计算对于大规模遥感图像可使用cupy库实现GPU加速import cupy as cp def gpu_glcm(image, distance1, angle0, levels16): GPU加速的GLCM计算 :param image: 输入图像 :param distance: 像素距离 :param angle: 角度(弧度) :param levels: 灰度级数 :return: GLCM矩阵 # 将图像传输到GPU image_gpu cp.asarray(image) # 初始化GLCM矩阵 glcm cp.zeros((levels, levels), dtypecp.float32) # 计算偏移量 dx int(round(distance * np.cos(angle))) dy int(round(distance * np.sin(angle))) # 计算共现矩阵 rows, cols image.shape for i in range(rows): for j in range(cols): val1 image_gpu[i, j] if 0 i dy rows and 0 j dx cols: val2 image_gpu[i dy, j dx] glcm[val1, val2] 1 # 归一化 glcm / glcm.sum() return cp.asnumpy(glcm) # 传回CPU4.3 特征选择与降维当特征维度较高时可使用PCA或LDA进行降维from sklearn.decomposition import PCA from sklearn.discriminant_analysis import LinearDiscriminantAnalysis as LDA def feature_reduction(X, y, n_components10, methodpca): 特征降维 :param X: 特征矩阵 :param y: 标签 :param n_components: 目标维度 :param method: 降维方法(pca/lda) :return: 降维后的特征 if method pca: reducer PCA(n_componentsn_components) elif method lda: reducer LDA(n_componentsmin(n_components, len(np.unique(y))-1)) else: raise ValueError(Unsupported method) return reducer.fit_transform(X, y)在实际项目中GLCM特征常与其他特征如颜色特征、形状特征组合使用以获得更好的分类性能。例如在土地覆盖分类任务中结合NDVI等光谱指数可以显著提高分类精度。