简介：面向土木工程与材料分析领域的Matlab图像处理脚本包，用于骨料颗粒轮廓形态的量化评估。通过提取图形边缘轮廓，将轮廓图形n等分并计算n边凹多边形的内角，进而利用骨料棱角计算公式得到棱角指数AI，可应用于集料棱角性试验、混凝土骨料形貌分析及科研绘图。压缩包共11个文件，核心为6个m脚本，覆盖轮廓获取、等分插值、凸性判断、角度求解与主演示流程；配3张示例图展示骨料轮廓处理前后对照，2个文本说明文件包含使用说明和许可信息，整包仅70KB，便于快速下载部署。已有649人浏览学习，适合具备一定Matlab基础、希望直接调用现成算法分析图像边缘特征的研究人员。借助示例图和模块化脚本，读者能一键复现从边缘检测到AI指数输出的完整流程，也能针对自身图像灵活调整n值等参数，减少重复编程，提升实验或工程分析效率。

1. 从边缘轮廓到凹多边形内角：一条能直接跑的 MATLAB 图像处理链路

拿到一张零件图，要做开口角度测量，第一反应是人工量。但图一多就得写程序。问题在于：图像处理拿到的“轮廓”是一串像素点，根本不存在所谓“n 个顶点”，你得先把轮廓按等间距拆成 n 份，再把相邻等分点连成一个 n 边多边形，最后才能求每个内角。麻烦的是这是个凹多边形，角度不能直接用 acos 一把梭，得先判断每个点是凸还是凹。这篇笔记把整条链路拆开：边缘提取、等弧长采样、凹点判定、内角计算，每一步都给了能跑的 MATLAB 代码和参数说明，适合做视觉检测和形状分析的一线工程师。

2. 边缘轮廓提取：先拿到有序且没毛刺的边界点集

这一章解决的是“轮廓怎么来的”和“拿到的轮廓能不能直接用”。常见的坑是：很多人直接对二值图像调 edge() 算子，得到一堆零散线段，既没有拓扑关系也不是闭合曲线。我推荐用 bwboundaries，它输出的边界是首尾相连的有序坐标序列，正好满足后续等分和角度计算的需求。

2.1 二值化与形态学清理：边缘干净比算法花哨更重要

先看第一段代码，做灰度图、自适应二值化和连通域清理：

img = imread('shape.png'); if size(img, 3) == 3 img = im2gray(img); % 彩色转灰度，老版本用 rgb2gray end bw = imbinarize(img, 'adaptive', 'ForegroundPolarity', 'bright'); bw = imfill(bw, 'holes'); % 孔洞填充，防止内边界被追踪 bw = bwareaopen(bw, 100); % 去掉小于100像素的小碎块

逻辑说明：im2gray 把 RGB 图压成单通道灰度图，后续所有操作都在这个灰度图上进行。imbinarize 的 adaptive 模式用的是局部自适应阈值，比全局 Otsu 阈值在光照不均匀时稳很多；ForegroundPolarity 设为 bright 表示亮部是目标，如果你的目标是暗色背景上的亮块，这行是对的，反过来要改成 dark。imfill 把目标内部的孔洞填平，否则后面提取边界时会把内孔也当成一条轮廓。bwareaopen 的参数 100 是面积阈值，单位是像素，需要根据图像分辨率调整，图很大时可以提到 500 或 1000。

提示：如果目标周围背景太复杂，自适应阈值也会失效。我一般会先用 imcrop 把 ROI 裁出来再跑二值化，后续所有处理都针对裁剪后的图，能少踩很多坑。

这一步做完后，bw 应该是一张干净的二值掩膜，目标区域为 1，背景为 0。我见过不少边缘断裂和毛刺问题，根源都在二值化这一步，而不是后面的算法选型。

2.2 边界追踪：bwboundaries 的输出没那么简单

边界追踪直接用 bwboundaries，它返回的是一个元胞数组，每个元素对应一个连通区域的边界坐标。代码：

boundaries = bwboundaries(bw, 'noholes'); [~, idx] = max(cellfun(@(B) size(B, 1), boundaries)); boundary = boundaries{idx}; boundary = double(fliplr(boundary)); % 把 [row, col] 换成 [x, y]

逻辑说明：bwboundaries 的第一个输出是边界坐标，格式是 [row, col]，注意不是 [x, y]。cellfun 遍历所有连通区域，取边界点数最多的那一条，这通常对应图像里最大的目标区域。fliplr 把列顺序从 [行, 列] 换成 [x, y]，这一步很关键，很多角度计算翻车都是因为 xy 没有对齐。double 转换是因为后续插值要求双精度输入。

参数说明：'noholes' 的作用是忽略孔洞边界，只返回每个连通区域的外边界。如果你的目标内部有孔洞且这些孔洞本身有分析价值，就不要加这个参数，但那样后续要自己按层级关系筛选边界。

2.3 边界平滑：等分前必须做的一步

bwboundaries 返回的是像素级边界，相邻点之间是锯齿状折线，直接拿去等弧长采样，等分点会沿着锯齿抖动，最终角度计算结果一塌糊涂。所以要先平滑：

k = 15; boundary(:, 1) = smoothdata(boundary(:, 1), 'gaussian', k); boundary(:, 2) = smoothdata(boundary(:, 2), 'gaussian', k);

逻辑说明：把边界点的 x 坐标和 y 坐标分别作为一维序列，用高斯窗口做平滑。窗口 k 的取值需要根据轮廓总点数调整，我一般取轮廓平均边长的 3～5 倍。如果目标周长是 5000 像素，k 取 100 左右比较合适；如果轮廓本身很平滑，k=15 也可能过大。这个参数后面还会再提到，它直接影响尖角处的内角精度。

3. 把轮廓按弧长 n 等分：累积长度插值实现等距采样

拿到有序边界点之后，下一步是把轮廓平均切成 n 段。这里有一个关键选择：按弧长等分，而不是按顶点序号等分。像素边界顶点疏密不均，直接每隔固定序号取一个点，结果会聚集在密的地方，凹特征点可能漏掉。

3.1 为什么不直接等分顶点序号

bwboundaries 返回的边界点数量是不均匀分布的，转弯处和尖锐区域的点更密，直线段上点更稀。如果按序号等分，等分点在稀疏区域间隔很大，在密集区域挤成一团，完全没有“均匀”的几何意义。按角度等分也不行，因为凹多边形轮廓不是绕形心均匀展开的。唯一可靠的标准是“等弧长”，也就是让每个等分点之间的曲线路径长度相等，这需要先算出每个边界点的累积长度，再做插值。

3.2 累积弧长线性插值代码

下面这个函数是整篇笔记的核心工具之一，完成任意闭合轮廓的等弧长采样：

function [P, arcStep] = resamplePolygonByArclength(boundary, n) % boundary: 有序边界点，[x, y]，首尾重合 % n: 等分份数 seg = diff(boundary, 1, 1); segLen = sqrt(sum(seg.^2, 2)); cumLen = [0; cumsum(segLen)]; totalLen = cumLen(end); s = (0 : n-1)' / n * totalLen; P = interp1(cumLen, boundary, s, 'linear'); P(end + 1, :) = P(1, :); % 闭合：末尾接回起点 arcStep = totalLen / n; end

逻辑说明：diff 对边界点做逐点差分，得到每一小段的位移向量；sqrt(sum(seg.^2, 2)) 计算每段长度；cumsum 得到从起点到每个边界点的累积弧长，这是一个单调递增的数组，可以作为插值的横轴。s 是期望的等分弧长位置序列，从 0 取到 totalLen 的 (n-1)/n，故意不取最后一点 totalLen，因为闭合曲线在 totalLen 处的坐标等于起点坐标，如果把它也取出来，会出现首尾重复点。interp1 在累积弧长轴上做线性插值，得到 n 个等分坐标。最后手动把 P(1,:) 追加到末尾，让多边形在几何上是闭合的。

参数说明：n 是等分份数，也是后续凹多边形顶点的个数。arcStep 返回每段弧长，可以用来人工检查采样是否均匀：理论上等分点之间的欧氏距离应该接近 arcStep，如果偏差过大，说明轮廓本身有严重的曲率突变。

3.3 等分点数 n 怎么选

n 的取值直接决定角度结果的可靠程度。下表是我在形状分析场景下常用的参考值：

n 的初始值带点“玄学”成分，我一般先用 32 跑一遍，再对比 64 的结果，如果角度分布差异很大，说明 n 取小了或平滑窗口不合适。还可以用逼近误差来辅助判断：

D = pdist2(boundary, P); distToPolyline = min(D(:, 1:end-1), [], 2); maxErr = max(distToPolyline);

逻辑说明：pdist2 计算原始边界上每个点到所有等分点的欧氏距离，min 取每个原始点到等分多边形的最短距离，max 得到最大逼近误差。如果 maxErr 超过 2 个像素，说明采样点太少，或者平滑把轮廓削掉太多；如果 maxErr 很小但角度结果乱跳，通常是 n 太大或噪声没滤干净。

4. 用 n 个等分点构造凹多边形并求内角：向量法实现凹点判定

这一章解决标题里的核心问题：n 个等分点组成一个凹多边形，怎么求每个顶点的内角。关键在于：多边形内角不一定是两条边的夹角。凸顶点的内角小于 180 度，凹顶点的内角大于 180 度，而 acos 函数只能算出 0 到 180 度之间的那个角，所以必须先用叉积判断每个点是凸还是凹。

4.1 凹多边形内角不是直接 acos 出来的“夹角”

对第 i 个顶点，取它前一个点和后一个点，用当前点构造两个向量 v1 和 v2，acos(dot(v1,v2) / (norm(v1)*norm(v2))) 得到的是两条线段之间的“小夹角”。这个夹角在凸顶点处正好等于内角，在凹顶点处，真实内角应该是 360 度减去这个小夹角。因此，判断凹点是整个算法成败的分水岭。

判断凹凸靠叉积 z 分量。在统一为逆时针（CCW）绕序的前提下，叉积大于 0 说明局部是凸的，叉积小于 0 说明是凹的。所以计算内角前，必须先用有向面积把多边形绕序统一，否则符号判断会整体反掉。

4.2 完整函数：统一绕序、逐点计算内角

下面这个函数接收等弧长采样后的闭合点集，输出每个顶点的内角（弧度）：

function angles = polygonInteriorAngles(P) % P: 闭合等分点，最后一行与第一行重合 P(end, :) = []; n = size(P, 1); if n < 3 error('至少需要3个点才能构成多边形'); end % 1) 用有向面积统一逆时针绕序 area2 = sum(P(:, 1) .* P([2:n 1], 2) - ... P([2:n 1], 1) .* P(:, 2)); if area2 < 0 P = flipud(P); end % 2) 逐点计算内角 angles = zeros(n, 1); for i = 1:n pPrev = P(mod(i - 2, n) + 1, :); pCur = P(i, :); pNext = P(mod(i, n) + 1, :); v1 = pPrev - pCur; v2 = pNext - pCur; cosA = dot(v1, v2) / (norm(v1) * norm(v2) + eps); theta = acos(max(-1, min(1, cosA))); crossK = v1(1) * v2(2) - v1(2) * v2(1); if crossK < 0 theta = 2 * pi - theta; end angles(i) = theta; end end

逻辑说明：第一步去掉闭合重复点，避免 v1 或 v2 变成零向量。area2 是二倍有向面积，大于 0 表示数学意义上的逆时针绕序，小于 0 就 flipud 把整个点列反过来。第二步逐点遍历，mod 索引实现首尾循环。cosA 计算时在分母加 eps 防止除零；acos 之前用 max/min 把 cosA 限制在 [-1, 1]，这一步不写的话，浮点误差会产生 NaN。crossK 是叉积的 z 分量，在 CCW 体系下为负代表该点是凹点，内角取 2*pi - theta。

参数说明：整个函数不需要手动指定内外方向，绕序统一是自动完成的。唯一要保证的是输入 P 是一个简单多边形，不能自交，否则内角定义本身就会失效。

4.3 自检方法：内角和必须等于 (n-2)*pi

对任意简单多边形，无论凹凸，内角和都等于 (n-2)*pi。这个公式是最便宜的自检手段：

P = resamplePolygonByArclength(boundary, 32); angles = polygonInteriorAngles(P); anglesDeg = rad2deg(angles); % 32个顶点，内角和应该等于30*pi assert(abs(sum(angles) - 30*pi) < 1e-9, '内角和自检失败');

说明：P 返回的是 33 行（32 个等分点加一个闭合重合点），polygonInteriorAngles 内部会去掉最后一行，所以实际顶点数是 32，内角和是 30*pi。如果断言失败，不要先怀疑公式，优先检查输入点集是否自交。自交检测可以用线段相交循环判断，但在大多数场景下，先确认二值化和等分参数更高效。

5. 凹多边形内角实战避坑：5 条高频问题排查记录

这套流程跑下来，最容易翻车的地方其实不在算法本身，而在数据准备和坐标系细节。下面 5 条来自我实际调试时的血泪经验，均按现象、原因、解决三部分列出。

5.1 角度列表出现 NaN

现象：输出 anglesDeg 里某些值是 NaN，另外还有几个角度接近 0 度。

原因：常见有两种。一是闭合点没处理，P 的最后一行和第一行相同，导致 pCur 和 pPrev 重合，v1 变成零向量；二是 cosA 因为浮点误差略大于 1，acos 函数直接返回 NaN。我调试时还遇到过边界坐标是 uint8 类型，差值计算溢出变成 0 的情况。

解决：在函数开头强制P = double(P)，去掉最后一行；acos 前用 max/min 做裁剪。调试时用find(isnan(anglesDeg))定位出错顶点，检查对应三点的坐标值是否合理。

5.2 凹点判定全部反了

现象：肉眼明显凹下去的位置，计算出来的角度小于 180 度；凸出来的位置反而出现大于 180 度的角。

原因：没有把 bwboundaries 返回的 [row, col] 转成 [x, y]，或者绕序没有被正确统一。图像坐标系里 y 轴向下，和数学坐标系相反，直接拿 [row, col] 当 [x, y] 用，叉积符号正好反过来。

解决：边界提取时强制boundary = fliplr(boundary)；内角函数里始终用有向面积统一绕序。只要这两点做到，凹点判定方向就不会错。

5.3 等分点丢了凹槽

现象：图像上有个很深的缺口，但内角列表里没有任何大于 180 度的角度。

原因：n 取值太小，整个凹槽落在两个等分点之间。等分点连线直接跨越凹槽，把凹陷截平了，多边形顶点里根本没有一个落在凹槽内部。

解决：调大 n，先按第 3 章的表格取 32 起步；再用第 3.3 节的逼近误差检查 maxErr。如果逼近误差突增，说明等分点没采到关键凹陷处。

5.4 尖角被平滑吞掉

现象：原本应该接近 90 度的尖角，算出来是 170 度。

原因：smoothdata 的窗口 k 取太大，高斯平滑把尖角磨圆了。等弧长插值再均匀采样时，被磨圆的区域会被多个等分点描述，角度自然被拉平。

解决：平滑窗口用轮廓平均边长的 3～5 倍，而不是固定值 15 或 100。对尖角特征特别重要的目标，我会同时跑原始边界和不平滑边界两版，对比角度差异，差异超过 5 度就说明平滑参数不合适。

5.5 边界毛刺造成角度抖动

现象：同一张图连续跑两次，等分点位置有 1～2 个像素的跳动，角度结果相差 3 度以上。

原因：二值化边缘落在噪声区域，边界点位置本身不稳定；自适应阈值在光照不均匀的边缘地带来回跳。

解决：增强预处理，adaptive 阈值后加一次 imclose 闭合间隙，再用 bwareafilt 保留最大连通域。批量处理时固定所有参数，不逐张调阈值。如果抖动还在，对同一张图重复采样多次，取角度中位数作为结果。

6. 把内角变成可用的形状特征：凹角标注与直方图校验

计算出一组内角后，直接看数字列表没有意义，我习惯把它转成可视化特征，这一步也是确认结果是否符合直觉的最快方式。

6.1 用凹角标注定位凹陷位置

凹多边形的内角大于 180 度，这些角度直接指向轮廓上的凹点。筛选出这些角度并画在原图上：

idx = find(anglesDeg > 180 + 5); % 加5度容差，避免180边缘抖动 imshow(shapeImg); hold on; plot(P(idx, 1), P(idx, 2), 'ro', 'MarkerSize', 8); for i = 1:numel(idx) text(P(idx(i), 1) + 5, P(idx(i), 2) + 5, ... sprintf('%.1f°', anglesDeg(idx(i))), 'Color', 'r'); end

逻辑说明：阈值取 185 度而不是 180 度，是因为数值计算时凸点内角可能因为采样误差跑到 182 度附近，太严的阈值会产生误报。画图时用 P(idx,1) 和 P(idx,2) 作为 x、y 坐标，注意 P 是 33 行闭合点集，index 是等分点序号，对应关系是对的。

6.2 用直方图快速评估角度分布

角度直方图能反映轮廓的几何特征：凹点数量、集中角度分布、是否有异常尖角。代码很简单：

histogram(anglesDeg, 0:15:360); xlabel('内角 / deg'); ylabel('顶点数');

配合第 3 章的逼近误差 maxErr 一起看：如果 maxErr 在 1 像素以内，且直方图分布符合预期，这组参数就可以固下来批量跑。一旦 maxErr 超过 2 像素，优先回查 n 和平滑窗口，不要调内角函数的逻辑。

从那以后，我每次拿到带凹口的轮廓任务，都会强制走一遍：二值化去噪、等弧长采样、绕序统一、内角和自检，最后再画凹角标注图。图形上标出来的角度和肉眼判断一致，这个流程才算真正跑通。希望这套方法能帮你在自己的图像处理任务里少走一段弯路。

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