做视觉检测、激光点云处理或者机器人导航这类活儿的朋友大概率都遇到过同一个需求一堆散乱的点里明明藏着一个圆可你用最小二乘去拟合结果被几个离群点带偏中心偏了半径也不对。这种场景最适合的方案就是 RANSAC。这篇文章把我自己用 C 实现 RANSAC 平面圆拟合的完整过程过一遍从数学原理到能跑的代码再到参数调节和踩过的坑全部写清楚希望能给正在做类似功能的同学省点时间。代码我会拆成小块讲每个函数都尽量说清楚为什么这么写你拿去改改也能直接用。1. 为什么最小二乘搞不定的场景RANSAC 反而能扛住1.1 离群点对最小二乘拟合的影响先看一个很典型的例子。假设你在图像里提取了 200 个边缘点其中 130 个点真的落在某个圆上剩下 70 个点是反光、遮挡、误检带来的杂散点。用传统最小二乘去拟合圆算法会把所有点的误差平方加起来然后求一组圆心和半径让这个总和最小。问题是离群点的误差通常很大哪怕只占一小部分它们的平方项也会把拟合结果狠狠拽走。最后拟合出来的圆可能半径偏大、中心偏移甚至完全偏离真实圆弧。这个现象说到底是因为最小二乘的损失函数是二次的对“大误差”极度敏感。你不是不知道圆在哪里是被那几个坏点带乱了节奏。RANSAC 的思路完全不同它不追求一次拟合所有点而是反复随机抽取最小数量的点去生成候选模型然后用“有多少点支持这个模型”来打分。这种方式天然对离群点不敏感只要真实圆上的点占比别太低它就能把正确的圆捞出来。1.2 RANSAC 的核心假设与适用边界RANSAC 的全称是 Random Sample Consensus核心假设就两条一是数据里存在一群“内点”能用一个模型统一描述二是剩下的“外点”跟模型没关系。对圆拟合来说模型就是圆心加半径内点就是到圆距离小于阈值的那些点外点就是其他杂散点。它的通用流程是随机取三个点拟合一个圆统计所有点到这个圆的距离把距离小于阈值的点记为内点保留内点数最多的那次采样作为最佳模型最后再用全部内点重新做一次精确拟合。适用范围我总结过如果内点占比低于 50%迭代次数会涨得很快但仍能跑如果内点占比低于 10%即使算法理论能跑实际也很容易收敛到错误模型最好换思路。另外如果数据里同时存在两个不同的圆RANSAC 只会锁定其中支撑点更多的一个想找另一个就得先把已找到圆的点删掉再跑一次。这些都是使用前要想清楚的边界条件。2. 平面圆拟合的数学底子从三点定圆到代数最小二乘2.1 三点定圆公式与退化判断三个不共线的点能唯一确定一个圆这个大家都知道。计算圆心坐标最常用的公式是求三角形外接圆设三点为 A、B、C圆心 O 满足 OAOBOC。解这个方程可以得到分母 D 2 * (Ax(By - Cy) Bx(Cy - Ay) Cx(Ay - By))圆心 x ( (Ax² Ay²)(By - Cy) (Bx² By²)(Cy - Ay) (Cx² Cy²)(Ay - By) ) / D圆心 y ( (Ax² Ay²)(Cx - Bx) (Bx² By²)(Ax - Cx) (Cx² Cy²)(Bx - Ax) ) / DD 的绝对值如果接近 0说明三点近似共线或者距离非常近这时候圆的半径会趋于无穷大数值上很不稳定。这个退化判断必须在代码里做掉不然你会在后面收到一个中心在十万八千里外的“圆”。2.2 Kasa 代数拟合法与正规方程三点定圆只是 RANSAC 里生成假设用的最终把所有内点统一拟合成一个更精细的圆时还需要支持任意数量点的拟合函数。最直接的是 Kasa 方法也叫代数最小二乘。它把圆方程改写成x² y² ax by c 0其中圆心是 (-a/2, -b/2)半径满足 r² (a² b²) / 4 - c。我们不去最小化点到圆的真实几何距离而是最小化每个点在代数方程里留下的残差也就是让 x² y² ax by c 尽量接近 0。这个问题是线性的构造正规方程组解一下就行矩阵未知量右端Σx² Σxy Σxa-Σ(x³ xy²)Σxy Σy² Σyb-Σ(x²y y³)Σx Σy nc-Σ(x² y²)如果你了解过平面拟合或者多项式拟合会发现套路一样把二次项当成已知量把一次项和常数项当成未知量用最小二乘解线性系统。这里的 n 是点数。3x3 的线性方程组直接高斯消元就行不需要引入 Eigen 这种重依赖代码我放在后面。2.3 为什么还要做几何距离精修代数拟合法虽然快但有个天生的毛病它对小半径圆更“偏心”拟合出的半径通常偏小一点。原因是它最小化的是代数残差而不是真实的点到圆距离。当内点噪声比较小、圆弧覆盖完整时这个偏差可以忽略但如果你做的是高精度测量或者数据只覆盖了一小段圆弧代数拟合的误差会明显变大。我常用的补救办法是在代数拟合结果上再做几轮 Gauss-Newton 迭代直接最小化目标函数 Σ(点到圆心距离 - r)²。这是个非线性最小化问题但一次迭代只需要解一个 3x3 正规方程代价非常小。用在 RANSAC 的内点精修阶段精度提升很直观。3. C 实现 RANSAC 圆拟合的每个关键步骤3.1 基础数据结构和距离函数先定义点、圆和距离计算。这里有个容易被忽略的小点算点到圆心距离要用std::hypot它的数值稳定性比sqrt(dx*dx dy*dy)好尤其当坐标值很大、差值很小的时候能少丢精度。#include algorithm #include cmath #include cstdio #include random #include vector struct Point2D { double x; double y; }; struct Circle2D { Point2D center; double radius; }; double PointToCircleDistance(const Point2D pt, const Circle2D circle) { return std::abs(std::hypot(pt.x - circle.center.x, pt.y - circle.center.y) - circle.radius); }PointToCircleDistance是后面所有评分的基础它返回点到圆周的绝对距离小于阈值就算内点。这个函数的性能也会直接影响 RANSAC 总耗时尽量保持精简。3.2 求解器3x3 线性方程组的高斯消元Kasa 拟合、Gauss-Newton 精修、三点定圆都不需要复杂的矩阵库一个通用的小求解器就够了。高斯消元加列主元代码不长但很靠谱bool SolveLinear3x3(const double A[3][3], const double b[3], double x[3]) { double M[3][4]; for (int i 0; i 3; i) { for (int j 0; j 3; j) M[i][j] A[i][j]; M[i][3] b[i]; } for (int col 0; col 3; col) { int pivot col; for (int r col 1; r 3; r) { if (std::abs(M[r][col]) std::abs(M[pivot][col])) pivot r; } if (std::abs(M[pivot][col]) 1e-12) return false; if (pivot ! col) { for (int j 0; j 4; j) std::swap(M[pivot][j], M[col][j]); } double div M[col][col]; for (int j 0; j 4; j) M[col][j] / div; for (int r 0; r 3; r) { if (r col) continue; double factor M[r][col]; for (int j 0; j 4; j) M[r][j] - factor * M[col][j]; } } for (int i 0; i 3; i) x[i] M[i][3]; return true; }选主元这步别省。圆拟合里如果近似共线的三点凑在一起矩阵很容易病态不做列主元解出来的结果会飞。3.3 三点定圆与代数拟合的具体代码三点定圆按照前面公式写bool FitCircle3Points(const Point2D p1, const Point2D p2, const Point2D p3, Circle2D circle) { const double ax p1.x, ay p1.y; const double bx p2.x, by p2.y; const double cx p3.x, cy p3.y; const double d 2.0 * (ax * (by - cy) bx * (cy - ay) cx * (ay - by)); if (std::abs(d) 1e-12) return false; const double a2 ax * ax ay * ay; const double b2 bx * bx by * by; const double c2 cx * cx cy * cy; const double ux (a2 * (by - cy) b2 * (cy - ay) c2 * (ay - by)) / d; const double uy (a2 * (cx - bx) b2 * (ax - cx) c2 * (bx - ax)) / d; circle.center.x ux; circle.center.y uy; circle.radius std::hypot(ax - ux, ay - uy); return true; }Kasa 代数拟合bool FitCircleAlgebraic(const std::vectorPoint2D pts, Circle2D circle) { if (pts.size() 3) return false; double sx 0, sy 0; double sxx 0, syy 0, sxy 0; double sx3 0, sy3 0, sx2y 0, sxy2 0; for (const auto p : pts) { double x p.x, y p.y; double x2 x * x, y2 y * y; sx x; sy y; sxx x2; syy y2; sxy x * y; sx3 x2 * x; sy3 y2 * y; sx2y x2 * y; sxy2 x * y2; } const double n static_castdouble(pts.size()); double A[3][3] { { sxx, sxy, sx }, { sxy, syy, sy }, { sx, sy, n } }; double b[3] { -(sx3 sxy2), -(sx2y sy3), -(sxx syy) }; double x[3]; if (!SolveLinear3x3(A, b, x)) return false; const double a x[0]; const double bcoef x[1]; const double c x[2]; double rr (a * a bcoef * bcoef) / 4.0 - c; if (rr 0) return false; circle.center.x -a / 2.0; circle.center.y -bcoef / 2.0; circle.radius std::sqrt(rr); return true; }注意这里右端项的符号别搞反。我一开始实现时把负号漏了拟合出来的圆心位置偏得一塌糊涂。如果你看到结果整体平移优先查这一步。3.4 RANSAC 主循环采样、评分、自适应终止主循环的逻辑可以浓缩成四个动作随机抽三个点、三点定圆、数内点、更新最优。我加了自适应迭代次数当找到足够好的模型后会按公式动态减小剩余迭代次数不用傻傻跑满初始值。Circle2D RansacCircleFit(const std::vectorPoint2D points, double threshold 0.1, int maxIterations 2000, int minInliers 3, unsigned int seed 12345) { const int n static_castint(points.size()); Circle2D best; int bestCount 0; if (n 3) return best; std::mt19937 rng(seed); std::uniform_int_distributionint dist(0, n - 1); int iterations maxIterations; for (int iter 0; iter iterations; iter) { int idx[3]; idx[0] dist(rng); do { idx[1] dist(rng); } while (idx[1] idx[0]); do { idx[2] dist(rng); } while (idx[2] idx[0] || idx[2] idx[1]); Circle2D candidate; if (!FitCircle3Points(points[idx[0]], points[idx[1]], points[idx[2]], candidate)) { continue; } int count 0; for (const auto p : points) { if (PointToCircleDistance(p, candidate) threshold) count; } if (count bestCount) { bestCount count; best candidate; double inlierRatio static_castdouble(count) / n; double outlierRatio 1.0 - inlierRatio; if (outlierRatio 1.0) { double p_success 0.99; int needed static_castint( std::log(1.0 - p_success) / std::log(1.0 - std::pow(1.0 - outlierRatio, 3))) 1; needed std::max(1, needed); if (needed iterations) iterations needed; } } } if (bestCount minInliers) return Circle2D{}; std::vectorPoint2D inliers; inliers.reserve(bestCount); // 两段式精修收集内点、代数拟合、几何精修再重复一次 for (int pass 0; pass 2; pass) { inliers.clear(); for (const auto p : points) { if (PointToCircleDistance(p, best) threshold) { inliers.push_back(p); } } if (static_castint(inliers.size()) 3) break; Circle2D refined; if (FitCircleAlgebraic(inliers, refined)) { RefineCircleGeometric(inliers, refined); best refined; } } return best; }为什么分两步精修第一步先用代数拟合把全体内点综合起来得到一个比三点定圆更稳的初值第二步再用几何距离精修消除代数拟合的半径偏小问题。重复两轮是因为精修后模型会微调内点集合可能略微变化再收集一次可以保证用上最完整的内点集。Gauss-Newton 几何精修实现如下void RefineCircleGeometric(const std::vectorPoint2D pts, Circle2D circle, int maxIter 10) { for (int iter 0; iter maxIter; iter) { double J11 0, J12 0, J13 0; double J22 0, J23 0, J33 0; double r1 0, r2 0, r3 0; for (const auto p : pts) { double dx p.x - circle.center.x; double dy p.y - circle.center.y; double dist std::hypot(dx, dy); if (dist 1e-12) continue; double residual dist - circle.radius; double d_cx -dx / dist; double d_cy -dy / dist; double d_r -1.0; J11 d_cx * d_cx; J12 d_cx * d_cy; J13 d_cx * d_r; J22 d_cy * d_cy; J23 d_cy * d_r; J33 d_r * d_r; r1 d_cx * residual; r2 d_cy * residual; r3 d_r * residual; } double A[3][3] { {J11, J12, J13}, {J12, J22, J23}, {J13, J23, J33} }; double b[3] { -r1, -r2, -r3 }; double delta[3]; if (!SolveLinear3x3(A, b, delta)) break; circle.center.x delta[0]; circle.center.y delta[1]; circle.radius delta[2]; double step std::sqrt(delta[0] * delta[0] delta[1] * delta[1] delta[2] * delta[2]); if (step 1e-10) break; } }这里我直接累加 J^T J 矩阵和 J^T r 向量避免显式构建雅可比矩阵再乘一遍C 里这样写更直观也更快。4. 可运行的完整 Demo 与参数调优4.1 生成带噪声和离群点的测试数据工具写好了总得验证一下。我用程序生成 200 个点真实圆是圆心 (0,0)、半径 5.0一部分点沿圆周分布并加高斯噪声另一部分点随机撒在周围矩形区域里模拟离群点。然后分别用普通最小二乘和 RANSAC 拟合对比结果。int main() { const int kNumPoints 200; const double kNoise 0.08; const double kOutlierRatio 0.4; std::mt19937 rng(42); std::normal_distributiondouble gauss(0.0, kNoise); std::uniform_real_distributiondouble unif(-8.0, 8.0); std::vectorPoint2D points; points.reserve(kNumPoints); for (int i 0; i kNumPoints; i) { double angle 2.0 * M_PI * i / kNumPoints; double x 5.0 * std::cos(angle) gauss(rng); double y 5.0 * std::sin(angle) gauss(rng); if (i static_castint(kNumPoints * kOutlierRatio)) { x unif(rng); y unif(rng); } points.push_back({x, y}); } Circle2D naive; FitCircleAlgebraic(points, naive); Circle2D result RansacCircleFit(points, 0.2, 5000, 20, 12345); std::printf(true : center(0,0) r5.0\n); std::printf(ls : center(%7.3f,%7.3f) r%.3f\n, naive.center.x, naive.center.y, naive.radius); std::printf(ransac : center(%7.3f,%7.3f) r%.3f\n, result.center.x, result.center.y, result.radius); return 0; }在我本机跑的结果基本是这样普通最小二乘被离群点拽得中心偏移超过两个单位半径直接跑到 11 以上RANSAC 恢复出的圆误差在零点几个像素以内。这个对比每次演示都很有说服力。4.2 阈值和迭代数怎么定一张经验表RANSAC 最影响结果的参数就是内点距离阈值 threshold。阈值太大离群点会被算成内点拟合精度下降阈值太小带噪声的内点数量减少模型可能找不到足够支撑。我一般用数据噪声标准差 sigma 的 2~3 倍作为初值。如果不知道 sigma可以先跑一次统计所有点与最优圆的距离分布再看 90% 的点集中在哪个区间。下面是不同场景的经验参考数据来源典型噪声阈值推荐图像像素坐标0.5~1 px1.5~3 px激光雷达点云0.01~0.05 m0.03~0.10 m机械接触测量0.005~0.01 mm0.01~0.03 mm迭代数可以通过公式算设内点率为 w三点全为内点的概率是 w³那么 N 次采样里至少有一次成功的概率 p 满足 p 1 - (1 - w³)^N。反过来就能解出 N。比如内点率 65%p 取 0.99N 大约只需要 15 次内点率降到 30%N 就要到 60 次左右内点率 10% 时 N 要超过 4000 次。这也是我代码里做自适应迭代的原因。maxIterations 设 2000 到 5000 是大多数二维场景的稳妥选择。如果数据已经预处理得很干净500 次也够如果离群率很高又对实时性有要求建议先做粗滤波降低离群率而不是无脑加大迭代数。4.3 从运行时间看 RANSAC 的工程取舍RANSAC 每轮都要把所有点扫一遍算距离数据量上到几十万点时5000 次迭代就会变成性能瓶颈。我在实际项目里做过两件事来提速。第一分块采样。如果点云是组织化的比如按行扫描得到的深度图可以先在若干局部区域随机选点避免每次都从全局抽但这要求你对数据分布有先验。第二提前终止加内存点剪枝。在数内点时计数器一旦超过当前最优值就立刻 break不用统计完整轮。这个优化在点很多且最优值明显时能省不少时间。另一个工程技巧是固定随机种子。研发阶段把 seed 固定住结果可复现方便对比参数上线阶段再用std::random_device生成种子让每次运行不依赖固定模式。这个看似不起眼却能避免很多“同一个输入两次结果不同”的排查时间。5. 工程化实操中会遇到的问题5.1 半径偏小、中心偏移是代数拟合的锅我第一次用纯代数最小二乘做精修时半径总是比真实值小那么零点几个百分点。开始以为是 RANSAC 阈值问题后来单独对一组干净内点做拟合才发现是 Kasa 方法的固有偏差。解决办法就是我前面加的RefineCircleGeometric。这里提醒一句几何精修也不是万能的。如果圆弧覆盖太少比如只有 20 度的短弧圆心在垂直于弧的方向上约束极弱几何精修可能迭代不收敛。面对短弧数据我更推荐在 RANSAC 验证阶段就要求采样点之间保持一定角度间隔避免三点聚在一小段弧上。5.2 三点共线导致的退化三点共线或距离极近时三点定圆公式的分母接近 0拟合出的圆半径会变得巨大。这种情况在真实数据里不少见尤其是点比较稀疏、一个局部区域只能采到三四个点的时候。处理办法有两层一是在FitCircle3Points里检查分母绝对值小于阈值直接返回失败二是在内点验证阶段加半径合理性检查超出预期的半径范围直接丢弃。我之前有个项目里的圆弧半径大约 20 像素结果 RANSAC 返回了一个半径 5000 像素的圆内点数量还不少。排查到最后就是采样到了近似共线的三个点加了一个“半径必须在 [minR, maxR] 之间”的约束后问题立刻消失。5.3 随机种子与结果复现RANSAC 是随机算法不固定种子时每次拟合结果可能会有微小波动。这个波动在低离群率下几乎看不出来但离群率一高就可能出现不同的局部最优。建议开发阶段固定种子调参时把算法里的 threshold、minInliers、maxIterations 逐个扫描你才能真正对比出哪个参数组合更稳。另外一个容易被忽略的问题是std::uniform_int_distribution每次构造也会消耗一点时间。循环外只构造一次分布对象循环内直接调用这个细节在高频调用时能省不少。5.4 常见问题速查表现象常见原因解决办法拟合半径明显偏小代数拟合固有偏差做几何距离精修每次结果不同未固定随机种子开发期固定 seed拟合出巨大半径采样到近似共线三点检查分母、限制半径范围内点数始终很少阈值太小或离群率过高调大阈值迭代数自适应只有一个短弧却拟合失败弧上内点约束弱限制采样点角度间隔大点云跑得慢每轮全量扫描提前终止、分块采样我还想补充一个习惯写完 RANSAC 之后建议把“最优那三点的索引、对应圆参数、内点数”打日志保存下来。真实项目里数据复杂没有这种中间信息出问题你会非常难定位。6. 完整代码清单前面所有片段组合起来就是下面这个完整文件按顺序粘贴到.cpp文件里就能编译运行。#include algorithm #include cmath #include cstdio #include random #include vector struct Point2D { double x; double y; }; struct Circle2D { Point2D center; double radius; }; double PointToCircleDistance(const Point2D pt, const Circle2D circle) { return std::abs(std::hypot(pt.x - circle.center.x, pt.y - circle.center.y) - circle.radius); } bool SolveLinear3x3(const double A[3][3], const double b[3], double x[3]) { double M[3][4]; for (int i 0; i 3; i) { for (int j 0; j 3; j) M[i][j] A[i][j]; M[i][3] b[i]; } for (int col 0; col 3; col) { int pivot col; for (int r col 1; r 3; r) { if (std::abs(M[r][col]) std::abs(M[pivot][col])) pivot r; } if (std::abs(M[pivot][col]) 1e-12) return false; if (pivot ! col) { for (int j 0; j 4; j) std::swap(M[pivot][j], M[col][j]); } double div M[col][col]; for (int j 0; j 4; j) M[col][j] / div; for (int r 0; r 3; r) { if (r col) continue; double factor M[r][col]; for (int j 0; j 4; j) M[r][j] - factor * M[col][j]; } } for (int i 0; i 3; i) x[i] M[i][3]; return true; } bool FitCircle3Points(const Point2D p1, const Point2D p2, const Point2D p3, Circle2D circle) { const double ax p1.x, ay p1.y; const double bx p2.x, by p2.y; const double cx p3.x, cy p3.y; const double d 2.0 * (ax * (by - cy) bx * (cy - ay) cx * (ay - by)); if (std::abs(d) 1e-12) return false; const double a2 ax * ax ay * ay; const double b2 bx * bx by * by; const double c2 cx * cx cy * cy; const double ux (a2 * (by - cy) b2 * (cy - ay) c2 * (ay - by)) / d; const double uy (a2 * (cx - bx) b2 * (ax - cx) c2 * (bx - ax)) / d; circle.center.x ux; circle.center.y uy; circle.radius std::hypot(ax - ux, ay - uy); return true; } bool FitCircleAlgebraic(const std::vectorPoint2D pts, Circle2D circle) { if (pts.size() 3) return false; double sx 0, sy 0; double sxx 0, syy 0, sxy 0; double sx3 0, sy3 0, sx2y 0, sxy2 0; for (const auto p : pts) { double x p.x, y p.y; double x2 x * x, y2 y * y; sx x; sy y; sxx x2; syy y2; sxy x * y; sx3 x2 * x; sy3 y2 * y; sx2y x2 * y; sxy2 x * y2; } const double n static_castdouble(pts.size()); double A[3][3] { { sxx, sxy, sx }, { sxy, syy, sy }, { sx, sy, n } }; double b[3] { -(sx3 sxy2), -(sx2y sy3), -(sxx syy) }; double x[3]; if (!SolveLinear3x3(A, b, x)) return false; const double a x[0]; const double bcoef x[1]; const double c x[2]; double rr (a * a bcoef * bcoef) / 4.0 - c; if (rr 0) return false; circle.center.x -a / 2.0; circle.center.y -bcoef / 2.0; circle.radius std::sqrt(rr); return true; } void RefineCircleGeometric(const std::vectorPoint2D pts, Circle2D circle, int maxIter 10) { for (int iter 0; iter maxIter; iter) { double J11 0, J12 0, J13 0; double J22 0, J23 0, J33 0; double r1 0, r2 0, r3 0; for (const auto p : pts) { double dx p.x - circle.center.x; double dy p.y - circle.center.y; double dist std::hypot(dx, dy); if (dist 1e-12) continue; double residual dist - circle.radius; double d_cx -dx / dist; double d_cy -dy / dist; double d_r -1.0; J11 d_cx * d_cx; J12 d_cx * d_cy; J13 d_cx * d_r; J22 d_cy * d_cy; J23 d_cy * d_r; J33 d_r * d_r; r1 d_cx * residual; r2 d_cy * residual; r3 d_r * residual; } double A[3][3] { {J11, J12, J13}, {J12, J22, J23}, {J13, J23, J33} }; double b[3] { -r1, -r2, -r3 }; double delta[3]; if (!SolveLinear3x3(A, b, delta)) break; circle.center.x delta[0]; circle.center.y delta[1]; circle.radius delta[2]; double step std::sqrt(delta[0] * delta[0] delta[1] * delta[1] delta[2] * delta[2]); if (step 1e-10) break; } } Circle2D RansacCircleFit(const std::vectorPoint2D points, double threshold 0.1, int maxIterations 2000, int minInliers 3, unsigned int seed 12345) { const int n static_castint(points.size()); Circle2D best; int bestCount 0; if (n 3) return best; std::mt19937 rng(seed); std::uniform_int_distributionint dist(0, n - 1); int iterations maxIterations; for (int iter 0; iter iterations; iter) { int idx[3]; idx[0] dist(rng); do { idx[1] dist(rng); } while (idx[1] idx[0]); do { idx[2] dist(rng); } while (idx[2] idx[0] || idx[2] idx[1]); Circle2D candidate; if (!FitCircle3Points(points[idx[0]], points[idx[1]], points[idx[2]], candidate)) { continue; } int count 0; for (const auto p : points) { if (PointToCircleDistance(p, candidate) threshold) count; } if (count bestCount) { bestCount count; best candidate; double inlierRatio static_castdouble(count) / n; double outlierRatio 1.0 - inlierRatio; if (outlierRatio 1.0) { double p_success 0.99; int needed static_castint( std::log(1.0 - p_success) / std::log(1.0 - std::pow(1.0 - outlierRatio, 3))) 1; needed std::max(1, needed); if (needed iterations) iterations needed; } } } if (bestCount minInliers) return Circle2D{}; std::vectorPoint2D inliers; inliers.reserve(bestCount); for (int pass 0; pass 2; pass) { inliers.clear(); for (const auto p : points) { if (PointToCircleDistance(p, best) threshold) { inliers.push_back(p); } } if (static_castint(inliers.size()) 3) break; Circle2D refined; if (FitCircleAlgebraic(inliers, refined)) { RefineCircleGeometric(inliers, refined); best refined; } } return best; } int main() { const int kNumPoints 200; const double kNoise 0.08; const double kOutlierRatio 0.4; std::mt19937 rng(42); std::normal_distributiondouble gauss(0.0, kNoise); std::uniform_real_distributiondouble unif(-8.0, 8.0); std::vectorPoint2D points; points.reserve(kNumPoints); for (int i 0; i kNumPoints; i) { double angle 2.0 * M_PI * i / kNumPoints; double x 5.0 * std::cos(angle) gauss(rng); double y 5.0 * std::sin(angle) gauss(rng); if (i static_castint(kNumPoints * kOutlierRatio)) { x unif(rng); y unif(rng); } points.push_back({x, y}); } Circle2D naive; FitCircleAlgebraic(points, naive); Circle2D result RansacCircleFit(points, 0.2, 5000, 20, 12345); std::printf(true : center(0,0) r5.0\n); std::printf(ls : center(%7.3f,%7.3f) r%.3f\n, naive.center.x, naive.center.y, naive.radius); std::printf(ransac : center(%7.3f,%7.3f) r%.3f\n, result.center.x, result.center.y, result.radius); return 0; }如果你处理的是三维点云里的圆记得先做一个平面拟合把所有点投影到该平面上再做圆拟合。RANSAC 主循环完全可以复用只需把输入从二维点换成投影后的二维坐标就行。这个思路在结构光测量、工业视觉定位、点云工件识别里都能直接用。