Mathematica--三维点云圆拟合(二)
·
接着上一节,这节介绍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在拟合圆中对异常点的鲁棒性。但是实际过程,如果数据存在一定程度缺失,或者理想状态,优先最小二乘拟合圆。
更多推荐




所有评论(0)