MATLAB纹理特征提取全攻略：GLCM、LBP、Gabor等六种方法实战对比

简介本资源是一套面向图像处理初学者与科研人员的MATLAB纹理特征提取工具集聚焦计算机视觉中的纹理分析核心任务适用于遥感图像分类、医学影像识别、工业缺陷检测等实际场景。压缩包共包含6类主流算法的完整可运行代码GLCM灰度共生矩阵、GLDS灰度差分统计、LBP局部二值模式、GMRF广义马尔可夫随机场、FD分形盒维数及Gabor滤波器每种方法均封装为独立函数并附带示例调用脚本便于理解原理、调试参数与对比性能。资源总计255KB以.m源码文件为主结构清晰、注释详尽无冗余依赖开箱即用。目前已有3558人学习下载适合课程设计、毕业课题或项目快速原型开发帮助用户系统掌握多种纹理建模思路并为后续特征融合与分类器构建提供可靠输入。 纹理特征提取这事我早期做图像分类时被坑过不少次。最典型的场景是用 MATLAB 跑完graycomatrix拿到一组统计量就往分类器里塞结果同一纹理换了个尺度准确率就崩了。后来才慢慢明白纹理不是一个单一概念——它既包含像素灰度在局部的统计分布也涉及相邻像素的空间关系、周期结构甚至表面粗糙度的自相似性。所以教材和论文里才会出现 GLCM、GLDS、LBP、GMRF、分形维数、Gabor 这一长串方法。这篇文章把我整理好的 MATLAB 代码和踩坑经验一次性放出来适合正在做图像检索、表面缺陷检测、遥感地物分类、医学图像分析的朋友直接参考。1. 纹理特征为什么需要六种方法先搞清楚每个特征在描述什么1.1 纹理特征的四个层次很多人第一次接触纹理特征时会觉得方法太多、不知道选哪个。我的经验是不要按方法去记而是按它在描述纹理的哪个侧面去理解。纹理特征大致分成四类。第一类是统计特征核心想法是纹理的差异可以通过灰度值的统计规律反映出来。GLCM灰度共生矩阵、GLDS灰度差分统计、LBP局部二值模式都属于这一类。它们的区别在于统计维度不同GLCM 统计的是固定位移下灰度对的联合分布GLDS 统计的是像素差分的概率分布LBP 统计的是局部邻域二值化后的结构模式。第二类是模型特征代表是 GMRF高斯马尔可夫随机场。它不直接统计纹理而是假设纹理图像是从某个随机过程生成的然后估计出这个过程的参数用参数表达纹理。可以理解为给图像做建模 参数估计。第三类是几何/分形特征代表是分形维数FD。它从尺度不变性的角度计算纹理表面的粗糙程度适合描述自然纹理、地形、裂纹这类具有自相似性的对象。第四类是频域特征代表是 Gabor 滤波器组。它把纹理看作不同频率、不同方向的周期信号组合利用一组带通滤波器提取特定频率和方向上的响应能量。1.2 六种方法各自擅长回答什么问题理解分类之后下一步是搞清楚每种方法的适用场景。我把它们总结成了一张表方便你做选型方法描述的是纹理的哪一面典型应用场景GLCM灰度值的空间共现关系通用图像分类、材料识别、医学图像GLDS局部灰度差的分布快速粗糙度评估、早期算法基线LBP局部二值结构模式人脸识别、光照变化明显的场景GMRF纹理的随机过程参数随机纹理分割、自然纹理建模FD表面粗糙度与自相似性地形分析、裂纹检测、材料表面评估Gabor频段与方向的响应强度周期纹理、方向性纹理织物、指纹打个比方描述同一块布料你可以说它纱线密度高统计也可以说它是斜纹织法结构还可以说它织造工艺参数是XX模型或者摸起来粗糙分形/频域。每一种说法都有道理但适用的描述边界完全不同。做特征提取时最忌讳的就是拿着一把锤子把所有钉子都当钉子砸。2. 统一的数据接口让六种方法在同一套代码风格下工作2.1 环境与工具箱准备这篇文章里的代码我默认你用的是 MATLAB R2018b 以上的版本。核心依赖是 Image Processing Toolbox因为graycomatrix、graycoprops、gabor、imgaborfilt都在这套工具箱里。LBP 部分如果想用现成函数extractLBPFeatures需要 Computer Vision Toolbox如果没有也没关系我在后面给出了一份手写实现教学和基本实验都够用。安装环境时最容易忽略的是版本差异。gabor和imgaborfilt是 R2015b 才加入 Image Processing Toolbox 的如果你的 MATLAB 比较老gabor函数会直接报未定义。老版本的话需要自己用fspecial加正弦波构造 Gabor 滤波器我会在 Gabor 那一节补充说明。2.2 统一函数入口避免六个脚本各写各的我早期写特征提取代码时习惯每个方法写一个独立脚本结果到了真正做实验时图像的预处理方式、灰度范围、返回格式都不一样每次都要重新改。后来我改成统一入口也就是把所有方法封装成函数每个函数接收灰度图像返回一个行向量特征。这样组合特征时只需要一个[]拼接。下面是我常用的包装函数框架function featureVector extractTextureFeatures(img, methods) % 统一入口img为灰度或彩色图methods为方法名元胞数组 if nargin 2 methods {glcm, glds, lbp, gmrf, fd, gabor}; end if size(img, 3) 3 img rgb2gray(img); end img im2double(img); % 统一转换到 double范围 [0,1] featureVector []; for i 1:length(methods) switch lower(methods{i}) case glcm featureVector [featureVector, glcmFeatures(img)]; case glds featureVector [featureVector, gldsFeatures(img)]; case lbp featureVector [featureVector, lbpFeatures(img)]; case gmrf featureVector [featureVector, gmrfFeatures(img)]; case fd featureVector [featureVector, fractalDimension(img)]; case gabor featureVector [featureVector, gaborFeatures(img)]; end end end这个入口最关键的一步是im2double(img)。很多人在不同方法里混用 uint8 和 double导致 GLCM 的灰度级范围、Gabor 的响应值、GMRF 的矩阵运算全都不在一个尺度上最后特征拼接出来自然不稳定。2.3 预处理里最容易埋雷的三个细节先说灰度化。如果输入是彩色图直接rgb2gray是按 ITU-R BT.601 公式做的加权平均这在多数场景没问题。但如果你的应用本身对颜色敏感比如皮肤检测、地物分类里的植被判别灰度化会把颜色信息丢掉那就不适合用纹理特征或者应该把颜色特征和纹理特征一起用。第二点是尺寸归一化。纹理特征对分辨率很敏感。同一块织物拍出来是 200×200 还是一张 1000×1000 的图GLCM 的对比度、Gabor 的响应均值和 LBP 直方图都会有差异。如果你的数据来自不同设备建议先统一到固定尺寸比如imresize(img, [256, 256])再做特征提取。第三点是关于灰度范围的统一。graycomatrix默认的GrayLimits是图像本身的最小值和最大值这意味着不同图像的灰度映射不一致。更稳妥的做法是把图像先归一化到 [0,1]然后显式指定GrayLimits, [0,1]。这样特征才在不同图像之间可比。3. GLCM灰度共生矩阵用的时候最顺手坑也最多3.1 核心原理一句话版本GLCM 要回答的问题是在图像里灰度值 i 和灰度值 j 在指定位移关系下共同出现的频次是多少。举个例子对于水平方向位移为 1 的像素对统计整幅图中左边是灰度 3、右边是灰度 5出现了多少次就得到共生矩阵中的一个元素。实际使用中一般不会只取一个方向而是取 0°、45°、90°、135° 四个方向每个方向算一个 GLCM再基于每个 GLCM 提取对比度、相关性、能量、同质性等统计量。3.2 MATLAB 代码实现function feat glcmFeatures(img, graylevels, dist) % img: double 灰度图范围 [0,1] % graylevels: 灰度级数默认 16 % dist: 步长默认 1 if nargin 2, graylevels 16; end if nargin 3, dist 1; end offsets [0 1; -1 1; -1 0; -1 -1] * dist; glcms graycomatrix(img, ... Offset, offsets, ... NumLevels, graylevels, ... GrayLimits, [0, 1], ... Symmetric, true); stats graycoprops(glcms, ... {Contrast, Correlation, Energy, Homogeneity}); if size(stats.Contrast, 1) 1 % 按方向拼接每个方向 4 个统计量 feat [stats.Contrast(:), stats.Correlation(:), ... stats.Energy(:), stats.Homogeneity(:)]; else % 自动退化为单方向 feat [stats.Contrast, stats.Correlation, ... stats.Energy, stats.Homogeneity]; end end几点说明。offsets里每一行是[row_offset, col_offset]第一行[0 1]表示水平方向第二行[-1 1]表示左上到右下对角线方向依次是垂直方向和另一条对角线。注意 MATLAB 的行列坐标是行号对应 y 方向、列号对应 x 方向所以[0 1]其实是水平位移不是垂直位移这个很容易记反。Symmetric, true 表示生成对称矩阵。简单说它会把每个方向的像素对再反过来统计一遍相当于把[i,j]和[j,i]都计入。这样做的好处是统计量对方向符号不敏感而且graycoprops提取的特征值更稳定。3.3 灰度级数和灰度范围怎么定灰度级数NumLevels是最容易被忽视的参数。默认是 8也就是把灰度压到 8 级。对于纹理差异细微的图像8 级往往糊掉了但改成 256 级也不是好事因为矩阵会非常稀疏统计量失去意义。我通常用 16 或 32先做实验对比再决定。GrayLimits也很关键。如果你不指定函数会取当前图像的实际灰度范围。但注意graycomatrix传入 double 图像时默认GrayLimits是[min(img(:)), max(img(:))]如果传入 uint8 图像默认是[0, 255]。这就是为什么我在统一入口里先im2double再在这里显式传[0,1]。踩坑提示如果你的图像是纯色或者近似纯色graycoprops的 Correlation 会返回 NaN因为分母是灰度方差方差为 0 时除数为 0。处理时可以做一次isnan检查或者强制把 Correlation 置 0。我在实际项目里处理过很多次了这个问题非常普遍。3.4 方向拼接还是方向平均GLCM 特征的方向处理有两种做法。一种是把四个方向的 16 个特征全部保留适合纹理方向性较强的场景比如织物瑕疵检测另一种是把四个方向的统计量取平均得到一个 4 维向量适合做旋转不敏感的分类。再补充一个技巧如果希望特征对旋转更鲁棒可以不取平均而是取四个方向的最小值和最大值组成 8 维特征。这样既保留了方向差异的强度信息又比直接拼接更紧凑。4. GLDS与LBP轻量级统计特征的两种典型姿势4.1 GLDS统计相邻像素的灰度差GLDS 是灰度差分统计的简称思路比 GLCM 更直接对于某个位移向量 d计算像素对之间的灰度差绝对值然后统计这些差值的概率分布。常用的统计量包括均值、对比度、角二阶矩能量和熵。function feat gldsFeatures(img, dist, numBins) if nargin 2, dist 1; end if nargin 3, numBins 256; end offsets [0 1; -1 1; -1 0; -1 -1] * dist; feat []; for k 1:size(offsets, 1) dy offsets(k, 1); dx offsets(k, 2); d shiftDiff(img, dy, dx); p histcounts(d(:), -0.5:1:(numBins - 0.5)) / numel(d); idx 0:(numBins - 1); meanV sum(idx .* p); contrast sum(idx.^2 .* p); asm sum(p.^2); ent -sum(p .* log(p eps)); feat [feat, meanV, contrast, asm, ent]; end end function d shiftDiff(img, dy, dx) % 只取两个像素都不越界的有效区域 [H, W] size(img); r1 max(1, 1 - dy):min(H, H - dy); c1 max(1, 1 - dx):min(W, W - dx); r2 r1 dy; c2 c1 dx; d abs(img(r1, c1) - img(r2, c2)); end这个shiftDiff子函数的技术要点是不采用circshift那种循环移位因为循环移位会把图像另一侧的像素拉过来导致边缘处产生虚假的大差值也不采用简单丢弃边缘的方式而是用索引裁剪出有效区域。这样既避免了边界污染又最大化利用了图像面积。GLDS 的缺点很明确它只统计差值大小不关心差值的空间位置所以对纹理的周期性、方向性描述能力弱。但它速度快、实现短、特征维度低特别适合作为全特征集里的基线特征。4.2 LBP从局部结构编码到直方图LBP 的基本操作是取一个中心像素比较它和周围 8 个邻域像素的灰度大小大于等于记 1小于记 0然后按固定顺序组成一个 8 位二进制数转成十进制就是该像素的 LBP 码。最后统计整幅图像的 LBP 码直方图就得到特征向量。LBP 最大的优点是灰度单调变化不敏感。比如整幅图像变暗了只要相对大小关系不变LBP 码就不变。所以它在人脸识别、光照变化大的室外场景里非常受欢迎。4.3 手写 LBP 与工具箱版本先放手写版便于理解原理也防止没有 Computer Vision Toolbox 时没法用。function feat lbpFeaturesManual(img) img im2double(img); [H, W] size(img); codeImg zeros(H - 2, W - 2); for i 2:(H - 1) for j 2:(W - 1) center img(i, j); nb [img(i-1, j-1), img(i-1, j), img(i-1, j1), ... img(i, j1), img(i1, j1), img(i1, j), ... img(i1, j-1), img(i, j-1)]; codeBits nb center; codeImg(i-1, j-1) sum(codeBits .* (2.^(7:-1:0))); end end feat histcounts(codeImg(:), -0.5:1:255.5) / numel(codeImg); end手写版本的问题在于速度。针对一张 512×512 的图双重循环在 MATLAB 里要跑好几秒只适合教学和小图验证。生产场景建议直接用工具箱函数feat extractLBPFeatures(I, ... NumNeighbors, 8, ... Radius, 1, ... Upright, false);这里Upright, false 表示使用旋转不变 LBP。注意这个参数名的理解容易反过来Upright为 true 时表示使用非旋转不变的直立 LBP为 false 时才做旋转不变处理。4.4 等价模式和旋转不变LBP 进阶必须知道的事原始 LBP 有 256 种模式但在实际纹理中大部分模式很少出现。统计发现有一类等价模式Uniform Pattern占了绝大多数它指二进制码中 0/1 跳变次数不超过 2 的模式。例如 10000000 跳变 2 次属于等价模式而 10101010 跳变 8 次不是。把等价模式单独编号其余模式合并成一类直方图就从 256 维降到了 59 维8 邻域情况。这一方面降低了特征维度、减少了噪声另一方面也保留了绝大多数纹理信息。如果再做旋转不变模式数可以进一步降到 10 维。我做实验时的习惯是如果数据量够大、分类器够强就用 256 维完整 LBP如果特征要和其他方法拼接或者样本量不大就用extractLBPFeatures默认的 uniform 版本。从结果看59 维版本在大多数场景下的分类精度并不比 256 维差训练还更稳定。5. GMRF与分形维数从模型和几何两个角度看纹理5.1 GMRF在做什么把纹理当作随机场GMRF高斯马尔可夫随机场的想法很不一样。它假设图像中每个像素的灰度值可以由它周围邻域像素的灰度值加权求和再加上一个高斯噪声来表示。写成公式就是I(s) Σ_{r∈N} θ_r (I(sr) I(s-r)) e(s)其中 N 是选定的邻域方向集合θ_r 是模型参数e(s) 是零均值高斯噪声。纹理不同邻域像素对中心像素的预测能力就不同θ 参数自然不同。因此一组 θ 参数和一个噪声方差就可以作为该纹理的特征向量。这里的关键是对称化。GMRF 一般要求 θ_r θ_{-r}也就是相对方向的系数相等这样模型才满足平稳性。所以代码里把正方向和反方向的邻域像素相加作为回归变量而不是分别设两个系数。5.2 最小二乘估计 GMRF 参数function feat gmrfFeatures(img, order) if nargin 2, order 2; end img im2double(img); [H, W] size(img); % 取内部区域避免边界 y 2:(H - 1); x 2:(W - 1); center img(y, x); Y center(:); % 二阶邻域的四组对称方向水平、垂直、两条对角线 nbrs {[1 0], [0 1], [1 1], [1 -1]}; nParams length(nbrs); X zeros(numel(Y), nParams); for k 1:nParams dy nbrs{k}(1); dx nbrs{k}(2); up img(y dy, x dx); down img(y - dy, x - dx); X(:, k) up(:) down(:); end theta X \ Y; resid Y - X * theta; sigma2 var(resid); feat [theta(:), sigma2]; end这段代码最需要注意的是索引范围。y 2:(H-1)、x 2:(W-1)是中心像素的有效区域这样ydy最大为 H最小为 1不会越界。MATLAB 的\运算符做最小二乘时比自己写pinv(X*X)*X*Y更稳定我建议直接用\。如果你的图像尺寸很小导致样本数少于参数个数X \ Y会给出警告甚至错误。这时候要么减小邻域阶数要么在原图上做块重叠采样扩大样本量。GMRF 特征的实际效果和应用场景密切相关。它估计的是空间依赖强度对自然纹理比如草地、云层、石材表面描述力不错但计算量比 GLCM 大特征解释性也差一些。如果项目里的纹理偏向规则、周期性强GMRF 不一定比 GLCM 有优势。5.3 分形维数的盒子计数法实现分形维数用在图像纹理上主要是描述灰度表面的粗糙程度。思想是把图像的灰度值看成三维空间里的一个曲面x、y 是像素坐标z 是灰度值然后用不同尺寸的立方体盒子去覆盖这个曲面统计需要的盒子数随尺度的变化。粗糙纹理占据的体积更大盒子数随尺度缩小的增长速度也更快分形维数就更高。function fd fractalDimension(img, scales) img im2double(img); [H, W] size(img); if nargin 2 scales [2, 4, 8, 16, 32]; end logInvS []; logN []; for s scales if s min(H, W) continue; end hCount floor(H / s); wCount floor(W / s); N 0; for i 1:hCount for j 1:wCount block img((i-1)*s 1 : i*s, (j-1)*s 1 : j*s); zRange max(block(:)) - min(block(:)); nBoxes ceil(zRange / (1 / s)) 1; N N nBoxes; end end logInvS(end1) log(1 / s); logN(end1) log(N); end if length(logInvS) 2 fd 1; return; end p polyfit(logInvS, logN, 1); fd p(1); end这个实现里nBoxes ceil(zRange / (1 / s)) 1的计算逻辑是因为灰度已经归一化到 [0,1]而盒子的空间边长是 s所以对应到 z 轴上的盒子高度取 1/s。每个块内灰度跨越几个盒子高度就用向上取整。1是为了保证最小也覆盖一个盒子避免灰度差为 0 时出现 0 个盒子的情况。有几个容易踩的坑。第一盒子的空间尺度和 z 轴尺度必须一致。如果图像灰度不归一化z 轴的范围可能是 0 到 255空间轴是像素数两者量纲不一致分形维数就会失真。第二分形维数依赖于尺度范围的选取一般从 2 开始翻倍取到图像短边的一半左右。第三polyfit的拟合要至少两个点最好四个点以上否则算出来的斜率毫无意义。5.4 分形维数的局限分形维数只输出一个标量信息量相对有限。它对纯色块图像会退化到 1 附近对高噪图像会虚高因为噪声本身在极小尺度上增大了粗糙度。我一般把它当辅助特征和 GLCM、Gabor 等特征拼接而不是单独作为分类依据。另外要提醒一句分形维数没有唯一的标准算法盒子计数法只是最经典的一种。除了它还有毯子法、功率谱法、差分计盒法等等。不同算法算出的数值不能直接横向比较同一算法下对比才有效。6. Gabor滤波器组频域视角的纹理特征构建6.1 Gabor 滤波器为什么对纹理有效Gabor 滤波器本质上是一个被高斯窗截断的正弦波。因为高斯窗的存在它只在局部区域内响应因为正弦波的存在它只对特定频率和特定方向的成分敏感。这种局部 频带 方向的组合让 Gabor 滤波器成为分析周期性纹理、方向性纹理的利器。一个更直观的理解纹理里经常有重复的图案比如织物的经纬线、指纹的脊线、木材的年轮。这些重复结构在频域里就是某个频率附近的能量峰。Gabor 滤波器就像一把尺子只测量特定频率和方向的能量强度。换不同的尺子就能拼出一幅完整的纹理频谱画像。6.2 使用 gabor 和 imgaborfilt 构建特征MATLAB 的 Image Processing Toolbox 从 R2015b 开始提供了gabor和imgaborfilt用法非常直接。function feat gaborFeatures(img, wavelengths, orientations) if nargin 2 wavelengths [3, 5, 8, 13, 21]; end if nargin 3 orientations [0, 45, 90, 135]; end img im2double(img); g gabor(wavelengths, orientations); mag imgaborfilt(img, g); feat []; for k 1:numel(g) m abs(mag(:,:,k)); feat [feat, mean(m(:)), std(m(:))]; end endimgaborfilt返回的mag是一个三维数组第三维长度等于滤波器个数。每个滤波器对应一幅响应图。特征提取时我通常取响应幅值的均值和标准差。均值反映该频率方向上的总体能量标准差反映该频带响应的波动幅度。有些论文还会加入能量、熵、局部极值密度等但对于大多数分类任务均值加标准差已经够用。6.3 参数调节的经验wavelength波长是最影响结果的参数。它和空间频率成反比波长越小滤波器越关注细小纹理。我给的[3, 5, 8, 13, 21]是一个比较通用的范围。如果纹理很细比如细砂纸表面可以往下调到 2如果纹理很粗比如大块岩石要往上调到 30 以上。orientation方向的选取要看你的任务。四方向 0、45、90、135 度覆盖了主方向足够大多数分类任务。如果是强方向性纹理的精确测量可以考虑每 30 度一个方向共 6 个方向。方向越多特征维度越高训练样本不够时容易过拟合。还有两个容易踩的坑。第一imgaborfilt的边缘填充方式由PadMode控制默认是symmetric镜像填充。如果你发现边缘响应异常高可以改成replicate或者none具体看你的纹理是否靠近边界。第二滤波器数量多、输入图像大时内存消耗会比较夸张。比如 5 个波长 × 4 个方向就是 20 个滤波器每张 1000×1000 的图要生成 20 幅 double 类型的响应图约 160MB。如果批量跑几千张建议分块处理或改用单精度存储。7. 六种特征放到同一张表里怎么选、怎么组合7.1 维度、速度和旋转敏感度对照把六个方法放在一起对比才能知道该往哪个方向投入时间。我根据自己的实现整理了一个典型对照表假设输入是 256×256 灰度图方法典型特征维度计算耗时参考旋转不变性GLCM164方向×4统计量快依赖方向处理方式GLDS164方向×4统计量很快依赖方向处理方式LBP256 或 59等价模式中本文还有配套的精品资源点击获取