Some time ago I wrote an R function to fit an ellipse to point data, using an algorithm developed by Radim Halíř and Jan Flusser¹ in Matlab, and posted it to the r-help list. The implementation was a bit hacky, returning odd results for some data.

A couple of days ago, an email arrived from John Minter asking for a pointer to the original code. I replied with a link and mentioned that I’d be interested to know if John made any improvements to the code. About ten minutes later, John emailed again with a much improved version ! Not only is it more reliable, but also more efficient. So with many thanks to John, here is the improved code:

<br />fit.ellipse <- function (x, y = NULL) {<br />  # from:<br />  # http://r.789695.n4.nabble.com/Fitting-a-half-ellipse-curve-tp2719037p2720560.html<br />  #<br />  # Least squares fitting of an ellipse to point data<br />  # using the algorithm described in: <br />  #   Radim Halir & Jan Flusser. 1998. <br />  #   Numerically stable direct least squares fitting of ellipses. <br />  #   Proceedings of the 6th International Conference in Central Europe <br />  #   on Computer Graphics and Visualization. WSCG '98, p. 125-132 <br />  #<br />  # Adapted from the original Matlab code by Michael Bedward (2010)<br />  # [email protected]<br />  #<br />  # Subsequently improved by John Minter (2012)<br />  # <br />  # Arguments: <br />  # x, y - x and y coordinates of the data points.<br />  #        If a single arg is provided it is assumed to be a<br />  #        two column matrix.<br />  #<br />  # Returns a list with the following elements: <br />  #<br />  # coef - coefficients of the ellipse as described by the general <br />  #        quadratic:  ax^2 + bxy + cy^2 + dx + ey + f = 0 <br />  #<br />  # center - center x and y<br />  #<br />  # major - major semi-axis length<br />  #<br />  # minor - minor semi-axis length<br />  #<br />  EPS <- 1.0e-8 <br />  dat <- xy.coords(x, y) <br />  <br />  D1 <- cbind(dat$x * dat$x, dat$x * dat$y, dat$y * dat$y) <br />  D2 <- cbind(dat$x, dat$y, 1) <br />  S1 <- t(D1) %*% D1 <br />  S2 <- t(D1) %*% D2 <br />  S3 <- t(D2) %*% D2 <br />  T <- -solve(S3) %*% t(S2) <br />  M <- S1 + S2 %*% T <br />  M <- rbind(M[3,] / 2, -M[2,], M[1,] / 2) <br />  evec <- eigen(M)$vec <br />  cond <- 4 * evec[1,] * evec[3,] - evec[2,]^2 <br />  a1 <- evec[, which(cond > 0)] <br />  f <- c(a1, T %*% a1) <br />  names(f) <- letters[1:6] <br />  <br />  # calculate the center and lengths of the semi-axes <br />  #<br />  # see http://www.ncbi.nlm.nih.gov/pmc/articles/PMC2288654/<br />  # J. R. Minter<br />  # for the center, linear algebra to the rescue<br />  # center is the solution to the pair of equations<br />  # 2ax +  by + d = 0<br />  # bx  + 2cy + e = 0<br />  # or<br />  # | 2a   b |   |x|   |-d|<br />  # |  b  2c | * |y| = |-e|<br />  # or<br />  # A x = b<br />  # or<br />  # x = Ainv b<br />  # or<br />  # x = solve(A) %*% b<br />  A <- matrix(c(2*f[1], f[2], f[2], 2*f[3]), nrow=2, ncol=2, byrow=T )<br />  b <- matrix(c(-f[4], -f[5]), nrow=2, ncol=1, byrow=T)<br />  soln <- solve(A) %*% b<br /><br />  b2 <- f[2]^2 / 4<br />  <br />  center <- c(soln[1], soln[2]) <br />  names(center) <- c("x", "y") <br />  <br />  num  <- 2 * (f[1] * f[5]^2 / 4 + f[3] * f[4]^2 / 4 + f[6] * b2 - f[2]*f[4]*f[5]/4 - f[1]*f[3]*f[6]) <br />  den1 <- (b2 - f[1]*f[3]) <br />  den2 <- sqrt((f[1] - f[3])^2 + 4*b2) <br />  den3 <- f[1] + f[3] <br />  <br />  semi.axes <- sqrt(c( num / (den1 * (den2 - den3)),  num / (den1 * (-den2 - den3)) )) <br />  <br />  # calculate the angle of rotation <br />  term <- (f[1] - f[3]) / f[2] <br />  angle <- atan(1 / term) / 2 <br />  <br />  list(coef=f, center = center, major = max(semi.axes), minor = min(semi.axes), angle = unname(angle)) <br />}<br /><br />

Next here is a utility function which takes a fitted ellipse and returns a matrix of vertices for plotting:

<br />get.ellipse <- function( fit, n=360 ) <br />{<br />  # Calculate points on an ellipse described by <br />  # the fit argument as returned by fit.ellipse <br />  # <br />  # n is the number of points to render <br />  <br />  tt <- seq(0, 2*pi, length=n) <br />  sa <- sin(fit$angle) <br />  ca <- cos(fit$angle) <br />  ct <- cos(tt) <br />  st <- sin(tt) <br />  <br />  x <- fit$center[1] + fit$maj * ct * ca - fit$min * st * sa <br />  y <- fit$center[2] + fit$maj * ct * sa + fit$min * st * ca <br />  <br />  cbind(x=x, y=y) <br />}<br /><br />

And finally, some demo code from John:

<br />create.test.ellipse <- function(Rx=300,         # X-radius<br />                                Ry=200,         # Y-radius<br />                                Cx=250,         # X-center<br />                                Cy=150,         # Y-center<br />                                Rotation=0.4,   # Radians<br />                                NoiseLevel=0.5) # Gaussian Noise level<br />{<br />  set.seed(42)<br />  t <- seq(0, 100, by=1)<br />  x <- Rx * cos(t)<br />  y <- Ry * sin(t)<br />  nx <- x*cos(Rotation)-y*sin(Rotation) + Cx<br />  nx <- nx + rnorm(length(t))*NoiseLevel <br />  ny <- x*sin(Rotation)+y*cos(Rotation) + Cy<br />  ny  <- ny + rnorm(length(t))*NoiseLevel<br />  cbind(x=nx, y=ny)<br />}<br /><br />X <- create.test.ellipse()<br />efit <- fit.ellipse(X)<br />e <- get.ellipse(efit)<br />plot(X) <br />lines(e, col="red") <br /><br />print(efit)<br /><br />

¹Halíř R., Flusser J.: Numerically stable direct least squares fitting of ellipses. In: Proceedings of the 6th International Conference in Central Europe on Computer Graphics and Visualization. WSCG ’98. CZ, Plzeň 1998, pp. 125-132. (postscript file)