多彩编程 多彩编程MZPH · CODE BLOG
ARTICLE DETAIL

文章详情

深耕前端与后端开发技术的一线实战笔记与踩坑复盘。

Matlab激光光斑定位:二值化与三点定圆拟合方法

Matlab激光光斑定位:二值化与三点定圆拟合方法 简介一份基于Matlab实现激光光斑中心定位与尺寸测量的完整课程设计讲解源自燕山大学电气工程学院课程设计项目适合学习数字图像处理、图像分割与圆拟合的本科生及科研入门者。内容从彩色图像二值化入手分别介绍基于max/min求阈值和graythreshim2bw两种分割路径随后利用bwlabel、regionprops等函数去除噪声通过提取与目标区域具有相同标准二阶中心矩的椭圆参数构造高边数正多边形逼近拟合圆并借助三点求圆心与半径。文档包含设计任务书、原理讲解、Matlab程序实现与二值化前后及去噪前后效果对比同时给出完整课程设计结构可作为相关课设、竞赛或论文的参考资料。整个包仅含1个PDF大小约917KB便于离线阅读。目前已有454人学习下载对光斑图像处理实验和课程设计答辩有务实帮助。1. 激光光斑定位不依赖Hough变换的Matlab拟合思路激光光斑定位最常见的做法是边缘检测加Hough圆变换但真正处理相机采集的激光光斑时这套组合经常翻车中心过曝导致边缘断裂、背景反光带来伪圆、批处理时Hough的参数还要反复调。这次拆解的课程设计走的是一条更稳的路线——先二值化把光斑从背景中切出来用bwlabel和regionprops锁定最大连通域并提取椭圆特征再借助linspace生成正多边形逼近圆最后用圆上任取三点求解圆心和半径。整个方案只依赖Matlab图像处理工具箱的原生函数代码量控制在百行以内不涉及循环拟合适合标定、光斑质心跟踪这类对速度和可复现性有要求的场景。下面按二值化、连通域去噪、椭圆拟合、三点定圆四个环节逐步展开。2. 图像二值化手动阈值与Otsu全局阈值的取舍2.1 为什么要先做二值化激光光斑在传感器上的成像理想情况是亮斑加暗背景实际还叠加环境光、传感器暗电流和随机噪声。如果直接对灰度图做边缘检测阈值和边缘位置的耦合关系会让人反复调参。二值化的本质是把每个像素判定为“目标”或“背景”灰度大于阈值的置逻辑1小于阈值的置0。这一步之后所有几何运算只关心像素的集合属性不再涉及灰度多级值数据量也大幅压缩后续bwlabel和regionprops才能高效工作。激光光斑这种高对比度目标非常适合二值化处理因为光斑灰度通常远高于背景。相比之下生物切片、遥感影像这类目标和背景灰度接近的图像直接二值化会损失大量信息更合适的做法是先做边缘增强或分水岭分割。所以二值化不是万能前置步骤它的适用前提是目标和背景在灰度上有明确的分界。2.2 方法一max/min均值阈值J imread(spot_01.jpg); P rgb2gray(J); % 光斑测量一般先转灰度再处理 ma max(P(:)); % 全图最大灰度 mi min(P(:)); % 全图最小灰度 T (ma mi) / 2; % 阈值取灰度极值的中点 I P T; % 逻辑矩阵1 为目标0 为背景 imshow(I);逻辑说明max和min各只依赖全图一个像素值T直接落在灰度范围的中点等于假设光斑和背景在直方图两侧各占一半。这个假设在背景干净、光斑灰度分布集中的图像上成立效果也直观。但一个亮点噪声或者一个过曝点就会同时抬高ma并可能拉低miT被带偏。我一般只在快速预览时用这个方法定量计算不推荐。2.3 方法二graythresh的Otsu全局阈值J imread(spot_01.jpg); P rgb2gray(J); level graythresh(P); % Otsu按类间方差最大化求归一化阈值 I im2bw(P, level); % level 范围 [0,1]直接作为二值化阈值 imshow(I);逻辑说明graythresh对灰度直方图做全局搜索遍历所有可能的阈值选择让目标和背景两类“类间方差最大、类内方差最小”的那个点。它统计的是整幅图的灰度分布不是单点极值所以抗噪性明显好于方法一。level返回的是归一化到[0,1]的阈值im2bw接收这个值完成分割。新版Matlab里更推荐用imbinarize替代im2bw参数语义保持一致迁移成本很低。2.4 两种方法的对比与失效边界方法阈值来源对噪声的敏感度适用场景max/min中点全图灰度极值的均值高背景干净、光斑灰度均匀graythresh (Otsu)类间方差最大化低光斑与背景灰度明显分离两个方法都是全局阈值光照不均时都会失效。如果光斑图像在暗场拍摄背景接近0Otsu效果很好如果背景有渐变光同一个阈值在不同位置含义不同这时需要分块阈值或先做背景减除。另一个实际经验当光斑占全图面积比例很小时Otsu的阈值会偏向背景一侧可以先用imhist观察直方图形状确认存在明显的双峰再使用。提示课程设计完整程序里直接用了im2bw(I)不传level时默认按0.5切这依赖光斑与背景反差极大的前提反差不足时建议显式传graythresh计算出的level。3. bwlabel连通域标记与regionprops去噪3.1 去噪思路按面积筛而不是按形态筛二值化之后光斑周围经常残留零星亮点来源包括灰尘反光、传感器坏点、光斑边缘的衍射环。这些噪声单点面积小但数量多直接用像素均值会把光斑中心坐标带偏。课程设计里的策略是bwlabel给每个连通域编号regionprops统计面积取面积最大的连通域作为目标其余位置全部置零。这个策略成立的前提是目标面积显著大于噪声面积实际处理激光光斑时这个前提很容易满足因为聚焦后的光斑至少占几十个像素而噪声通常只有一两个像素。如果光斑边缘有破损、和噪声块粘连面积最大区域可能包含噪声区域形态就不准。批量处理前先抽两到三张图观察连通域分布白色像素如果是分散的几十个小块说明阈值取高了背景一大片连通说明阈值取低了。需要做形态学处理时用imopen先开运算断开粘连再继续。3.2 bwlabel连通域编号L bwlabel(I); % 对逻辑图做连通域标记返回与I同尺寸的标签矩阵参数说明L中每个独立连通域对应一个正整数编号从1递增。bwlabel默认按8邻域判断连通即上下左右加四个对角都算“相连”。需要留意的是如果光斑中心过曝形成黑洞外圈会被隔成多个区域这时最大面积区域的面积会明显偏小。常见做法是先用imfill(I, holes)填充空洞再做连通域标记。3.3 regionprops一次取齐后续要用的字段stats regionprops(L, {Area, ConvexHull, MajorAxisLength, ... MinorAxisLength, Eccentricity, Centroid});参数说明regionprops第二个参数用元胞数组一次性请求多个属性比逐个调用省时间。重点看这几个字段的含义属性字段类型含义Area标量区域像素总数Centroid1×2向量质心坐标 (x, y)作为光斑中心初值MajorAxisLength标量与区域同标准二阶中心矩椭圆的长轴长度MinorAxisLength标量同一椭圆的短轴长度Eccentricity标量离心率越接近0说明区域越圆ConvexHullp×2矩阵包含区域的最小凸多边形顶点坐标MajorAxisLength和MinorAxisLength的计算基准是区域像素的二阶中心矩把每个像素当作等质量质点求协方差矩阵主轴方向对应特征向量轴长由特征值换算。所以这里的椭圆描述的是像素分布不是边缘拟合结果。区域内部有缺洞时轴长会偏小这也是上一小节建议先填洞的原因。ConvexHull返回的凸包顶点在本方案里不直接参与定圆但可以用来绘制区域边界或判断区域是否外凸。3.4 提取最大面积区域并剔除噪声A [stats.Area]; % 把所有区域的面积收集成数组 [~, idx] max(A); % 取最大面积的索引即目标区域的标签 I1 I; I1(L ~ idx) 0; % 非目标区域全部置0 figure; imshow(I1);逻辑说明max返回两个值第一个是最大面积第二个是它对应的标签idx。I1先复制二值图再用(L~idx)作为逻辑索引把所有不属于目标区域的像素清零。剩下的就是单一连通域。批量处理时可以在这一步把idx保存下来避免对每张图重复显示。提示老版本bwlabel是兼容性最好的连通域接口但处理上千张大图时内存占用偏高。追求性能可以用bwconncomp替代得到PixelIdxList后按numel排序取最大连通域效果一致速度更快。4. 椭圆参数与正多边形逼近的光斑圆拟合4.1 为什么用二阶矩椭圆描述光斑regionprops的MajorAxisLength、MinorAxisLength基于标准二阶中心矩计算个别坏像素对它的影响远小于对边缘检测的影响。对光斑这个像素集合求协方差矩阵等价于把每个像素当作等质量质点算出两个主轴方向的标准差再换算成等效椭圆的长短轴。只要光斑的像素分布在统计上稳定长短轴就不会因为边缘有个缺口而跳变。这个稳定性正是后续拟合圆所需要的。需要留意的是二阶矩椭圆描述的是像素集合的“离散程度”不是光斑的物理边缘。光斑被图像边界截断时质心统计中心仍然有效但轴长会明显偏小这种情况不能直接用这套流程需要先判断光斑是否完整位于视场内。4.2 linspace生成正多边形并逼近圆圆的半径在所有方向相等而二值化后得到的是一个不规则连通域。需要从椭圆参数推导出“光斑等效圆”。实现方式是用linspace生成等角度间隔的点序列再把每个点画到以质心为中心、以半长轴和半短轴为幅度分量的椭圆轨迹上t linspace(0, 2*pi, 500); % 0到2π取500个点实际构成499条边的正多边形 c1 stats(idx).Centroid; a1 stats(idx).MajorAxisLength; b1 stats(idx).MinorAxisLength; x1 c1(1) a1/2 * cos(t); % x方向幅度用半长轴 y1 c1(2) b1/2 * sin(t); % y方向幅度用半短轴 plot(x1, y1, g-, LineWidth, 1.2);参数说明linspace(0, 2*pi, N)生成N个等间距点首尾不重合因此曲线实际是N-1条线段连接的闭合多边形。N越大多边形与圆的偏差越小。regionprops返回的长短轴是完整轴长画轨迹时要除以2。如果只关心等效圆半径可以直接用r (a1 b1) / 4作为近似。4.3 N的取值对拟合精度的影响N实际边数正多边形与圆的最大偏差典型用途76约6.7% r演示算法思想10099约0.05% r常规显示与粗测500499约0.002% r定量计算推荐偏差来源是弦到弧的最大垂直距离数学表达式为r*(1 - cos(pi/N))。N7时偏差为0.067r视觉上明显有棱角N500时偏差不足0.00002r肉眼和程序都分辨不出与圆的差别。课程设计里为了展示效果贴了N7和N500的对比图定量计算时N取500是稳妥值。继续增大到5000偏差降到0.0000002r量级但三点定圆只取轨迹上的3个点轨迹点数对计算结果的提升已经可以忽略反而增加内存开销。4.4 原始拟合公式的适用边界课程设计原文的轨迹公式是x1 c1(1) d1 * b1 * cos(t); y1 c1(2) d1 * a1 * sin(t);这里d1是Eccentricity离心率。问题在于离心率e sqrt(a^2 - b^2) / a光斑越接近正圆e越趋近0整个轨迹会塌缩成质心附近的一个小点半径不再是光斑的真实尺度。只有当光斑本身呈明显椭圆、e足够大时这个式子才有近似意义。原文能出拟合效果是因为当时的光斑长短轴差异较大不代表这个公式通用。实现等效圆时应该用4.2里的半长轴/半短轴方式或在区域接近正圆时直接用等效半径。5. 三点定圆心与半径不依赖拟合误差的精确求解5.1 圆方程做差消去二次项圆的标准方程(x - x0)^2 (y - y0)^2 r^2展开后x0^2、y0^2、r^2都是常数项。圆上任意三点P2、P3、P4分别代入后两两做差这些常数项全部消掉剩下一个关于x0、y0的二元一次方程组。以P2和P3为例2(x3-x2)x0 2(y3-y2)y0 x3^2 y3^2 - x2^2 - y2^2P3和P4同样处理得到第二条方程。这就是“三点定圆”的数学本质两个未知数两条线性方程。5.2 代码实现与参数说明x2x1(1,1); y2y1(1,1); % 取拟合轨迹上的第1个点 x3x1(1,30); y3y1(1,30); % 第30个点 x4x1(1,80); y4y1(1,80); % 第80个点 a2*(x3-x2); b2*(y3-y2); nx3^2y3^2-x2^2-y2^2; d2*(x4-x3); e2*(y4-y3); fx4^2y4^2-x3^2-y3^2; det_val b*d - e*a; x0 (b*f - e*n) / det_val; % 圆心横坐标 y0 (d*n - a*f) / det_val; % 圆心纵坐标 r0 sqrt((x0-x2)^2 (y0-y2)^2); % 半径参数说明x1是linspace(0, 2*pi, 500)产生的1×500行向量索引1、30、80对应的角度分别约为0度、20.9度、57度。三点角度间隔在20度到36度之间不会接近共线。如果取索引1、2、3三个点挤在一起行列式det_val趋近0x0和y0会漂移出图像范围结果完全不可用。这里没有加eps因为取点策略本身已经保证行列式远离0如果从任意连通域边缘取点则需要用det_val eps防止除零。注意取点索引要避开x1首尾附近两端点处角度步长不完整容易让三点接近共线。通用规则是三个点的角度间隔保持在20度以上。5.3 用全部轨迹点回代验证拟合轨迹有500个点求圆只用其中3个剩下497个点都可以用来检验。回代全部轨迹点计算每个点到圆心的距离看方差是否足够小dist hypot(x1 - x0, y1 - y0); r_mean mean(dist); r_std std(dist); fprintf(R%.2f, std%.3f, std/mean%.2f%%\n, ... r_mean, r_std, r_std/r_mean*100);逻辑说明如果std/mean小于1%说明500个点与圆心的距离高度一致三点求圆的结果可信如果超过5%说明4.2里的轨迹本身不是圆回去检查长短轴幅度的写法不要继续用三点定圆。另外一个更快的思路如果只需要光斑半径这个数值不需要整条拟合曲线可以跳过三点定圆直接取sqrt(MajorAxisLength*MinorAxisLength)/2作为等效半径regionprops返回值单位是像素标定过像素尺寸后可以直接换算物理长度。本文还有配套的精品资源点击获取
返回列表