做图像处理的朋友应该都有过这种体验图像二值化之后画面里一堆白色目标想对它们做分类识别第一反应是把像素全存下来那数据量太大了而且对位置、旋转都很敏感。我处理这类问题有个固定套路先把目标边界提出来用边界跟踪得到闭合轮廓再用链代码把轮廓编码成一串数字最后用主成分描述概括轮廓的方向和主轴长度。这套组合拳在零件分拣、字符识别、叶片形态测量里都特别实用而且Matlab代码量不大很适合作为图像表示的基础工具箱。这篇文章我会把三个环节拆开讲边界跟踪怎么从二值图里提取轮廓坐标链代码怎么把坐标序列变成方向序列主成分描述怎么用几个数字概括整个形状。每部分都会给出可直接复现的Matlab代码最后再聊聊实际调试中容易翻车的细节。1. 边界跟踪到底解决了什么问题从像素矩阵到闭合轮廓点集1.1 图像表示的本质是降维一张二值图像本身就是个矩阵但直接把像素矩阵当特征用问题很多。一是维度太高一个200x200的目标就有四万个像素二是对平移、旋转、缩放都太敏感同一个目标稍微转个角度像素矩阵就完全变了。所以需要找到一种更紧凑、更稳定的表示方式。边界跟踪做的就是这一步把目标的外轮廓提取成一个闭合的点集每个点只有两个坐标。这样数据量瞬间从几万像素变成几百个点并且保留了形状最关键的几何信息。轮廓内部的纹理、灰度变化如果不需要可以被完全丢掉。1.2 边界表示和区域表示怎么选图像表示大致分两类基于区域和基于边界。区域表示关心的是目标覆盖了哪些像素常用面积、质心、区域矩、连通域标记来描述边界表示关心的是目标边缘长什么样常用周长、边界点坐标、链代码、傅里叶描述子来描述。两者的应用场景很不一样。如果你做的是目标存在性判断比如检测画面里有没有瑕疵区域表示更直接如果你做的是形状区分比如圆形和方形、不同型号的零件、不同品种的叶片边界表示往往更有效因为人眼识别形状时依赖的也是轮廓。边界跟踪得到的轮廓点集是后续链代码和主成分描述的共同输入。所以这一步的质量直接决定后面所有特征的好坏。轮廓漏了一个关键凹点或者多绕了一圈噪声链代码和主成分特征都会跟着变。2. 手写边界跟踪代码Moore邻域搜索的完整实现2.1 Moore邻域和搜索方向边界跟踪最常见的思路是Moore邻域搜索也就是当前像素外围的8个邻居。从起点出发沿边界走每一步都检查当前点周围的像素找到下一个边界点。关键在搜索顺序。假设当前边界点是P前一个边界点是C需要从C的邻域位置开始按固定顺序检查P的8个邻居找到第一个前景像素就把它当作下一个边界点。这样做能保证始终贴着目标走不会跳到内部去。方向顺序我习惯用逆时针从正右方开始右、右上、上、左上、左、左下、下、右下对应的方向码就是0到7。这个顺序如果写错跟踪会直接乱掉输出一堆乱七八糟的点。2.2 一个适合理解的简化Matlab实现下面这段代码是一个适合理解原理的简化版本只对单连通、无空洞、且不与图像边界粘连的简单形状有效。实际项目里我一般直接调bwboundaries但手写一遍能让你彻底搞懂搜索原则。function B traceBoundary(bw) % 简化版 Moore 边界跟踪 % bw: 二值图像目标为1背景为0 % B : N*2矩阵[行, 列]闭合轮廓点集不含重复起点 bw bw 0; [H, W] size(bw); % 找到第一个非零像素作为起点 idx find(bw, 1, first); if isempty(idx) B []; return; end [r0, c0] ind2sub([H, W], idx); % 方向定义从正右开始逆时针索引0~7 dr [0 -1 -1 -1 0 1 1 1]; dc [1 1 0 -1 -1 -1 0 1]; B [r0 c0]; cur [r0 c0]; % 起点是最左上角的前景像素左侧、上方一定是背景 % 初始搜索方向从正右开始 dirIdx 0; maxIter H * W * 4; % 防呆上限避免死循环 iter 0; while iter maxIter iter iter 1; found false; for step 1:8 d mod(dirIdx step - 1, 8) 1; % 从当前方向开始顺时针探测 nr cur(1) dr(d); nc cur(2) dc(d); if nr 1 nr H nc 1 nc W bw(nr, nc) B [B; nr nc]; %#okAGROW % 搜索方向回退一格避免漏掉内部转折点 dirIdx mod(d - 2, 8); cur [nr nc]; found true; break; end end if ~found break; end % 回到起点且已经走了3步以上认为闭合 if size(B, 1) 3 isequal(cur, B(1, :)) break; end end end这段代码的核心就一个循环每次在8个方向里找第一个前景点找到就加入点集并更新方向找不到就说明边界断裂。dirIdx mod(d - 2, 8)这行如果你直接抄别改掉它决定了跟踪是贴着边界走而不是朝远处乱跑。2.3 为什么我实际项目中会用内置函数手写版本的价值在于教学但它对复杂形状很不友好。目标如果和图像边界粘连起点选择就不对了目标如果有洞外边界和内边界要分别处理而且这个简化版没有做逆时针/顺时针一致性判断不同方向的跟踪结果可能会反向。所以我平时更推荐用bwboundaries它对8连通、4连通、是否有孔洞都能正确处理返回的边界点集顺序稳定。手写代码的意义是让你在调试链码和PCA时万一结果不对能自己推哪儿出问题而不是对着内置函数黑盒发呆。3. 链代码不只是数字串编码、起点归一化和旋转不变量3.1 8方向链码的编码规则拿到闭合轮廓点集之后下一步是把坐标序列变成方向序列这就是链代码。8方向链码用0到7这8个数字表示每一步的移动方向编码规则和上一节的方向定义完全一致0表示右1表示右上2表示上3表示左上4表示左5表示左下6表示下7表示右下。给你一个直观例子一个3x3的方形轮廓从左上角开始顺时针走方向序列大概是2、2、4、4、6、6、0、0这样的循环。注意这里的4不是数值大小而是“左移”这个方向别搞混。链代码的优点很明显它直接把二维坐标转换成一维序列非常适合作为匹配和压缩的输入。但缺点同样突出对起点敏感对旋转敏感。同一个轮廓从不同的边界点开始编码得到的数字串完全不同旋转一个角度编码序列也会大变。3.2 起点归一化和差分链码起点归一化的目标是无论从哪个边界点开始编码都能得到同一条数字串。方法很简单把所有可能的循环移位都生成取字典序最小的那个序列作为标准形态。差分链码的目标则是消除旋转影响。它的思路是不记录每一步的绝对方向而是记录相邻两步之间的方向变化量。比如上一步方向是0下一步方向是1差分就是1上一步是7下一步是0差分是mod(0 - 7, 8) 1表示只转了45度。差分链码对旋转不敏感是因为整体旋转会让所有绝对方向码同时偏移但相邻方向之间的相对差值保持不变。这一点在形状匹配里非常关键也是很多基于链码的识别算法能稳定工作的原因。3.3 链码生成的Matlab实现function code chainCode(B, directionType) % 输入: B 为闭合轮廓点集N*2矩阵[行, 列] % 输入: directionType eight 或 four % 输出: code 为链码序列 if nargin 2 directionType eight; end % 确保闭合如果首尾不同补上起点 if ~isequal(B(1, :), B(end, :)) B [B; B(1, :)]; end if strcmpi(directionType, eight) dr [0 -1 -1 -1 0 1 1 1]; dc [1 1 0 -1 -1 -1 0 1]; else dr [0 -1 0 1]; dc [1 0 -1 0]; end n size(B, 1) - 1; code zeros(n, 1); for i 1:n drr B(i1, 1) - B(i, 1); dcc B(i1, 2) - B(i, 2); pos find(dr drr dc dcc, 1); if isempty(pos) error(轮廓点集存在非相邻跳变请检查边界跟踪结果); end code(i) pos - 1; end end function codeN normalizeStart(code) % 起点归一化取所有循环移位中字典序最小的序列 n numel(code); code2 [code(:); code(:)]; candidates zeros(n, n); for i 1:n candidates(i, :) code2(i:in-1); end [~, idx] sortrows(candidates); codeN candidates(idx(1), :); end function diffCode diffChainCode(code, modBase) % 差分链码相邻方向变化量 % 返回循环差分长度与原链码相同 if nargin 2 modBase 8; end n numel(code); diffCode mod(diff([code(:); code(1)]), modBase); end实际使用里要注意bwboundaries返回的坐标顺序是按行、列排列的和图像坐标的x、y差一个转置关系。如果你直接把系数矩阵丢进链码函数先确认一下方向对应关系最快捷的验证方式是画一个已知形状打印前几个链码值人工检查对不对。3.4 链码的小型实战验证我用一个任意形状做过测试先用边界跟踪得到轮廓再生成8方向链码然后把轮廓旋转30度、平移若干像素重新算链码。绝对链码确实变了但差分链码几乎完全一致。差别只出现在图像离散化带来的锯齿边缘上不影响整体匹配。这个性质在字符识别里很有用。手写数字A和B即使位置不同、倾角不同差分链码的统计特征依然能区分开不同类别。配合简单的归一化处理就能作为分类器的输入特征。4. 主成分描述从二维点集到方向与主轴长4.1 边界点的协方差矩阵在讲什么链代码描述的是轮廓的局部走向主成分描述则负责整体形态。做法是把轮廓点集看作平面上的一个二维点云对这个点云做主成分分析找出方差最大的方向以及与之垂直的次大方向。想象一堆轮廓点铺在平面上PCA要回答的问题就是这些点整体朝哪个方向拉得最长垂直于这个方向又有多宽这个信息对判断细长物体的朝向特别有用比如螺丝刀、叶片、条形零件。数学上很简单先把坐标中心化减去均值构造2x2的协方差矩阵然后求特征值和特征向量。最大特征值对应的特征向量就是主方向特征值的平方根近似描述点集沿该方向的标准差。4.2 协方差矩阵、特征向量和椭圆近似这套东西可以类比成“给形状套一个椭圆”中心是点集的均值主轴方向是第一个特征向量长短轴长度是特征值的平方根的两倍。椭圆虽然丢掉了细节但抓住了形状最基本的姿态和比例。如果轮廓接近圆两个特征值会非常接近主方向就不稳定随便一点噪声都能让角度跳变。遇到这种情况别慌这不是算法错了而是形状本身就没有明确主方向。在做朝向估计时要提前判断特征值比例只有比例明显大于1时主方向才有意义。4.3 主成分描述的Matlab实现function [theta, axesLen, meanPt] pcaShape(B) % 输入: B 为轮廓点集N*2矩阵[行, 列] % 输出: theta 主方向角度单位度 % axesLen 长短轴长度[长轴; 短轴] % meanPt 轮廓中心坐标[行均值, 列均值] pts double(B); meanPt mean(pts, 1); X pts - meanPt; C (X * X) / (size(X, 1) - 1); % eig返回的特征值可能不按大小排序需要重排 [V, D] eig(C); [vals, idx] sort(diag(D), descend); V V(:, idx); theta atan2(V(2, 1), V(1, 1)) * 180 / pi; axesLen 2 * sqrt(max(vals, 0)); end这里有个容易被忽略的细节eig求出的特征向量方向不唯一可能是正向也可能是反向所以theta可能会差180度。这在实际使用中通常不影响因为方向角本身是模180度的衡量朝向不需要区分头尾。但如果要旋转校正图像需要额外判断方向向量符号让物体统一朝上。4.4 主成分描述能做什么主成分描述最常见的用途有三个。第一个是朝向估计。将主方向角度作为旋转校正的基准可以把不同角度的物体统一转正方便后续模板匹配。第二个是形态分类。长短轴比例就是最简单的形状特征细长物体接近0圆形物体接近1。在零件分拣场景里这个比例能直接区分铁钉和垫圈。第三个是降维重建。轮廓点投影到主轴上理论上前两个主成分就能描述点云的主体结构。如果边界本身接近椭圆用20个点重建和用500个点重建的视觉区别很小这对节省存储和加速计算很有价值。5. 一套能直接跑的Matlab实测流程及输出解读5.1 从图像到特征的三步主流程前面三部分讲的是原理和单步函数现在把它们串起来跑一遍完整流程。假设你有一张包含单个目标的二值图像预处理之后你想提取边界链码和PCA特征主脚本可以这样写% 读图并预处理 rgb imread(shape.png); bw imbinarize(im2gray(rgb)); bw imclose(bw, strel(disk, 3)); % 闭合小孔平滑锯齿 % 第一步边界跟踪 boundaries bwboundaries(bw, 8, noholes); if isempty(boundaries) error(没有找到目标); end % 取面积最大的目标避免取到噪声 B boundaries{1}; for k 2:numel(boundaries) if size(boundaries{k}, 1) size(B, 1) B boundaries{k}; end end % 第二步链代码 code chainCode(B, eight); codeN normalizeStart(code); diffCode diffChainCode(codeN, 8); % 第三步主成分描述 [theta, axesLen, meanPt] pcaShape(B); % 打印特征 fprintf(主方向角度: %.2f°\n, theta); fprintf(长轴长度: %.2f, 短轴长度: %.2f\n, axesLen(1), axesLen(2)); fprintf(长短轴比: %.3f\n, axesLen(2) / axesLen(1)); fprintf(链码长度: %d\n, numel(codeN));这个脚本里的预处理顺序值得说一句先转灰度再二值化再闭运算。闭运算的半径不要太大否则会把真正的轮廓细节抹掉。我一般先用半径为1或2的小核效果不够再慢慢加。5.2 输出结果怎么看我用一个旋转了的矩形轮廓做过测试输出的大致形态是特征数值说明主方向角度约45°对应矩形长边方向长轴长度约120和矩形长边接近短轴长度约60和矩形短边接近长短轴比约0.5明显小于1说明呈条状链码起点最小字典序序列与起点选择无关看到这个结果第一件事不是高兴而是检查数值是否合理。长轴长度应该和目标在图像中的实际尺寸差不远长短轴比应该符合肉眼对形状的判断。如果角度和预期差太多优先怀疑坐标行列顺序这是Matlab图像处理最常见的坑。5.3 旋转前后的归一化验证我多次测试后发现旋转前后差分链码的相似度很高但如果没有先做起点归一化差分链码和绝对链码都会乱。正确的顺序是绝对链码 - 起点归一化 - 差分链码。不要先差分再归一化那样归一化的对象已经不是原始方向序列语义会发生偏移。具体做法可以这样验证将同一目标旋转到30度和60度分别跑一遍特征提取。旋转后的主方向角度应该差30度长短轴比应该几乎不变差分链码的差异控制在锯齿噪声允许的范围内。如果这些条件都满足说明你的整套流程是可靠的。6. 真实调试中容易翻车的几个细节6.1 边界跟踪的死循环和漏点问题手写Moore边界跟踪最容易出现的情况就是死循环或者在边界上来回重复。原因基本都是方向回退逻辑写错。我建议在调试时加一个迭代上限并且每走一步打印当前点和方向观察是否出现“原地打转”的模式。更隐蔽的问题是漏点。当边界有凹槽时搜索方向如果回退过头会跳过凹槽底部的像素导致轮廓被拉直。这类问题在简单矩形上看不出来但一到不规则形状就会现原形。处理办法是先对二值图像做形态学平滑尽量减少亚像素级锯齿再跑跟踪。6.2 链码对噪声极其敏感链码是逐像素的方向序列最怕的就是边界上的小凸起和凹陷。原本平滑的轮廓因为二值化噪声多出两三个像素的毛刺链码里就会多出一串高频变化的方向码严重干扰形状匹配。我的经验是三步走第一步是形态学闭运算先把小孔和小凸起抹平第二步是链码层面做差分用差分序列替代绝对序列这样大部分高频噪声会被自然抑制第三步如果还需要更强平滑可以对轮廓点坐标做三次样条重采样每隔一个固定间距取一个点再重新编码。6.3 PCA的180度方向模糊和退化情形PCA求出的主方向在0度和180度上是等价的这在很多应用里不是问题但在做旋转校正时就是大麻烦。比如希望物体始终朝上PCA给出的主方向可能是朝下图像翻转之后纹理会完全不一样。解决方法是额外判断一个参考点比如质心到第一个边界点的大致方向用它来锁定主方向的符号。另一个需要警惕的是退化情形。圆形、正方形这类各向同性形状第二特征值和第一特征值几乎相等主方向角度对噪声极度敏感。我处理这类目标时会放弃PCA角度特征只用长短轴比和面积再加链码的统计直方图来区分形状。6.4 坐标体系的统一问题最后一类坑也是最容易忽视的Matlab中bwboundaries返回的坐标是[行, 列]而很多数学概念里的二维坐标是[x, y]即[列, 行]。如果你不统一链码方向会被旋转90度PCA主方向也会差一个角度。我个人的习惯是在所有函数入口统一转换成[行, 列]在输出特征值时再根据需要转为[x, y]。转换只需一行代码pts [pts(:,2), pts(:,1)]但必须在最前面就完成否则后面所有中间结果都要跟着错。最后分享一个非常实用的习惯我会把链码和PCA封装成两个独立的小函数输入统一是轮廓点集输出统一是标准结构体。这样换数据时只需要改预处理那一段后面的特征提取逻辑永远不用动。你先按这套流程跑通一个简单形状再慢慢增加复杂度后面遇到真实项目里的脏图也不会太慌。