How to find a best fit circle/ellipse using R?

Question

I've been reading about a few methods to fit a circle to data (like this). I would like to see how the methods work on real data and thought of using R for this. I tried searching rseek for packages that can help with this but came up with nothing useful.

So, are there packages that help to easily compute the best fit circle for a given data set (similar to how lm() will fit a linear model to a data set)? Otherwise, how might one perform such a task in R?

Answer 1

Here's a fairly naive implementation of a function that minimises SS(a,b,r) from that paper:

fitSS <- function(xy,
                  a0=mean(xy[,1]),
                  b0=mean(xy[,2]),
                  r0 = mean(sqrt((xy[,1]-a0)^2 + (xy[,2]-b0)^2)),
                  ...){
    SS <- function(abr){
        sum((abr[3] - sqrt((xy[,1]-abr[1])^2 + (xy[,2]-abr[2])^2))^2)
    }
    optim(c(a0,b0,r0), SS, ...)
}

I've written a couple of supporting functions to generate random data on circles and to plot circles. Hence:

> xy = sim_circles(10)
> f = fitSS(xy)

The fit$par value is a vector of xcenter, ycenter, radius.

> plot(xy,asp=1,xlim=c(-2,2),ylim=c(-2,2))
> lines(circlexy(f$par))

Note it doesn't use the gradients nor does it check the error code for convergence. You can supply it with initial values or it can have a guess.

Code for plotting and generating circles follows:

circlexy <- function(xyr, n=180){
    theta = seq(0,2*pi,len=n)
    cbind(xyr[1] + xyr[3]*cos(theta),
          xyr[2] + xyr[3]*sin(theta)
          )
}
sim_circles <- function(n,x=0,y=0,r=1,sd=0.05){
    theta = runif(n, 0, 2*pi)
    r = r + rnorm(n, mean=0, sd=sd)
    cbind(x + r*cos(theta),
          y + r*sin(theta)
      )
}