Circle CircleFitByHyper (Data& data)
/*
Circle fit to a given set of data points (in 2D)
This is an algebraic fit based on the journal article
A. Al-Sharadqah and N. Chernov, "Error analysis for circle fitting algorithms",
Electronic Journal of Statistics, Vol. 3, pages 886-911, (2009)
It is an algebraic circle fit with "hyperaccuracy" (with zero essential bias).
The term "hyperaccuracy" first appeared in papers by Kenichi Kanatani around 2006
Input: data - the class of data (contains the given points):
data.n - the number of data points
data.X[] - the array of X-coordinates
data.Y[] - the array of Y-coordinates
Output:
circle - parameters of the fitting circle:
circle.a - the X-coordinate of the center of the fitting circle
circle.b - the Y-coordinate of the center of the fitting circle
circle.r - the radius of the fitting circle
circle.s - the root mean square error (the estimate of sigma)
circle.j - the total number of iterations
This method combines the Pratt and Taubin fits to eliminate the essential bias.
It works well whether data points are sampled along an entire circle or
along a small arc.
Its statistical accuracy is theoretically higher than that of the Pratt fit
and Taubin fit, but practically they all return almost identical circles
(unlike the Kasa fit that may be grossly inaccurate).
It provides a very good initial guess for a subsequent geometric fit.
Nikolai Chernov (September 2012)
*/
{
int i,iter,IterMAX=99;
reals Xi,Yi,Zi;
reals Mz,Mxy,Mxx,Myy,Mxz,Myz,Mzz,Cov_xy,Var_z;
reals A0,A1,A2,A22;
reals Dy,xnew,x,ynew,y;
reals DET,Xcenter,Ycenter;
Circle circle;
data.means(); // Compute x- and y- sample means (via a function in the class "data")
// computing moments
Mxx=Myy=Mxy=Mxz=Myz=Mzz=0.;
for (i=0; i=abs(y)) break;
x = xnew; y = ynew;
}
// computing paramters of the fitting circle
DET = x*x - x*Mz + Cov_xy;
Xcenter = (Mxz*(Myy - x) - Myz*Mxy)/DET/Two;
Ycenter = (Myz*(Mxx - x) - Mxz*Mxy)/DET/Two;
// assembling the output
circle.a = Xcenter + data.meanX;
circle.b = Ycenter + data.meanY;
circle.r = sqrt(Xcenter*Xcenter + Ycenter*Ycenter + Mz - x - x);
circle.s = Sigma(data,circle);
circle.i = 0;
circle.j = iter; // return the number of iterations, too
return circle;
}