先放一段能直接复制的代码。下面这个函数，输入一张图像，不管是灰度图还是 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);

我实测得到的一组典型结果为：

坦白说，不同版本的砖墙图像因为分辨率、光照、砖缝宽窄不同，数值会有浮动。但整体趋势是一致的：方向度高、线像度高、规整度高，这三项让砖墙纹理和木材纹理、布料纹理等明显区分开。

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 这种几十年前的特征，真正落地时拼的不是公式背得有多熟，而是对边界条件、数值稳定性和工程细节的控制。你不需要完全照搬论文的原始公式，抓住物理含义，用工程化的手段把特征算稳、算快，才是项目里真正有价值的部分。如果你用它跑自己的纹理图时发现某些数值不符合直觉，别急着怀疑代码，先把自己的图放大看看纹理基元大小，往往答案就在图像本身。