2-53 一张灰度纹理图，怎样一次提取 15 个灰度—梯度统计特征？

MATLAB 灰度梯度共生矩阵 + 15 维图像纹理特征提取程序

一句话介绍： 这套 MATLAB 程序读取一张灰度纹理图像，先计算水平与垂直方向的局部梯度并压缩成 32 个梯度等级，再建立 256×32 的灰度—梯度共生矩阵，最终一次输出小梯度优势、大梯度优势、能量、均值、熵、惯性、逆差矩等 15 个纹理特征。

先看它解决什么问题

只看一张灰度图的平均亮度，很难说明它的纹理究竟是平滑、粗糙、均匀，还是包含大量强边缘。

而纹理真正有价值的信息往往同时藏在两个维度里：

像素本身有多亮；

像素附近变化有多剧烈。

当前程序没有只统计灰度，也没有直接把整张图送入分类器，而是先把这两个维度组合起来：

灰度图像 → 局部梯度 → 梯度离散化 → 灰度—梯度共生矩阵 → 15 个统计纹理特征 → 1×15 特征向量

所以它更准确地说是一套：面向灰度图像离线分析的灰度—梯度联合纹理特征提取程序。

1. 当前工程真正由哪些文件组成？

压缩包中的工程比较精简，核心文件如下。

当前真正被 dutu.m 调用的处理链只有：

D68.gif→double转换 →GrayGradinet→OUT

也就是说，项目说明中的 15 种纹理特征，实际对应的是 GrayGradinet.m 输出的 T1～T15。

GrayExitMonentR.m 和 texture.m 虽然也存在于工程中，但当前主脚本并没有调用它们，因此不能把它们的 6 个灰度共生矩阵统计量混入当前 15 维输出。

2. 当前程序分析的是怎样一张图？

工程自带输入图像尺寸为：

640×640 像素

它是一张灰度纹理图，像素范围覆盖 0～255，当前文件中实际出现 88 个不同灰度值。

主脚本首先执行：

img=imread('D68.gif');

随后把图像转换为双精度数据：

img=double(img);

这个转换并不只是格式变化。

GrayGradinet.m 后面需要直接做相邻像素相减和平方。如果按照当前主入口先转换为 double，梯度差值就可以正常保留正负和幅值，不会在整数运算阶段丢失需要的差分信息。

因此当前输入链可以概括成：

8 位灰度图 → 双精度灰度矩阵 → 梯度与纹理统计

3. 为什么先计算梯度，而不是直接统计灰度？

同样的灰度直方图，可以对应完全不同的空间纹理。

例如两张图可能拥有相似的黑、白、灰像素比例，但一张变化缓慢，另一张包含密集条纹。

仅统计灰度时，这种差别可能体现得并不充分。

当前程序因此对每个有效位置计算一个局部梯度：

GM(i,j) = √[(I(i,j+1) − I(i,j))² + (I(i+1,j) − I(i,j))²]

这里同时比较：

当前像素与右侧像素；

当前像素与下方像素。

然后对两个方向的差值平方求和并开平方。

简单理解：

灰度告诉程序“这里是什么亮度”，梯度告诉程序“这里变化得有多快”。

对于 640×640 的 D68.gif，最后一行和最后一列无法同时取得右侧、下方邻点，因此得到的梯度矩阵尺寸为：639×639

一共包含：408321 个有效灰度—梯度位置。

按附件图像复算，当前梯度范围约为：0～200.09

4. 为什么还要把梯度压缩成 32 个等级？

原始梯度是连续数值，直接与 256 个灰度等级组合，会使统计空间比较分散。

当前源码首先找到梯度矩阵中的：Gmin 和 Gmax

然后把整个梯度范围线性映射到：0～31

也就是固定的：new_gray=32

其关系可以写成：

Gq = uint8[(GM − Gmin) / (Gmax − Gmin) × 31]

因此后续程序不再直接使用原始梯度，而使用 32 级离散梯度。

这样一来，每个有效像素位置都可以表示成一对状态：

当前灰度等级 i + 当前梯度等级 j

其中：灰度 i：0～255；梯度 j：0～31。

这一步也是当前程序中最值得关注的尺度参数之一：梯度等级越多，梯度描述越细；等级越少，统计会更加集中。

5. 256×32 的灰度—梯度共生矩阵怎样建立？

有了灰度和离散梯度之后，程序建立：

256×32 的矩阵 H

矩阵中的每个元素 H(i,j) 表示：

灰度等级 i 与梯度等级 j 在图像中共同出现了多少次。

程序遍历 639×639 个有效位置，每遇到一次对应组合，就把该位置的计数加 1。

因此：

横向 32 列描述梯度强弱，纵向 256 行描述灰度高低。

接着使用全部有效位置数量 N 对矩阵归一化：

P(i,j) = H(i,j) / N

当前示例中：N = 639×639 = 408321

这样就得到后续 15 个特征真正使用的归一化灰度—梯度矩阵。

如果再定义：

pᵢ = Σⱼ P(i,j)

qⱼ = Σᵢ P(i,j)

那么：

pᵢ 表示灰度边缘分布；

qⱼ 表示梯度边缘分布。

后面的均值、方差和熵等统计量，基本都围绕 P、pᵢ 和 qⱼ 展开。

6. 这 15 个纹理特征分别描述什么？

这里必须以当前 GrayGradinet.m 的实际实现为准，而不是直接套用其他文献中同名特征的公式。

几个最核心的关系如下。

小梯度优势：

T1 = Σⱼ qⱼ / (j+1)²

当前代码给低梯度等级更高权重，因此低梯度区域越集中，T1 越容易增大。

大梯度优势：

T2 = Σⱼ j·qⱼ

这里直接对梯度等级做线性加权，所以高梯度区域增多会推动 T2 增大。

能量：

T5 = ΣᵢΣⱼ P(i,j)²

联合分布越集中在少数灰度—梯度组合上，能量通常越高。

灰度平均与梯度平均：

T6 = Σᵢ i·pᵢ

T7 = Σⱼ j·qⱼ

因此按当前代码，T2 与 T7 在数值上实际上采用了相同的梯度加权平均关系。

灰度和梯度离散程度：

T8 = √[Σᵢ (i−T6)²pᵢ]

T9 = √[Σⱼ (j−T7)²qⱼ]

虽然源码注释写的是“均方差”，但因为最后执行了平方根，当前实现形式更接近加权标准差。

7. 熵、相关、惯性和逆差矩又在看什么？

前三个熵特征分别从不同范围描述分布的不确定程度。

灰度熵：

T11 = −Σᵢ pᵢ·log₁₀(pᵢ+eps)

梯度熵：

T12 = −Σⱼ qⱼ·log₁₀(qⱼ+eps)

混合熵：

T13 = −ΣᵢΣⱼ P(i,j)·log₁₀[P(i,j)+eps]

这里有一个复现时很重要的细节：

当前 15 特征函数使用的是 log10，不是 log2。

因此如果换成其他教材或代码中的二进制熵公式，最终数值不会与当前程序完全一致。

T10 在源码中标注为“相关”，当前计算关系为：

T10 = ΣᵢΣⱼ (i−T6)(j−T7)P(i,j)

它没有继续除以灰度和梯度标准差，所以严格按当前实现理解，它更接近一个协方差型的灰度—梯度相关量，而不是标准化到固定范围内的相关系数。

惯性 T14 为：

T14 = ΣᵢΣⱼ (i−j)²P(i,j)

逆差矩 T15 为：

T15 = ΣᵢΣⱼ P(i,j) / [1+(i−j)²]

两者都在利用灰度索引与梯度索引之间的距离，只是：

T14 对距离大的组合赋予更高权重；

T15 对距离大的组合进行衰减。

当前实现中灰度索引范围是 0～255，而梯度索引范围是 0～31，因此复现现有程序时应保留这一尺度关系，不能自行把两者重新归一化后再计算。

8. 自带 D68.gif 最终得到怎样的 15 维特征？

依据当前 dutu.m、GrayGradinet.m、附件中的 D68.gif 以及源码使用的梯度量化规则复算，当前示例得到的 15 维结果约为：

最终变量：OUT就是一个1×15 特征向量。

这里尤其需要注意 T3 和 T4。

当前源码计算的是“边缘计数平方后只除一次样本总数”，等价于：

T3 = N·Σᵢ pᵢ²

T4 = N·Σⱼ qⱼ²

所以它们会随着有效统计位置 N 的变化而改变尺度，并不是简单的归一化概率平方和。

因此如果拿不同尺寸图像直接比较 T3、T4，应当确认图像尺寸和统计口径是否一致。

这些结果最适合作为：核对当前程序是否执行到同一算法实现的参考值，而不是所有纹理图像都应该接近的固定标准。

9. 为什么工程里还有传统灰度共生矩阵函数？

除了当前主函数，工程还保存了：

GrayExitMonentR.m以及：texture.m

其中 texture.m 可以按照距离 d 和方向 theta 建立灰度共生计数，代码支持：

0°、45°、90°、135°

GrayExitMonentR.m 则包含二阶矩、相关度、熵、对比度、逆差矩和和方差等 6 个统计量。

不过当前 dutu.m 并没有把这两个函数连接到主流程。

因此当前项目不能描述成：

“15 个灰度—梯度特征 + 6 个灰度共生矩阵特征全部自动计算”

更准确的工程事实是：

主入口自动完成的是 GrayGradinet.m 的 15 维灰度—梯度特征；另外两个文件为工程中保留的灰度共生矩阵相关算法代码，可供后续单独研究或扩展。

10. 这套程序适合用在哪里？

灰度纹理图像的离线特征提取

程序已经包含从原始灰度图到 15 维数值特征的完整链路，适合对表面纹理、材料纹理以及其他灰度纹理图进行离线统计分析。

它的输出不是新的图像，而是一个固定长度的数值向量，因此特别适合作为后续数据分析的输入。

图像分类或识别前的特征构造

如果后续还有分类器，可以把每张图像得到的：T1～T15

作为一组人工纹理特征。

当前工程本身没有包含分类器、训练集划分或识别模型，所以它负责的是：

图像 → 特征

而不是：

图像 → 自动类别判定。

纹理算法教学与实验对比

程序把梯度计算、梯度量化、共生统计和 15 个指标全部直接写在 MATLAB 源码中，没有封装成不可见的黑盒函数。

因此很适合观察：

图像纹理变化 → 灰度/梯度联合分布变化 → 特征值变化这一过程。

小规模图像特征研究

当前实现大量使用 MATLAB 双层循环，结构直观，但并不是面向高吞吐实时处理设计的。

因此它更适合：

单图实验、小批量离线分析、课程设计和算法验证

而不是直接作为工业在线视觉平台。

11. 换成自己的图像时，哪些地方最值得注意？

真正影响当前算法结果的关键设置并不多。

如果切换自己的图片，最简单的方式就是修改 dutu.m 中的图像文件名，并继续保持：

单通道灰度图 + 0～255 整数灰度值

这一输入形式。

当前 GrayGradinet 直接使用：

IN(i,j)+1

作为共生矩阵行索引，因此不要直接把已经归一化到 0～1 的非整数浮点图像送入函数。

彩色 RGB 图也不能直接按当前函数处理，因为 GrayGradinet 使用的是二维灰度矩阵逻辑。需要分析彩色图片时，应先转换成符合当前索引要求的灰度数据。

另外，当前梯度量化使用：

Gmax − Gmin

作为分母。如果输入是一张完全恒定的纯色图，可能出现 Gmax = Gmin，此时原源码没有另外设置零分母保护。

这类输入在使用前需要额外处理。

12. 怎么运行这套程序？

运行dutu

程序会依次执行：

读取D68.gif→ 转为double→ 调用GrayGradinet(img)→ 返回OUT

计算完成后，在 MATLAB 工作区查看：OUT

即可得到 1×15 的纹理特征向量。

如果需要处理自己的灰度图片，只需要把：

imread('D68.gif')

中的文件名替换成目标图像，并保证输入满足当前函数的灰度范围和二维矩阵要求。

当前这条 15 特征主链使用的是 MATLAB 基础矩阵、循环、数学运算和图像读取函数，源码中没有训练模型、外部数据接口或额外模型权重需要先加载。

13. 一句话看懂这个项目

这是一个灰度—梯度共生统计的 15 维图像纹理特征提取程序：主脚本读取工程自带的 640×640 灰度图D68.gif，先把图像转换为双精度矩阵，再通过水平与垂直相邻像素差构造 639×639 梯度图，把梯度线性离散成 32 级，与 256 级灰度共同形成 256×32 共生矩阵，随后计算小梯度优势 T1、大梯度优势 T2、灰度/梯度不均匀性 T3/T4、能量 T5、灰度/梯度平均 T6/T7、灰度/梯度离散程度 T8/T9、相关量 T10、三类熵 T11～T13、惯性 T14 和逆差矩 T15，最终输出 1×15 数值特征向量；它适合灰度纹理图像的离线特征提取、算法教学和后续分类识别前的人工特征构造，但当前主程序本身不包含批量训练、分类器或实时视觉系统。