接着上一节,这节介绍RANSAC拟合圆。

RANSAC拟合圆:

原理其实与RANSAC拟合直线一样,随机选3个点拟合圆,并统计其余点到这个圆的残差,当残差小于设定的阀值,就记为内点,最后统计内点数量,满足条件就停止,找到最优的圆。

MMA代码:

ransacCircleFit[pts_, nIter_ : 1000, threshold_ : 0.01] := 
  Module[{n, bestCx, bestCy, bestR, bestInliers, sample, cx, cy, r, 
    residuals, inliers, nInliers}, n = Length[pts];
   bestInliers = 0;

   Do[sample = RandomSample[pts, 3];
   {cx, cy, r} = threePointCircle[sample];

    If[NumericQ[r] && r > 0,(*计算内点数*)
     residuals = Abs[Norm[# - {cx, cy}] - r] & /@ pts;
     nInliers = Count[residuals, x_ /; x < threshold];
     If[nInliers > bestInliers, bestInliers = nInliers;
      {bestCx, bestCy, bestR} = {cx, cy, r};];], {nIter}];

   (*用所有内点重新拟合*)
   inliers = 
    Select[pts, Abs[Norm[# - {bestCx, bestCy}] - bestR] < threshold &];

   fitCircleNonlinearLeastSquares[inliers]];(*最小二乘拟合圆*,具体代码见上一节)

(*三点确定圆的辅助函数*)
threePointCircle[{{x1_, y1_}, {x2_, y2_}, {x3_, y3_}}] := 
  Module[{a, b, c, d, cx, cy, r}, 
   a = x1*(y2 - y3) - y1*(x2 - x3) + x2*y3 - x3*y2;
   If[Abs[a] < 10^-10, Return[{0, 0, -1}]];
   b = (x1^2 + y1^2)*(y3 - y2) + (x2^2 + y2^2)*(y1 - y3) + (x3^2 + 
        y3^2)*(y2 - y1);
   c = (x1^2 + y1^2)*(x2 - x3) + (x2^2 + y2^2)*(x3 - x1) + (x3^2 + 
        y3^2)*(x1 - x2);
   cx = -b/(2 a);
   cy = -c/(2 a);
   r = Sqrt[(x1 - cx)^2 + (y1 - cy)^2];
   {cx, cy, r}];

代码中需要注意,fitCircleNonlinearLeastSquares就是上一节的线性最小二乘拟合圆算法。

重点说一下,这个三点确定圆的辅助函数,代码中利用了圆心坐标公式,并结合行列式中的特性解出来,可以说是最简单的实现。如果要去直接解方程组,要显式考虑共线性。

可视化:

与上一节的最小二乘法相比,拟合的效果最好,体现了RANSAC在拟合圆中对异常点的鲁棒性。但是实际过程,如果数据存在一定程度缺失,或者理想状态,优先最小二乘拟合圆。

Logo

汇聚全球AI编程工具,助力开发者即刻编程。

更多推荐