简介:面向Matlab图像处理学习者与研究人员,这份压缩包提供了一套完整的轮廓几何分析工具,可获取图形边缘轮廓并将轮廓n等分,进而求出n边凹多边形各内角,同时结合骨料计算公式得出棱角指数AI,适用于骨料形态定量评价、颗粒形状分析等场景。包内共11个文件,以6个m函数脚本为核心,覆盖轮廓提取、等分插值、内角计算、指数求算等关键步骤,另有3张png示例图用于效果对照,2个txt说明文件辅助解读代码逻辑与运行方式,整包仅70KB,轻量易用。已有649人学习下载,适合需要快速复现骨料棱角指数算法或研究凹多边形几何特征的读者。通过参考配套博客,可清晰理解各函数调用关系与参数含义,借助示例图直接运行验证,各脚本拆分细致,便于在此基础上扩展或迁移到其他图像处理任务。
1. 从图像到凹多边形内角:这事比想象中麻烦,但值得做透
做 matlab 图像处理的人多半都遇到过这个需求:给一张零件的照片,要自动把外形轮廓提出来,然后按弧长等分成若干段,最后把每个折角的内角量出来。听起来像常规操作,但一旦轮廓里有凹部,事情就变味了——用常见的多边角测量方法,凹进去那个角算出来往往只有几十度,和实际完全对不上。这份资源的核心,就是把「轮廓提取 → 弧长 n 等分 → 凹多边形内角求解」整条链路打通,连凹凸点判断和角度修正这种细节都给了可跑通的代码。适合做视觉测量、形状分析或 CAD 逆向的工程师,也适合被凹角问题卡住的学生。下面按我自己的拆解顺序来复现。
2. 轮廓提取与坐标预处理:从像素矩阵到有序边界点
2.1 提取轮廓:为什么用 bwboundaries 而不是找边缘再加连通域
拿到一张图形,第一步是把边缘轮廓变成一组有序的坐标点。常见做法是边缘检测后找连通域,但这一步有个隐患:边缘检测出来的点是二值化的细线,可能断裂或不连续,后续做等分时会出现莫名其妙的跳变。我一般直接对二值图用 bwboundaries,它返回的是按连通域边界顺序排列的像素坐标,天然有序,省掉自己排序的麻烦。
% 读取图像并转为二值图 img = imread('part.png'); gray = rgb2gray(img); % 如果本身就是灰度图可跳过 bw = imbinarize(gray, 'adaptive'); % 自适应阈值,应对光照不均 % 提取边界坐标,'noholes' 表示只取外轮廓不要内部孔洞 B = bwboundaries(bw, 'noholes'); boundary = B{1}; % 取第一个连通域的边界多数场景下直接读图、二值化、取边界就够用了。注意bwboundaries返回的坐标是[row, col],也就是[y, x],这和常规的笛卡尔坐标系相反,如果后面要做向量计算,第一步就得先统一为[x, y]。另外boundary是闭合的,首尾是同一个点,等分前要把尾点去掉。
2.2 坐标预处理:去重复、统一方向、转成双精度
边界坐标直接来自像素网格,采样步长不均匀——斜线处相邻点间距是 1.414,横平竖直处是 1。这本身不影响角度计算,但会影响弧长等分的精度,所以要先做一次轻度平滑和去重。像素边界还有一个特点:可能存在相邻点完全重复(某些边界提取算法会产生),不处理的话累计弧长会出现零长度段,插值的时候会出 NaN。
% 转换为双精度并统一为 [x, y] 顺序 boundary = double(boundary(:, [2 1])); % 去重:删除连续重复点 diff_idx = find(any(diff(boundary), 2)); boundary = boundary([1; diff_idx+1], :); % 闭合处理:如果首尾相同则去掉末尾点 if isequal(boundary(1,:), boundary(end,:)) boundary(end, :) = []; end % 轻度平滑:用滑动平均降低像素锯齿,窗口取 5 k = 5; boundary_smooth = smoothdata(boundary, 'movmean', k);这里smoothdata的窗口参数值得说一句:窗口太小锯齿还在,窗口太大拐角会被削圆,凹进去的角会被磨平,导致后续算出来的内角偏大。对于一般图片轮廓,我试下来窗口 5 到 7 比较稳妥,具体可以在调试时画出来对比一下。等分前把坐标统一成双精度,能避免后面插值时整数矩阵带来的精度丢失。
3. 按弧长把轮廓 n 等分:核心算法与实现细节
3.1 为什么不能直接按像素索引等分
很多人拿到轮廓点序列后,直接boundary(1:k:end)来取等分点,这个做法在轮廓接近圆形、点密度均匀时勉强能用,但对真实图形基本翻车。原因在于边界点的密度不均匀:斜率越大的位置相邻像素距离越远(斜着走一步是 1.414),平直区域距离近(一步是 1),凹进去的曲线部分密度更是忽高忽低。按索引等分等于把物理长度不等的段当成了等长段,后续内角会系统性偏差。
正确的做法是先把轮廓看成一条参数曲线,计算每个点处的累计弧长,再在累计弧长上做等间距采样。这个思路和 CAD 里曲线的等分完全一致,只是落到 MATLAB 上需要手写一小段逻辑。
% 计算相邻点之间的欧氏距离 seg_len = sqrt(sum(diff(boundary_smooth).^2, 2)); % 累计弧长 cum_len = [0; cumsum(seg_len)]; % 总弧长 total_len = cum_len(end); % 生成 n 个等距弧长位置(含首尾,n 段对应 n+1 个点) n = 8; target_pos = linspace(0, total_len, n+1); % 在累计弧长上做插值,重新采样得到等分点坐标 point_x = interp1(cum_len, boundary_smooth(:,1), target_pos, 'spline'); point_y = interp1(cum_len, boundary_smooth(:,2), target_pos, 'spline'); equipoints = [point_x', point_y'];这段代码的逻辑一句话就能讲清:先算每段微元长度,累加得到每个原始点的弧长位置,然后构造等间距的目标弧长序列,用插值函数反查坐标。interp1前两个参数分别是「原采样位置」和「原采样值」,第三个参数是「目标位置」。这里务必使用spline或pchip,不要用默认的linear。线性插值在轮廓弯曲处会把角切成直线段,凹角直接变没。spline光滑但容易在突变处出现过冲,pchip保形更好,我通常优先pchip,视觉上更接近原始轮廓。
3.2 n 取多少合适:分辨率与稳定性的平衡
等分数 n 是这个流程里最敏感的参数。n 太小,轮廓的弯曲特征被强行拉直,凹角位置偏移;n 太大,等分点密集后对噪声和锯齿敏感,个别点落在轮廓的微小波动上,角度计算结果抖动剧烈。经验法则是先看轮廓的总像素周长,n 取周长的 1/50 到 1/80,比如周长 400 像素,n 取 6 到 8。这里给的n = 8适合中等复杂度轮廓。
还有一点容易被忽略:插值后等分点的首尾要闭合。linspace生成的序列天然包含起点和终点,而终点恰好是起点本身(因为轮廓闭合),所以计算内角时最后一个角由equipoints(end)和equipoints(1)共同构成,不要丢掉这个角。很多新手算完 n 个点只算了 n-1 个角,就是没意识到闭合链要多算一次。
4. 凹多边形内角计算:向量夹角法与凹凸点判定
4.1 基础内角:三个相邻等分点的向量夹角
拿到 n+1 个等分点后,内角计算就可以开始了。基本思路是对每个顶点P_i,用它的前一个点P_{i-1}和后一个点P_{i+1}构造两个向量,然后求夹角。夹角公式是atan2加余弦定理,不能用简单的acos(dot(...)),后者数值不稳定,角度接近 0 度或 180 度时误差会放大。
% 计算多边形各顶点的内角,返回角度制(0~360) function angles = polygon_interior_angles(pts) m = size(pts, 1); angles = zeros(m, 1); for i = 1:m i_prev = mod(i-2, m) + 1; % 前一个点(环形索引) i_next = mod(i, m) + 1; % 后一个点 v1 = pts(i_prev, :) - pts(i, :); v2 = pts(i_next, :) - pts(i, :); % 用 atan2 求 v1 到 v2 的有向夹角,范围 (-pi, pi] ang = atan2(det([v1; v2]), dot(v1, v2)); % 转换为 0 ~ 360 度 if ang < 0 ang = ang + 2*pi; end angles(i) = ang * 180 / pi; end end这里关键是理解有向夹角。atan2(det, dot)求得的是 v1 绕顶点旋转到 v2 的有向角度,范围是 -180 到 180 度。对于凸多边形,每个顶点的有向角方向一致(全正或全负),修正后直接就是内角;而凹多边形内部会出现反方向顶点。上面的代码把所有负角都加 360 度,得到的是「转向角」,在凹点处会给出小于 180 度的角——这本身没错,但它不是内角,是外角的补关系,需要结合凹凸性做修正。
4.2 凹凸点判定与内角修正
凹多边形和凸多边形的核心区别:凸多边形所有顶点转向方向一致;凹多边形存在个别顶点转向相反。上一节的代码如果直接用,凹点的内角会被算成小角,而实际凹进去的内角是大于 180 度的。
判断凹凸点的方法很多,常用的是多边形有向面积符号法或局部叉积法。局部叉积的做法是:对顶点P_i,计算边向量P_i-P_{i-1}和P_{i+1}-P_i的叉积,符号正负代表在闭合路径上该点是顺时针方向还是逆时针方向的转向。如果整个轮廓统一按逆时针(或顺时针)排列,那么凸点的叉积符号相同,凹点的符号相反。
% 统一轮廓方向为逆时针(通过有向面积判断并翻转) area = sum(boundary_smooth(:,1) .* ... boundary_smooth([2:end 1],2) - ... boundary_smooth([2:end 1],1) .* boundary_smooth(:,2)) / 2; if area < 0 boundary_smooth = flipud(boundary_smooth); end % 等分点同样按逆时针顺序排序 equipoints = sort_points_ccw(equipoints); % 计算内角(修正凹点) angles = zeros(n+1, 1); for i = 1:n+1 i_prev = mod(i-2, n+1) + 1; i_next = mod(i, n+1) + 1; v1 = equipoints(i_prev, :) - equipoints(i, :); v2 = equipoints(i_next, :) - equipoints(i, :); ang = atan2(det([v1; v2]), dot(v1, v2)) * 180 / pi; % 逆时针轮廓下,叉积符号为负说明是凹点 cross_val = v1(1)*v2(2) - v1(2)*v2(1); if cross_val > 0 ang = 360 - abs(ang); % 凹点的内角修正 else ang = abs(ang); end angles(i) = ang; end修正逻辑背后是几何直觉:逆时针方向走轮廓时,凸点的转向角是正的(右转,小于 180 度),凹点的转向角是负的(左转),其内角应为 360 度减去转向角的绝对值。sort_points_ccw是对等分点按极角排序的辅助函数,确保计算顺序和轮廓走向一致。
这个修正策略在一个前提下才能用:轮廓方向统一。如果原图轮廓是顺时针,所有符号都会反过来,凹点会被误判成凸点,内角整体变成「镜像」。所以前面用有向面积做方向统一的步骤不能省。
5. 避坑与常见问题:这条链路里最容易翻车的四个细节
5.1 现象:计算出来的凹角全部小于 180 度,凸角反而大于 180 度
原因:轮廓方向没有统一,叉积符号判断与预期相反,凹凸点完全对调。解决:在等分前检查有向面积,强制把轮廓和等分点都调整为逆时针;如果图像坐标系 Y 轴向下(图像常见坐标系),做笛卡尔转换时记得翻转 Y 轴,否则方向判断也会反。
5.2 现象:n 等分点落在轮廓附近但不在轮廓上,画出来是锯齿折线
原因:插值方式用了默认的线性,或者等分数太大,spline 过冲把点推到轮廓外。解决:先画出等分点和原始边界对比,检查是否贴合;优先改用pchip插值,它对尖锐转角的过冲控制更好;同时把等分点数量降到周长 1/60 左右再试。
5.3 现象:角度结果与预期相差 90 度左右,而且方向性混乱
原因:图片坐标系下行列坐标未转换为 XY 笛卡尔坐标,导致向量运算时把行当成了 y 但方向向下未翻转。解决:bwboundaries返回的[row, col]在转[x, y]时,要把 row 取负后才符合数学坐标系;或者全程都在图像坐标系里计算,只要向量运算前后统一,角度结果是可以保持不变的,但方向判断必须始终一致。
5.4 现象:同一个轮廓换一张图片,角度结果波动很大
原因:分割阈值或图像尺寸不同导致边界位置偏移,等分点对轻微变化过度敏感。解决:先对二值图做形态学闭运算,填补边界缺口;再做边界平滑;最后输出角度前打印等分点坐标,肉眼核对凹点位置是否命中真实凹部,定位到个别等分点漂移时,可以微调插值方法而不是整体改参数。
6. 进一步验证:用合成凹多边形校验整套算法
实际项目里最怕的是「算完了不知道对不对」。我的习惯是用合成多边形做 ground truth 验证:人为构造一个顶点坐标已知的凹多边形,跑完整流程,看程序输出的角度能否回推出原有的顶点角度。如果吻合,说明整条链路可靠,再切换到真实图像。
% 构造一个已知的凹五边形(含一个凹点) theta = linspace(0, 2*pi, 6)'; r = 100 * ones(6, 1); r(4) = 40; % 第四个点径向缩小,形成一个凹部 xs = r .* cos(theta); ys = r .* sin(theta); % 在轮廓上密集采样模拟图像的边缘点 t_interp = linspace(0, 1, 1000)'; cx = interp1(linspace(0,1,6), xs, t_interp, 'spline'); cy = interp1(linspace(0,1,6), ys, t_interp, 'spline'); boundary_smooth = [cx cy]; % 然后走第 3、4 章相同的等分与角度计算流程r(4) = 40就是在构造一个可预期的凹点。跑完对比输出的角度序列,凹点那个角应该接近 240 度多一点,其余凸点接近标准五边形的内角。这里提一个我个人踩过坑后的固定动作:每次跑完等分,我都会额外输出一张等分点和原轮廓叠加的图,等分点有「悬空」立刻能发现;从那以后,所有角度结果我都会先做一次方向性检查,再进数据表,这已经成为惯性。希望你也能把这套验证流程写进自己的代码里,能省下不少调试的时间——希望帮到你。
本文还有配套的精品资源,点击获取