R help - Jul 2012 - Fit circle with R

Dear Researchers,

I wrote two function to fit a circle using noisy data.

1- the fitCircle() is derived from MATLAB code of * zhak Bucher* from the
link
http://www.mathworks.com/matlabcentral/fileexchange/5557-circle-fit/content/circfit.m
2- the CircleFitByPratt() from MATLAB code of *Nikolai Chernov *from the
link
http://www.mathworks.com/matlabcentral/fileexchange/22643-circle-fit-pratt-method/content/CircleFitByPratt.m,
based on:
*V. Pratt, "Direct least-squares fitting of algebraic surfaces",
Computer
Graphics, Vol. 21, pages 145-152 (1987)*

 I am looking for new methods to compare and improve my analysis because
the error increase with decreasing of points used in the functions.

Thanks for all suggestions
Gianni


Here the funtions with example

# fitCircle, returns:
# xf,yf = centre of the fitted circle
# Rf = radius of the fitted circle
# Cf = circumference of the fitted circle
# Af = Area of the fitted circle

fitCircle <- function(x,y){
     a = qr.solve(cbind(x,y,rep(1,length(x))),cbind(-(x^2+y^2)))
     xf = -.5*a[1]
     yf = -.5*a[2]
     Rf  =  sqrt((a[1]^2+a[2]^2)/4-a[3])
    Cf = 2*pi*Rf
    Af = pi*(Rf^2)
    m <- cbind(xf,yf,Rf,Cf,Af)
    return(m)}

# CircleFitByPratt, returns:
#   [,1] and [,2]  = centre of the fitted circle
# [,3] = radius of fitted cirlce


CircleFitByPratt <- function(x,y){

    n <- length(x)
    centroid <- cbind(mean(x),mean(y))

    Mxx=0; Myy=0; Mxy=0; Mxz=0; Myz=0; Mzz=0;

    for(i in 1:n){
        Xi <- x[[i]] - centroid[1]
        Yi <- y[[i]] - centroid[2]
        Zi <- (Xi*Xi) + (Yi*Yi)
        Mxy = Mxy + Xi*Yi;
        Mxx = Mxx + Xi*Xi;
        Myy = Myy + Yi*Yi;
        Mxz = Mxz + Xi*Zi;
        Myz = Myz + Yi*Zi;
        Mzz = Mzz + Zi*Zi;
    }

    Mxx = Mxx/n
    Myy = Myy/n
    Mxy = Mxy/n
    Mxz = Mxz/n
    Myz = Myz/n
    Mzz = Mzz/n

    # computing the coefficients of the characteristic polynomial
    Mz = Mxx + Myy;
    Cov_xy = Mxx*Myy - Mxy*Mxy;
    Mxz2 = Mxz*Mxz;
    Myz2 = Myz*Myz;

    A2 = 4*Cov_xy - 3*Mz*Mz - Mzz;
    A1 = Mzz*Mz + 4*Cov_xy*Mz - Mxz2 - Myz2 - Mz*Mz*Mz;
    A0 = Mxz2*Myy + Myz2*Mxx - Mzz*Cov_xy - 2*Mxz*Myz*Mxy + Mz*Mz*Cov_xy;
    A22 = A2 + A2;

    epsilon=1e-12;
    ynew=1e+20;
    IterMax=20;
    xnew = 0;

    # Newton's method starting at x=0
    epsilon=1e-12;
    ynew=1e+20;
    IterMax=20;
    xnew = 0;
    iter=1:IterMax

    for (i in 1:IterMax){
        yold = ynew;
        ynew = A0 + xnew*(A1 + xnew*(A2 + 4.*xnew*xnew));
        if (abs(ynew) > abs(yold)){
            print('Newton-Pratt goes wrong direction: |ynew| >
|yold|')
            xnew = 0;
            break
        }
        Dy = A1 + xnew*(A22 + 16*xnew*xnew);
        xold = xnew;
        xnew = xold - ynew/Dy;
        if (abs((xnew-xold)/xnew) < epsilon) {break}
        if(iter[[i]] >= IterMax){
            print('Newton-Pratt will not converge');
            xnew = 0;
        }
        if(xnew < 0.){
            print('Newton-Pratt negative root:  x=',xnew);
        }
    }

    DET = xnew*xnew - xnew*Mz + Cov_xy;
    Center = cbind(Mxz*(Myy-xnew)-Myz*Mxy , Myz*(Mxx-xnew)-Mxz*Mxy)/DET/2;

    #    computing the circle parameters

    DET = xnew*xnew - xnew*Mz + Cov_xy;
    Center = cbind(Mxz*(Myy-xnew)-Myz*Mxy , Myz*(Mxx-xnew)-Mxz*Mxy)/DET/2;

    Par = cbind(Center+centroid , sqrt(Center[2]*Center[2]+Mz+2*xnew));
    return(Par)
}

#EXAMPLE
library(plotrix)

# Create a Circle of radius=10, centre=5,5
R = 10; x_c = 5; y_c = 5;
thetas = seq(0,pi,(pi/64));
xs = x_c + R*cos(thetas)
ys = y_c + R*sin(thetas)
# Now add some random noise
mult = 0.5;
xs = xs+mult*rnorm(rnorm(xs));
ys = ys+mult*rnorm(rnorm(ys));
plot(xs,ys,pch=19,cex=0.5,col="red",xlim=c(-10,20),ylim=c(-10,20),asp=1)
# real circle
draw.circle(x_c,y_c,radius=10,border="black")
points(x_c,y_c,,pch=4,col="black")

CPrat <- CircleFitByPratt(xs,ys)
draw.circle(CPrat[1],CPrat[2],radius=CPrat[3],border="blue")
points(CPrat[1],CPrat[2],pch=4,col="blue")

MyC <- fitCircle(xs,ys)
draw.circle(MyC[1],MyC[2],radius=MyC[3],border="green")
points(MyC[1],MyC[2],pch=4,col="green")

# Select less points
points(xs[20:49],ys[20:49])


MyC1 <- fitCircle(xs[20:49],ys[20:49])
draw.circle(MyC1[1],MyC1[2],radius=MyC1[3],border="blue",lty=2,lwd=2)

CPrat1 <- CircleFitByPratt(xs[20:49],ys[20:49])
draw.circle(CPrat1[1],CPrat1[2],radius=CPrat1[3],border="green",lty=2,lwd=2)
points(CPrat[1],CPrat[2],pch=4,col="red")

	[[alternative HTML version deleted]]