做视觉检测、激光点云处理或者机器人导航这类活儿的朋友，大概率都遇到过同一个需求：一堆散乱的点里，明明藏着一个圆，可你用最小二乘去拟合，结果被几个离群点带偏，中心偏了，半径也不对。这种场景最适合的方案就是 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 满足 OA=OB=OC。解这个方程可以得到：
分母 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) ) / D
D 的绝对值如果接近 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 Σx | a | -Σ(x³ + xy²) |
| Σxy Σy² Σy | b | -Σ(x²y + y³) |
| Σx Σy n | c | -Σ(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::vector<Point2D>& 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_cast<double>(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::vector<Point2D>& points, double threshold = 0.1, int maxIterations = 2000, int minInliers = 3, unsigned int seed = 12345) { const int n = static_cast<int>(points.size()); Circle2D best; int bestCount = 0; if (n < 3) return best; std::mt19937 rng(seed); std::uniform_int_distribution<int> 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_cast<double>(count) / n; double outlierRatio = 1.0 - inlierRatio; if (outlierRatio < 1.0) { double p_success = 0.99; int needed = static_cast<int>( 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::vector<Point2D> 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_cast<int>(inliers.size()) < 3) break; Circle2D refined; if (FitCircleAlgebraic(inliers, refined)) { RefineCircleGeometric(inliers, refined); best = refined; } } return best; }为什么分两步精修？第一步先用代数拟合把全体内点综合起来，得到一个比三点定圆更稳的初值；第二步再用几何距离精修，消除代数拟合的半径偏小问题。重复两轮是因为精修后模型会微调，内点集合可能略微变化，再收集一次可以保证用上最完整的内点集。
Gauss-Newton 几何精修实现如下：
void RefineCircleGeometric(const std::vector<Point2D>& 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_distribution<double> gauss(0.0, kNoise); std::uniform_real_distribution<double> unif(-8.0, 8.0); std::vector<Point2D> 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_cast<int>(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) r=5.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 px | 1.5~3 px |
| 激光雷达点云 | 0.01~0.05 m | 0.03~0.10 m |
| 机械接触测量 | 0.005~0.01 mm | 0.01~0.03 mm |
迭代数可以通过公式算：设内点率为 w，三点全为内点的概率是 w³，那么 N 次采样里至少有一次成功的概率 p 满足 p = 1 - (1 - w³)^N。反过来就能解出 N。比如内点率 65%，p 取 0.99，N 大约只需要 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::vector<Point2D>& 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_cast<double>(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::vector<Point2D>& 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::vector<Point2D>& points, double threshold = 0.1, int maxIterations = 2000, int minInliers = 3, unsigned int seed = 12345) { const int n = static_cast<int>(points.size()); Circle2D best; int bestCount = 0; if (n < 3) return best; std::mt19937 rng(seed); std::uniform_int_distribution<int> 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_cast<double>(count) / n; double outlierRatio = 1.0 - inlierRatio; if (outlierRatio < 1.0) { double p_success = 0.99; int needed = static_cast<int>( 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::vector<Point2D> 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_cast<int>(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_distribution<double> gauss(0.0, kNoise); std::uniform_real_distribution<double> unif(-8.0, 8.0); std::vector<Point2D> 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_cast<int>(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) r=5.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 主循环完全可以复用，只需把输入从二维点换成投影后的二维坐标就行。这个思路在结构光测量、工业视觉定位、点云工件识别里都能直接用。
