 实战：Python 提取 14 种纹理特征与遥感图像分类)
灰度共生矩阵(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, distances=distances, angles=angles, levels=16, symmetric=True, normed=True) # 计算各统计特征 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, desc=class_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], levels=16, symmetric=True, normed=True) # 计算各方向特征均值 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, index=False)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'], axis=1) y = df['class_idx'] # 划分训练测试集 X_train, X_test, y_train, y_test = train_test_split( X, y, test_size=0.2, random_state=42) # 训练随机森林 clf = RandomForestClassifier(n_estimators=100, random_state=42) clf.fit(X_train, y_train) # 评估模型 y_pred = clf.predict(X_test) print(classification_report(y_test, y_pred, target_names=df['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, distances=distances, angles=angles, levels=levels, symmetric=True, normed=True) # 计算特征区分度指标(示例) 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], levels=16, symmetric=True, normed=True) # 计算各方向特征均值 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, distance=1, angle=0, levels=16): """ GPU加速的GLCM计算 :param image: 输入图像 :param distance: 像素距离 :param angle: 角度(弧度) :param levels: 灰度级数 :return: GLCM矩阵 """ # 将图像传输到GPU image_gpu = cp.asarray(image) # 初始化GLCM矩阵 glcm = cp.zeros((levels, levels), dtype=cp.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_components=10, method='pca'): """ 特征降维 :param X: 特征矩阵 :param y: 标签 :param n_components: 目标维度 :param method: 降维方法(pca/lda) :return: 降维后的特征 """ if method == 'pca': reducer = PCA(n_components=n_components) elif method == 'lda': reducer = LDA(n_components=min(n_components, len(np.unique(y))-1)) else: raise ValueError("Unsupported method") return reducer.fit_transform(X, y)在实际项目中，GLCM特征常与其他特征（如颜色特征、形状特征）组合使用，以获得更好的分类性能。例如，在土地覆盖分类任务中，结合NDVI等光谱指数可以显著提高分类精度。
