简介面向图像处理与骨料形态分析研究者的MATLAB实用代码工具包用于提取图形边缘轮廓、将轮廓n等分并求出凹多边形内角进而按骨料棱角计算公式获得AI指数。压缩包共11个文件、70KB其中6个m脚本涵盖主程序、轮廓采样、内角求解等核心环节3张png图片提供前后对照示例2个txt文档补充说明与参考信息整体结构紧凑便于直接运行。该工具已有649人学习下载适合从事粗骨料形状表征、混凝土微观分析或数字图像处理的学生与工程师使用。借助其中代码读者可掌握轮廓离散化与凹多边形内角计算的完整流程同时示例图像与说明文档可帮助快速上手并根据自身图像数据迁移应用为骨料棱角性评价提供稳定可复用的MATLAB计算分析基础。1. 从像素图到凹多边形内角轮廓等分与角度计算的一次完整实现做图像几何测量的同学大都碰到过这个需求输入一张带轮廓的图形把边缘等分成 n 份再求出这 n 个顶点连成的多边形内角最好还能把凹点标出来。听起来就是“边缘提取 余弦定理”的事可真上手跑一遍就会发现直接用acos算内角凹点角度永远被算成锐角内角和怎么都对不上。这份 Matlab 资源的核心就两件事按弧长把轮廓 n 等分再用有符号转向角求凹多边形内角。前者解决采样偏差后者解决凹凸方向中间还埋着边界顺序、闭合点、分辨率这几个坑。适合做尺寸测量、形状判异、医学图像轮廓分析的人。新手能照着脚本跑通熟手可以直接拿走核心函数改参数。下面我从bwboundaries提取轮廓开始逐步拆到弧长重采样、转向角求内角最后用内角和定理做自检。2. 边缘轮廓提取bwboundaries 顺序追踪与预处理三件套2.1 为什么顺着 find 找点走不通很多第一次做边缘提取的人习惯这么写用edge拿二值边缘图再用find把非零像素坐标全部取出来。img imread(part.png); E edge(img, canny); [row, col] find(E); plot(col, row, .); axis equal;这段代码能画出轮廓形状但find返回的坐标是按行扫描顺序排列的和轮廓线本身的走向没有任何关系。也就是说你拿到的是“一堆散点”不是“一条有顺序的路径”。如果直接用这些点连多边形结果往往是折线来回交叉像蜘蛛网一样后续根本没法算角度。我一般会绕开find改用bwboundaries。它沿着二值区域的邻接关系做边界追踪输出的坐标天然就是有序的一层轮廓对应一个 cell边界点按行进方向排列拿过来就能当多边形的顶点序列用。这是整条链路的第一个关键选型。2.2 预处理三步阈值、闭运算、孔洞填充与面积筛选直接对原始灰度图跑bwboundaries通常不干净因为边缘毛刺、断线、小噪点都会干扰轮廓顺序。我会先做一套固定预处理再提取轮廓。img imread(part.png); if size(img, 3) 3 img rgb2gray(img); end bw img graythresh(img); % 全局阈值分割 bw imclose(bw, strel(disk, 3)); % 闭运算填补断线和小窄缝 bw bwareaopen(bw, 80); % 删除面积小于80像素的连通块 bw imfill(bw, holes); % 填充内部孔洞保证外轮廓封闭参数说明graythresh用的是 OTSU 全局阈值光照不均匀的图可以改成局部阈值但那不是这份资源的主线。strel(disk, 3)的半径我一般取 2 到 5。半径太小堵不上断口太大则会把相邻结构黏连在一起轮廓形状变形。bwareaopen(bw, 80)的 80 是噪点面积阈值按你图上最小真实目标的像素数来设。目标本身很小就调低背景杂点多就调高。imfill是否使用要看目标求外轮廓内角孔洞一定得填否则bwboundaries会把内孔洞也当作边界轮廓输出。预处理完提取轮廓boundaries bwboundaries(bw, noholes); B boundaries{1}; % B: Nx2第一列是行坐标x第二列是列坐标y P [B(:, 2), B(:, 1)]; % 转成 (x, y)每行一个轮廓点noholes是第二个关键参数。它告诉函数“只取外边界不要内孔边界”否则你拿到的第一层边界可能只是某个内部空洞。2.3 多个连通域时怎么选目标边界如果图里有多个物体boundaries的 cell 数量会大于 1。我习惯按边界点数筛选目标轮廓通常是最长的那一圈。lenList zeros(numel(boundaries), 1); for k 1:numel(boundaries) lenList(k) size(boundaries{k}, 1); end [~, idx] max(lenList); B boundaries{idx}; P [B(:, 2), B(:, 1)];这样做的前提是目标在画面里占比最大。如果目标不是最大的就手动指定idx或者加一个“边界点数量落在某个区间”的筛选条件比盲目取第一个 cell 稳得多。3. 轮廓 n 等分按弧长重采样而不是按像素下标等距取点3.1 按像素下标等距取点错在哪拿到有序轮廓点后第一个直觉是“等分嘛每隔 m 个点取一个坐标”。step floor(size(P, 1) / n); Pq P(1:step:end, :);这种做法在轮廓点分布均匀时勉强能用但真实图像轮廓几乎不可能均匀圆弧部分曲率大边缘检测出来的像素点密度高长直线段部分点很稀疏。按下标取点结果就是圆弧上等分点扎堆直线段上点间距又太大连出来的多边形严重偏离原始轮廓。另一个隐蔽问题是凹角附近。轮廓像素在转角内侧密集、外侧稀疏按索引取点可能把凹角顶点直接跳过去凹特征被抹平。正确做法是按“累计弧长”做参数化先算每两个相邻轮廓点的距离累加得到总弧长然后在总弧长上均匀取 n 个位置用插值求每个位置对应的坐标。这个思路等价于把轮廓曲线拉直成一条线段再在线段上等分。3.2 弧长重采样函数与参数说明下面这个函数是我在多个项目里一直沿用的版本function Pq resampleContour(P, n) % P: m x 2有序轮廓点第一列x第二列y % n: 等分数 if size(P, 2) ~ 2 error(输入必须是 Nx2 的坐标数组); end % 如果首尾没闭合先补上起点 if norm(P(1,:) - P(end,:)) 1e-6 P [P; P(1,:)]; end d diff(P, 1, 1); segLen sqrt(sum(d.^2, 2)); cumLen [0; cumsum(segLen)]; % 累计弧长 total cumLen(end); s (0:n-1) * total / n; % 等分点所在弧长位置 Pq interp1(cumLen, P, s, linear); Pq(end1, :) Pq(1, :); % 补闭合点方便画图和多边形运算 end几个容易忽略的点cumLen必须严格单调递增interp1才不出错。轮廓如果存在连续重复像素要先unique(P, rows, stable)去重。s从 0 取到total*(n-1)/n不取total本身。这样最后一个等分点和起点之间仍然隔着一段total/n不会和起点重合。线性插值出来的等分点是亚像素浮点坐标。它可能落在两个原始像素点之间这是特性不是错误角度计算完全不受影响。调用方式n 10; Pq resampleContour(P, n); figure; plot(P(:,1), P(:,2), -, Color, [0.6 0.6 0.6]); hold on; plot(Pq(:,1), Pq(:,2), r-o, MarkerSize, 6); axis equal;画出来之后应该看到等分点在原轮廓线上均匀分布圆弧密集区域点距也保持一致。如果发现某个角点上挤了两个点大概率是原始轮廓该处有重复像素回去去重再跑。3.3 等分数 n 怎么选n 的选择决定后续角度质量。我给两条实用经验第一n 必须小于原始轮廓点数 m实用上建议n round(m / 2)。等分点间距小于 1 个像素时相邻顶点会落在这段直线上角度趋近 180 度凹点被稀释成平角。第二n 根据轮廓上的角点数来定。先肉眼看轮廓有几个明显角点记为 kn 取2*k到4*k之间。角点数少比如三角形n 取 6 到 12轮廓带锯齿或小弧度特征多n 就得往大取。也可以按弧长密度折算n round(total / 5);意思是每隔 5 个像素取一个等分点。像素单位下这个经验值比较稳但最终还是要看等分结果图来调。4. 凹多边形内角用有符号转向角区分凹点与凸点4.1 为什么代码里不写 acos常规向量夹角公式是这样的对顶点 i取来向向量v1 V(i) - V(i-1)去向向量v2 V(i1) - V(i)那么夹角余弦等于dot(v1, v2)除以两个向量的模长乘积再用acos反解角度。cosTheta dot(v1, v2) / (norm(v1) * norm(v2)); theta acos(cosTheta) * 180 / pi;这个公式用在凸多边形上没问题但一遇到凹多边形就翻车。原因在于acos的返回值被限制在 0 到 180 度之间它天然丢失了“转向方向”。一个真正的凹点内角是 270 度但来向向量和去向向量之间的夹角只有 90 度acos算出来只有 90 度凹点被直接当成凸点。还有一个连带问题凹点被误判为凸点之后内角和定理(n-2)*180怎么都校验不过整个计算过程变成黑匣子。4.2 有符号转向角与多边形方向统一要区分凹凸必须引入有符号角度。核心思路是在顶点 i 处求来向向量到去向向量的转向角用atan2同时取出方向和大小。turn atan2(cross2, dot2);其中cross2 v1(1)*v2(2) - v1(2)*v2(1)是二维向量叉积的标量dot2是点积。这个turn落在[-pi, pi]正负代表左转、右转。接下来的问题是内角和转向角怎么换算这里依赖多边形顶点顺序的方向。如果多边形整体是逆时针走向左转对应凸点右转对应凹点内角公式是interiorAngle pi - turn;如果顶点顺序是顺时针同样的公式会把凸点算成凹点。所以计算之前必须先统一方向。注意不要用polyarea判断方向它返回的是面积绝对值拿不到正负。要自己写带符号面积。function sA signedArea2(P) % 带符号面积逆时针为正顺时针为负 x P(:,1); y P(:,2); sA 0.5 * sum(x .* circshift(y, -1) - circshift(x, -1) .* y); end完整内角计算函数function angDeg interiorAngles(V) % V: 多边形顶点坐标首尾可以重复也可以不重复 if norm(V(1,:) - V(end,:)) 1e-6 V V(1:end-1, :); % 去掉闭合点 end if signedArea2(V) 0 V flipud(V); % 统一为逆时针方向 end n size(V, 1); angDeg zeros(n, 1); for i 1:n pPrev V(mod(i-2, n) 1, :); pCur V(i, :); pNext V(mod(i, n) 1, :); v1 pCur - pPrev; v2 pNext - pCur; cross2 v1(1)*v2(2) - v1(2)*v2(1); dot2 v1(1)*v2(1) v1(2)*v2(2); turn atan2(cross2, dot2); angDeg(i) (pi - turn) * 180 / pi; end end说明几点flipud是上下翻转顶点数组能把顺时针序列变成逆时针序列左右顺序也跟着反但多边形形状不变。atan2对零向量会返回 0所以输入顶点不能有重复点。等分时 n 取太大就会出现这种情况等分完先unique去重。pi - turn当turn接近 0 时得到 180 度对应共线顶点这在 n 偏大时可能出现不影响整体判断。4.3 用 L 形轮廓验证角度分类我自己验证算法惯用一个 L 形六边形凹点、凸点、直角都包含Vtest [0 0; 2 0; 2 1; 1 1; 1 2; 0 2]; angDeg interiorAngles(Vtest)输出90.0000 90.0000 90.0000 270.0000 90.0000 90.0000第 4 个点是凹点内角 270 度其余五个是 90 度。内角和正好 720 度而六边形理论内角和(6-2)*180 720对得上。凹点判定规则就很直观了concaveIdx find(angDeg 180 1e-6);容差1e-6是为了防止浮点误差把理论直角算成 179.9999 度。大于 180 度的顶点就是凹点这正是“n 边凹多边形内角”这个需求里要输出的核心信息。5. 避坑与常见问题边界方向、重复点与退化角度的排查记录5.1 角度全变成外角凸点凹点完全反了现象同一组轮廓点有时算出来的角度像“外角”凸点显示 120 度凹点显示 60 度整体感觉是内外角互换。原因顶点顺序是顺时针而公式默认按逆时针处理。pi - turn在顺时针序列下得到的是外角不是内角。解决在interiorAngles函数入口处强制统一方向。我习惯先用signedArea2判断面积为负就flipud。所有调用这个函数的场景都走同一套入口不要在外面手动翻转。5.2 n 取太大等分点出现重叠角度乱跳现象n 从 8 改成 30等分点一眼看过去好几个叠在一起角度序列里出现 0 度或接近 0 度的异常值。原因轮廓像素分辨率有限插值点间距小于 1 个像素时多个等分点落在同一条直边上甚至同一像素附近。转向角在零向量附近失去意义。解决限制n size(P, 1) - 1实用上限是round(size(P, 1) / 2)。等分生成后加一步去重[Pq, ia] unique(Pq, rows, stable); angDeg interiorAngles(Pq);注意unique行数变了n 也要同步更新否则后面标注顶点编号会错位。5.3 从 edge 结果拿点云轮廓自交现象不用bwboundaries直接用edge加find等分和角度计算全乱连多边形都画不出来。原因find输出的坐标按行扫描排列和几何走向无关。散点之间没有顺序关系连线必然交叉。解决不要从边缘图直接提取坐标。先得到二值区域再用bwboundaries或bwtraceboundary做边界追踪拿到的点序列天然按轮廓走向排列。这是我反复提醒自己的一条边缘检测给的是“哪些像素是边缘”边界追踪给的是“边缘像素以什么顺序连接”后者才是几何计算需要的。5.4 在仿真坐标里对的代码换到图像坐标全错现象同一套内角代码跑人工构造的坐标数组正常跑bwboundaries输出的图像坐标就反了。原因图像的行坐标向下增长本质是左手坐标系。同样一段路径在图像坐标里算出的转向角符号和数学坐标系正好相反。解决不要试图混两个坐标系。要么全程用图像坐标[col, row]要么转换成[x, y]后再算。我的习惯是进入处理链路第一步就统一成P [B(:,2), B(:,1)]后续所有函数只认(x, y)。方向正负不用管signedArea2统一处理即可。5.5 多目标图像选错边界角度结果没有规律现象一张图里有好几个工件boundaries返回多个 cell。取{1}直接算角度序列毫无规律内角和也不对。原因bwboundaries对每个连通域都返回一条边界第一个 cell 不一定是目标轮廓还有可能把背景中的某个亮斑当成目标。解决用边界点数筛选或者人工指定目标索引。lenList cellfun((b) size(b, 1), boundaries); [~, idx] max(lenList); B boundaries{idx};如果目标不是最长边界就画出来确认目标在boundaries中的序号再写死索引。这个步骤属于跑批量任务前必须做的人眼核对跳过它迟早要翻车。6. 验证技巧内角和自检与凹点自动标注的完整脚本6.1 一条命令验证内角和定理凹多边形也好凸多边形也好只要是简单多边形内角和恒等于(n-2)*180。这条定理是最便宜的质量门禁。sum(angDeg) % 实际内角和 expected (numel(angDeg) - 2) * 180; % 理论内角和 assert(abs(sum(angDeg) - expected) 1, 内角和校验失败);我一般把偏差阈值设为 1 度。超过 1 度优先查三件事等分点有没有重复点、顶点顺序有没有统一方向、轮廓预处理是否把两个物体黏连在了一起。这个断言能在批量处理几十张图时自动报警替我挡住大部分脏数据。6.2 把凹点自动标注出来的完整函数把前面所有代码合并成一个可复用的函数输入图片路径和等分数输出角度和凹点位置并自动画图。function [angDeg, Pq, concaveIdx] analyzeContourAngle(imgPath, n) img imread(imgPath); if size(img, 3) 3 img rgb2gray(img); end bw img graythresh(img); bw imclose(bw, strel(disk, 3)); bw bwareaopen(bw, 80); bw imfill(bw, holes); boundaries bwboundaries(bw, noholes); lenList cellfun((b) size(b, 1), boundaries); [~, idx] max(lenList); B boundaries{idx}; P [B(:, 2), B(:, 1)]; Pq resampleContour(P, n); angDeg interiorAngles(Pq); concaveIdx find(angDeg 180 1e-6); figure; plot(Pq(:,1), Pq(:,2), b-o, MarkerSize, 5); hold on; plot(Pq(concaveIdx, 1), Pq(concaveIdx, 2), r*, MarkerSize, 14); for k 1:n text(Pq(k,1), Pq(k,2), sprintf(%d: %.1f deg, k, angDeg(k)), ... FontSize, 8); end axis equal; title(sprintf(n%d, sumA%.1f, expected%.1f, ... n, sum(angDeg), (n-2)*180)); end调用示例[angDeg, Pq, concaveIdx] analyzeContourAngle(part.png, 12);图上蓝色圆点是等分顶点红色星号是凹点每个顶点旁标注序号和角度值。标题里的实际内角和与理论值一眼就能对上不用单独打印验证。从那以后我每次拿图算轮廓角度都会强制走一遍“方向统一 → 弧长等分 → 内角和自检”这三步。尤其是sum(angDeg)和(n-2)*180偏差超过 1 度时一定回去查输入而不是硬调公式。这套流程帮我挡住了不少脏数据希望也帮到你。本文还有配套的精品资源点击获取