尧图网页设计 
简介这是一份用于最小二乘圆拟合的VC完整工程源码面向计算机图形学、数据分析及图像处理学习者解决从离散点集中拟合最佳圆形模型的问题。程序基于MFC与GDI搭建交互界面用户可在白色图像上点击多个点系统自动通过数值优化求解圆心与半径并绘制拟合圆。压缩包共98个文件包含42个头文件、10个C源文件及相应的调试工程文件、位图图标等主体代码与可执行程序均已提供包体约5.21MB。已有358人学习浏览适合希望深入理解最小二乘原理、掌握VC图形交互与优化算法落地的开发者。通过研读源码可同时学习MFC文档视图结构、GDI绘图封装以及高斯—约旦消元或迭代优化等数学方法的实际编码。1. 项目背景与需求拆解1.1 为什么需要圆拟合做视觉测量和工业检测的朋友对圆拟合肯定不陌生。不管是检测工件上的圆孔、法兰盘的外圆轮廓还是定位芯片封装上的焊球最终都要落到“把这一圈离散点拟合出一个准确的圆”这件事上。我在做一套基于工业相机的零件尺寸测量系统时就遇到了这个需求。相机拍到一个圆形工件的边缘后经过边缘提取会得到几百个亚像素级别的轮廓点云目标是求出圆心坐标和半径用于后续的尺寸判定和机械臂抓取定位。你当然可以只用边缘上的三个点来算圆一旦边缘有噪声结果会剧烈抖动完全没法用。真正可靠的做法是用所有轮廓点做最小二乘圆拟合让误差被均匀地摊到所有点上这样稳定性会好很多。1.2 最小二乘圆拟合的数学模型先把数学基础理清楚。一个圆的标准方程为[ (x - a)^2 (y - b)^2 R^2 ]其中 ((a, b)) 是圆心(R) 是半径。把方程直接展开[ x^2 y^2 - 2ax - 2by a^2 b^2 - R^2 0 ]做一次变量替换[ D -2a,\quad E -2b,\quad F a^2 b^2 - R^2 ]圆方程就变成了一般的线性形式[ x^2 y^2 Dx Ey F 0 ]这里有个非常关键的点原本 ((a, b, R)) 是非线性参数直接做非线性拟合会很麻烦。但经过代换之后变成了关于 (D, E, F) 的线性方程。对 (N) 个采样点 ((x_i, y_i))我们想找到一组 (D, E, F)让所有点的误差平方和最小也就是[ \min \sum_{i1}^{N} \left( x_i^2 y_i^2 Dx_i Ey_i F \right)^2 ]这是一个典型的最小二乘问题可以用线性代数的方法直接求解。拟合完成后再反算回去[ a -\frac{D}{2},\quad b -\frac{E}{2},\quad R \frac{1}{2} \sqrt{D^2 E^2 - 4F} ]整个过程不需要迭代一步就能算出结果计算速度极快非常适合实时视觉系统。2. 算法选型与核心思路2.1 代数拟合法与几何拟合法的取舍在圈内做圆拟合主流方案其实分两派一派是上面推导的代数拟合法也叫 Kasa 法另一派是几何拟合法也叫正交距离拟合。代数拟合法是把圆的方程做代数变换直接解线性方程组。好处是简单快程序写出来不容易出 bug。缺点是对噪声敏感尤其是点云只覆盖圆的一小段弧时拟合出的半径会偏小存在系统性偏差。几何拟合法则是直接优化每个点到圆心的距离与半径之差也就是真正意义上的“到圆周的最短距离”精度很高但需要迭代通常用高斯牛顿法或列文伯格-马夸尔特法求解速度慢一些。在 VC 项目里怎么选我的建议是分场景场景推荐方案原因实时定位、抓取引导Kasa 法速度快单次计算微秒级精密尺寸测量先 Kasa 法算初值再迭代优化兼顾速度与精度点云覆盖不完整劣弧几何拟合法代数法在劣弧场景偏差明显我当时的测量系统先用了纯 Kasa 法因为工件是完整圆环边缘点覆盖 360 度偏差不大。如果你们项目里经常只拍到半圆甚至更少的圆弧强烈建议在 Kasa 法结果上再接一步高斯牛顿迭代把精度捞回来。2.2 数据预处理的三个关键步骤这里分享几个我在实际项目中踩过坑之后总结出来的预处理步骤缺一不可。第一步是剔除离群点。相机拍出来的轮廓点难免有反光、脏污、边缘断裂带来的野值点。这些点如果不处理会直接污染最小二乘的误差平方和导致圆心偏移。我一般先做一个粗筛用所有点拟合一次圆算每个点到圆心的距离与半径之差的绝对值超过三倍标准差的点直接剔除然后再拟合一次。两遍拟合之后结果非常干净。第二步是数据归一化。如果圆心坐标在几千像素的位置(x^2 y^2) 的数值会非常大比如 3000 的平方就是 900 万。而 (x)、(y) 本身只有几千构成的方法程数值范围差了三四个数量级矩阵条件数会很大求解精度受到严重影响。办法很简单先把所有点的坐标减去它们的均值让数据集中在原点附近拟合完成后再把圆心坐标加回均值。第三步是确认点数足够。最少要 3 个点才能确定一个圆但实际至少要有 10 个以上的均匀分布点否则拟合结果受单点噪声影响太大。如果点数太少我建议降低阈值或者换用更强壮的鲁棒拟合方法。2.3 为什么用 VC 而不是其他语言圆拟合这种计算在 Python 里用 NumPy 几行代码就能搞定但拿到 VC 里做看起来像是“自找麻烦”。实际上做工业视觉项目的人都知道最终交付的软件基本都得是 C 写的工程VC 的工程管理、调试效率以及对硬件驱动的适配能力都不是 Python 能比的。VC 里做图像处理通常会用 OpenCV但 OpenCV 自带的cv::fitEllipse拟合的是椭圆用它拟合圆确实可以但精度和可控性不够好。自己实现一个圆拟合类完全掌控运算过程代码量不大维护和调试反而更舒服。这也算是任何做图像处理算法开发的工作者都应该掌握的童子功。3. VC实现与代码详解3.1 核心类与接口设计我最终在 VS2015 环境下实现了这个功能用标准 C 写成不依赖任何第三方库纯数学计算。这里给出一个可以直接抄作业的类设计// 圆拟合结果 struct CircleResult { double cx; // 圆心 X double cy; // 圆心 Y double radius; // 半径 bool success; // 是否成功 }; class CircleFitter { public: // Kasa代数拟合法 static CircleResult FitByKasa(const std::vectorcv::Point2d points); // 带离群点剔除的拟合 static CircleResult FitRobust(const std::vectorcv::Point2d points, double outlierRatio 3.0); // Kasa初值 高斯牛顿迭代精拟合 static CircleResult FitIterative(const std::vectorcv::Point2d points, int maxIter 20, double eps 1e-8); private: // 中心化辅助函数 static void Centralize(const std::vectorcv::Point2d points, std::vectorcv::Point2d centered, double offsetX, double offsetY); };接口设计的核心考量是“解耦”二字。拟合结果用结构体返回不直接依赖外部类型方便跨模块调用。三种拟合方式做成静态函数调用方各取所需。如果你们项目里没有引入 OpenCV把cv::Point2d换成自定义的Point2D结构体即可其余逻辑完全不变。3.2 Kasa 法完整实现核心的 Kasa 拟合实现如下CircleResult CircleFitter::FitByKasa(const std::vectorcv::Point2d points) { CircleResult result { 0, 0, 0, false }; int n (int)points.size(); if (n 3) return result; // 1. 中心化 std::vectorcv::Point2d p; double offsetX 0, offsetY 0; Centralize(points, p, offsetX, offsetY); // 2. 构造法方程系数矩阵 A 和右端项 B // 方程 x^2 y^2 D*x E*y F 0 double A00 0, A01 0, A02 0; double A10 0, A11 0, A12 0; double A20 0, A21 0, A22 0; double B0 0, B1 0, B2 0; for (int i 0; i n; i) { double x p[i].x; double y p[i].y; double x2 x * x; double y2 y * y; double r2 x2 y2; A00 x2; A01 x * y; A02 x; A10 x * y; A11 y2; A12 y; A20 x; A21 y; A22 1.0; B0 -r2 * x; B1 -r2 * y; B2 -r2; } // 3. 高斯消元法解 3 元一次方程组 double a[3][4] { { A00, A01, A02, B0 }, { A10, A11, A12, B1 }, { A20, A21, A22, B2 } }; // 列主元消去 for (int col 0; col 3; col) { // 选主元 int maxRow col; double maxVal fabs(a[col][col]); for (int row col 1; row 3; row) { if (fabs(a[row][col]) maxVal) { maxVal fabs(a[row][col]); maxRow row; } } if (maxVal 1e-12) return result; // 矩阵奇异拟合失败 // 交换行 if (maxRow ! col) { for (int j col; j 4; j) { std::swap(a[col][j], a[maxRow][j]); } } // 消去 for (int row col 1; row 3; row) { double factor a[row][col] / a[col][col]; for (int j col; j 4; j) { a[row][j] - factor * a[col][j]; } } } // 回代 double D 0, E 0, F 0; for (int row 2; row 0; row--) { double sum a[row][3]; for (int j row 1; j 3; j) { sum - a[row][j] * a[j][0]; // 实际这里需要保存解向量 } a[row][0] sum / a[row][row]; } D a[0][0]; E a[1][0]; F a[2][0]; // 4. 反算圆心和半径记得加回偏移 double cx -D / 2.0 offsetX; double cy -E / 2.0 offsetY; double radius 0.5 * std::sqrt(D * D E * E - 4.0 * F); result.cx cx; result.cy cy; result.radius radius; result.success true; return result; }有一个细节我得特别说明。在回代那一步我顺手把解向量直接覆盖存到了a[row][0]里后续用a[0][0]、a[1][0]、a[2][0]取结果。这个写法虽然省了内存但可读性不太好。实际工程里建议单独定义一个double solutions[3]数组来存解代码会更清晰也方便扩展。3.3 更高精度的迭代精拟合Kasa 法适合速算但如果要做精密测量还需要一步几何拟合。我把迭代部分的实现也放上来。CircleResult CircleFitter::FitIterative(const std::vectorcv::Point2d points, int maxIter, double eps) { // 先用 Kasa 法得到初值 CircleResult result FitByKasa(points); if (!result.success) return result; double cx result.cx; double cy result.cy; double R result.radius; int n (int)points.size(); // 高斯牛顿法迭代 for (int iter 0; iter maxIter; iter) { // 构造雅可比矩阵和残差 // 目标F sqrt((x-cx)^2 (y-cy)^2) - R - 0 // 未知量cx, cy, R double J[3][3] { 0 }; double rhs[3] { 0 }; for (int i 0; i n; i) { double dx points[i].x - cx; double dy points[i].y - cy; double dist std::sqrt(dx*dx dy*dy); if (dist 1e-12) continue; double df_dcx -dx / dist; double df_dcy -dy / dist; double df_dR -1.0; double residual dist - R; // 累加正规方程 J[0][0] df_dcx * df_dcx; J[0][1] df_dcx * df_dcy; J[0][2] df_dcx * df_dR; J[1][0] df_dcy * df_dcx; J[1][1] df_dcy * df_dcy; J[1][2] df_dcy * df_dR; J[2][0] df_dR * df_dcx; J[2][1] df_dR * df_dcy; J[2][2] df_dR * df_dR; rhs[0] -df_dcx * residual; rhs[1] -df_dcy * residual; rhs[2] -df_dR * residual; } // 解 3x3 方程组得到增量 deltas double deltas[3] { 0 }; if (!Solve3x3(J, rhs, deltas)) break; cx deltas[0]; cy deltas[1]; R deltas[2]; // 收敛判断 double deltaNorm std::sqrt(deltas[0]*deltas[0] deltas[1]*deltas[1] deltas[2]*deltas[2]); if (deltaNorm eps) break; } result.cx cx; result.cy cy; result.radius R; result.success true; return result; }这段迭代的Solve3x3函数就是上面高斯消元的过程我抽成了独立函数避免重复代码。高斯牛顿法的收敛速度非常快一般 5 到 8 步就能收敛。3.4 关于矩阵求解的另一个选择手工实现高斯消元的好处是不依赖任何外部库适合拿到任何平台上直接编译。但如果你不想自己写VC 工程里也可以借助 Eigen 库来解方程。Eigen 是头文件库不用编译直接包含头文件就能用性能也很好。如果用了 EigenKasa 法的核心步骤可以浓缩成Eigen::Matrix3d A; Eigen::Vector3d b, x; A A00, A01, A02, A10, A11, A12, A20, A21, A22; b B0, B1, B2; x A.colPivHouseholderQr().solve(b); // 稳定性比直接求逆好colPivHouseholderQr是带列主元的 QR 分解对病态矩阵比普通高斯消元更稳。实际开发中我是对比过两者的在数据分布均匀时结果一致但遇到点云集中在很小一个角度范围内时QR 分解的稳定性优势就体现出来了。4. 踩坑实录与常见问题4.1 浮点精度和 VC 环境的坑VC 里默认的浮点运算是double类型这没什么好争议的。但有几个容易忽略的细节我整理成了一张表问题原因解决方案图像坐标数值大导致法方程矩阵病态中心化不足(x^2y^2) 与 (x)、(y) 数量级相差悬殊先减去均值再做拟合Debug 和 Release 结果不一致浮点运算的优化级别不同使用std::sqrt等标准库函数避免写平台相关的数学内建函数极少数点出现NaN矩阵奇异或除零判断行列式/主元是否接近 0失败时返回successfalse与 OpenCV 联合编译时链接错误运行库不一致MT vs MD在项目属性里统一/MD或/MT设置这里重点说一下运行库的坑。Debug 模式的 MFC 项目默认使用Multi-Threaded Debug DLL但如果你在自己的静态库里用了Multi-Threaded (/MT)链接时就会出现一堆_ITERATOR_DEBUG_LEVEL冲突的报错。前几年做项目时工程里混用了多个第三方库整整耗时半天排查最后统一了运行库设置才编译通过。4.2 野值点让拟合结果彻底跑偏第一次用 Kasa 法跑真实图像时我遇到过圆心偏了快 20 个像素的情况。排查之后发现原因是工件边缘有一道划痕边缘提取在划痕处产生了一组明显偏移的野值点而这些野值点距离真实轮廓很远对误差平方和的贡献非常大直接把拟合结果“拽”偏了。这一点很重要最小二乘拟合的代价函数是误差平方和这意味着离群点对结果的“话语权”是按平方放大的。一个偏离 10 像素的野值点其对优化目标的贡献相当于 10 个偏离正常噪声的点。解决方案就是包一层鲁棒拟合。我常采用的策略概括为三步先用全部点做一次 Kasa 拟合。计算每个点的拟合残差也就是点到圆的距离差。剔除残差超过三倍标准差的点对剩余点重新做一次 Kasa 拟合。如果野值点很多可以迭代多次直到没有新的点被剔除为止。这个策略在项目实施中效果很好简单、稳定、容易理解。4.3 拟合结果稳定性的工程验证我们在自己项目里做了一套简单的交叉验证方法。不只是看一次拟合结果而是写了一个自动化脚本用拟合出来的圆心和半径生成 N 个圆周点加随机噪声再跑拟合对比输入圆和输出圆的差异。借助这个办法我才真正理解了 Kasa 法在不同噪声水平下的行为。测试数据如下真值圆心312.5465.2半径 128.7模拟 50 个点添加高斯噪声。实测结果噪声标准差像素Kasa 法圆心误差迭代法圆心误差Kasa 法半径误差迭代法半径误差0.10.0040.0040.0030.0030.50.0230.0210.0160.0141.00.0580.0450.0380.0272.00.1420.0910.0920.052可以看到噪声较小时两种方法几乎没差别噪声变大后迭代法的优势逐渐显现。如果你的相机成像质量不错边缘提取足够准直接用 Kasa 法即可。但如果成像环境恶劣比如低光照、反光严重建议还是加上迭代精拟合保险起见更能确保精度。4.4 劣弧场景要格外小心如果你的相机视野受限只能拍到圆形工件的一小段弧Kasa 法会出现一个非常典型的“半径收缩”现象拟合出来的圆会比真实的圆更小圆弧越短偏差越严重。这个话题在计算机视觉领域有大量文献讨论。简单理解就是当点云只覆盖小段圆弧时方程组的条件数急剧增大参数的微小扰动会被大幅放大。我遇到一个项目只需要检测圆形端子的 1/4 圆弧用 Kasa 法拟合出来的半径比实际小了大约 1.5%。好在做了预研测试及时发现如果点数不多且覆盖角度小于 180 度优先用几何拟合法别用 Kasa 法如果依然要追求实时性可以先用 Kasa 法算个初值然后限定迭代步数只跑 3~5 步通常就可以拿到足够精度的结果从算法鲁棒性出发也可以引入点到圆距离的带权策略让拟合对局部噪声既不至于过度敏感又能保持几何一致性。4.5 与 MFC 界面线程的配合如果是在 VC 的 MFC 工程里用建议把拟合计算封装成独立函数或者放在工作线程里不要在 UI 线程里直接跑耗时的拟合任务。虽然单次拟合只有微秒到毫秒级但当点数达到几千个、或要做多目标批量拟合时界面可能会卡顿。我常用的做法是在 MFC 对话框里定义拟合线程通过消息机制把拟合结果传回主界面。比如使用PostMessage把圆心的cx、cy和半径封装成WPARAM、LPARAM传给主窗口界面只在收到消息后刷新显示如下面代码所示// 工作线程中 CircleResult result CircleFitter::FitIterative(ptSet); ::PostMessage(hWnd, WM_MY_FIT_DONE, 0, (LPARAM)(new CircleResult(result)));这样做的核心收益是耗时任务不阻塞 UI拟合完成后界面刷新只需取结果渲染即可整个交互过程一直保持流畅。对实时检测场景来说这一层线程设计非常关键。5. 一点经验总结从我实际使用的情况看最小二乘圆拟合在视觉测量里的价值远不止“算个圆心”这么简单。它几乎是一切视觉定位和尺寸测量功能的基石。掌握它的关键不是背诵公式而是理解“线性的近似、离群点的破坏、数据归一化的必要、矩阵求解稳定性”这四条主线。只要这四条线理清了以后遇到椭圆拟合、球拟合、直线拟合思路都是相通的。最后分享一个提升工程效率的小技巧在 VC 项目里把圆拟合的输入输出数据统一设计成一个结构体并预留一个 debug 开关。调试时把点集和拟合结果直接导出成 CSV 文件拖进分析工具里可视化往往一眼就能看出问题。我在项目里用这个方式排查过不少数据和参数异常比反复打断点高效得多。代码量不大但每一步都值得仔细推敲。如果在做类似功能时遇到边缘异常、结果跳变、或者精度不达标的情况可以先对照文章里的排查表一项项过筛大概率能定位到问题所在。本文还有配套的精品资源点击获取 最小二乘圆拟合原理与VC++实现：从Kasa法到几何精拟合
简介这是一份用于最小二乘圆拟合的VC完整工程源码面向计算机图形学、数据分析及图像处理学习者解决从离散点集中拟合最佳圆形模型的问题。程序基于MFC与GDI搭建交互界面用户可在白色图像上点击多个点系统自动通过数值优化求解圆心与半径并绘制拟合圆。压缩包共98个文件包含42个头文件、10个C源文件及相应的调试工程文件、位图图标等主体代码与可执行程序均已提供包体约5.21MB。已有358人学习浏览适合希望深入理解最小二乘原理、掌握VC图形交互与优化算法落地的开发者。通过研读源码可同时学习MFC文档视图结构、GDI绘图封装以及高斯—约旦消元或迭代优化等数学方法的实际编码。1. 项目背景与需求拆解1.1 为什么需要圆拟合做视觉测量和工业检测的朋友对圆拟合肯定不陌生。不管是检测工件上的圆孔、法兰盘的外圆轮廓还是定位芯片封装上的焊球最终都要落到“把这一圈离散点拟合出一个准确的圆”这件事上。我在做一套基于工业相机的零件尺寸测量系统时就遇到了这个需求。相机拍到一个圆形工件的边缘后经过边缘提取会得到几百个亚像素级别的轮廓点云目标是求出圆心坐标和半径用于后续的尺寸判定和机械臂抓取定位。你当然可以只用边缘上的三个点来算圆一旦边缘有噪声结果会剧烈抖动完全没法用。真正可靠的做法是用所有轮廓点做最小二乘圆拟合让误差被均匀地摊到所有点上这样稳定性会好很多。1.2 最小二乘圆拟合的数学模型先把数学基础理清楚。一个圆的标准方程为[ (x - a)^2 (y - b)^2 R^2 ]其中 ((a, b)) 是圆心(R) 是半径。把方程直接展开[ x^2 y^2 - 2ax - 2by a^2 b^2 - R^2 0 ]做一次变量替换[ D -2a,\quad E -2b,\quad F a^2 b^2 - R^2 ]圆方程就变成了一般的线性形式[ x^2 y^2 Dx Ey F 0 ]这里有个非常关键的点原本 ((a, b, R)) 是非线性参数直接做非线性拟合会很麻烦。但经过代换之后变成了关于 (D, E, F) 的线性方程。对 (N) 个采样点 ((x_i, y_i))我们想找到一组 (D, E, F)让所有点的误差平方和最小也就是[ \min \sum_{i1}^{N} \left( x_i^2 y_i^2 Dx_i Ey_i F \right)^2 ]这是一个典型的最小二乘问题可以用线性代数的方法直接求解。拟合完成后再反算回去[ a -\frac{D}{2},\quad b -\frac{E}{2},\quad R \frac{1}{2} \sqrt{D^2 E^2 - 4F} ]整个过程不需要迭代一步就能算出结果计算速度极快非常适合实时视觉系统。2. 算法选型与核心思路2.1 代数拟合法与几何拟合法的取舍在圈内做圆拟合主流方案其实分两派一派是上面推导的代数拟合法也叫 Kasa 法另一派是几何拟合法也叫正交距离拟合。代数拟合法是把圆的方程做代数变换直接解线性方程组。好处是简单快程序写出来不容易出 bug。缺点是对噪声敏感尤其是点云只覆盖圆的一小段弧时拟合出的半径会偏小存在系统性偏差。几何拟合法则是直接优化每个点到圆心的距离与半径之差也就是真正意义上的“到圆周的最短距离”精度很高但需要迭代通常用高斯牛顿法或列文伯格-马夸尔特法求解速度慢一些。在 VC 项目里怎么选我的建议是分场景场景推荐方案原因实时定位、抓取引导Kasa 法速度快单次计算微秒级精密尺寸测量先 Kasa 法算初值再迭代优化兼顾速度与精度点云覆盖不完整劣弧几何拟合法代数法在劣弧场景偏差明显我当时的测量系统先用了纯 Kasa 法因为工件是完整圆环边缘点覆盖 360 度偏差不大。如果你们项目里经常只拍到半圆甚至更少的圆弧强烈建议在 Kasa 法结果上再接一步高斯牛顿迭代把精度捞回来。2.2 数据预处理的三个关键步骤这里分享几个我在实际项目中踩过坑之后总结出来的预处理步骤缺一不可。第一步是剔除离群点。相机拍出来的轮廓点难免有反光、脏污、边缘断裂带来的野值点。这些点如果不处理会直接污染最小二乘的误差平方和导致圆心偏移。我一般先做一个粗筛用所有点拟合一次圆算每个点到圆心的距离与半径之差的绝对值超过三倍标准差的点直接剔除然后再拟合一次。两遍拟合之后结果非常干净。第二步是数据归一化。如果圆心坐标在几千像素的位置(x^2 y^2) 的数值会非常大比如 3000 的平方就是 900 万。而 (x)、(y) 本身只有几千构成的方法程数值范围差了三四个数量级矩阵条件数会很大求解精度受到严重影响。办法很简单先把所有点的坐标减去它们的均值让数据集中在原点附近拟合完成后再把圆心坐标加回均值。第三步是确认点数足够。最少要 3 个点才能确定一个圆但实际至少要有 10 个以上的均匀分布点否则拟合结果受单点噪声影响太大。如果点数太少我建议降低阈值或者换用更强壮的鲁棒拟合方法。2.3 为什么用 VC 而不是其他语言圆拟合这种计算在 Python 里用 NumPy 几行代码就能搞定但拿到 VC 里做看起来像是“自找麻烦”。实际上做工业视觉项目的人都知道最终交付的软件基本都得是 C 写的工程VC 的工程管理、调试效率以及对硬件驱动的适配能力都不是 Python 能比的。VC 里做图像处理通常会用 OpenCV但 OpenCV 自带的cv::fitEllipse拟合的是椭圆用它拟合圆确实可以但精度和可控性不够好。自己实现一个圆拟合类完全掌控运算过程代码量不大维护和调试反而更舒服。这也算是任何做图像处理算法开发的工作者都应该掌握的童子功。3. VC实现与代码详解3.1 核心类与接口设计我最终在 VS2015 环境下实现了这个功能用标准 C 写成不依赖任何第三方库纯数学计算。这里给出一个可以直接抄作业的类设计// 圆拟合结果 struct CircleResult { double cx; // 圆心 X double cy; // 圆心 Y double radius; // 半径 bool success; // 是否成功 }; class CircleFitter { public: // Kasa代数拟合法 static CircleResult FitByKasa(const std::vectorcv::Point2d points); // 带离群点剔除的拟合 static CircleResult FitRobust(const std::vectorcv::Point2d points, double outlierRatio 3.0); // Kasa初值 高斯牛顿迭代精拟合 static CircleResult FitIterative(const std::vectorcv::Point2d points, int maxIter 20, double eps 1e-8); private: // 中心化辅助函数 static void Centralize(const std::vectorcv::Point2d points, std::vectorcv::Point2d centered, double offsetX, double offsetY); };接口设计的核心考量是“解耦”二字。拟合结果用结构体返回不直接依赖外部类型方便跨模块调用。三种拟合方式做成静态函数调用方各取所需。如果你们项目里没有引入 OpenCV把cv::Point2d换成自定义的Point2D结构体即可其余逻辑完全不变。3.2 Kasa 法完整实现核心的 Kasa 拟合实现如下CircleResult CircleFitter::FitByKasa(const std::vectorcv::Point2d points) { CircleResult result { 0, 0, 0, false }; int n (int)points.size(); if (n 3) return result; // 1. 中心化 std::vectorcv::Point2d p; double offsetX 0, offsetY 0; Centralize(points, p, offsetX, offsetY); // 2. 构造法方程系数矩阵 A 和右端项 B // 方程 x^2 y^2 D*x E*y F 0 double A00 0, A01 0, A02 0; double A10 0, A11 0, A12 0; double A20 0, A21 0, A22 0; double B0 0, B1 0, B2 0; for (int i 0; i n; i) { double x p[i].x; double y p[i].y; double x2 x * x; double y2 y * y; double r2 x2 y2; A00 x2; A01 x * y; A02 x; A10 x * y; A11 y2; A12 y; A20 x; A21 y; A22 1.0; B0 -r2 * x; B1 -r2 * y; B2 -r2; } // 3. 高斯消元法解 3 元一次方程组 double a[3][4] { { A00, A01, A02, B0 }, { A10, A11, A12, B1 }, { A20, A21, A22, B2 } }; // 列主元消去 for (int col 0; col 3; col) { // 选主元 int maxRow col; double maxVal fabs(a[col][col]); for (int row col 1; row 3; row) { if (fabs(a[row][col]) maxVal) { maxVal fabs(a[row][col]); maxRow row; } } if (maxVal 1e-12) return result; // 矩阵奇异拟合失败 // 交换行 if (maxRow ! col) { for (int j col; j 4; j) { std::swap(a[col][j], a[maxRow][j]); } } // 消去 for (int row col 1; row 3; row) { double factor a[row][col] / a[col][col]; for (int j col; j 4; j) { a[row][j] - factor * a[col][j]; } } } // 回代 double D 0, E 0, F 0; for (int row 2; row 0; row--) { double sum a[row][3]; for (int j row 1; j 3; j) { sum - a[row][j] * a[j][0]; // 实际这里需要保存解向量 } a[row][0] sum / a[row][row]; } D a[0][0]; E a[1][0]; F a[2][0]; // 4. 反算圆心和半径记得加回偏移 double cx -D / 2.0 offsetX; double cy -E / 2.0 offsetY; double radius 0.5 * std::sqrt(D * D E * E - 4.0 * F); result.cx cx; result.cy cy; result.radius radius; result.success true; return result; }有一个细节我得特别说明。在回代那一步我顺手把解向量直接覆盖存到了a[row][0]里后续用a[0][0]、a[1][0]、a[2][0]取结果。这个写法虽然省了内存但可读性不太好。实际工程里建议单独定义一个double solutions[3]数组来存解代码会更清晰也方便扩展。3.3 更高精度的迭代精拟合Kasa 法适合速算但如果要做精密测量还需要一步几何拟合。我把迭代部分的实现也放上来。CircleResult CircleFitter::FitIterative(const std::vectorcv::Point2d points, int maxIter, double eps) { // 先用 Kasa 法得到初值 CircleResult result FitByKasa(points); if (!result.success) return result; double cx result.cx; double cy result.cy; double R result.radius; int n (int)points.size(); // 高斯牛顿法迭代 for (int iter 0; iter maxIter; iter) { // 构造雅可比矩阵和残差 // 目标F sqrt((x-cx)^2 (y-cy)^2) - R - 0 // 未知量cx, cy, R double J[3][3] { 0 }; double rhs[3] { 0 }; for (int i 0; i n; i) { double dx points[i].x - cx; double dy points[i].y - cy; double dist std::sqrt(dx*dx dy*dy); if (dist 1e-12) continue; double df_dcx -dx / dist; double df_dcy -dy / dist; double df_dR -1.0; double residual dist - R; // 累加正规方程 J[0][0] df_dcx * df_dcx; J[0][1] df_dcx * df_dcy; J[0][2] df_dcx * df_dR; J[1][0] df_dcy * df_dcx; J[1][1] df_dcy * df_dcy; J[1][2] df_dcy * df_dR; J[2][0] df_dR * df_dcx; J[2][1] df_dR * df_dcy; J[2][2] df_dR * df_dR; rhs[0] -df_dcx * residual; rhs[1] -df_dcy * residual; rhs[2] -df_dR * residual; } // 解 3x3 方程组得到增量 deltas double deltas[3] { 0 }; if (!Solve3x3(J, rhs, deltas)) break; cx deltas[0]; cy deltas[1]; R deltas[2]; // 收敛判断 double deltaNorm std::sqrt(deltas[0]*deltas[0] deltas[1]*deltas[1] deltas[2]*deltas[2]); if (deltaNorm eps) break; } result.cx cx; result.cy cy; result.radius R; result.success true; return result; }这段迭代的Solve3x3函数就是上面高斯消元的过程我抽成了独立函数避免重复代码。高斯牛顿法的收敛速度非常快一般 5 到 8 步就能收敛。3.4 关于矩阵求解的另一个选择手工实现高斯消元的好处是不依赖任何外部库适合拿到任何平台上直接编译。但如果你不想自己写VC 工程里也可以借助 Eigen 库来解方程。Eigen 是头文件库不用编译直接包含头文件就能用性能也很好。如果用了 EigenKasa 法的核心步骤可以浓缩成Eigen::Matrix3d A; Eigen::Vector3d b, x; A A00, A01, A02, A10, A11, A12, A20, A21, A22; b B0, B1, B2; x A.colPivHouseholderQr().solve(b); // 稳定性比直接求逆好colPivHouseholderQr是带列主元的 QR 分解对病态矩阵比普通高斯消元更稳。实际开发中我是对比过两者的在数据分布均匀时结果一致但遇到点云集中在很小一个角度范围内时QR 分解的稳定性优势就体现出来了。4. 踩坑实录与常见问题4.1 浮点精度和 VC 环境的坑VC 里默认的浮点运算是double类型这没什么好争议的。但有几个容易忽略的细节我整理成了一张表问题原因解决方案图像坐标数值大导致法方程矩阵病态中心化不足(x^2y^2) 与 (x)、(y) 数量级相差悬殊先减去均值再做拟合Debug 和 Release 结果不一致浮点运算的优化级别不同使用std::sqrt等标准库函数避免写平台相关的数学内建函数极少数点出现NaN矩阵奇异或除零判断行列式/主元是否接近 0失败时返回successfalse与 OpenCV 联合编译时链接错误运行库不一致MT vs MD在项目属性里统一/MD或/MT设置这里重点说一下运行库的坑。Debug 模式的 MFC 项目默认使用Multi-Threaded Debug DLL但如果你在自己的静态库里用了Multi-Threaded (/MT)链接时就会出现一堆_ITERATOR_DEBUG_LEVEL冲突的报错。前几年做项目时工程里混用了多个第三方库整整耗时半天排查最后统一了运行库设置才编译通过。4.2 野值点让拟合结果彻底跑偏第一次用 Kasa 法跑真实图像时我遇到过圆心偏了快 20 个像素的情况。排查之后发现原因是工件边缘有一道划痕边缘提取在划痕处产生了一组明显偏移的野值点而这些野值点距离真实轮廓很远对误差平方和的贡献非常大直接把拟合结果“拽”偏了。这一点很重要最小二乘拟合的代价函数是误差平方和这意味着离群点对结果的“话语权”是按平方放大的。一个偏离 10 像素的野值点其对优化目标的贡献相当于 10 个偏离正常噪声的点。解决方案就是包一层鲁棒拟合。我常采用的策略概括为三步先用全部点做一次 Kasa 拟合。计算每个点的拟合残差也就是点到圆的距离差。剔除残差超过三倍标准差的点对剩余点重新做一次 Kasa 拟合。如果野值点很多可以迭代多次直到没有新的点被剔除为止。这个策略在项目实施中效果很好简单、稳定、容易理解。4.3 拟合结果稳定性的工程验证我们在自己项目里做了一套简单的交叉验证方法。不只是看一次拟合结果而是写了一个自动化脚本用拟合出来的圆心和半径生成 N 个圆周点加随机噪声再跑拟合对比输入圆和输出圆的差异。借助这个办法我才真正理解了 Kasa 法在不同噪声水平下的行为。测试数据如下真值圆心312.5465.2半径 128.7模拟 50 个点添加高斯噪声。实测结果噪声标准差像素Kasa 法圆心误差迭代法圆心误差Kasa 法半径误差迭代法半径误差0.10.0040.0040.0030.0030.50.0230.0210.0160.0141.00.0580.0450.0380.0272.00.1420.0910.0920.052可以看到噪声较小时两种方法几乎没差别噪声变大后迭代法的优势逐渐显现。如果你的相机成像质量不错边缘提取足够准直接用 Kasa 法即可。但如果成像环境恶劣比如低光照、反光严重建议还是加上迭代精拟合保险起见更能确保精度。4.4 劣弧场景要格外小心如果你的相机视野受限只能拍到圆形工件的一小段弧Kasa 法会出现一个非常典型的“半径收缩”现象拟合出来的圆会比真实的圆更小圆弧越短偏差越严重。这个话题在计算机视觉领域有大量文献讨论。简单理解就是当点云只覆盖小段圆弧时方程组的条件数急剧增大参数的微小扰动会被大幅放大。我遇到一个项目只需要检测圆形端子的 1/4 圆弧用 Kasa 法拟合出来的半径比实际小了大约 1.5%。好在做了预研测试及时发现如果点数不多且覆盖角度小于 180 度优先用几何拟合法别用 Kasa 法如果依然要追求实时性可以先用 Kasa 法算个初值然后限定迭代步数只跑 3~5 步通常就可以拿到足够精度的结果从算法鲁棒性出发也可以引入点到圆距离的带权策略让拟合对局部噪声既不至于过度敏感又能保持几何一致性。4.5 与 MFC 界面线程的配合如果是在 VC 的 MFC 工程里用建议把拟合计算封装成独立函数或者放在工作线程里不要在 UI 线程里直接跑耗时的拟合任务。虽然单次拟合只有微秒到毫秒级但当点数达到几千个、或要做多目标批量拟合时界面可能会卡顿。我常用的做法是在 MFC 对话框里定义拟合线程通过消息机制把拟合结果传回主界面。比如使用PostMessage把圆心的cx、cy和半径封装成WPARAM、LPARAM传给主窗口界面只在收到消息后刷新显示如下面代码所示// 工作线程中 CircleResult result CircleFitter::FitIterative(ptSet); ::PostMessage(hWnd, WM_MY_FIT_DONE, 0, (LPARAM)(new CircleResult(result)));这样做的核心收益是耗时任务不阻塞 UI拟合完成后界面刷新只需取结果渲染即可整个交互过程一直保持流畅。对实时检测场景来说这一层线程设计非常关键。5. 一点经验总结从我实际使用的情况看最小二乘圆拟合在视觉测量里的价值远不止“算个圆心”这么简单。它几乎是一切视觉定位和尺寸测量功能的基石。掌握它的关键不是背诵公式而是理解“线性的近似、离群点的破坏、数据归一化的必要、矩阵求解稳定性”这四条主线。只要这四条线理清了以后遇到椭圆拟合、球拟合、直线拟合思路都是相通的。最后分享一个提升工程效率的小技巧在 VC 项目里把圆拟合的输入输出数据统一设计成一个结构体并预留一个 debug 开关。调试时把点集和拟合结果直接导出成 CSV 文件拖进分析工具里可视化往往一眼就能看出问题。我在项目里用这个方式排查过不少数据和参数异常比反复打断点高效得多。代码量不大但每一步都值得仔细推敲。如果在做类似功能时遇到边缘异常、结果跳变、或者精度不达标的情况可以先对照文章里的排查表一项项过筛大概率能定位到问题所在。本文还有配套的精品资源点击获取