基于MATLAB图像处理的大型几何锻件尺寸与圆度测量

简介聚焦大型几何锻件的非接触式测量面向机械制造、质量检测及数字图像处理方向的学习者提供一套基于 MATLAB 的实践案例用于理解从图像采集到尺寸与形状参数计算的完整思路。压缩包共 3 个文件包含 MATLAB 源代码、BMP 原始图像和 JPG 效果图整体约 77KB体积轻量便于直接查阅算法实现、对照输入图像与最终可视化结果。借助其中的 MATLAB 脚本读者可梳理灰度转换、滤波、边缘检测、特征点定位等数字图像处理环节并观察霍夫变换、模板匹配一类测量算法如何服务于锻件边界与关键点提取。目前已有 144 人学习适合作为课程设计、项目复现或质量检测方案调研的参考样例帮助快速建立 MATLAB 图像测量流程的整体认知。1. 从一张锻件照片到几个尺寸数字这套 MATLAB 测量方案解决什么问题车间里刚下线的轴类锻件还带着余温卡尺伸不进去三坐标排期又要等上半天。这种活儿最磨人的地方往往不是测不准而是等不起。这套基于 MATLAB 数字图像处理的大型几何锻件测量方案走的就是「拍一张、算一遍」的路子机位固定好脚本自动完成去背景、分割、亚像素边缘定位和几何量计算直接给出长度、直径、圆度误差。资源里源代码、图面、效果图三样齐全效果图尤其关键它让你看见灰度图怎么一步步变成数字而不是只甩给你一堆函数名。它适合两类人手上已有相机和工装、想快速搭一套在机测量原型的设备工程师以及拿数字图像处理当课题、需要一条完整可复现链路的同学。它不追求微米级精度追求的是把流程跑通、把误差来源交代清楚。2. 技术栈与目录结构一条 MATLAB 图像处理管线怎么搭起来整套东西没有用深度学习也没有上 OpenCV 混编就是纯 MATLAB 加图像处理工具箱。这个选型乍看有点「保守」但对现场调试反而是好事每一步中间结果都能imshow出来看参数改一个数字立刻见效不用重新编译、不用配环境变量。锻件测量的难点从来不在算法多花哨而在于光照、氧化皮、油污这些脏东西会不断把你的分割结果带偏能快速定位是哪一步出了问题比换一个更强的模型实在得多。2.1 三份资源怎么配合源代码、图面、效果图三份东西各自的位置要分清不然拿到手会不知道从哪看起。源代码是主链路负责从读图到出数。它是可以改的部分参数、阈值、标定系数都在这边。图面解决的是「测哪里」的问题。锻件的被测要素不止一个总长、法兰外径、台阶直径、圆度每一个都要在图面上标出测量部位和基准。图面还顺带告诉你公差带是多少这样脚本算出来的数才有判据——没有公差数字就只是数字。效果图是最容易被忽略、却最该先看的一份。它记录的是每一步处理的中间结果原始图、去背景后的图、二值掩膜、边缘叠加图、拟合圆叠加图。调试的时候如果最终尺寸不对正确的顺序不是去改参数而是拿自己的中间图和效果图逐张比对看是从哪一步开始分叉的。我一般会把这个比对过程叫作「对齐中间态」比盲调阈值快十倍。提示先跑通一张图把它的五个中间态和自己的结果并排看一遍再谈批量。跳过这一步直接批量跑出来的 CSV 大概率是一堆没法解释的数。2.2 目录结构与文件职责一份能长期维护的 MATLAB 工程最怕所有逻辑堆在一个.m文件里。这套资源按处理阶段拆函数主入口只负责编排我按自己的习惯整理成下面这样路径职责改动频率src/main_measure.m主入口串起读图→分割→边缘→计算→出图低src/preprocess.m灰度化、背景均匀化、保边去噪中src/segment_part.m阈值分割 形态学清理输出工件掩膜高src/subpix_edge.m沿灰度剖面做亚像素边缘定位低src/calib_px2mm.m标定块求像素当量低src/geom_feature.m外接矩形、圆拟合、圆度误差中src/visualize_step.m中间结果出图供与效果图比对低src/batch_run.m批量遍历目录结果写 CSV低config/params.m所有可调参数集中返回一个 struct高data/calib/标定块图像低data/raw/待测原始图低doc/drawing/图面测量部位与基准标注低doc/result/效果图每步一张低把参数抽到params.m是个小动作收益很大不同批次、不同光照条件下你只需要换一个参数文件主链路代码一行不动。% config/params.m —— 所有可调项集中在这里主流程只读不改 function P params() P.bgDisk 60; % 背景估计的结构元半径要大于工件最宽处 P.medFilt [5 5]; % 中值滤波窗口 P.threshBias 0.35; % Otsu 阈值的下限兜底 P.minArea 500; % 小于这个像素面积的连通域当碎屑丢掉 P.openDisk 3; % 开运算半径断细连接 P.closeDisk 8; % 闭运算半径补边缘缺口 P.cannyThresh [0.06 0.18]; P.px2mm []; % 空表示现场标定填数值则跳过标定 P.outDir fullfile(doc, result); endbgDisk这个值要跟你的视野匹配结构元半径必须大于工件在图像里的最大宽度否则开运算会把工件本身也当成背景估进去减完以后工件中间直接塌一块。这是新手最常踩的一个坑后面第 5 章还会细说。2.3 运行环境与依赖自检工具箱依赖只有 Image Processing Toolbox 一个硬性的。imbilatfilt保边双边滤波是 R2018a 之后才有的im2gray是 R2018b 之后bwareafilt更早一些。我给的建议是 R2019b 及以上能省掉一堆「未定义函数」的报错。真要跑在老版本上imbilatfilt换成imgaussfilt也能凑合只是边缘会糊一点亚像素定位的收益就打折了。跑主流程之前先自检一次比跑到一半崩掉体面得多% check_env.m —— 开局自检缺工具箱直接断言退出 function ok check_env() need {Image Processing Toolbox}; v ver; have {v.Name}; fprintf(MATLAB 版本: %s\n, version); ok true; for k 1:numel(need) has any(strcmp(have, need{k})); fprintf( %-30s %d\n, need{k}, has); ok ok has; end assert(ok, 缺少 Image Processing Toolbox装完再跑主流程); endver返回的是结构体数组{v.Name}把名字抽成 cell再用strcmp逐项比对。这段代码本身没什么技术含量但它在团队协作里非常值钱别人拿到你的包报错第一句话往往是「你这跑不起来」有了自检输出你一眼就能回他「你工具箱没装」。还有一点MATLAB 的当前路径很关键。main_measure.m里我习惯用mfilename(fullpath)反推工程根目录再addpath(genpath(root))这样不管用户从哪个目录启动相对路径都不会错位。% main_measure.m 开头几行 root fileparts(mfilename(fullpath)); % src 的上一级 root fileparts(root); addpath(genpath(root)); P params(); check_env();3. 核心链路拆解一张锻件照片到边缘坐标的六步处理这条链路我拆成六步读图、预处理、分割、形态学清理、边缘提取、标定换算。每一步的输出都是下一步的输入中间任何一步歪了后面的数字全是白算。所以这一章我按顺序把代码摊开讲重点放在「为什么这么写」而不是「这函数什么意思」。3.1 预处理先把照明不均匀干掉车间打光很难做到均匀尤其锻件尺寸大视野一宽中间亮四角暗几乎是必然的。直接对原图做 Otsu 阈值分割结果往往是中间被切成工件、边缘被切成背景或者反过来。常见的做法是先估计背景再相减。用一个大半径的结构元对灰度图做形态学开运算得到的就是「缓慢变化的照明曲面」原图减掉它剩下的就是相对均匀的工件。% preprocess.m —— 灰度化 背景均匀化 保边去噪 function Ig preprocess(imgPath, P) I imread(imgPath); if size(I, 3) 3 Ig rgb2gray(I); % 锻件表面色差不大灰度信息足够 else Ig I; end Ig im2double(Ig); % 大结构元开运算估计照明背景半径要大于工件在图像中的最大宽度 seBg strel(disk, P.bgDisk); bg imopen(Ig, seBg); Ic Ig - bg; % 减去照明曲面 Ic mat2gray(Ic); % 拉回 0~1避免负值 % 中值滤波先干掉氧化皮造成的椒盐点再用双边滤波保边 Ic medfilt2(Ic, P.medFilt); Ig imbilatfilt(Ic, 0.08, 3); endimbilatfilt的三个参数分别是邻域空间标准差、值域标准差和邻域大小。值域标准差越小保边越狠但噪声也留得越多0.08 是个偏保守的起点实际调的时候看边缘叠加图边缘发毛就往上加一点边缘被磨圆就往下减。mat2gray这一步不能省。减法之后必然有负值im2double不管这个imbinarize遇到负值的行为也不直观先归一化到[0,1]后面所有阈值参数才有可比性。3.2 阈值分割与形态学把锻件从背景里抠出来分割的目标不是「好看」是给后面的边缘搜索划一个可信范围。所以这一步的产物是一个二值掩膜它允许边缘位置粗糙一点但不允许出现大片空洞或粘在画幅边缘的杂块。% segment_part.m —— Otsu 初分割 形态学清理输出单一工件掩膜 function mask segment_part(Ig, P) level graythresh(Ig); % Otsu 自适应阈值 level max(level, P.threshBias); % 兜底阈值过低会把阴影吃进来 mask imbinarize(Ig, level); mask imfill(mask, holes); % 填掉氧化皮剥落造成的小孔洞 mask bwareaopen(mask, P.minArea); % 去掉飞边、碎屑这类小连通域 mask imopen(mask, strel(disk, P.openDisk)); % 开断开细连接 mask imclose(mask, strel(disk, P.closeDisk)); % 闭补边缘缺口 mask imclearborder(mask); % 清掉贴着画幅边缘的干扰 mask bwareafilt(mask, 1); % 只保留最大连通域 工件本体 end顺序有讲究。imfill必须放在bwareaopen前面否则小孔洞会被当成独立连通域先删掉——虽然结果看着差不多但bwareaopen的计算量白费了一部分。imopen放在imclose前面是因为开运算容易在工件腰部切出一道口子闭运算正好补回来反过来做的话闭运算会先把噪声连成块开运算就啃不动了。imclearborder和bwareafilt这两个是「保险丝」。现场拍照很难保证工件完全在画幅内一旦有别的金属件从边角入镜不做这两步掩膜就会变成两个大块后面的外接矩形直接翻倍。注意如果工件本身就必须贴着画幅边缘视野受限的情况imclearborder会连工件一起清掉这时候要注释掉它改成按面积和长宽比筛选。3.3 亚像素边缘为什么要做抛物线插值像素级边缘的精度上限就是 1 个像素换算到实物取决于你的像素当量。假设视野 600 mm 对应 2400 像素一个像素就是 0.25 mm那像素级边缘的量化误差直接就是 ±0.25 mm对于要求 0.1 mm 的场合完全不够看。亚像素的常规做法是在粗边缘附近沿法向取一维灰度剖面对梯度做极值点插值。梯度极值点落在两个像素之间时用三点抛物线拟合求偏移量公式是Δ 0.5*(g1 - g3)/(g1 - 2*g2 g3)g1 g2 g3是相邻三个点的梯度值。这个公式假设局部是二次曲面在边缘过渡带只有两三个像素宽的时候精度最好。% subpix_edge.m —— 沿列扫描灰度剖面三点抛物线插值到亚像素 function [xs, ys] subpix_edge(Ig, yRange, xRange, polarity) xs zeros(0,1); ys zeros(0,1); for x xRange prof Ig(yRange, x); % 取一列灰度剖面 g gradient(prof); % 数值梯度中心差分 if polarity 0 [~, k] max(g); % 由暗到亮找梯度极大 else [~, k] min(g); % 由亮到暗找梯度极小 end if k 1 || k numel(prof) continue; % 极值贴着剖面两端插值不可靠跳过 end g1 g(k-1); g2 g(k); g3 g(k1); den g1 - 2*g2 g3; if abs(den) 1e-9 continue; % 三点共线没有极值 end d 0.5 * (g1 - g3) / den; % 亚像素偏移量纲是像素 xs(end1, 1) x; %#okAGROW ys(end1, 1) yRange(1) k - 1 d; %#okAGROW end endpolarity这个参数表示边缘的明暗方向工件比背景亮就是1工件比背景暗就是-1。它决定你找的是梯度极大还是极小方向搞反的直接后果是边缘整体偏移一个像素左右而且因为偏移方向一致误差不会自己抵消量出来的直径会系统性偏大或偏小。这种系统性偏差比随机噪声讨厌得多因为它看着很「稳定」。yRange限定扫描的行带这个不是可选项。整列全扫的话任何一个反光点都可能抢走极值位置。我的做法是用上一步的掩膜做一次腐蚀取掩膜的上下边界各向内缩 10 个像素作为yRange把搜索范围锁死在工件附近。3.4 标定与像素当量mm/pixel 从哪来所有像素单位的结果都必须乘一个像素当量才能变成毫米。这个系数不能靠「相机分辨率除以视野」硬算镜头畸变、安装高度误差、甚至是工件不在焦平面上都会让实际值偏离理论值。标准做法是拍一个已知尺寸的标定块在同一机位、同一焦距下算出来% calib_px2mm.m —— 用已知尺寸的标定块求像素当量 function k calib_px2mm(calPath, refMM) I imread(calPath); Ig im2double(im2gray(I)); bw imbinarize(Ig, graythresh(Ig)); bw bwareafilt(bw, 1); % 只留标定块本体 st regionprops(bw, BoundingBox); wPix st(1).BoundingBox(3); % 标定块在图像上的像素宽度 k refMM / wPix; % mm / pixel fprintf(像素当量 %.6f mm/px, 标定块像素宽 %.2f\n, k, wPix); end注意BoundingBox给的是[x y width height]宽度的边界是像素的外边界和像素中心之间差半个像素。当wPix有好几百的时候半个像素的误差可以忽略但如果标定块只占几十个像素这个偏差就能到 1%所以标定块要尽量占满视野。P.px2mm留空表示现场标定填死数值表示跳过——有些场合机位是焊死的标一次能用几个月那就把值写进参数文件省一次读图。但每次换镜头、动过相机支架这个值必须重新标这是我见过最多人偷懒的地方。4. 几何量计算外接矩形、圆拟合与圆度误差的代码实现边缘点云拿到了接下来是把点云变成图面上那些带公差的尺寸。锻件最常测的就三类轴向长度、外径、圆度。前两个用最小外接矩形第三个用最小二乘圆代码都不长但每一行的物理意义要清楚。4.1 主轴方向与最小外接矩形长度宽度怎么取锻件在图像里基本不会是水平摆放的直接取BoundingBox得到的宽高会明显偏大——斜着放的长方体其水平外接矩形比实际长度大得多。正确做法是先求主轴方向把工件旋转到水平再取外接矩形。regionprops的Orientation给的是等效椭圆的倾角对细长形工件足够准。旋转用imrotate配合crop旋转后重新取BoundingBox% geom_feature.m 片段主轴对齐后取外接矩形 function [L, W, ang] part_length_width(mask, k) st regionprops(mask, Orientation, Centroid); ang st(1).Orientation; % 主轴倾角单位度 mr imrotate(mask, -ang, bilinear, crop); % 转正 st2 regionprops(mr, BoundingBox); bb st2(1).BoundingBox; L bb(4) * k; % 高度方向 轴向长度 W bb(3) * k; % 宽度方向 径向尺寸 fprintf(长度 %.2f mm, 宽度 %.2f mm, 主轴 %.2f°\n, L, W, ang); end有两个细节。imrotate的角度要取负Orientation是工件相对水平线的倾角转正就得反着转。旋转用bilinear而不是nearest最近邻会让边缘产生阶梯状锯齿外接矩形的边界就跳变一个像素。crop参数会让旋转后的画幅跟原图一样大工件转出边界就被裁掉了。所以对超大工件我一般改成loose再imclearborder一次虽然画幅变大但不会丢数据。提示Orientation对接近正方形的工件会不稳定因为等效椭圆的长短轴差不多长倾角就失去了意义。遇到这种情况改用最小面积外接矩形对轮廓做convhull后旋转卡壳别硬套Orientation。4.2 圆拟合与圆度误差最小二乘圆不是万能的圆度误差的标准定义是同心圆中半径差最小的两圆半径之差。工程上常用最小二乘圆做基准把轮廓点到拟合圆心的距离极差作为圆度误差的近似。这个近似在轮廓接近理想圆的时候很准一旦轮廓有明显的椭圆度或者局部塌陷最小二乘圆会「和稀泥」算出来的圆度偏乐观。% geom_feature.m 片段最小二乘圆拟合 圆度误差 function r fit_circle(mask, k) b bwboundaries(mask, noholes); xy fliplr(b{1}); % bwboundaries 返回 (row,col)翻成 (x,y) % 代数法Kasa拟合x^2 y^2 D*x E*y F 0 A [xy(:,1), xy(:,2), ones(size(xy,1),1)]; rhs -(xy(:,1).^2 xy(:,2).^2); p A \ rhs; % 最小二乘解 cx -p(1)/2; cy -p(2)/2; R sqrt(cx^2 cy^2 - p(3)); % 半径 d hypot(xy(:,1)-cx, xy(:,2)-cy); r struct(cx, cx, cy, cy, R, R*k, ... roundErr, (max(d)-min(d))*k, ... formErr, max(abs(d-R))*k); endrhs取负号是因为把方程整理成了D*x E*y F -(x^2y^2)的形式A\rhs求的是超定方程的最小二乘解。R的表达式来自圆心半径与D E F的关系R sqrt(cx^2 cy^2 - F)代码里-p(3)就是-F。formErr那一项是最大偏差比roundErr更严格——它只允许单向偏差适合评估「只允许材料多、不允许少」的场合。两个都输出判读的时候按图面上的公差要求选。bwboundaries给的是闭轮廓点数可能有几千个A\rhs解一个 3 列的方程组耗时可以忽略。但如果做批量还是建议先xy xy(1:5:end,:)抽稀精度损失极小速度提升明显。4.3 批量跑与结果落盘一张 CSV 交差单张图调通之后真正的价值在批量。目录里几十上百张图跑完直接出一张带文件名的结果表质检那边拿去就能对。% batch_run.m —— 遍历目录逐张测量结果写 CSV function batch_run(inDir, outDir, k) files dir(fullfile(inDir, *.jpg)); T table(); for i 1:numel(files) p fullfile(inDir, files(i).name); P params(); Ig preprocess(p, P); mask segment_part(Ig, P); [L, W] part_length_width(mask, k); r fit_circle(mask, k); % 单张失败不能拖垮整批用 NaN 占位 if isnan(L) || isnan(W) warning(%s 测量失败结果记 NaN, files(i).name); end T [T; table({files(i).name}, L, W, r.R, r.roundErr, ... VariableNames, {file,L_mm,W_mm,D_mm,roundErr_mm})]; end stamp datestr(now, yyyymmdd_HHMM); writetable(T, fullfile(outDir, [result_ stamp .csv])); fprintf(共处理 %d 张结果已写入 CSV\n, height(T)); endparams()放在循环里面是有意的如果哪天你在params.m里加了「根据图像尺寸动态调整bgDisk」的逻辑放外面就永远只有第一张图的参数生效。这个坑很隐蔽表现出来是「第一张准后面全偏」。datestr在 R2022a 之后被标记为不推荐新版本可以换char(datetime(now,Format,yyyyMMdd_HHmm))。老版本继续用datestr也没问题。5. 避坑与排查这套流程最容易翻车的五个地方代码跑起来只是第一步能不能出可信的数取决于你对这几个坑的警觉程度。下面五条都是我在实际调试里反复遇到的按「现象 → 原因 → 解决」写遇到问题可以对着查。5.1 采集与预处理阶段的坑第一坑工件中间出现一个黑窟窿掩膜变成环形。现象是二值图上工件中心区域被判成背景bwareafilt之后只剩一个圈。原因是bgDisk设得太小imopen的结构元装不下工件本体把工件当成照明背景估进去了相减之后工件中心自然被减没了。解决很直接把bgDisk调到大于工件在图像中的最大宽度宁可大一点——结构元过大只是让背景估计变粗糙结构元过小是直接把工件吃掉。第二坑边缘在某一小段被「咬掉」一块长度总是偏短。现象是边缘叠加图上有一段边缘明显内缩不是噪声那种毛刺而是一段平滑的内凹。原因通常是那一小段反光灰度接近背景Canny 在那里找不到梯度。解决办法是预处理阶段改用双边滤波替代高斯或者把cannyThresh的低阈值往下压到 0.04——但要配合掩膜做限定否则噪声边缘会全涌出来。更稳的做法是在边缘点集上做一次 RANSAC 直线拟合把偏离主线的点当外点剔掉那几个被咬掉的点自然就被过滤了。第三坑同一批图前面的准后面的整体偏大。现象是批处理结果里越靠后的图尺寸越大偏差几十微米级别看着像温漂。实际原因是下一批图换了一批料工件离相机更近了像素当量变了。P.px2mm写死在参数文件里机位动过不知道误差就全部转移到尺寸上。解决办法是把标定做成流程的一环每批图至少带一张标定图batch_run开头先跑calib_px2mm跟参数文件里的值比一下偏差超过 0.5% 就告警。5.2 算法与标定阶段的坑第四坑圆度误差算出来是 0 或者接近 0明显不合理。现象是roundErr小数点后三四位才有数而视觉上工件明显是椭的。原因是fit_circle拿到的轮廓点被bwareafilt处理得太干净了轮廓本身就是个完美的圆——这通常意味着你拟合的不是工件外轮廓而是某个形态学操作之后的圆滑边界。检查方法很简单把bwboundaries的点叠加到原图上出一张图看它是不是真的贴着工件边缘。如果叠在掩膜边界而不是原始边缘上就说明前一步的边缘提取没生效几何计算用的是分割结果的边界精度自然只剩像素级。第五坑imclearborder把工件的一部分吃掉了。现象是细长工件的一端没了长度少了一大截。原因就是字面意思工件有部分接触画幅边缘imclearborder分不清是干扰还是主体。这个函数的行为是按连通域整体删除不看你面积大小。解决办法是在调用之前先判断mask(1,:)、mask(end,:)、mask(:,1)、mask(:,end)四条边上是否有大面积非零占比超过 5% 就跳过imclearborder改成用长宽比加面积来筛连通域。这个判断不写换一批视野受限的图就得重新翻车一次。6. 二次开发把它接到自己的测量工位上这套东西真正要落地改的不是算法是接口。我的经验是算法层一行不动在它外面套一层「现场适配层」具体做三件事。第一件是把参数从代码里彻底拔出来。现在的params.m还是 MATLAB 函数改成读json或者ini现场工程师不用开 MATLAB 就能调% load_params.m —— 现场参数外置改 json 不用碰代码 function P load_params(jsonPath) txt fileread(jsonPath); P jsondecode(txt); % 需要 R2016b 以上 def params(); % 代码里的默认值兜底 f fieldnames(def); for i 1:numel(f) if ~isfield(P, f{i}) P.(f{i}) def.(f{i}); % 缺字段就用默认值不致于报错 end end end这个兜底循环很关键。现场配置文件往往是复制上一版的新版本加了字段就会缺不做兜底直接P.newField报错排查起来要翻半天。第二件是换掉硬编码的测量部位。现在part_length_width一把抓总长和总宽实际图面上被测要素可能有三四个台阶直径。可行的做法是把测量部位定义成一份配置表每行一个部位写明边缘的polarity、扫描的yRange、对应的标称值和公差循环跑一遍就出全套结果。字段含义示例name部位名对应图面标注flange_dpolarity边缘明暗方向1 暗到亮1yRange扫描行带限定搜索范围[420 560]nominal标称尺寸mm180.00tol公差带mm0.15第三件是验证。改完任何一处我的习惯是拿三张图做回归一张标定图、一张标准件图、一张最脏的现场图。标定图验证像素当量没漂标准件图验证尺寸在公差内现场图验证鲁棒性没退化。三张都过再上线少一张都不敢放。从那以后我每次动这套脚本的参数哪怕只改一个阈值都强制走一遍这三张图的回归——有一次图省事跳过了结果bgDisk调大之后标准件的孔洞填充逻辑失效跑了一整批才发现所有内径都多算了 0.3 mm返工重测花的时间远超省下的那十分钟。希望这些能帮到你少走一遍我走过的弯路。本文还有配套的精品资源点击获取