简介面向激光光斑定位与尺寸测量需求的Matlab讲解PDF适合数字图像处理课程设计或相关项目入门者。文档以燕山大学课程设计为背景系统讲解如何对激光光斑图像进行二值化再利用bwlabel、regionprops等函数去噪并通过获取区域的标准二阶中心矩椭圆参数与最小凸多边形顶点最终拟合出接近圆形区域并求出圆心和半径。资源为1个PDF文件压缩包大小917KB完整呈现了从图像预处理到圆拟合的Matlab实现思路与代码。已有454人学习内容包含彩色图像二值化原理与程序、去除噪声方法、圆拟合原理、完整Matlab程序及课程设计任务书等模块可直接作为数字图像处理与Matlab编程的参考资料。1. 激光光斑中心定位与尺寸测量的Matlab实现思路在激光器装调、光束质量分析和光路对准这类场景里光斑分析的核心就两个数中心坐标和光斑直径。很多人的第一版代码会直接找全图最大值作为中心这在噪声很小的干净画面上勉强能用一旦出现饱和像素、漫反射光晕或者激光器本身的衍射条纹结果就会偏移好几个像素。而光学系统里一个像素往往对应几微米中心偏移加上尺寸测量偏差足以让耦合效率计算结果失效。Matlab做这件事的优势在于它把常用图像处理函数和曲线拟合工具箱放在同一个环境里从读图、预处理到亚像素定位和批量脚本都能一次跑通不需要跨语言处理。这篇文章就按实际做光束质量软件时的流程把光斑中心定位和尺寸测量拆开讲覆盖从算法选型、参数设置到真实光斑排错的完整路径。适合要用CCD截屏图像做光斑判读的设备工程师也适合配合自动对准平台做数据处理控件的上位机开发者。2. 光斑图像预处理从背景扣除到阈值分割的Matlab代码2.1 读取图像与灰度化fspecial、imfilter的基础用法实际拿到的光斑图片通常是相机直接保存的BMP或PNG这里面有一个容易踩的坑相机位深可能不是8位直接读进来以后用imshow看不出问题但数值范围可能只有0到255也可能是0到4095。但无论如何在Matlab里先统一读取再转灰度是第一步。I imread(spot.bmp); if size(I, 3) 3 I rgb2gray(I); end I double(I);转double不是多此一举。后面要做灰度加权质心、做背景统计这些全是浮点运算如果停留在uint8求和过程中超过255就自动截断中心坐标直接算错。显示的时候再转回uint8即可运算过程保持double。接下来是滤波。光斑图片最典型的噪声来源是传感器暗电流和振铃噪声它们表现为孤立亮点或细小颗粒。这个阶段用高斯平滑或中值滤波都可以区别在于中值滤波对脉冲噪声更有效但对边缘保持不好高斯滤波会稍微抹平光斑边缘。我一般倾向用高斯核做轻平滑因为后面要做亚像素拟合平滑带来的边缘偏移可以通过去卷积方式补偿一部分而中值滤波引入的非线性很难建模到估计误差里。kernelSize 5; sigma 1; ker fspecial(gaussian, [kernelSize, kernelSize], sigma); Ismooth imfilter(I, ker, replicate);fspecial(gaussian, [5, 5], 1)生成的核尺寸5像素sigma为1这组参数对大多数光斑图像足够。核尺寸如果取到9以上光斑边缘会被明显拉伸实测直径会偏大约0.5到1个像素。imfilter里的replicate表示边界像素做复制填充避免边缘出现黑边。fspecial在新版本里虽然被官方标注为可能移除但在教学中仍然最直观因为它把核生成和滤波分离方便你检查核的形状。如果追求性能可以替换成imgaussfilt(I, 1)两者的结果在sigma相同的情况下非常接近。2.2 背景扣除与阈值分割graythresh、adaptthresh的取舍光斑图片的背景不是零。传感器有暗电流偏置光学系统有杂散光如果直接对大图做阈值分割背景的凸起有可能被错当成光斑的一部分最后拟合出的sigma会偏大。因此必须先做背景扣除。背景的一个稳健估计是图像四周边缘区域的中值前提是光斑没有占满整个视场。bgRegion I(1:20, :); bg median(bgRegion(:)); Ibg Ismooth - bg; Ibg(Ibg 0) 0;取上边界20行作为背景区中值比均值更能抵抗偶发亮斑干扰。如果光斑位置漂移导致压住上边界这个估计就会失效所以更通用的做法是取四边各取一定宽度取所有边的中值。光斑占满全图的情况后面5.3节再单独讲。背景扣除后进入阈值分割阶段。Matlab里两个典型选项是graythresh配合imbinarize做全局阈值以及adaptthresh做局部自适应阈值。Tg graythresh(uint8(Ibg)); maskG imbinarize(uint8(Ibg), Tg); Ts adaptthresh(uint8(Ibg), 0.45, ForegroundPolarity, bright); maskA imbinarize(uint8(Ibg), Ts);graythresh实现的是Otsu方法它假设图像灰度直方图有双峰在光斑亮区与背景暗区分明时效果稳定。adaptthresh里的第二个参数是灵敏度范围0到10.45表示前景判定比0.5更严格只有相对周围更亮的区域才被保留。灵敏度调得越低保留的前景就越少。局部阈值适合照明不均、背景有渐晕的光斑图而全局阈值在背景干净时更快且不容易把噪声切成细碎区域。下面是两种方法的选择参考阈值方法适用情况常见问题graythresh全局阈值背景均匀、光照稳定、光斑区域明显背景混入渐晕时阈值偏大或偏小adaptthresh局部阈值背景有渐变、照明不均、光斑边缘弱灵敏度参数需要反复试敏感度高手动常数值相机和光源位置固定、口径固定换一次环境就要重新标定分割完的mask常含有孤立的噪点像素或者光斑内部由于过度饱和出现空洞这时要做一次形态学清理maskClean imopen(maskG, strel(disk, 3)); maskClean imfill(maskClean, holes);strel(disk, 3)的半径是光斑有效半径的十分之一以下时不会损伤光斑边缘。imfill填充饱和区域造成的空洞这在中心光强超过相机量程时非常常见。2.3 连通域标记与眩光抑制regionprops的妙用预处理最后一步是把光斑从背景和其他杂物里分离出来。用bwconncomp标记连通域再用regionprops读取属性。cc bwconncomp(maskClean, 8); props regionprops(cc, Ibg, Area, Centroid, BoundingBox, MeanIntensity);第二个参数传图是让regionprops在后续能算出灰度加权属性。这里要把面积筛选写进去否则一两个像素的亮点也可能被当成光斑区域。areas [props.Area]; validIdx find(areas 50 areas size(Ibg,1)*size(Ibg,2)*0.5);50像素的限制用来排除孤立噪点0.5的上限排除全图大亮斑。筛选后再从props(validIdx)取对应区域。眩光是实际测试里最容易影响中心定位的因素。高功率光斑穿过镜头时在传感器上形成一圈环形伪影灰度不高但范围大。质心法对这种环形结构敏感光斑中心会被拉向亮环一侧。一个简单办法是计算每个区域的圆度perims regionprops(cc(maskClean), Perimeter); % 补取周长 circularity 4 * pi * areas ./ (perims^2);规则光斑的圆度接近1眩光形成的环带圆度明显小于0.6。把circularity 0.7作为筛选条件能挡掉大多数低角度入射造成的弧形眩光。这一步做完得到的还是像素级的区域信息真正的中心定位和尺寸测量放到后面章节。3. 光斑中心定位的两种经典算法质心法与高斯拟合3.1 质心法Centroid的Matlab实现与加权方式质心法的理论依据是图像灰度函数的一阶矩公式为x0 sum(x * I(x,y)) / sum(I(x,y)) y0 sum(y * I(x,y)) / sum(I(x,y))这个计算给出的是亚像素位置不需要额外插值。但直接把全图代入会引入背景噪声背景像素虽然灰度低数量庞大累积起来足以把质心拉偏。所以必须对背景做截断处理。maxI max(Ibg(:)); threshold 0.2 * maxI; maskT Ibg threshold; [yy, xx] ndgrid(1:size(Ibg,1), 1:size(Ibg,2)); totalIntensity sum(Ibg(maskT)); cx sum(sum(xx .* Ibg .* maskT)) / totalIntensity; cy sum(sum(yy .* Ibg .* maskT)) / totalIntensity;0.2 * maxI这个阈值是经验值作用是把背景和光斑底部切掉。这里的xx和yy坐标是通过ndgrid生成的注意它不是meshgrid。两者生成的维度方向不同ndgrid第一个输出按行变化第二个按列变化与图像的行列坐标自然对应不容易把xy写反。regionprops里直接有WeightedCentroid属性它的底层就是灰度加权质心。区别在于regionprops是对整个连通域做计算没有先做阈值截断。如果mask本身是从graythresh得到的直接调用WeightedCentroid通常就够了但遇到双峰光斑或强度分布很平的光斑自己加的阈值截断更可控。质心法在光斑近圆形、背景干净时误差能控制在0.05像素以内但对背景非对称、光斑被透镜遮挡形成半圆环等场景会系统性偏移。这是质心法的固有问题换任何加权方式都没法根治只能靠高斯拟合兜底。3.2 二维高斯拟合lsqcurvefit、fitnlm的参数设置高斯拟合的逻辑是把光斑灰度分布建模成二维椭圆高斯函数I(x,y) A * exp( -((x-x0)^2/(2*sx^2) (y-y0)^2/(2*sy^2)) ) B这里的A是峰值幅度x0/y0是中心坐标sx/sy是光斑在x和y方向的标准差B是残留背景。之所以用椭圆高斯而不是圆形高斯是因为激光二极管和光纤输出光斑普遍存在像散圆高斯拟合会把椭圆的长短轴折中中心虽然没有偏差但尺寸参数失真。Matlab实现用优化工具箱的lsqcurvefit它采用最小二乘迭代不需要额外安装统计类工具箱兼容性最好。[rows, cols] size(Ibg); [yy, xx] ndgrid(1:rows, 1:cols); model (p, xy) p(1) * exp( -((xy(:,1)-p(3)).^2 / (2*p(5)^2) (xy(:,2)-p(4)).^2 / (2*p(6)^2)) ) p(7); xyData [xx(:), yy(:)]; data Ibg(:); p0 [max(Ibg(:)), cx, cy, 10, 10, median(Ibg(:))]; lb [0, cx-100, cy-100, 0.5, 0.5, 0]; ub [inf, cx100, cy100, min(rows, cols), min(rows, cols), mean(Ibg(:))]; opts optimoptions(lsqcurvefit, Display, off, MaxFunctionEvaluations, 5000); pfit lsqcurvefit(model, p0, xyData, data, lb, ub, opts);注意model的参数形式第一个输入是参数向量p第二个是自变量坐标矩阵Matlab要求solve函数的回调必须这样排列。p0里的中心位置用了第3.1节质心法的结果这是很关键的联动。高斯拟合是迭代算法初值离真实值太远会收敛到局部极小用质心结果做初值能保证稳定收敛。sx/sy的初值取10像素如果光斑实际明显更大可以改为max(rows,cols)/10。边界条件中lb的0.5限制了高斯宽度不为零ub中min(rows, cols)防止sigma超过图像尺寸。fitnlm方法来自统计和机器学习工具箱代码行数更少p0cell [max(Ibg(:)), cx, cy, 10, 10, median(Ibg(:))]; mdl fitnlm(xyData, data, model, p0cell, Lower, lb, Upper, ub, ... Options, statset(Display, off)); pfit mdl.Coefficients.Estimate;如果两个工具箱都有fitnlm可以顺带输出参数置信区间对写测试报告很有用。但论速度lsqcurvefit略快大数据量时建议用前者。3.3 中心算法的误差对比与适用场景在做批量处理前先把三种中心定位算法的适用边界说清楚算法精度对噪点敏感度对饱和敏感度建议场景最大值像素定位1像素高高粗调对心灰度加权质心0.05像素中中圆对称光斑、实时性要求高二维椭圆高斯拟合0.01像素低高光束质量分析、椭圆光斑、弱背景干扰饱和像素对高斯拟合影响最大。光斑中心一旦饱和灰度顶部被削平拟合算法会试图用更宽的高斯去匹配这个平顶导致sigma偏大、中心偏移可能不大但仍然可达0.3像素。处理方式是在拟合前做一次饱和度检测记录max(I(:)) 255的像素比例超过1%时要么建议重新调曝光要么改用第4.2节的截断法测直径中心仍然用质心法。4. 光斑大小计算直径、束腰与椭圆度的Matlab测量4.1 等效直径与二阶矩的计算激光行业对“光斑大小”的定义不统一常见的是D4σ、FWHM和1/e²。D4σ由ISO 11146定义公式借助广义二阶矩D4σx 4 * sqrt(σx²) σx² sum((x - x̄)² * I(x,y)) / sum(I(x,y))这里的σx²就是灰度能量在x方向上的空间方差D4σ因此也叫二阶矩直径。它跟高斯拟合得到的sigma关系是如果光斑严格高斯基模D4σ 4 * sigma。这是D4σ和拟合width最直接的换算关系。Matlab计算D4σ的代码Ienergy Ibg / sum(Ibg(:)); sigmaX2 sum(sum((xx - cx).^2 .* Ienergy)); sigmaY2 sum(sum((yy - cy).^2 .* Ienergy)); D4x 4 * sqrt(sigmaX2); D4y 4 * sqrt(sigmaY2);Ibg不需要做阈值截断否则二阶矩被低估。但因为背景噪声分布在离中心很远的大范围区域sigmaX2对远处背景敏感所以必须先做背景扣除并且保证背景均值接近0。如果担心边缘噪点的贡献可以先加一个半径约为3倍预估sigma的圆形窗把窗外像素置零。4.2 按峰值强度定义尺寸1/e²的截断法工程上更直观的测法是直接从中心出发沿若干方向找到光强下降到峰值1/e²约13.5%的位置取这些点到中心的平均距离作为半径。angles linspace(-pi, pi, 24); maxI max(Ibg(:)); thresholdI maxI * exp(-2); radius zeros(size(angles)); t 0:0.5:sqrt(rows^2 cols^2); for k 1:length(angles) xs cx t * cos(angles(k)); ys cy t * sin(angles(k)); prof interp2(xx, yy, Ibg, xs, ys, linear, 0); idx find(prof thresholdI, 1, first); if ~isempty(idx) radius(k) t(idx); else radius(k) t(end); end end D 2 * mean(radius);interp2的第7个参数填0表示边界外的插值结果按0处理否则返回NaN。linspace(-pi, pi, 24)生成24条射线已经能覆盖大部分椭圆光斑的形状。T方向步长0.5像素这决定半径测量值的离散误差如果要求更高可以设0.1计算量会上升。1/e²和FWHM的换算关系为D_FWHM D_1e2 / sqrt(2) * sqrt(ln2) ≈ 0.588 * D_1e2这个换算只在高斯光斑时严格成立。如果你的光斑形态偏离高斯较多用D4σ更可靠。4.3 椭圆拟合与倾斜角检查激光器输出光斑往往不是正圆而是带有椭圆度。椭圆的主轴方向对应光束的像散方向这个参数在光纤耦合调整中很关键。常见做法是取掩模的边缘点做特征分解避开第三方工具函数。boundaries bwboundaries(maskClean); B boundaries{1}; bx B(:, 2) - cx; by B(:, 1) - cy; scatterMtx [bx, by]; covMtx scatterMtx * scatterMtx / size(scatterMtx, 1); [V, D] eig(covMtx); theta atan2(V(2, 1), V(1, 1)) * 180 / pi; ellipticity sqrt(D(2,2)) / sqrt(D(1,1));bwboundaries输出的是行列坐标第一列是y第二列是x所以构造点集时要反着写。theta是椭圆长轴相对图像x轴的夹角单位度。ellipticity 1表示长轴在x方向小于1则相反。注意这个特征值法对mask边缘锯齿敏感可以先对边缘点做5点滑动平均平滑再算协方差。5. 真实光斑图像的实战MATLAB批量处理脚本与参数调优5.1 批量处理序列图片的核心脚本把前面所有步骤封装成一个函数输入图片路径输出结构体然后在主脚本里用dir遍历整个文件夹。这样做的意义在于参数一旦调好几百张图不需要再手动干预。function r processSpot(fname, opts) I imread(fname); if size(I, 3) 3 I rgb2gray(I); end I double(I); Ismooth imgaussfilt(I, opts.sigma); bg median(Ismooth(1:20, :), all); Ibg max(Ismooth - bg, 0); mask imbinarize(uint8(Ibg), graythresh(uint8(Ibg))); mask imopen(mask, strel(disk, 3)); mask imfill(mask, holes); cc bwconncomp(mask, 8); props regionprops(cc, Ibg, Area, Centroid); areas [props.Area]; [~, maxIdx] max(areas); mask false(size(Ibg)); mask(cc.PixelIdxList{maxIdx}) true; % 质心初步定位 [yys, xxs] ndgrid(1:size(Ibg,1), 1:size(Ibg,2)); totalI sum(Ibg(mask)); cx0 sum(sum(xxs .* Ibg .* mask)) / totalI; cy0 sum(sum(yys .* Ibg .* mask)) / totalI; % 子窗口高斯拟合 win opts.winWidth; cxMin max(1, round(cx0 - win)); cxMax min(size(Ibg, 2), round(cx0 win)); cyMin max(1, round(cy0 - win)); cyMax min(size(Ibg, 1), round(cy0 win)); crop Ibg(cyMin:cyMax, cxMin:cxMax); [y2, x2] ndgrid(1:size(crop,1), 1:size(crop,2)); xyData [x2(:), y2(:)] [cxMin-1, cyMin-1]; data crop(:); % lsqcurvefit 拟合... r.cx pfit(3); r.cy pfit(4); r.sx pfit(5); r.sy pfit(6); r.D4x 4 * pfit(5); r.D4y 4 * pfit(6); r.theta 0; % 椭圆角度在需要时再算 endopts.winWidth的取值要大于预估光斑半径否则拟合窗口直接把光斑切一半。一般参考值是max(rows, cols) / 10。如果拟合窗口过大把杂散光源圈进来结果同样会偏所以只增加计算量并不会提高精度。主脚本的批量循环files dir(fullfile(data, *.bmp)); opts.sigma 1.2; opts.winWidth 80; resultTable table; for i 1:numel(files) r processSpot(fullfile(data, files(i).name), opts); resultTable(i).fileName files(i).name; resultTable(i).cx r.cx; resultTable(i).cy r.cy; resultTable(i).D4x r.D4x; resultTable(i).D4y r.D4y; end writetable(struct2table(resultTable), result.csv);dir在Windows下不会带路径拼路径时用fullfile保证跨平台。struct2table把结构体数组直接转成table正好能被writetable写成CSV。5.2 实战中常见的三个误用信号与排错第一个现象是中心坐标在图像序列里来回跳变幅度超过0.5像素。原因多数是阈值把光斑边缘切成碎块导致最大连通域交替选取。处理办法是保留面积最大的连通域并且把mask重建时的操作从bwconncomp改用bwareafilt(mask, 1, largest)一步到位。第二个现象是高斯拟合得到的sigma明显小于从图像上肉眼估计的半径。这通常是光斑过饱和顶部被削平后拟合算法只能收缩sigma去匹配陡峭的腰。此时看max(I(:))是否等于255饱和面积比例是否超过1%。如果超了改用D4σ方法或者直接调低相机增益重拍。第三个现象是质心总是偏向某个高亮灰尘点。灰尘在CCD上形成的亮斑通常只有几个像素但灰度极高对加权质心的影响大。排错时先看最大强度像素是否总在固定位置如果是则在连通域筛选阶段加入MeanIntensity检查灰尘点区域面积小但平均强度高可以按面积与峰值的比值过滤。下面是参数调节的最小影响表改参数前先记录这几个状态现象优先检查项建议修正中心抖动掩模最大连通域是否稳定用bwareafilt只保留最大区域直径偏大阈值分割是否把背景照亮提高阈值截断比例直径偏小饱和像素比例改用D4σ或降低曝光中心偏移视野边缘是否有眩光增大滤波核或缩小ROI窗口5.3 批量脚本的参数调优建议参数不建议用全局变量保存因为调参过程会来回覆盖最终无法复现。较稳妥的做法是定义一个opts结构体每次调完参数把opts连同结果表格一起保存为mat文件这样后面发现问题时能追溯当时用了哪组参数。批量执行的耗时主要集中在高斯拟合的迭代上。两千张512x512灰度图每张都跑全图拟合大约需要半小时。如果只需要中心位置可以先只跑质心法只对质心偏移超过阈值的图像做高斯拟合复核这样用20%的计算量覆盖90%的场景。Matlab R2024a等新版本中lsqcurvefit会自动使用多线程如果处理器核心数不多可以设置opts.UseParallel false避免频繁创建线程池带来的额外开销。6. 验证与进阶仿真图像自测与亚像素精度提升6.1 用已知参数的仿真光斑做自测真实图像永远没有标准答案算法对不对只能靠仿真验证。生成一张已知中心、尺寸和噪声的高斯光斑然后把它输入给第5节的processSpot函数对比算法输出与真实值之间的偏差。[X, Y] ndgrid(1:512, 1:512); trueCx 245.3; trueCy 268.7; trueSx 18.5; trueSy 15.2; I_gt 200 * exp(-(((X - trueCx).^2) / (2 * trueSx^2) ... ((Y - trueCy).^2) / (2 * trueSy^2))) 5 2 * randn(512); imwrite(uint8(I_gt), sim_spot.bmp);仿真时把背景设成5噪声标准差2这组参数的信号噪声比和真实相机相当。处理完以后统计偏差如果中心误差超过0.05像素说明质心初值或者拟合下限设置有问题优先检查背景扣除是否把真实光斑底部切了一部分。这个测试脚本每次改完算法都应该重跑一遍并且在算法参数改动后保留结果文件。调参时你会发现有时候把winWidth减少20像素中心误差反而降低原因是窗口外的微弱杂散光被排除在拟合之外干扰减少。6.2 亚像素精度的增益手段质心法本身具备亚像素精度但受限于离散采样对光斑的对称性和灰度分布要求严格。抛物线插值是一个互补手段对沿一维剖面搜索峰值的位置很有效[peakVal, idx] max(profile); if idx 1 idx length(profile) x1 idx - 1; x2 idx; x3 idx 1; delta 0.5 * (profile(x1) - profile(x3)) / ... (profile(x1) - 2 * peakVal profile(x3)); subPixelIdx x2 delta; else subPixelIdx idx; end分母是标准抛物线的二阶差分正常情况下不会为零但光斑被严重削平使三点灰度相等时会出现除零因此要保留idx1的判断路径。这个方法只适合中心附近灰度变化近似抛物线的情况用在光斑中心定位时通常是先取穿过峰值围的一维灰度剖线再做抛物线插值最终中心取x和y两个方向的插值结果。精度在0.2像素左右不如高斯拟合优点是速度极快适合实时粗定位。6.3 把ROI粗剪到光斑附近再拟合多数真实光斑图像的有效区域只占全图25%以下。把整幅图送入高斯拟合有两个坏处远处噪声进入拟合域压低sigma估计以及求解时间变长。常规做法是先粗定位再裁剪从第3.1节的质心粗结果出发裁剪边长为8倍粗估sigma的正方形区域。sigmaGuess max(sx, sy); % 从前面任何方法得到的粗估 win max(ceil(4 * sigmaGuess), 40); cxSpot round(cx0); cySpot round(cy0); roi Ibg(max(1, cySpot - win):min(rows, cySpot win), ... max(1, cxSpot - win):min(cols, cxSpot win));裁剪后重新跑lsqcurvefit时拟合的坐标原点要加上裁剪偏移否则出的中心坐标是窗口内的相对坐标第5.1节里的代码已经通过给xyData加偏移处理了这一点。实际操作中这个裁剪窗口直接决定拟合稳定性窗口过大会把光斑旁边的二级衍射环圈进来窗口过小则高斯模型的边缘被截断sigma被低估。每次调整之后可以同时记录fitquality残差model(pfit, xyData) - data的RMS残差突然变大就要回看ROI边界是否切到有效信号。把上面这套仿真加裁剪的检查流程固化到批处理脚本开头每次跑真实数据前自动执行一次输出的结果表格里如果没有明显的中心突变这批数据就可以直接用了。本文还有配套的精品资源点击获取