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)
      )
}

Answer 2

Well, looky here: an R-blogger column has written some code to fit to ellipses and circles. His code, which I won't repost here, is based on previous work done by Radim Halíř and Jan Flusser in Matlab. His code includes (commented) the original Matlab lines for comparison.

I've peeked at a number of papers on this topic, and can only say that I'm not qualified to determine which algorithms are the most robust. For those interested, take a look at these papers:

http://www.emis.de/journals/BBMS/Bulletin/sup962/gander.pdf

http://ralph.cs.cf.ac.uk/papers/Geometry/fit.pdf

http://autotrace.sourceforge.net/WSCG98.pdf

Followup edit: I ran Spacedman's code against the linked R-code for fitting ellipses, using the same "noisy" set of 1e5 points on a circle as input. The results are:

testcircle<-create.test.ellipse(Rx=200,Ry=200,Rot=.56,Noise=5.5,leng=100000)
 dim(testcircle)
[1] 100000      2

microbenchmark(fitSS(testcircle),fit.ellipse(testcircle))
Unit: milliseconds
                    expr       min        lq    median        uq       max
       fitSS(testcircle) 649.98245 704.05751 731.61282 787.84212 2053.7096
 fit.ellipse(testcircle)  25.74518  33.87718  38.87143  95.23499  256.2475
 neval
   100
   100

For reference, the output of the two fitting functions were:

From SSfit, the list

ssfit
$par
[1] 249.9530 149.9927 200.0512

$value
[1] 185.8195

$counts
function gradient 
     134       NA 

$convergence
[1] 0

$message
NULL

From fit.ellipse, we get

ellfit
$coef
            a             b             c             d             e 
-7.121109e-01 -1.095501e-02 -7.019815e-01  3.563866e+02  2.136497e+02 
            f 
-3.195427e+04 

$center
       x        y 
249.0769 150.2326 

$major
[1] 201.7601

$minor
[1] 199.6424

$angle
[1] 0.412268

You can see that the elliptic equation's coefficients are near-zero for terms which "deviate" from a circle; plotting the two results yields almost indistinguishable curves.

Answer 3

To fit an ellipse, there is the fitEllipse function in the PlaneGeometry package. It uses the fitConic package.

library(PlaneGeometry)

library(PlaneGeometry)
# the "true" ellipse: 
ell <- Ellipse$new(center = c(1, 1), rmajor = 3, rminor = 2, alpha = 25)
# We add some noise to 30 points on this ellipse:
set.seed(666L)
points <- ell$randomPoints(30, "on") + matrix(rnorm(30*2, sd = 0.2), ncol = 2)
# Now we fit an ellipse to these points:
ellFitted <- fitEllipse(points)
# let's draw all this stuff, true ellipse in blue, fitted ellipse in green:
box <- ell$boundingbox()
plot(NULL, asp = 1, xlim = box$x, ylim = box$y, xlab = NA, ylab = NA)
draw(ell, border = "blue", lwd = 2)
points(points, pch = 19)
draw(ellFitted, border = "green", lwd = 2)