Circle CircleFitByTaubin (Data& data)
/*
Circle fit to a given set of data points (in 2D)
This is an algebraic fit, due to Taubin, based on the journal article
G. Taubin, "Estimation Of Planar Curves, Surfaces And Nonplanar
Space Curves Defined By Implicit Equations, With
Applications To Edge And Range Image Segmentation",
IEEE Trans. PAMI, Vol. 13, pages 1115-1138, (1991)
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
The method is based on the minimization of the function
sum [(x-a)^2 + (y-b)^2 - R^2]^2
F = -------------------------------
sum [(x-a)^2 + (y-b)^2]
This method is more balanced than the simple Kasa fit.
It works well whether data points are sampled along an entire circle or
along a small arc.
It still has a small bias and its statistical accuracy is slightly
lower than that of the geometric fit (minimizing geometric distances),
but slightly higher than that of the very similar Pratt fit.
Besides, the Taubin fit is slightly simpler than the Pratt fit
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,A3,A33;
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);
circle.s = Sigma(data,circle);
circle.i = 0;
circle.j = iter; // return the number of iterations, too
return circle;
}