简介面向土木工程与材料分析领域的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 的取值需要根据轮廓总点数调整我一般取轮廓平均边长的 35 倍。如果目标周长是 5000 像素k 取 100 左右比较合适如果轮廓本身很平滑k15 也可能过大。这个参数后面还会再提到它直接影响尖角处的内角精度。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典型用途特点8粗粒度形状对比抗噪强但细节丢得多凹槽容易消失16-24零件开口、缺口角度测量平衡细节和稳健性32-64凹槽数量多或角度要求高接近真实轮廓但噪声放大明显大于100轮廓描述子计算等价于原始轮廓角度意义不大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 和 v2acos(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 因为浮点误差略大于 1acos 函数直接返回 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 取太大高斯平滑把尖角磨圆了。等弧长插值再均匀采样时被磨圆的区域会被多个等分点描述角度自然被拉平。解决平滑窗口用轮廓平均边长的 35 倍而不是固定值 15 或 100。对尖角特征特别重要的目标我会同时跑原始边界和不平滑边界两版对比角度差异差异超过 5 度就说明平滑参数不合适。5.5 边界毛刺造成角度抖动现象同一张图连续跑两次等分点位置有 12 个像素的跳动角度结果相差 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 和平滑窗口不要调内角函数的逻辑。从那以后我每次拿到带凹口的轮廓任务都会强制走一遍二值化去噪、等弧长采样、绕序统一、内角和自检最后再画凹角标注图。图形上标出来的角度和肉眼判断一致这个流程才算真正跑通。希望这套方法能帮你在自己的图像处理任务里少走一段弯路。本文还有配套的精品资源点击获取