图像纹理特征提取：GLCM、Tamura与Gabor的MATLAB实现全解析

简介面向MATLAB环境下的图像纹理特征提取这份资源适用于图像处理入门者、计算机视觉研究者及需要快速上手纹理分析算法的工程师。内容围绕统计特征、结构特征和模型化特征展开涉及灰度共生矩阵GLCM、局部二值模式LBP、Gabor滤波器等经典方法并提供基于MATLAB的实现思路与评估手段可辅助理解纹理识别的完整流程。资源包共4个文件含1个Matlab脚本与3张PNG效果示意图压缩包仅375KB轻量便于下载与二次开发。目前已有2264人学习下载适合用于课堂实验、论文复现或算法对比测试。脚本与图示结合可直观观察不同纹理特征提取的效果差异有助于快速搭建自己的纹理分析原型。1. 图像纹理特征提取为什么 GLCM 不是唯一答案以及 MATLAB 里最高效的落地方式做过图像识别的人应该都有过这种体验一张布匹的照片颜色相近、形状没有显著差异但人眼一看就知道是棉还是麻。这种区分感知靠的是纹理。图像纹理特征提取算法要解决的就是把这种“手感”变成数值。在工业缺陷检测、遥感地物分类、医学影像辅助诊断里纹理特征往往是比颜色和形状更稳定的判据。MATLAB 在这条技术路径上有天然优势底层矩阵运算快Image Processing Toolbox 和 Statistics Toolbox 里的函数基本覆盖了主流特征算法的标准实现从读图到特征向量输出代码量比 OpenCV Python 短一半以上。这篇文章从像素统计到滤波响应从 GLCM 到 Tamura 和 Gabor把纹理特征提取的选型逻辑、MATLAB 实现步骤、参数调法和坑一次讲清。2. 灰度共生矩阵GLCM纹理特征的基座算法与 MATLAB 参数含义2.1 GLCM 的数学逻辑二阶统计特征为什么能描述纹理纹理本质上是灰度在空间上的重复排列。一阶统计量均值、方差只描述灰度分布丢失了空间位置关系所以无法区分“棋盘格”和“随机噪点”。灰度共生矩阵Gray-Level Co-occurrence Matrix统计的是满足特定空间关系的像素对灰度组合出现的频率。比如“向右偏移 1 个像素、灰度 5 和灰度 8 相邻出现”这个事件在整幅图中发生了多少次。这个思路之所以能在几十种纹理算法里一直占据核心位置是因为它把“排列”拆成了方向和距离两个可量化的维度。方向通常取 0°、45°、90°、135°距离 d 一般取 1 到几个像素。距离选 1 能捕捉最精细的微观纹理距离选大值能捕捉周期性的粗纹理两者物理含义完全不同。合成孔径雷达图像里的纹理通常用较大距离而医学病理切片的细胞纹理距离 1 就足够。MATLAB 里graycomatrix函数默认将灰度级压缩到 8 级偏移方向默认是 [0 1]即水平向右一个像素。这里注意一个反直觉的细节灰度级压缩得越狠从 256 级到 8 级对灰度微小变化的敏感度越低但对光照渐变的鲁棒性越高。这在实际项目中是一个直接要找平衡的参数。2.2 graycomatrix graycoprops 最小可运行代码% 读取灰度图像缩短变量名便于后续循环 I imread(texture_sample.png); if size(I, 3) 3 I rgb2gray(I); end % 核心步骤1计算GLCMNumLevels设64级Offset四方向 % Symmetric, true 让矩阵对称这样0°和180°的特征值一致计算量减半 offsets [0 1; -1 1; -1 0; -1 -1]; glcm graycomatrix(I, Offset, offsets, NumLevels, 64, Symmetric, true); % 核心步骤2从GLCM提取6个标量统计量 stats graycoprops(glcm, {Contrast, Correlation, Energy, Homogeneity}); % 展示第一个方向的结果 disp(方向 [0 1] 的4个特征值:); disp([stats.Contrast(1), stats.Correlation(1), stats.Energy(1), stats.Homogeneity(1)]); % 把4个方向的特征拼接成一行向量作为分类器输入 % 注意转置graycoprops返回行向量需要按列拼接 feat_vec [stats.Contrast(:), stats.Correlation(:), stats.Energy(:), stats.Homogeneity(:)]; fprintf(特征向量长度: %d\n, length(feat_vec));代码里offsets每行是一个[行偏移 列偏移]注意[-1 1]是左上到右下方向即 135°。这里有个常见隐患graycoprops返回的四个统计量各是一个 1x4 的行向量想拼成特征向量要统一转置否则顺序会乱。如果图像较暗或过曝先做histeq(I)再做 GLCM纹理特征会更稳定这与光照不均条件下的实测结果有关后面排错章节还会提到。NumLevels的值每降低一半GLCM 计算耗时大约减少四分之一但Contrast的区分度会下降。在我的实测里64 级是工业图像的最优值256 级容易把噪声当作纹理8 级会把细节抹平。2.3 四个标准统计量的物理意义与适用场景差异统计量公式方向表征含义适合场景对噪声敏感度Contrast对比度灰度差加权局部强度变化幅度裂纹检测、砂眼检测高椒盐噪声会拉高值Correlation相关性行/列方向线性依赖纹理的方向性强度木材纹理方向判定中Energy能量/ASM灰度对出现概率平方和纹理的均匀性和稠密性布料均匀度检测低Homogeneity逆差矩灰度差倒数加权局部同质性细腻纹理区域分割较低需要特别指出的是Energy高不代表图像“能量高”而是纹理规则、重复性强。纯色图像的 Energy 接近 1而随机噪声图像的 Energy 趋近 0。这个命名误导了很多人实际看分类结果时建议先单独对每个特征做归一化后再进入分类器因为Contrast的量纲和Correlation差两个数量级。3. Tamura 特征与边缘方向直方图GLCM 之外的算法互补与 MATLAB 实现3.1 Tamura 六特征里的粗糙度和方向度为什么要自己写GLCM 有一个结构性短板它描述的是统计分布不能直接回答“这幅图的纹理颗粒多大”这种人类视觉层面的问题。Tamura 等人发现人眼对纹理的感知可以归纳为六个维度粗糙度、对比度、方向度、线像度、规整度、粗略度。其中粗糙度和方向度与 GLCM 特征相关系数低能补充 GLCM 缺失的语义层次。MATLAB 里没有内置 Tamura 函数这是多数教程不写它的原因。但粗糙度计算并不复杂对每个像素计算 2^k × 2^k 邻域的平均灰度差找到使差最大的 k然后对全图取平均。这个操作用imgaussfilt加循环就能完成比想象中简单得多。方向度的计算稍微复杂一些。先要用fspecial(sobel)或imgradient求梯度幅值和方向然后统计方向直方图用直方图的峰值尖锐程度表示方向度。3.2 一个能直接跑的 Tamura 粗糙度与方向度计算函数function [coarseness, directionality] tamura_texture(I) % I: 灰度图像double类型范围[0,1] % 输出: coarseness标量, directionality标量 I im2double(I); [h, w] size(I); % ---- 粗糙度计算 ---- % 对6个尺度分别做均值滤波窗口尺寸为2^kk1到6 max_k 6; A cell(1, max_k); for k 1:max_k win 2^k; % 用same填充保证输出尺寸不变边界用replicate镜像避免边缘伪差 A{k} imfilter(I, ones(win)/win^2, replicate); end % 计算每个像素在各尺度上的邻域平均差 E zeros(h, w); for k 1:max_k-1 diffmap abs(A{k1} - A{k}); % 逐像素取多尺度中差异最大的那个尺度 E max(E, diffmap); end coarseness mean(E(:)); % ---- 方向度计算 ---- % 用Sobel梯度计算梯度方向量化到18个区间每20度一个bin [gx, gy] imgradientxy(I); gmag hypot(gx, gy); gdir atan2(gy, gx) * 180 / pi; gdir(gdir 0) gdir(gdir 0) 180; % 方向角转到[0,180) % 只取梯度幅值大于阈值的像素参与统计减少平坦区域干扰 thr mean(gmag(:)) * 1.2; valid gmag thr; if sum(valid(:)) 10 directionality 0; return; end bins 18; hist_counts histcounts(gdir(valid), bins, BinLimits, [0 180]); hist_counts hist_counts / sum(hist_counts); % 归一化为概率 p hist_counts; % 方向度的经典定义直方图所有波峰的二阶矩之和 directionality 1 - bins * sum(p.^2); % 这个值越大说明方向分布越均匀方向性弱越小说明有明显主方向 end这个函数里值得说明两部分。粗糙度部分用imfilter的replicate填充模式比默认的零填充更适合纹理分析——零填充会在图像四条边产生假的强边缘让粗糙度偏大。方向度部分没有直接输出“最主方向角度”而是输出方向直方图分布的集中度。实际项目里如果需要“纹理偏转了多少度”这个量直接从gdir直方图的峰值位置读取。Tamura 对比度与 GLCM 的 Contrast 不同Tamura 版本是std(I) / (mean(I))^(1/4)对整体亮度不敏感适合待检测目标存在整体亮度漂移的场景比如生产线上不同批次产品的打光亮度不一致。3.3 边缘方向直方图的 MATLAB 实现与参数选择方向直方图Edge Direction HistogramEDH在 MPEG-7 标准里是一个独立的纹理描述子。它把边缘方向量化为 16 个方向区间外加一个“无边缘”区间共 17 个 bin。MATLAB 里用imgradient即可在几行内完成完整实现% 边缘方向直方图 function feat_edh edge_direction_histogram(I, num_bins) if size(I, 3) 3 I rgb2gray(I); end % 高斯预平滑sigma1.5能有效过滤噪声但保留主要边缘 I imgaussfilt(I, 1.5); [gmag, gdir] imgradient(I); % 边缘响应阈值取Otsu阈值的0.5倍 level graythresh(uint8(I)) * 255; edge_mask gmag level * 0.5; % 将方向范围[-180,180)偏移到[0,360)再量化 gdir mod(gdir, 360); bin_width 360 / num_bins; bin_idx floor(gdir / bin_width) 1; bin_idx(~edge_mask) 0; % 非边缘点归零 feat_edh histcounts(bin_idx, 0:num_bins1); feat_edh feat_edh(2:end) / sum(feat_edh(2:end)); % 去掉下标0所在bin end边缘方向直方图对旋转非常敏感。图像旋转 10°直方图的峰整体位移 10°如果直接把这个特征喂给分类器分类精度会断崖式下降。解决办法是每张图像都找到主峰位置然后把直方图按主峰位置对齐这一步在特征提取前先做否则同一类别因为摆放角度不同会被分到不同类。4. 从单张图到数据集批量纹理特征提取的 MATLAB 工程化流程与分类验证4.1 批处理脚本结构文件夹遍历、特征存储与标签映射纹理特征提取进了实际项目就是批处理。假设数据集目录结构是dataset/class_name/*.png每个子文件夹名就是类别标签。用dir递归遍历后按类别编码为整数标签。特征矩阵存成.mat文件比逐行写 CSV 要快一个数量级尤其特征向量长度达到几十到几百时CSV 的文本 I/O 会成为瓶颈。% 批量提取纹理特征 data_dir dataset; class_names dir(data_dir); class_names class_names([class_names.isdir] ~ismember({class_names.name}, {., ..})); all_feats []; all_labels []; for i 1:length(class_names) class_path fullfile(data_dir, class_names(i).name); img_files dir(fullfile(class_path, *.png)); for j 1:length(img_files) img imread(fullfile(class_path, img_files(j).name)); if size(img, 3) 3, img rgb2gray(img); end % 组合三类特征GLCM 16维 Tamura 2维 EDH 16维 LBP直方图 glcm_feat glcm_feature_vector(img); tamura_feat tamura_texture(double(img)); edh_feat edge_direction_histogram(img, 16); lbp_feat extractLBPFeatures(img, Upright, false, NumNeighbors, 8); % 每个特征源的取值范围差异较大先各自归一化到[0,1] glcm_feat rescale(glcm_feat); tamura_feat rescale(tamura_feat); edh_feat rescale(edh_feat); lbp_feat rescale(lbp_feat); feat [glcm_feat, tamura_feat, edh_feat, lbp_feat]; all_feats [all_feats; feat]; all_labels [all_labels; i]; end fprintf(完成类别 %s (%d 张)\n, class_names(i).name, length(img_files)); end save(texture_features.mat, all_feats, all_labels);这段脚本里extractLBPFeatures是 Computer Vision Toolbox 里的函数默认提取带旋转不变性的 LBP 直方图。LBP 与 GLCM 配合使用时LBP 擅长捕捉微观纹理的局部模式差异GLCM 擅长捕捉全局的灰度分布关系组合后对很多自然纹理数据集的分类准确率能提升 3~6 个百分点。特征融合最忌讳的是直接拼在一起就进分类器因为NumNeighbors8时 LBP 特征维度是 59它会把 GLCM 的 16 维特征在距离计算中淹没掉归一化这一步不能省。4.2 特征选择与可视化验证PCA 降维观察可分性特征拼完不要急着训分类器先用 PCA 降到二维或三维看各类别的散点是否可分。这一步在 MATLAB 里只需要三行load(texture_features.mat); [coeff, score, latent] pca(all_feats); % score是降维后的坐标 gscatter(score(:,1), score(:,2), all_labels); % 可视化前两个主成分如果 PCA 散点图中各类别混叠严重说明提取的特征没有捕捉到类别间的核心差异。这时有两条路一是回看 GLCM 的Offset和NumLevels参数二是考虑图像没有做预处理。预处理最常见的是光照归一化——每张图减去自身均值再除以标准差。纹理特征里Correlation对光照不敏感但Contrast和 LBP 直方图对光照非常敏感。如果样本来自不同时间、不同设备采集一定要做光照归一化这不是一个可以不做的优化项。分类器选型上特征维度不高比如 50~200 维时线性 SVM 或随机森林都合适。实测中纹理特征在随机森林上的表现通常比 SVM 好一圈原因是部分纹理特征与类别之间不是线性可分关系比如同一类别内存在不同朝向的纹理单元。% 用随机森林评估特征有效性 rng(42); cv cvpartition(all_labels, HoldOut, 0.3); train_feats all_feats(training(cv), :); train_labels all_labels(training(cv), :); test_feats all_feats(test(cv), :); test_labels all_labels(test(cv), :); mdl TreeBagger(200, train_feats, train_labels, Method, classification, ... MinLeafSize, 1, NumPredictorsToSample, all); pred predict(mdl, test_feats); pred_numeric str2double(pred); acc mean(pred_numeric test_labels); fprintf(随机森林分类准确率: %.1f%%\n, acc * 100);MinLeafSize1在纹理特征上一般能取到最高准确率因为纹理特征维度有限树的深度不需要太大就能过拟合但随机森林的随机性天然抵御过度拟合。NumPredictorsToSample默认是特征数的三分之一在特征维度低于 50 时改为all效果更好。4.3 特征维度与计算时间的权衡表特征类别提取函数特征维度单张 512×512 图像耗时备注GLCM 四方向×4统计量graycomatrix graycoprops16约 0.12sNumLevels64LBP 直方图extractLBPFeatures59约 0.08s旋转不变模式Tamura自定义函数2约 0.35s主要耗时在 6 次 imfilterGabor 滤波响应imgaborfilt244方向×6尺度约 0.40s每个滤波核一次卷积EDH自定义函数16约 0.05s依赖 imgaussfilt如果数据量达到几万张特征提取时间就会相当可观。常见做法是把批处理脚本改成parfor并行循环用parpool开启本地多核并行速度提升接近核心数倍。Gabor 滤波因为每个核独立是并行效率最高的环节。这里要提醒的是graycomatrix在并行循环里会创建多个 OpenMP 线程池与 MATLAB 的parfor叠加反而会变慢需要外层留一个核给主进程。5. 高频坑点与少有人讲的进阶技巧5.1 GLCM 遇到直方图均衡化特征失效的根源很多教程会让图像先做histeq增强对比度再计算 GLCM这在学术示例里有合理性但在工程上要警惕。histeq是把灰度分布拉平它的直接后果是让原本灰度集中的纹理单元变得层次分明这对于人眼观察是增益但对于 GLCM 来说则可能破坏原有纹理的相对灰度差异。比如区别两种材质时原本材质 A 的灰度集中分布在 30~50材质 B 分布在 120~140直方图均衡化会把两者的灰度区间整体拉伸导致原本可以稳定区分的分布被改变。实际经验是如果所有样本的直方图形状接近比如都是暗背景不需要均衡化只有样本间存在整体亮度差异时才做而且事后要在验证集上看特征分布的稳定度。5.2 Gabor 滤波器的参数微调与特征降维组合Gabor 滤波对纹理周期性和方向性的表征能力比 GLCM 更强但需要调好方向数和尺度数。方向取 40°、45°、90°、135°是基础配置如果纹理本身有更强的方向性可以考虑 6 个方向每 30° 一个。尺度的选择更关键imgaborfilt(I, wavelength, orientation)的波长参数决定能捕获的纹理周期。对于布匹类细纹理波长 2~4 个像素对于木材类粗纹理波长 8~16 个像素。最优做法是用gabor函数自动生成滤波组% 生成 4 方向 × 6 尺度的 Gabor 滤波组 g gabor(2:2:12, [0 45 90 135]); [mag, phase] imgaborfilt(I, g); % 每张滤波响应图的均值、方差作为特征 gabor_feat zeros(1, length(g) * 2); for k 1:length(g) gabor_feat(2*k-1) mean(mag(:,:,k), all); gabor_feat(2*k) std(mag(:,:,k), 0, all); end只取均值和方差其实是高度压缩的做法会丢失空间分布信息。更优的做法是把每张滤波响应图分成 4×4 的小块对每个小块取能量均值这样特征维度从 48 涨到 768配合 PCA 降到 100 维左右分类效果有明显提升。这种“滤波分块能量降维”的组合是工程中真正常用的。5.3 验证特征提取有效性的一个快速方法调完参数不知道该往哪个方向优化时有一个极快的诊断方法在训练集内做 KNN 分类K1看留一法准确率。如果留一法准确率不高调 SVM 或随机森林的收益非常有限——说明特征本身区分度不够的问题不在分类器。反过来如果 KNN 留一法率高但测试集准确率低问题出在特征过拟合或训练集过小要回到“样本多样性”而不是继续加特征。这个思路在纹理特征工程里比任何高级分类器调参都先走一步。本文还有配套的精品资源点击获取