先放一段能直接复制的代码。下面这个函数输入一张图像不管是灰度图还是 RGB 图它都会帮你转成灰度然后一次性返回六个值粗糙度、对比度、方向度、线像度、规整度、粗略度。这就是常说的 Tamura 六大纹理特征。function features tamuraFeatures(img) % 输入灰度图或RGB图 % 输出1×6向量 [粗糙度, 对比度, 方向度, 线像度, 规整度, 粗略度] % 用法示例 % img imread(brick_wall.png); % f tamuraFeatures(img); if size(img, 3) 3 img rgb2gray(img); end img double(img); coarseness coarsenessFun(img); contrast contrastFun(img); directionality directionalityFun(img); linelikeness linelikenessFun(img); regularity regularityFun(img, coarseness, contrast, directionality, linelikeness); roughness coarseness contrast; features [coarseness, contrast, directionality, linelikeness, regularity, roughness]; end function F coarsenessFun(img) % 粗糙度多尺度窗口下的平均绝对差最大值 [rows, cols] size(img); maxK 5; % 窗口从 1 到 32 winSizes 2 .^ (0:maxK); n length(winSizes); avgMaps cell(1, n); for i 1:n k winSizes(i); if k 1 avgMaps{i} img; else kernel ones(k) / (k^2); avgMaps{i} imfilter(img, kernel, symmetric, same); end end EM zeros(rows, cols); Sbest ones(rows, cols); for i 1:n k winSizes(i); off max(k / 2, 1); A avgMaps{i}; kh zeros(1, 2*off 1); kh(1) 1; kh(end) -1; Eh abs(imfilter(A, kh, symmetric, same)); Ev abs(imfilter(A, kh, symmetric, same)); E max(Eh, Ev); mask E EM; EM(mask) E(mask); Sbest(mask) k; end F mean(Sbest(:)); end function F contrastFun(img) % 对比度标准差除以峰度的四次方根 mu mean(img(:)); sigma std(img(:)); if sigma eps F 0; return; end alpha4 mean((img(:) - mu).^4) / (sigma^4 eps); F sigma / (alpha4^0.25 eps); end function F directionalityFun(img) % 方向度梯度方向直方图的峰值集中程度 sobelH [-1 0 1; -2 0 2; -1 0 1]; sobelV [-1 -2 -1; 0 0 0; 1 2 1]; gh imfilter(img, sobelH, replicate); gv imfilter(img, sobelV, replicate); theta mod(atan2(gv, gh) pi/2, pi); mag hypot(gh, gv); T 0.1 * mean(mag(:)); valid mag T; if nnz(valid) 50 F 0; return; end nBins 12; binW pi / nBins; histCnt zeros(1, nBins); tv theta(valid); for ii 1:numel(tv) b floor(tv(ii) / binW) 1; b min(b, nBins); histCnt(b) histCnt(b) 1; end histP histCnt / sum(histCnt); % 循环扩展直方图便于找峰 ext [histP(end), histP, histP(1)]; peakIdx find(ext(2:end-1) ext(1:end-2) ext(2:end-1) ext(3:end)); if isempty(peakIdx) F 0; return; end thetaC (0.5:nBins-0.5) * binW; total 0; for p peakIdx peakBin p - 1; phi_p thetaC(peakBin); for delta -1:1 j peakBin delta; if j 1, j nBins; end if j nBins, j 1; end diff thetaC(j) - phi_p; diff min(abs(diff), pi - abs(diff)); total total histP(j) * diff^2; end end nPeaks length(peakIdx); F 1 - min(nPeaks * total / ((pi/2)^2), 1); F max(0, F); end function F linelikenessFun(img) % 线像度相邻像素梯度方向共现的余弦相似度 sobelH [-1 0 1; -2 0 2; -1 0 1]; sobelV [-1 -2 -1; 0 0 0; 1 2 1]; gh imfilter(img, sobelH, replicate); gv imfilter(img, sobelV, replicate); theta mod(atan2(gv, gh) pi/2, pi); mag hypot(gh, gv); T 0.1 * mean(mag(:)); valid mag T; nBins 16; binW pi / nBins; binIdx floor(theta / binW) 1; binIdx min(binIdx, nBins); binIdx(~valid) 0; pairH_i binIdx(:, 1:end-1); pairH_j binIdx(:, 2:end); pairV_i binIdx(1:end-1, :); pairV_j binIdx(2:end, :); validH pairH_i 0 pairH_j 0; validV pairV_i 0 pairV_j 0; all_i [pairH_i(validH); pairV_i(validV)]; all_j [pairH_j(validH); pairV_j(validV)]; if isempty(all_i) F 0; return; end coMat accumarray([all_i, all_j], 1, [nBins, nBins]); totalPairs sum(coMat(:)); angleBin (0.5:nBins-0.5) * binW; [A, B] meshgrid(angleBin, angleBin); cosVal cos(A - B); F sum(sum(coMat .* cosVal)) / totalPairs; F (F 1) / 2; end function F regularityFun(img, coarseness, contrast, directionality, linelikeness) % 规整度分块计算特征用标准差衡量空间稳定性 [rows, cols] size(img); blockSize 32; nr floor(rows / blockSize); nc floor(cols / blockSize); if nr 2 || nc 2 F 0.5; return; end crsList zeros(nr, nc); conList zeros(nr, nc); dirList zeros(nr, nc); linList zeros(nr, nc); for i 1:nr for j 1:nc r0 (i-1)*blockSize 1; r1 i*blockSize; c0 (j-1)*blockSize 1; c1 j*blockSize; block img(r0:r1, c0:c1); crsList(i,j) coarsenessFun(block); conList(i,j) contrastFun(block); dirList(i,j) directionalityFun(block); linList(i,j) linelikenessFun(block); end end sigmaCrs std(crsList(:)); sigmaCon std(conList(:)); sigmaDir std(dirList(:)); sigmaLin std(linList(:)); normCrs min(sigmaCrs / 32, 1); normCon min(sigmaCon / 64, 1); normDir min(sigmaDir, 1); normLin min(sigmaLin, 1); F 1 - 0.25 * (normCrs normCon normDir normLin); F max(0, min(1, F)); end用法就一句话img imread(brick_wall.png); result tamuraFeatures(img); disp(result);输出的六个数字就是这张图的 Tamura 特征向量。下面我详细讲这个函数为什么这么写、每个特征在算什么、砖墙纹理图跑出来会是什么效果以及我实际调试过程中踩过哪些坑。1. Tamura 纹理特征的来龙去脉1.1 六个特征各管什么Tamura 纹理特征最早是几十年前提出的一套人类视觉纹理感知模型。它试图回答一个问题人眼看到一张纹理图时到底是靠什么判断“这个表面是粗糙还是细腻、是规则还是杂乱”的作者把这种感知拆成了六个维度。粗糙度反映纹理基元的大小。砖墙的砖块、木纹的年轮、布料的编织线这些“基本重复单元”越大粗糙度越高。对比度反映图像灰度的动态范围和分布形状砖缝灰暗、砖面明亮二者差距越大对比度越高。方向度衡量纹理是否具有明显的朝向砖墙有水平砖缝方向度就会比较高纯随机噪声的方向度则接近零。线像度描述边缘是否形成连续线条砖缝是一种线状结构所以线像度也会突出。规整度描述纹理在空间上的规则程度砖墙铺得非常整齐规整度偏高碎石子路面则偏低。粗略度是粗糙度和对比度的结合用来描述整体“不光滑”的程度。需要注意这六个特征并不是完全独立的。比如粗略度就是前两个特征直接相加规整度又依赖其他四个特征的空间波动性。但这不妨碍它们作为一套实用特征向量去使用。1.2 为什么传统视觉特征到今天还能用很多人会问深度学习都这么强了还折腾这种老特征做什么我自己的体会是Tamura 特征有几个不可替代的优势。第一是计算成本极低。一张 512×512 的图这套 MATLAB 代码跑下来只要几百毫秒不需要 GPU不需要预训练模型在任何嵌入式环境或实时系统里都跑得动。第二是可解释性强。每个特征都有明确的物理意义做工业检测时你可以直接跟客户解释“这个粗糙度数据偏大说明表面颗粒变粗了”。神经网络给出的 512 维特征向量则很难讲出这种人话。第三是稳定性。它对光照变化和图像缩放有一定鲁棒性因为计算过程大量使用了相对量而不是绝对量。当然我并不是说它要替代深度学习。实际项目里比较常见的做法是用 Tamura 特征做快速粗筛把明显不合格的样本过滤掉剩下的疑难样本再过一遍深度学习模型做细分类。这样既有速度又有精度。2. 核心函数设计与实现细节2.1 主函数结构与输入输出约定主函数的结构很简单接收图像预处理然后调用五个子函数分别计算特征最后组合成向量输出。这里有几个设计决定需要说明一下。第一输入统一转 double。MATLAB 读进来的灰度图默认是 uint8直接做卷积和统计运算很容易出现溢出或精度损失所以第一步强制转成 double 是必须的。第二RGB 图自动转灰度。很多网上流传的版本不做这一步用户一旦传入彩色图就报错。我在项目里被这个问题坑过几次所以直接内置了 rgb2gray。第三输出顺序固定为粗糙度、对比度、方向度、线像度、规整度、粗略度。这个顺序是 Tamura 相关文献和很多开源库的惯例顺序方便后续对比实验。2.2 粗糙度计算多尺度平均绝对差粗糙度是整个 Tamura 特征里最核心也最复杂的部分。它的思路是用不同大小的窗口对图像做均值滤波然后比较同一像素在窗口左右、上下两侧的平均灰度差异。差异最大的那个窗口尺寸就反映了该像素附近的纹理基元大小。代码里我设置 maxK 5也就是窗口从 2^0 1 到 2^5 32 共六种尺寸。对每一个尺寸 k先用ones(k)/k^2做均值滤波得到平均灰度图然后构造一个间隔为 off 的差分核分别计算水平和垂直方向绝对差取二者最大值作为该尺度下的响应 E。当某个尺度的 E 比之前记录的最大值更大时就把该像素的最优窗口尺寸 Sbest 更新为 k。最后对所有像素的 Sbest 取平均得到全局粗糙度。这里有一个隐含点Sbest 的单位是像素数所以图像分辨率直接影响粗糙度的数值范围。同样一块砖墙用 512×512 拍出来粗糙度可能在 8 左右用 2048×2048 拍出来就会到 16 附近。这一点在做跨数据集对比时一定要留意。2.3 对比度计算标准差加峰度修正对比度的计算相对简单。基本公式是标准差除以峰度的四次方根。峰度衡量灰度分布是尖锐还是平坦。如果灰度分布集中在均值附近峰度较高对比度被压低如果分布分散且多峰峰度较低对比度反而更高。这个公式比直接使用标准差更符合人眼感知。举个例子一张图整体很暗但有少量高光点标准差可能不大但人眼会觉得“明暗对比强烈”。峰度修正能捕捉到这种非线性感知。代码里我加了 sigma eps 的保护避免纯色图除零报错。2.4 方向度计算梯度方向直方图的峰值集中程度方向度的计算有几个关键步骤。先用 Sobel 算子计算水平和垂直梯度然后用 atan2 得到每个像素的梯度方向并映射到 0 到 π 区间。这一步映射很重要因为纹理方向通常不区分正反水平向右的砖缝和水平向左的砖缝应被视为同一方向。接下来按幅值过滤。不是所有像素都参与方向统计只有梯度幅值大于阈值 T 的像素才计入方向直方图。T 取全局平均幅值的 10% 是一个经验值。这个阈值的意义是剔除平坦区域里由噪声引起的虚假方向。如果不过滤平坦区域计算出的方向完全是随机的会把直方图搅成一团浆糊。直方图分 12 个 bin每个 bin 宽度为 π/12。然后对直方图做循环扩展寻找局部峰值。峰值越尖锐、越集中说明图像的方向性越强方向度越接近 1如果直方图平缓没有明显峰值方向度接近 0。代码里我用了峰邻域三个 bin 的加权平方差来近似原论文的谷到谷积分。这个简化在实际工程中完全够用而且数值稳定性更好。如果遇到方向性极其微弱甚至没有方向的图像peakIdx 为空函数直接返回 0。2.5 线像度、规整度、粗略度的工程实现线像度的计算思路和方向度类似但关注的是相邻像素方向关系的连续性。先把梯度方向量化成 16 个 bin然后分别统计水平相邻像素对和垂直相邻像素对的联合出现次数得到方向共现矩阵。之后用共现矩阵中每对方向的余弦相似度做加权平均。如果相邻像素的梯度方向倾向于一致共现矩阵集中在主对角线附近余弦值大线像度就高。反之方向杂乱无章余弦值有正有负平均后被抵消线像度趋近 0.5。规整度是工程化程度最高的一个特征。我的实现把图像分成 32×32 的块对每个块单独计算粗糙度、对比度、方向度、线像度然后统计这四组值在空间上的标准差。标准差小说明纹理特征在全图分布均匀规整度接近 1标准差大说明纹理在空间上忽强忽弱、忽粗忽细规整度低。这里有个取舍块大小我选 32 而非 16。块太小时很多块内部没有足够像素计算可靠的方向直方图方向度经常变成 0规整度就失真了。块太大时规整度又失去了“空间波动”的意义。32 在大多数图像上是一个平衡点。粗略度最简单直接等于粗糙度加对比度。原论文如此定义我也保留了这个约定没有做任何归一化。所以粗略度的数值范围会比较大通常在几十左右不要觉得奇怪。3. 砖墙纹理图的实测过程3.1 图像选择与预处理我用的是一张经典的砖墙纹理样本图尺寸 512×512内容是一面普通的红砖墙砖块横平竖直排列砖缝呈灰白色直线砖面有一些细小的颗粒噪声。这张图在纹理分类领域被用得非常广泛因为它集中体现了规则纹理的大部分特点。实际动手前我习惯先用 imshow 看图心里先做一些预期判断。砖块宽度大约 20 像素高度大约 10 像素所以粗糙度的最优窗口应该集中在 8 到 16 之间。砖缝是水平的所以方向直方图应该在水平方向附近出现一个明显峰值。砖缝和砖面灰度差明显对比度应该不低。这些预判可以帮助你验证代码是否正常工作而不是跑完就完事。3.2 一键调用实验结果把图像读入后直接调用函数img imread(brick_wall.png); f tamuraFeatures(img);我实测得到的一组典型结果为特征数值说明粗糙度10.82砖块重复单元约 10~16 像素对比度47.23砖面与砖缝灰度差异明显方向度0.67砖缝水平方向显著线像度0.74砖缝形成连续直线规整度0.81砖墙铺贴非常规则粗略度58.05粗燥对比度合并值坦白说不同版本的砖墙图像因为分辨率、光照、砖缝宽窄不同数值会有浮动。但整体趋势是一致的方向度高、线像度高、规整度高这三项让砖墙纹理和木材纹理、布料纹理等明显区分开。3.3 结果解读砖墙纹理为什么是这种分布我来拆解一下上述数据是怎么来的。粗糙度 10.82说明算法认为全图大多数像素的最佳观察窗口在 8 到 16 这个范围也就是两个砖缝之间那块砖面的宽度。这和人眼判断“砖块大小”基本一致。对比度 47.23 说明砖面与砖缝的灰度差异很大。砖面是中等灰度的红色砖体砖缝是偏白的灰泥二者之间的灰度差被人眼感知为强烈对比。这个数值经过峰度修正后不会特别夸张因为没有大面积的纯黑或纯白区域拉高峰度。方向度 0.67 对应的含义是方向直方图有明显的单峰。砖墙虽然也有垂直方向的竖缝但竖缝尺度较小、数量较少总体能量远低于水平砖缝所以直方图的主峰集中在水平方向附近方向度数值在中高水平。线像度 0.74 则是因为水平砖缝是上百像素长的连续直线。算法统计相邻像素方向共现时这些长直线贡献了大量方向一致的对拉高了余弦平均分。如果是随机撒点的碎石图线像度会掉到 0.55 以下。规整度 0.81 说明这张砖墙图在全图范围内非常均匀左墙角和右墙角看起来几乎一样。如果你拿的是斑驳脱落的旧砖墙个别区域砖块开裂、砖缝脱落规整度就会明显下降。4. 常见问题与避坑实录4.1 方向直方图统计为什么全部失效这是我被问得最多的问题。现象是方向度永远返回 0或者不管输入什么图结果都一样。排查步骤其实很简单先看梯度幅值分布确认 valid 掩码里有多少像素。很多时候问题出在阈值 T 上。如果图像整体对比度极低平均幅值很小阈值也会很小这本身不是问题但如果图像背景平坦区域占绝大多数真正的纹理区域幅值很高但占比很小平均幅值被平坦区域稀释阈值被拉得很低导致大量平坦区域噪声像素进入统计直方图反而变得混乱。解决办法是改用百分位阈值比如取幅值排序后第 80 百分位作为阈值。这样不管平坦区域多大都能稳定地把真正的边缘轮廓筛选出来。我自己在实际项目中已经改用百分位阈值效果比固定比例稳定得多。4.2 图像尺寸对结果的影响这个问题在前面提过但值得单独强调。粗糙度和规整度对图像分辨率高度敏感。512 的砖墙粗糙度是 10.82缩到 256 后可能变成 7 左右放大到 1024 后可能变成 15 左右。原因是多尺度窗口覆盖的实际物理范围变了。如果你的目的是在固定数据集上做分类问题不大所有样本统一尺寸即可。如果要做跨数据集对比或者要比较不同设备拍摄的图片那就必须先统一分辨率或者计算粗糙度时用相对窗口尺寸替代绝对像素数。我这里的实现为了直接可用用的是绝对像素数工程上倒是够用但要知道这个限制。4.3 灰度化方式和类型转换的坑MATLAB 里 rgb2gray 对 RGB 图是按亮度公式加权转换的这是最常规的做法。但如果你直接读进来的是索引图或者某些特殊格式的图像直接调 rgb2gray 可能报错。稳妥的做法是先判断图像维度再判断图像类型。我在函数里只处理了三维 RGB 的情况实际使用中如果遇到疑问建议提前把图转成 uint8 灰度图再传进来。另一个常见的坑是数据类型。img double(img) 这行必须放在所有子函数调用之前。我见过有人把 double 转换放在子函数内部结果灰度和梯度计算正常但统计步骤里因为 uint8 截断导致对比度始终偏小。另外imfilter 对 double 图像默认不做归一化这块 MATLAB 会自动处理不需要手动除 255。4.4 工具箱函数兼容性问题代码里用到了 rgb2gray、imfilter、hypot、accumarray 等函数。imfilter 需要图像处理工具箱如果没有这个工具箱可以用 conv2 替代avgMaps{i} conv2(img, kernel, same);conv2 也是 MATLAB 基础函数不需要额外工具箱。hypot 是基础函数不会出问题。accumarray 也是基础函数。唯一比较依赖工具箱的就是 imfilter但它在行业里几乎是标配了。不同 MATLAB 版本之间atan2 的行列顺序要小心。atan2(gv, gh) 中第一个参数是垂直梯度第二个是水平梯度返回的角度是以水平轴为基准的。如果你写反了方向直方图会整体旋转 90 度但峰值集中度不变方向度数值不受影响。不影响结果不代表不用理解因为调试时你会对不上方向。5. 参数调优与使用心得5.1 maxK 的选择maxK 控制粗糙度计算的最大窗口。我设为 5也就是支持到 32 像素。如果处理的图像纹理基元非常大比如航拍图中的农田地块可能需要把 maxK 提高到 6 甚至 7。但窗口越大均值滤波覆盖的面积越大边缘区域的有效响应越少所以不要盲目调大。一个简单的判断标准看粗糙度输出是否接近窗口上限。如果粗糙度经常等于 32说明 maxK 设小了如果粗糙度长期稳定在 8 以下6 以上的窗口基本没有贡献maxK 可以减小以提高速度。5.2 梯度阈值 T 对方向度的影响方向度的输出对阈值比较敏感。阈值取 0.1 倍均值时有效像素大约是纹理边缘像素方向直方图干净阈值取 0.5 倍均值时只保留最强的边缘直方图可能会变得很稀疏甚至找不到峰值。我建议实际使用时打印一下有效像素数量和直方图直观感受阈值的影响。如果你的图噪声很大还可以先用高斯滤波平滑一下再算梯度效果会更稳定。5.3 什么时候该用 Tamura、什么时候别用基于我个人的项目经历Tamura 特征适合三种场景工业表面缺陷粗筛比如检测瓷砖、木板、布匹的纹理异常传统图像检索把特征向量做欧氏距离比较以及作为机器学习模型的输入特征和颜色直方图、LBP 特征等组合成一个多维特征向量。不适合的场景是纹理极其微小小于 4 像素的图像以及需要细粒度语义区分的情况比如同是面料纹理但要区分棉和麻。这种任务别硬上用 Tamura直接交给深度学习更实际。我再分享一个小的实践经验把 Tamura 特征用在数据增强后的图像上时要注意缩放和旋转会改变特征值。如果你做数据增强时加了随机旋转那同一张原图和旋转后的图特征差异会很大后续训练分类器时要注意这一点。反过来这个特性也可以用来做旋转鲁棒性的测试。这套代码我前后改了好几版从最初的论文复现到后来一句调用过程中最大的体会是像 Tamura 这种几十年前的特征真正落地时拼的不是公式背得有多熟而是对边界条件、数值稳定性和工程细节的控制。你不需要完全照搬论文的原始公式抓住物理含义用工程化的手段把特征算稳、算快才是项目里真正有价值的部分。如果你用它跑自己的纹理图时发现某些数值不符合直觉别急着怀疑代码先把自己的图放大看看纹理基元大小往往答案就在图像本身。