跳转至

第10章 图像特征提取

一、特征提取基础(Fundamentals of Feature Extraction)

1.1 什么是特征

特征(Feature)是图像中我们感兴趣的、可以定量描述的信息片段。特征提取将分割后的 ROI 转化为可分析的数值,是图像分析的关键步骤。

特征提取的两个环节:

环节 英文 含义
特征检测 Feature Detection 在边界/区域/图像中找到特征
特征描述 Feature Description 为检测到的特征赋予定量属性

1.2 特征的性质

良好的特征应当具备以下性质:

性质 含义 示例
独立性(Independent) 各特征之间互不冗余
不变性(Invariant) 在平移/旋转/反射下保持不变 面积、圆度
协变性(Covariant) 在缩放下按比例变化 周长(缩放 \(s\) → 周长 \(\times s\)

1.3 特征分类总览

特征类型 考虑的信息 计算方法 典型特征
形状特征(Shape) 仅空间信息(mask) 基于 ROI 的二值掩膜 面积、周长、傅里叶描述子
强度特征(Intensity) 仅灰度信息(histogram) 基于直方图统计 均值、方差、偏度、峰度
纹理特征(Texture) 灰度 + 空间信息 GLCM、频谱分析 对比度、同质性、能量

二、形状特征(Shape Features)

形状特征仅依赖 ROI 的空间几何信息,与灰度分布无关。分为边界描述子区域描述子两大类。

2.1 基本边界描述子(Basic Boundary Descriptors)

描述子 定义 公式/说明
边界长度(Length) 边界上的像素数目 边界曲线的近似长度
直径(Diameter) 边界上相距最远的两点距离 \(\text{Diam}(B) = \max_{i,j}[D(p_i, p_j)]\)
主轴(Major Axis) 直径对应的线段 连接两个最远边界点
次轴(Minor Axis) 垂直于主轴的线段
离心率(Eccentricity) 焦距与主轴长度之比 描述形状的扁长程度
曲率(Curvature) 边界斜率的变化率 描述边界的弯曲程度

2.2 弯曲度(Tortuosity)

弯曲度衡量血管等曲线结构的蜿蜒程度,通常定义为曲线实际长度与端点直线距离之比。

2.3 标记(Signatures)

标记(Signature) 是将二维边界降维为一维函数的表示方法。最简单的方式是绘制质心到边界距离 \(r(\theta)\) 随角度 \(\theta\) 变化的函数

标记图中显著峰的数量足以区分不同形状的物体。

2.4 傅里叶描述子(Fourier Descriptors)

傅里叶描述子是边界形状特征中最重要的一类,其核心思想是将二维空间边界转化为一维复信号,再利用 DFT 转换到频域进行分析。

2.4.1 边界的复数表达

设边界由 \(K\) 个坐标点组成,按逆时针遍历:

\[ (x_0, y_0), (x_1, y_1), \ldots, (x_{K-1}, y_{K-1}) \]

将每个坐标对表示为复数:

\[ s(k) = x(k) + j \cdot y(k), \quad k = 0, 1, \ldots, K-1 \]

2.4.2 离散傅里叶变换(DFT)获取描述子

\(s(k)\) 做 DFT,得到傅里叶描述子 \(a(u)\)

\[ a(u) = \frac{1}{K} \sum_{k=0}^{K-1} s(k) \cdot e^{-j \frac{2\pi u k}{K}}, \quad u = 0, 1, \ldots, K-1 \]

2.4.3 逆变换重建边界

通过 IDFT,可以从傅里叶描述子重建边界:

\[ s(k) = \sum_{u=0}^{K-1} a(u) \cdot e^{j \frac{2\pi u k}{K}} \]

2.4.4 截断近似

仅用前 \(P\) 个低频系数来近似重建边界:

\[ \hat{s}(k) = \sum_{u=0}^{P-1} a(u) \cdot e^{j \frac{2\pi u k}{K}} \]
  • 低频分量(小 \(u\)):决定边界的整体宏观形状(大小、位置、大体轮廓)
  • 高频分量(大 \(u\)):刻画边界的微观细节(锯齿、凹陷、突起、噪声)

2.4.5 傅里叶描述子的性质

傅里叶描述子对几何变换具有良好的响应特性:

变换 \(a(u)\) 的影响
平移 \(a(0)\) 改变(加入平移量)
旋转 所有 \(a(u)\) 乘以 \(e^{j\theta}\)(相位旋转)
缩放 所有 \(a(u)\) 乘以缩放因子 \(s\)
起点变化 所有 \(a(u)\) 乘以 \(e^{-j\frac{2\pi u k_0}{K}}\)

利用这些性质可以构造平移、旋转、缩放不变的形状特征(如取 \(|a(u)|\)\(|a(u)|/|a(1)|\))。

2.5 基本区域描述子(Basic Region Descriptors)

描述子 公式 说明
周长(Perimeter) \(p = \text{边界像素数}\) ROI 边界的长度
面积(Area) \(A = \text{区域内像素数}\) ROI 的大小
紧凑度(Compactness) \(p^2 / A\) 值越大越不紧凑
圆度(Circularity) \(4\pi A / p^2\) 完美圆 \(=1\),与方向和位置无关
离心率(Eccentricity) \(\sqrt{1 - (b/a)^2}, \ (a \ge b)\) 基于拟合椭圆,两端:0=正圆,1=极扁

离心率计算(基于协方差矩阵)

\(K\) 个边界点 \(\boldsymbol{z}_k\),求均值 \(\bar{\boldsymbol{z}}\) 和协方差矩阵 \(\boldsymbol{C}\)

\[ \boldsymbol{C} = \frac{1}{K-1} \sum_{k=1}^{K} (\boldsymbol{z}_k - \bar{\boldsymbol{z}})(\boldsymbol{z}_k - \bar{\boldsymbol{z}})^T \]

协方差矩阵的特征向量给出主轴方向,特征值之比决定离心率。

2.6 拓扑描述子(Topological Descriptors)

拓扑性质不依赖于距离概念,在拉伸或旋转变换下保持不变。

欧拉数(Euler Number)

\[ E = C - H \]

其中 \(C\) 为连通分量数(Connected Components),\(H\) 为孔洞数(Holes)。

2.7 特征向量(Feature Vectors)

将多个形状特征组合为特征向量,可以在特征空间中比较不同形状的相似性。

形状 紧凑度 圆度 离心率
圆形 ≈1 ≈0
正方形 中等 <1
星形
水滴形 中等 中等

三、强度特征(Intensity Features)

强度特征仅考虑 ROI 内像素/体素的灰度分布不考虑空间位置信息

3.1 基本统计量

设 ROI 区域为 \(N\),共 \(|N|\) 个像素。

特征 公式 说明
能量(Energy) \(\sum_{(x,y)\in N} f(x,y)^2\) 灰度值平方和
均值(Mean) \(\bar{f} = \frac{1}{\lvert N \rvert} \sum_{(x,y)\in N} f(x,y)\) 整体平均亮度
中位数(Median) \(\text{median}(f)\) ROI 内灰度中值
最大值/最小值 \(\max/\min\) 灰度范围端点
极差(Range) \(\max(f) - \min(f)\) 灰度分布宽度
标准差(Std) \(\sigma = \sqrt{\frac{1}{\lvert N \rvert}\sum(f - \bar{f})^2}\) 对比度/异质性
方差(Variance) \(\sigma^2 = \frac{1}{\lvert N \rvert}\sum(f - \bar{f})^2\) 标准差的平方

3.2 直方图特征

\(z\) 为灰度变量,\(p(z_i)\) 为归一化直方图(\(i = 0, 1, \ldots, L-1\))。

特征 公式 取值范围 说明
均匀度(Uniformity) \(U = \sum_{i=0}^{L-1} p^2(z_i)\) \((0, 1]\) 恒值图像取最大值 1
平均熵(Average Entropy) \(e = -\sum_{i=0}^{L-1} p(z_i) \log_2 p(z_i)\) \([0, +\infty)\) 恒值图像取最小值 0

均匀度高 → 灰度集中分布;熵高 → 灰度分布分散、随机性强。

3.3 \(n\) 阶中心矩(\(n\)-th Moment of Mean)

\[ \mu_n = \frac{1}{|N|} \sum_{(x,y)\in N} \left( f(x,y) - \bar{f} \right)^n \]
阶数 名称 物理含义
\(\mu_2\) 方差(Variance) 对比度 / 灰度分散程度
\(\mu_3\) 偏度(Skewness) 灰度分布关于均值的不对称性
\(\mu_4\) 峰度(Kurtosis) 灰度分布在尾部的权重(高阶矩)

标准化偏度与峰度

\[ \text{Skewness} = \frac{\mu_3}{\sigma^3}, \quad \text{Kurtosis} = \frac{\mu_4}{\sigma^4} \]
  • 偏度 \(>0\):分布右偏(长尾在右侧,大多数像素偏暗但少数高亮)
  • 偏度 \(<0\):分布左偏(长尾在左侧)
  • 峰度大:数据在尾部/离群值较多;峰度小:数据分布较为平缓

3.4 平滑度(Smoothness)

\[ R(z) = 1 - \frac{1}{1 + \sigma^2(z)} \]

对于均匀区域(\(\sigma^2 = 0\)),\(R = 0\);对于高对比度区域(\(\sigma^2\) 大),\(R \to 1\)


四、纹理特征(Texture Features)

纹理特征同时考虑灰度值空间位置关系,是三阶特征中信息量最丰富的一类。

4.1 灰度共生矩阵(GLCM, Gray-Level Co-occurrence Matrix)

GLCM 是描述纹理的二阶统计量:统计图像中成对像素灰度的联合分布。

4.1.1 GLCM 的构建步骤

  1. 确定灰度级数 \(K\)(通常量化\(8 \sim 64\) 级以减小矩阵尺寸)
  2. 确定位置算子 \(Q\)(方向 \(\theta\) + 步长 \(d\)
  3. 初始化 \(K \times K\) 矩阵 \(\boldsymbol{G}\)
  4. 遍历所有像素对:若像素值为 \(z_i\),其 \(Q\) 方向邻居为 \(z_j\),则 \(\boldsymbol{G}[z_i, z_j] \mathrel{+}= 1\)
  5. (可选)归一化:\(p_{ij} = g_{ij} / \sum_{i,j} g_{ij}\)

4.1.2 位置算子 \(Q\) 的定义

常见的 \(Q\) 定义:

方向 \(Q\) 描述 典型应用
水平(0°) 右侧相邻像素 横向纹理检测
垂直(90°) 下方相邻像素 纵向纹理检测
对角(45°) 右下方相邻像素 对角纹理检测
对角(135°) 右上方相邻像素 反斜对角纹理检测

步长 \(d\) 可以是 1、2、3…,步长越大,捕获的纹理尺度越大。

4.1.3 GLCM 的直观理解

  • 随机纹理:GLCM 元素分散分布
  • 周期纹理:GLCM 呈现规律性结构
  • 平滑区域:GLCM 集中在对角线附近(相邻像素灰度相近)

4.2 GLCM 的描述子

基于归一化 GLCM \(p_{ij}\) 可计算以下二阶纹理特征:

描述子 公式 含义 取值范围
最大概率(Max Probability) \(\max_{i,j}(p_{ij})\) GLCM 中最强响应 \([0, 1]\)
相关性(Correlation) \(\frac{\sum(i-m_r)(j-m_c)p_{ij}}{\sigma_r \sigma_c}\) 像素灰度的线性相关程度 \([-1, 1]\)
对比度(Contrast) \(\sum_{i,j} (i-j)^2 p_{ij}\) 局部灰度变化剧烈程度 \([0, (K-1)^2]\)
能量/均匀度(Energy) \(\sum_{i,j} p_{ij}^2\) 纹理的规则性和稳定性 \([0, 1]\)
同质性(Homogeneity) \(\sum_{i,j} \frac{p_{ij}}{1+\lvert i-j \rvert}\) 局部均匀程度 \([0, 1]\)
熵(Entropy) \(-\sum_{i,j} p_{ij} \log_2 p_{ij}\) GLCM 元素的随机性 \([0, 2\log_2 K]\)

其中 \(m_r, m_c\) 是沿行和沿列的均值,\(\sigma_r, \sigma_c\) 是对应的标准差。

各描述子的典型关系

  • 平滑图像 → GLCM 集中对角 → 对比度小、同质性大、能量大
  • 粗糙/随机纹理 → GLCM 分散 → 对比度大、同质性小、能量小
  • 周期性纹理 → GLCM 呈规则结构 → 相关性高

4.3 频谱纹理(Spectral Texture)

从傅里叶频谱也可提取纹理特征。将频谱 \(S(u,v)\) 转换为极坐标 \(S(r, \theta)\)

\[ S(r) = \sum_{\theta=0}^{\pi} S_\theta(r), \quad S(\theta) = \sum_{r=1}^{R} S_r(\theta) \]
  • \(S(r)\):反映不同空间尺度的能量分布 → 纹理粗细
  • \(S(\theta)\):反映不同方向的能量分布 → 纹理方向性

在处理含周期噪声的图像时,频谱中会出现对应频率的尖峰,可通过陷波滤波器消除周期性成分后,再用统计方法描述剩余的非周期纹理。


五、影像组学(Radiomics)

5.1 什么是影像组学

影像组学(Radiomics)是将医学数字图像转化为可挖掘的高维数据的方法,通过定量分析反映潜在的病理生理特征。

5.2 影像组学工作流程

步骤 内容
1. 图像获取 MR / PET / CT,含分期、示踪剂、对比相信息
2. 图像分割 手动/自动分割 ROI(由放射科医师或分割模型完成)
3. 插值(重采样) 最近邻 / 双线性 / 三线性插值,统一体素尺寸
4. 灰度离散化 确定 bin width 或 bin number
5. 特征提取 形状 + 强度 + 纹理 + 滤波特征
6. 分析与建模 分类 / 回归 / 生存分析

5.3 滤波增强特征

在提取特征前,可对图像施加多种滤波,每个滤波输出都单独计算特征,极大扩充特征池:

滤波类型 示例
统计滤波 中值滤波
边缘滤波 LoG(高斯拉普拉斯)滤波
特殊滤波 分形维数(Fractal Dimension)滤波

5.4 应用场景

应用 方法 任务类型
组织状态评估 SVM 分类
肿瘤类型预测 随机森林(Random Forest) 分类
无进展生存期(PFS)预测 LASSO 回归 回归
生存风险分层 生存分析 预后分析

5.5 工具与局限

常用工具

  • MATLAB:传统影像组学工具
  • PyRadiomics:Python 开源影像组学库,标准化特征提取

局限性

  • 可重复性差(缺乏标准化、报告不充分、开源代码和数据有限)
  • 缺乏适当的外部验证,假阳性结果风险高
  • 特征的可解释性不足
  • 多为回顾性研究,证据等级较低

六、本章总结

特征提取三大类别

类别 核心工具 考虑的信息 典型应用
形状特征 边界/区域描述子、傅里叶描述子 空间(mask) 肿瘤形态学分析
强度特征 直方图统计、\(n\) 阶矩 灰度(histogram) 组织密度评估
纹理特征 GLCM、频谱分析 灰度 + 空间 组织异质性量化

傅里叶描述子总结

属性 说明
核心变换 边界坐标 \((x,y)\) → 复数序列 \(s(k)\) → DFT → \(a(u)\)
低频作用 \(a(0)\) = 质心位置;\(a(1)\) = 等效半径;\(a(0)\) + \(a(1)\) = 圆
高频作用 刻画凹凸细节,极高频率 = 噪声
截断策略 \(P\) 小 → 平滑宏观形状;\(P\) 大 → 精确但含噪声
不变性构造 $

GLCM 核心指标速查

指标 物理含义 高值 → 低值 →
Contrast 局部灰度变化剧烈程度 粗糙纹理、边界多 平滑区域
Homogeneity 局部均匀度 平滑、同质区域 粗糙纹理
Energy 纹理规则度/稳定性 规则重复纹理 随机纹理
Entropy GLCM 随机性 复杂随机纹理 规则结构
Correlation 灰度线性相关 线性纹理结构 无关纹理

七、课后作业(HW6)

题目 1:Signature 与傅里叶描述子

\[ r(\theta) = 2 - 2\sin\theta + \frac{\sin\theta \sqrt{|\cos\theta|}}{\sin\theta + 1.5}, \quad \theta = \frac{2\pi n}{1280}, \quad n = 0, 1, \ldots, 1279 \]

解答

(1) 真实 Signature 判断

答案为 A。需要从肿瘤轮廓图中识别对应的 \(r(\theta)\) 曲线。

(2) 绘制边界

使用 Python 计算 \(x = r(\theta)\cos\theta\)\(y = r(\theta)\sin\theta\),直接绘制即可。

(3) 傅里叶描述子重建——第一部分

将边界坐标表示为复数序列 \(s(n) = x(n) + j y(n)\),做 DFT 得到 \(a(u)\)

\[ a(u) = \frac{1}{N} \sum_{n=0}^{N-1} s(n) e^{-j\frac{2\pi u n}{N}} \]

通过 Python 计算得:

\[ a(0) = -0.7j, \quad a(1) = 1.8 \]

仅保留 \(a(0)\)\(a(1)\) 后,逆变换重建:

\[ \begin{aligned} s'(n) &= a(0) + a(1) e^{j\theta_n}, \quad \theta_n = \frac{2\pi n}{1280} \\ &= -0.7j + 1.8(\cos\theta_n + j\sin\theta_n) \\ &= 1.8\cos\theta_n + j(-0.7 + 1.8\sin\theta_n) \end{aligned} \]

拆分为直角坐标:

\[ \begin{cases} x'(n) = 1.8 \cos\theta_n \\ y'(n) = -0.7 + 1.8 \sin\theta_n \end{cases} \]

消去 \(\theta_n\)(即 \(x'^2 + (y' + 0.7)^2 = 1.8^2\)):

\[ \boxed{x^2 + (y + 0.7)^2 = 1.8^2} \]

结果:仅用 \(a(0)\)\(a(1)\) 两个系数,重建出的边界是一个圆心 \((0, -0.7)\)、半径 \(1.8\) 的圆——低频系数捕捉了边界的全局位置和大小,但丢失了所有细节。

(4) 傅里叶描述子重建——第二部分

新的极坐标方程 \(\rho(\theta) = \frac{1}{|\cos\theta| + |\sin\theta|}\) 在极坐标下是一个正方形。同样构建复数序列、DFT、仅保留 \(a(0)\)\(a(1)\)

通过 Python 计算得:

\[ a(0) = 0, \quad a(1) = 0.79 \]

重建:

\[ \begin{cases} x'(n) = 0.79 \cos\theta_n \\ y'(n) = 0.79 \sin\theta_n \end{cases} \]

即:

\[ \boxed{x^2 + y^2 = 0.79^2} \]

同样是一个圆。这展示了傅里叶描述子的核心特性:仅用最低频的两个系数 \(a(0)\)\(a(1)\) 只能恢复一个圆——\(a(0)\) 决定圆心(此处为原点),\(a(1)\) 决定半径。

傅里叶描述子的频率层级直觉

  • \(a(0)\) 仅描述边界的质心位置
  • \(a(1)\) 描述大致尺寸(等效半径)
  • \(a(0)\)\(a(1)\) 合在一起 = 一个
  • 需要更高阶的 \(a(2), a(3), \ldots\) 才能描述椭圆变形、不规则突触等复杂形状

题目 2:一阶统计特征

题目(PPT HW2):

给定图像 \(I_0\)

\[ I_0 = \begin{bmatrix} 3 & 1 & 0 & 0 & 2 \\ 1 & 0 & 1 & 1 & 0 \\ 0 & 2 & 0 & 3 & 1 \\ 1 & 2 & 1 & 0 & 0 \\ 0 & 0 & 2 & 3 & 1 \end{bmatrix} \]
  1. 计算 \(I_0\)energy、range、\(\mu_1\)\(\mu_2\)\(\mu_3\)、uniformity 和 entropy

  2. 假设图像采用循环填充(circular padding),每一步做 \(3 \times 3\) 均值滤波的循环卷积。描述该扩散过程的最终状态,并选择最恰当的描述

  3. 填写 \(I_\infty\) 的各项指标

解答

(1) 各指标的逐项计算

先统计灰度频数:

灰度 \(z\) 0 1 2 3
频数 10 8 4 3
概率 \(p(z)\) 0.4 0.32 0.16 0.12
  • Energy\(E = \sum f^2 = 3^2 \times 3 + 2^2 \times 4 + 1^2 \times 8 + 0^2 \times 10 = 27 + 16 + 8 + 0 = 51\)

  • Range\(\max - \min = 3 - 0 = 3\)

  • \(\mu_1\)(均值)

\[ \mu_1 = \sum z_i p(z_i) = 0 \times 0.4 + 1 \times 0.32 + 2 \times 0.16 + 3 \times 0.12 = 1.0 \]
  • \(\mu_2\)(方差)
\[ \mu_2 = \sum (z_i - 1)^2 p(z_i) = 1 \times 0.4 + 0 \times 0.32 + 1 \times 0.16 + 4 \times 0.12 = 0.4 + 0.16 + 0.48 = 1.04 \]
  • \(\mu_3\)(三阶中心矩)
\[ \mu_3 = \sum (z_i - 1)^3 p(z_i) = (-1) \times 0.4 + 0 \times 0.32 + 1 \times 0.16 + 8 \times 0.12 = -0.4 + 0.16 + 0.96 = 0.72 \]

因此偏度(Skewness)= \(\mu_3 / \sigma^3 = 0.72 / (1.04)^{3/2} \approx 0.68\) > 0,表明灰度分布右偏

  • Uniformity
\[ U = \sum p(z_i)^2 = 0.4^2 + 0.32^2 + 0.16^2 + 0.12^2 = 0.16 + 0.1024 + 0.0256 + 0.0144 = 0.3024 \]
  • Entropy
\[ \begin{aligned} e &= -(0.4\log_2 0.4 + 0.32\log_2 0.32 + 0.16\log_2 0.16 + 0.12\log_2 0.12) \\ &\approx 1.8449 \end{aligned} \]

(2) 扩散过程的最终状态

每次 \(3 \times 3\) 均值滤波(循环卷积)将每个像素替换为周围 9 个像素的平均值。由于使用循环填充,总亮度守恒。

随着 \(t \to \infty\),所有像素趋近于同一个值——原图像的全局均值 \(\mu_1 = 1.0\)。这符合选项 C 的描述。

(3) \(I_\infty\) 的指标

Image energy range \(\mu_1\) \(\mu_2\) \(\mu_3\) uniformity entropy
\(I_\infty\) 25 0 1.0 0 0 1 0

解释\(I_\infty\) 是恒值图像(所有像素 \(=1.0\)),因此:

  • energy = \(25 \times 1^2 = 25\)
  • range = 0(无灰度变化)
  • \(\mu_2, \mu_3\) = 0(无偏离均值的像素)
  • uniformity = 1(仅一个灰度级)
  • entropy = 0(无不确定性)

题目 3:GLCM(灰度共生矩阵)

题目(PPT HW3):

给定矩阵:

\[ M = \begin{bmatrix} 0 & 1 & 2 & 1 & 0 \\ 1 & 2 & 1 & 2 & 1 \\ 0 & 1 & 2 & 1 & 0 \\ 1 & 2 & 1 & 2 & 1 \\ 0 & 1 & 2 & 1 & 0 \end{bmatrix} \]
  1. 位置算子 \(Q\) = "一个像素右侧" → 计算 GLCM
  2. 位置算子 \(Q\) = "两个像素右侧" → 计算 GLCM
  3. 分别计算两者的 ContrastHomogeneity

解答

(1) \(Q\) = "一个像素右侧"(\(d=1\)

逐行统计相邻像素对 \((i, j)\),共 5 行 \(\times\) 4 对 = 20 对。

统计计数矩阵:

\[ M_1 = \begin{bmatrix} 0 & 3 & 0 \\ 3 & 0 & 7 \\ 0 & 7 & 0 \end{bmatrix} \]

矩阵大小是 \(3 \times 3\)(因灰度值仅含 0, 1, 2)。归一化:

\[ P_1 = \frac{M_1}{20} = \begin{bmatrix} 0 & 0.15 & 0 \\ 0.15 & 0 & 0.35 \\ 0 & 0.35 & 0 \end{bmatrix} \]

(2) \(Q\) = "两个像素右侧"(\(d=2\)

逐行统计间隔一个像素的组合对 \((i, j)\),共 5 行 \(\times\) 3 对 = 15 对。

统计计数矩阵:

\[ M_2 = \begin{bmatrix} 0 & 0 & 3 \\ 0 & 7 & 0 \\ 3 & 0 & 2 \end{bmatrix} \]

归一化:

\[ P_2 = \frac{M_2}{15} = \begin{bmatrix} 0 & 0 & 0.2 \\ 0 & 7/15 & 0 \\ 0.2 & 0 & 2/15 \end{bmatrix} \]

(3) Contrast 和 Homogeneity

使用公式:

\[ \text{Contrast} = \sum_{i,j} (i-j)^2 P(i,j), \quad \text{Homogeneity} = \sum_{i,j} \frac{P(i,j)}{1 + |i-j|} \]

GLCM 1(\(P_1\):所有非零项均满足 \(|i-j| = 1\)

  • Contrast = \(1^2 \times (0.15 + 0.15 + 0.35 + 0.35) = 1.0\)
  • Homogeneity = \(\frac{0.15}{2} + \frac{0.15}{2} + \frac{0.35}{2} + \frac{0.35}{2} = 0.5\)

GLCM 2(\(P_2\)\((0,2)\)\((2,0)\)\(|i-j| = 2\)\((1,1)\)\((2,2)\)\(|i-j| = 0\)

  • Contrast = \(4 \times (0.2 + 0.2) + 0 \times (7/15 + 2/15) = 1.6\)
  • Homogeneity = \(\frac{0.2}{3} + \frac{0.2}{3} + \frac{7/15}{1} + \frac{2/15}{1} \approx 0.13 + 0.47 + 0.13 \approx 0.73\)

对比分析

指标 \(d=1\) \(d=2\)
Contrast 1.0 1.6
Homogeneity 0.5 0.73

步长增大 → Contrast 增大、Homogeneity 也增大?这里的 Homogeneity 增大是因为矩阵 \(M\) 具有 2 像素周期的结构。 具体来说,\(d=1\) 时相邻像素对的灰度差始终为 1;\(d=2\) 时出现了 \((1,1)\)\((2,2)\) 的对角线项(\(|i-j|=0\)),显著提升了 Homogeneity。这个现象说明采样步长与纹理周期会产生共振


八、实验(LAB7)

实验来源:LAB7(Project 1-3)

Project 1:傅里叶描述子与形状保真度

实验内容:提取真实脑肿瘤 MRI 的轮廓,计算傅里叶描述子,逐步增加保留的描述子数量 \(P\)\(P = 2, 4, 8, 16, 32, 64, 128, N/2\)),重建肿瘤轮廓并观察形状保真度的变化。

实验原理

  1. 轮廓复数化\(s(n) = x(n) + j \cdot y(n)\)\(n = 0, 1, \ldots, N-1\)
  2. DFT 获取描述子\(a(k) = \frac{1}{N} \sum_{n=0}^{N-1} s(n) e^{-j \frac{2\pi k n}{N}}\)
  3. 截断与重建:通过 fftshift 将低频分量集中到频谱中心,仅保留中心 \(P\) 个系数(其余置零),再做 ifftshift + ifft 重建轮廓

关键发现—\(P\) 与保真度的三个阶段

阶段 \(P\) 范围 重建效果 物理含义
低频粗糙 \(2 \sim 8\) \(P=2\) 为圆,\(P=4/8\) 为粗糙块状多边形 仅捕捉位置、大小和大致方向
核心特征 \(16 \sim 64\) 快速捕获凹陷、突起和不规则拐角,保真度迅速提升 临床评估肿瘤体积和侵袭方向的最佳平衡区
细节饱和 \(\ge 128\) 肉眼与原始轮廓几乎无法区分,增加 \(P\) 后收益递减 剩余高频 = 像素级噪声

核心结论\(P\) 与形状保真度存在非线性收益递减关系——中等 \(P\) 是数据压缩与几何精度的最优权衡点。


Project 2:一阶统计特征与迭代模糊

实验内容

  1. 基于肿瘤掩膜计算质心,以质心为中心提取 \(120 \times 120\) 的 ROI
  2. 计算 ROI 与全图的均值、方差、均匀度、熵
  3. 执行 30 次迭代 \(3 \times 3\) 均值滤波(模拟渐进高斯模糊),记录四项指标的动态演变
  4. 对比 ROI 与全图的指标演变差异

实验原理

  1. 质心计算(基于图像矩):
\[ \bar{x} = \frac{M_{10}}{M_{00}}, \quad \bar{y} = \frac{M_{01}}{M_{00}} \]
  1. 一阶统计特征:基于直方图概率 \(p(z_i)\) 计算(参见第三章公式)

  2. 迭代滤波模拟高斯模糊:基于中心极限定理,多次 \(3 \times 3\) 均值滤波等效于逐渐增大高斯核的标准差

关键发现

(I) ROI 指标随滤波迭代的单调演变

指标 趋势 物理解释
均值 稳定不变 平滑是能量重新分配,总亮度守恒
方差 持续下降(1408 → 1023) 像素差异被拉近,对比度减弱
均匀度 持续上升(0.0081 → 0.0094) 直方图向均值聚拢
持续下降(7.12 → 6.89) 细节丢失,纹理复杂度降低

(II) ROI vs 全图的显著差异

对比维度 ROI(肿瘤区域) 全图
均值 76.24(高) 33.09(低,因大量黑色背景)
方差 1407(中等) 2177(高,因颅骨极亮与背景极暗的反差)
均匀度 0.008(低) 0.109(高,背景峰主导)
7.12(高) 5.17(低,背景峰高同质化)

启示:精准的 ROI 提取是防止背景干扰、获取有效特征的绝对前提。全图在模糊过程中还会发生组织与背景之间的灰度扩散渗透,导致一致性不升反降——这是宏观边界效应。


Project 3:GLCM 纹理分析与采样步长

实验内容

  1. 将肿瘤 ROI 量化到 16 个灰度级
  2. 计算 10 个不同采样步长(\(d = 1, 2, 3, 5, 8, 12, 15, 20, 25, 30\))下水平方向的 GLCM
  3. 提取对比度(Contrast)、同质性(Homogeneity)、能量(Energy)随 \(d\) 的变化曲线
  4. 可视化 \(d=1\)\(d=15\) 的 256 级原图 GLCM 矩阵(对数变换增强显示)

实验原理

  1. 灰度量化:将 256 级降为 16 级,使统计更集中、特征更具代表性
  2. GLCM 计算skimage.feature.graycomatrix,对称+归一化
  3. GLCM 可视化np.log1p 对数变换突出暗部细节,对角线附近亮度 = 像素灰度相近的频率

关键发现

(I) 纹理特征随步长的单调演变

步长 \(d\) Contrast Homogeneity Energy
1 0.30 0.87 0.27
30 6.50 0.42 0.14
  • Contrast 单调上升:步长越大,采样的两个像素跨越越大的组织差异
  • Homogeneity 单调下降:非对角线元素比例增加
  • Energy 单调下降:GLCM 集中度降低,纹理越「不规则」

(II) 步长与周期性纹理的「共振交互」

如果 ROI 内存在周期性结构纹理(空间周期 \(T\)):

  • \(d = T\)(同相采样)→ 两个像素灰度一致 → Contrast 骤降至局部极小 → Homogeneity 骤升至局部极大
  • \(d = 0.5T\)(反相采样)→ 波峰对波谷 → Contrast 极大

本实验肿瘤数据的曲线平滑无振荡 → 脑肿瘤内部纹理是高度无序的,不存在宏观周期性结构。

(III) GLCM 可视化的空间直觉

  • \(d=1\):GLCM 呈现极其明亮的细对角线 → 能量高度集中 → 近邻像素灰度几乎一致
  • \(d=15\):对角线「散焦」弥漫到非对角线区域 → 能量发散 → 远距离像素灰度关系弱化

九、历年卷解答

一、三阶矩与偏度(2021 选择)

知识点定位:强度特征 — \(n\) 阶中心矩

题目:"三阶矩计算的是什么特征?" 选项:均值 / 方差 / 偏度 / 峰度

答案偏度(Skewness)

解释

阶数 特征
一阶中心矩 \(\mu_1\) 恒为 0
二阶中心矩 \(\mu_2\) 方差(Variance)/ 对比度
三阶中心矩 \(\mu_3\) 偏度(Skewness)——分布的不对称性
四阶中心矩 \(\mu_4\) 峰度(Kurtosis)——尾部的权重

二、峰度(Kurtosis)与离群噪声(2022 判断)

知识点定位:强度特征 — 四阶矩 / 峰度

题目:"低信号的肿瘤区域中,存在零散的高信号噪声点,会让峰度(Kurtosis,四阶矩)变大。"

答案正确

解释:峰度(\(\mu_4 / \sigma^4\))对远离均值的离群点(Outlier)极度敏感——因为它以 \((f - \bar{f})^4\) 加权。低信号区域中零星出现的高信号噪声点,在四次方放大下会极大拉升 \(\mu_4\),从而使峰度显著增大。

峰度的直觉

峰度大 \(\neq\) "尖峰",而是意味着分布尾部厚重,即存在大量/大幅度的离群值。在医学图像中,异常的峰度往往是噪声或者是病理异常(如钙化点)的指示。


三、GLCM 与周期噪声的关系(2022 选择)

知识点定位:纹理特征 — GLCM / 频谱纹理

题目:"给一个存在周期噪声的 MRI 图像,问哪个是可能的 GLCM。"

答案:选择形如细长叶子状的 GLCM(集中在一条对角线 + 向上/向下弯曲的曲线)。

解释

  • 周期噪声在图像中表现为规律性的明暗条纹
  • 在 GLCM 中,周期噪声的规律性结构产生非对角线的规则分支
  • 纯随机纹理 → GLCM 分散;平滑图像 → GLCM 集中在对角线
  • 周期结构 → GLCM 呈现对角线 + 规则侧支的叶状/羽状结构

自编计算题

一、傅里叶描述子的系数含义(自编计算题)

题目:某闭合边界的傅里叶描述子中仅 \(a(0) = 2 + 3j\)\(a(1) = 5\),其余系数均为 0。问:

  1. 该边界在复平面上是什么形状?
  2. 该形状的质心位置和等效半径是多少?

解答

仅保留 \(a(0)\)\(a(1)\) 的重建公式:

\[ s'(k) = a(0) + a(1) e^{j \frac{2\pi k}{K}} \]

代入 \(a(0) = 2 + 3j\)\(a(1) = 5\)

\[ \begin{aligned} s'(k) &= (2 + 3j) + 5(\cos\theta_k + j\sin\theta_k) \\ &= (2 + 5\cos\theta_k) + j(3 + 5\sin\theta_k) \end{aligned} \]

即:

\[ \begin{cases} x(k) = 2 + 5\cos\theta_k \\ y(k) = 3 + 5\sin\theta_k \end{cases} \]

答案

  1. 这是一个(当仅保留 \(a(0)\)\(a(1)\) 时,重建结果恒为圆)
  2. 圆心在 \((2, 3)\),半径 \(= |a(1)| = 5\)

二、GLCM 计算练习(自编计算题)

题目:给定 \(4 \times 4\) 图像,灰度级为 \(\{0, 1\}\)(二值):

\[ I = \begin{bmatrix} 0 & 1 & 0 & 1 \\ 1 & 0 & 1 & 0 \\ 0 & 1 & 0 & 1 \\ 1 & 0 & 1 & 0 \end{bmatrix} \]

\(Q\) = "右侧相邻像素" 计算 GLCM,并求 Contrast 和 Energy。

解答

逐行统计:

  • 行 1:\((0,1), (1,0), (0,1)\)\((0,1) \times 2\), \((1,0) \times 1\)
  • 行 2:\((1,0), (0,1), (1,0)\)\((1,0) \times 2\), \((0,1) \times 1\)
  • 行 3:同第 1 行
  • 行 4:同第 2 行

总计:\((0,1)\) = \(2+1+2+1 = 6\)\((1,0)\) = \(1+2+1+2 = 6\),共 12 对。

计数矩阵:

\[ G = \begin{bmatrix} 0 & 6 \\ 6 & 0 \end{bmatrix}, \quad P = \begin{bmatrix} 0 & 0.5 \\ 0.5 & 0 \end{bmatrix} \]
  • Contrast = \((0-1)^2 \times 0.5 + (1-0)^2 \times 0.5 = 0.5 + 0.5 = 1.0\)
  • Energy = \(0^2 + 0.5^2 + 0.5^2 + 0^2 = 0.5\)

二值棋盘格纹理的 GLCM 特征

此例中的图像是完美的二值棋盘格,相邻像素永远不同\(0 \to 1\)\(1 \to 0\))。因此 GLCM 对角线为 0(没有相同值的相邻对),所有非零元素集中在反对角线上。这是检测棋盘格/条纹纹理的极具判别力的特征。