简介：面向图像处理与骨料形态分析研究者的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(end+1, :) = 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 根据轮廓上的角点数来定。先肉眼看轮廓有几个明显角点记为 k，n 取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(i+1) - 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 度时，一定回去查输入，而不是硬调公式。这套流程帮我挡住了不少脏数据，希望也帮到你。

本文还有配套的精品资源，点击获取