Python+OpenCV实战：用灰度共生矩阵(GLCM)分析医学影像纹理特征

Python+OpenCV实战：用灰度共生矩阵(GLCM)深度解析医学影像纹理

在医学影像分析领域，我们常常需要超越人眼，去量化那些看似“感觉”到的组织特性——比如肝脏的粗糙度、肺结节的异质性，或是肿瘤边界的模糊程度。这些特性，在图像处理中统称为“纹理”。对于CT、MRI等影像，纹理特征不仅是像素灰度的简单统计，更是像素间空间关系的复杂表达，它蕴含着组织微观结构、生理状态乃至病理变化的宝贵信息。作为一名医学图像处理工程师或AI医疗开发者，掌握一种能够精准捕捉并量化这些空间关系的方法，是从“看图像”走向“读数据”的关键一步。

灰度共生矩阵（GLCM）正是这样一把钥匙。它不像深度学习那样是个“黑箱”，而是基于严谨的统计原理，将图像的纹理模式转化为可计算的矩阵，进而衍生出一系列具有明确物理意义的特征值。本文不会重复那些教科书上的理论推导，我们将直接切入实战。我将结合Python和OpenCV，手把手带你构建一套从原始DICOM图像到GLCM特征向量的完整分析流程。重点会放在那些真正困扰实践者的核心问题上：面对一张具体的肝脏CT图像，我该如何选择那个神秘的距离d和角度θ？计算出来的角二阶矩（ASM） 是0.15，这个数字到底意味着组织是均匀的还是粗糙的？熵（ENT） 的大小与肿瘤的恶性程度有怎样的潜在关联？我们将通过具体的代码、真实的医学影像案例（已脱敏）和深入的解读，让GLCM从数学公式变成你诊断辅助或研究分析中的得力工具。

1. 理解GLCM：超越像素的统计视角

在开始写代码之前，我们需要在概念层面重新认识一下GLCM。很多人把它理解为一个复杂的数学变换，但实际上，它的核心思想非常直观：描述一对像素“结伴出现”的规律。

想象一下，你正在观察一片森林的航拍图。如果这片森林是人工种植的整齐杉树林，那么你在一个像素点看到“深绿色”（代表树冠）时，它右边紧邻的像素点极有可能也是“深绿色”。这种“深绿-深绿”的组合会频繁出现。反之，如果是一片原始杂木林，树冠大小、间距不一，“深绿-浅绿”（树冠与间隙）或“深绿-中绿”的组合会更加多样。GLCM做的就是这件事——它统计整幅图像中，所有满足特定空间关系（比如，右边1个像素）的像素对，其灰度值组合（i, j）出现的频率。

这个“特定的空间关系”，就是由参数d（距离）和θ（角度）定义的。d=1, θ=0° 表示统计每个像素与其正右方相邻像素的灰度对。通过改变d和θ，我们可以捕捉不同尺度、不同方向上的纹理模式。例如，分析肌肉纤维的走向可能需要关注特定角度；而分析脂肪浸润的粗糙度，可能需要一个较小的d来捕捉细微变化。

生成的GLCM本身是一个方阵，其行和列索引分别代表像素对中第一个像素（参考像素）和第二个像素（邻居像素）的灰度级。如果我们将图像灰度级量化为8级（0-7），那么GLCM就是一个8x8的矩阵。矩阵中元素P(i, j)的值，就代表了灰度组合(i, j)出现的概率（归一化后）。

注意：在医学影像中，直接使用原始HU值（CT）或信号强度（MRI）作为灰度级通常范围太大（如-1000到+3000），直接计算GLCM会导致矩阵稀疏且计算量大。因此，灰度级量化是必不可少的前置步骤，通常将灰度范围重新映射到8、16、32或64级。这个量化级数也是一个重要参数，会影响特征的敏感性。

为了更直观地理解不同纹理对应的GLCM差异，我们可以看一个简单的对比：

这个表格告诉我们，GLCM特征不是孤立的数字，它们与视觉感知的纹理属性有着直接的对应关系。接下来，我们就用代码来生成并分析这些矩阵。

2. 实战环境搭建与数据预处理

工欲善其事，必先利其器。我们的实战将基于Python生态中强大的科学计算和图像处理库。首先确保你的环境已安装以下核心包：

pip install opencv-python-headless numpy scikit-image pydicom matplotlib

• OpenCV (cv2): 用于基础的图像读写、灰度转换和矩阵操作。
• NumPy: 所有数值计算的基石，GLCM本质上就是NumPy数组的运算。
• scikit-image (skimage): 它提供了skimage.feature.greycomatrix和greycoprops函数，这是计算GLCM及其特征的官方且高效实现，强烈推荐使用，避免重复造轮子。
• PyDicom: 用于读取医学影像标准格式DICOM文件。
• Matplotlib: 用于可视化结果。

假设我们有一张腹部的CT扫描DICOM文件ct_abdomen.dcm。第一步是将其加载并转换为适合纹理分析的灰度图像数组。

import pydicom
import numpy as np
import cv2
from skimage.feature import greycomatrix, greycoprops
import matplotlib.pyplot as plt

# 1. 读取DICOM文件
ds = pydicom.dcmread('ct_abdomen.dcm')
# 获取像素数组，CT值通常以HU为单位
ct_array = ds.pixel_array.astype(np.float32)

# 2. 应用可能的窗宽窗位（例如肝脏窗）
# 假设窗位为40，窗宽为400
window_center = 40
window_width = 400
img_min = window_center - window_width // 2
img_max = window_center + window_width // 2
# 将HU值裁剪并缩放到0-255范围
ct_windowed = np.clip(ct_array, img_min, img_max)
ct_normalized = ((ct_windowed - img_min) / (img_max - img_min) * 255).astype(np.uint8)

# 3. 感兴趣区域(ROI)提取
# 假设我们通过手动或分割算法得到了一个肝脏区域的二值掩膜 `liver_mask`
# 这里为了演示，我们简单地从图像中心截取一个128x128的区域作为ROI
height, width = ct_normalized.shape
roi = ct_normalized[height//2-64:height//2+64, width//2-64:width//2+64]

# 可视化原始ROI
plt.figure(figsize=(10, 4))
plt.subplot(1, 2, 1)
plt.imshow(roi, cmap='gray')
plt.title('肝脏CT ROI (灰度图)')
plt.axis('off')
plt.subplot(1, 2, 2)
plt.hist(roi.ravel(), bins=50, range=[0, 256])
plt.title('ROI灰度直方图')
plt.xlabel('灰度值')
plt.ylabel('频数')
plt.tight_layout()
plt.show()

预处理的关键步骤包括窗宽窗位调整和ROI提取。医学影像的原始值范围很广，直接分析会引入大量无关噪声（如空气、骨骼的极端值）。根据目标组织（如肝脏、肺结节）调整窗宽窗位，能显著增强相关纹理的对比度。ROI提取则确保了我们的纹理分析聚焦于特定的解剖结构或病灶，这是获得有临床意义结果的前提。

3. 核心参数选择：距离、角度与灰度级

这是GLCM应用中最具技巧性的一环。参数的选择没有绝对的“黄金标准”，必须结合具体的图像内容、纹理尺度和你关心的生物学问题。

3.1 距离 d 的选择

d 代表了分析纹理的“尺度”。较小的d（如1, 2）对像素间的细微变化敏感，适合捕捉细密纹理（如早期纤维化）。较大的d（如8, 16）则反映更大范围内的灰度关系，适合分析粗糙或周期性较强的纹理（如某些肿瘤内部的坏死区域）。

提示：一个实用的策略是进行多尺度分析。即计算一组不同d值下的特征，观察特征值随d变化的趋势。这种多尺度特征本身可能就是强有力的鉴别因子。例如，恶性肿瘤的熵值可能随d增大而持续升高，而良性病变的熵值变化可能较平缓。

3.2 角度 θ 的选择

θ 定义了像素对的方向。通常我们会计算四个方向（0°, 45°, 90°, 135°）的GLCM，然后对每个特征取四个方向上的平均值或最大值。这样做有两个好处：1) 获得旋转不变的特征，减少因患者体位或扫描方向带来的影响；2) 如果纹理本身具有方向性（如肌肉），分析各个方向特征的差异（即各向异性）本身就是一个重要信息。

3.3 灰度级 levels 的选择

这是greycomatrix函数中的levels参数。它将输入图像的灰度值重新量化到[0, levels-1]的范围内。级数太少（如4）会丢失大量纹理细节，导致特征区分度下降；级数太多（如256）会使GLCM矩阵非常稀疏，计算量大且容易受噪声干扰。在医学影像分析中，16级或32级是一个常见的折中选择。

让我们用代码演示如何计算多参数GLCM并提取特征：

def compute_glcm_features(image, distances, angles, levels=16):
    """
    计算图像在多距离、多角度下的GLCM特征。
    返回一个特征字典。
    """
    # 量化灰度级
    # 将图像灰度范围线性映射到 [0, levels-1]
    max_val = image.max()
    min_val = image.min()
    if max_val == min_val:
        quantized = np.zeros_like(image, dtype=np.uint8)
    else:
        quantized = ((image - min_val) / (max_val - min_val) * (levels - 1)).astype(np.uint8)

    # 计算GLCM
    # skimage的greycomatrix要求图像为整数类型，且灰度值在[0, levels-1]
    glcm = greycomatrix(quantized,
                        distances=distances,
                        angles=angles,
                        levels=levels,
                        symmetric=True, # 通常使用对称GLCM，即同时考虑(i,j)和(j,i)
                        normed=True)    # 归一化为概率

    # 定义要计算的特征
    features = {}
    feature_names = ['contrast', 'dissimilarity', 'homogeneity', 'energy', 'correlation', 'ASM']
    # 注意：skimage中的'energy'即ASM的平方根，'ASM'才是角二阶矩。

    for prop in feature_names:
        # greycoprops返回的形状为 (len(distances), len(angles))
        feat_val = greycoprops(glcm, prop)
        # 我们取所有距离和角度上的平均值作为该特征的最终标量值
        features[prop] = np.mean(feat_val)

    # 单独计算熵 (skimage未直接提供，需从GLCM计算)
    # 熵 ENT = -sum_i sum_j P(i,j) * log(P(i,j) + eps)
    eps = 1e-10
    entropy = -np.sum(glcm[:, :, 0, 0] * np.log(glcm[:, :, 0, 0] + eps)) # 这里先取第一个距离和角度为例
    # 更严谨的做法是遍历所有距离和角度的GLCM片计算平均熵
    entropy_list = []
    for d in range(len(distances)):
        for a in range(len(angles)):
            p = glcm[:, :, d, a]
            entropy_list.append(-np.sum(p * np.log(p + eps)))
    features['entropy'] = np.mean(entropy_list)

    return features, glcm

# 定义参数
distances = [1, 3, 5]  # 多尺度分析
angles = [0, np.pi/4, np.pi/2, 3*np.pi/4]  # 0, 45, 90, 135度
levels = 16

# 计算ROI的特征
features, glcm_matrix = compute_glcm_features(roi, distances, angles, levels)

print("提取的GLCM特征值：")
for key, value in features.items():
    print(f"{key}: {value:.4f}")

这段代码的核心是compute_glcm_features函数。它完成了灰度量化、GLCM计算和特征提取的全过程。我们选择了三个距离和四个角度，并对结果取平均，从而得到一组相对稳健的旋转和尺度平均的特征。

4. 特征解读与临床意义关联

拿到一堆特征值后，如何解读？这才是连接图像处理与临床应用的桥梁。下面我们逐一拆解：

• 角二阶矩 (ASM) / 能量 (Energy)：它衡量图像灰度分布的均匀性。ASM值高（接近1），说明GLCM中的元素集中分布，图像纹理均匀、规则。在医学影像中，均匀的液体（如囊肿、单纯性肾囊肿）、均匀的软组织（如正常的脾脏）通常表现出较高的ASM。反之，ASM值低，表明纹理杂乱、异质，常见于恶性肿瘤、复杂囊肿或炎症区域。

• 熵 (Entropy, ENT)：它衡量图像纹理的随机性或复杂程度。熵值越大，说明灰度分布越随机，纹理越复杂。高度异质性的肿瘤，由于其内部包含坏死、出血、钙化等多种成分，灰度变化无常，因此熵值通常较高。而均匀的组织熵值较低。有研究表明，在肺结节、乳腺肿块的分级中，熵是一个重要的鉴别特征。

• 对比度 (Contrast)：它反映图像的局部变化强度，即纹理的清晰度。对比度值高，意味着像素对之间的灰度差异大，图像看起来“棱角分明”。例如，骨骼与软组织的交界处、钙化点边缘的对比度会很高。在肿瘤分析中，边界清晰的肿瘤可能比浸润性、边界模糊的肿瘤有更高的对比度。

• 逆差分矩 (Homogeneity, IDM)：它衡量图像纹理的局部均匀性。与ASM类似，但IDM对主对角线附近的元素赋予更高的权重。因此，它对靠近对角线的灰度变化更敏感。IDM值高，表示纹理局部均匀。在肝纤维化评估中，随着纤维化程度加重，肝脏纹理从均匀变得粗糙，IDM值可能会下降。

• 相关性 (Correlation)：它度量图像中像素灰度的线性依赖关系。相关性高，说明在某个方向上，像素灰度值以线性方式相关（例如，沿某个方向灰度缓慢递增或递减）。这可能对应有方向性的纹理，如肌肉纤维或血管束。

为了更系统地理解这些特征在鉴别诊断中的潜在作用，可以参考以下对比分析：

注意：绝对不要孤立地看待任何一个特征。单个特征的诊断价值有限，且容易受成像参数（如CT层厚、MRI序列）影响。在实际应用中，通常将多个GLCM特征（常结合其他形状、强度特征）组成一个高维特征向量，输入到机器学习分类器（如SVM、随机森林）中，进行综合判断。此外，建立自己科室或研究项目内部的、基于特定设备和协议的“正常参考范围”至关重要。

5. 高级技巧与性能优化

掌握了基础流程后，我们可以探讨一些提升分析质量和效率的高级技巧。

技巧一：动态灰度级与ROI自适应 对于不同对比度的图像或不同组织的ROI，使用固定的灰度级范围（如0-255）可能不合适。可以采用基于ROI灰度直方图的分位数进行自适应映射。例如，将ROI灰度值的2%和98%分位数分别映射到0和levels-1，这样可以增强ROI内部纹理的对比度，抑制极端离群值的影响。

def adaptive_quantize(roi, levels=16):
    """基于ROI直方图分位数的自适应灰度量化"""
    p_low, p_high = np.percentile(roi, (2, 98))
    roi_clipped = np.clip(roi, p_low, p_high)
    quantized = ((roi_clipped - p_low) / (p_high - p_low + 1e-10) * (levels - 1)).astype(np.uint8)
    return quantized

技巧二：多尺度与多方向特征的融合 如前所述，计算多个d和θ下的特征后，除了取平均，还可以构造更丰富的特征描述符。例如，可以计算每个特征在不同d下的变化率（纹理尺度特性），或计算不同θ下特征的标准差（纹理各向异性）。这些衍生特征有时比原始平均值更具鉴别力。

技巧三：与深度学习特征融合 GLCM特征属于手工设计的“浅层”特征。在现代AI医疗研究中，常将其与深度学习模型（如CNN）提取的“深层”特征进行融合。思路是：用预训练的CNN提取高层语义特征，同时用GLCM提取中层的纹理统计特征，将两者拼接后输入分类器。这种方法结合了深度学习的表示能力和传统纹理特征的可解释性。

性能优化建议：

• ROI尺寸：ROI不宜过小，否则统计量不可靠。通常建议至少包含32x32个像素。对于特别异质的区域，可能需要更大的ROI。
• 计算加速：对于需要处理大量图像的研究，skimage的GLCM计算可以并行化。或者，可以考虑使用更高效的C++实现进行封装，并通过Python调用。
• 特征选择：当计算了大量GLCM特征（多参数、多尺度）后，可能会面临特征维度过高的问题。使用递归特征消除（RFE）、基于树模型的特征重要性排序或LASSO等方法进行特征选择，剔除冗余特征，能提升后续模型的性能和泛化能力。

在我处理一批肝脏CT图像用于脂肪肝分级时，最初使用了固定的d=1和levels=32，结果发现轻度与中度脂肪肝的特征区分不明显。后来改为使用distances=[1,3,5]并计算了熵在不同距离下的斜率（(熵_d5 - 熵_d1)/4），发现这个“熵-距离斜率”在两组间有显著差异，成为了一个有效的补充特征。这提醒我们，灵活地创造性地运用GLCM参数和特征，往往比死板地套用标准流程更能解决实际问题。纹理分析的世界没有标准答案，最好的参数组合永远依赖于你的数据和你要回答的具体问题。