Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Some time ago I wrote an R function to fit an ellipse to point data, using an algorithm developed by Radim Halíř and Jan Flusser1 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:

```fit.ellipse <- function (x, y = NULL) {
# from:
# http://r.789695.n4.nabble.com/Fitting-a-half-ellipse-curve-tp2719037p2720560.html
#
# Least squares fitting of an ellipse to point data
# using the algorithm described in:
#   Radim Halir & Jan Flusser. 1998.
#   Numerically stable direct least squares fitting of ellipses.
#   Proceedings of the 6th International Conference in Central Europe
#   on Computer Graphics and Visualization. WSCG '98, p. 125-132
#
# Adapted from the original Matlab code by Michael Bedward (2010)
# [email protected]
#
# Subsequently improved by John Minter (2012)
#
# Arguments:
# x, y - x and y coordinates of the data points.
#        If a single arg is provided it is assumed to be a
#        two column matrix.
#
# Returns a list with the following elements:
#
# coef - coefficients of the ellipse as described by the general
#        quadratic:  ax^2 + bxy + cy^2 + dx + ey + f = 0
#
# center - center x and y
#
# major - major semi-axis length
#
# minor - minor semi-axis length
#
EPS <- 1.0e-8
dat <- xy.coords(x, y)

D1 <- cbind(dat\$x * dat\$x, dat\$x * dat\$y, dat\$y * dat\$y)
D2 <- cbind(dat\$x, dat\$y, 1)
S1 <- t(D1) %*% D1
S2 <- t(D1) %*% D2
S3 <- t(D2) %*% D2
T <- -solve(S3) %*% t(S2)
M <- S1 + S2 %*% T
M <- rbind(M[3,] / 2, -M[2,], M[1,] / 2)
evec <- eigen(M)\$vec
cond <- 4 * evec[1,] * evec[3,] - evec[2,]^2
a1 <- evec[, which(cond > 0)]
f <- c(a1, T %*% a1)
names(f) <- letters[1:6]

# calculate the center and lengths of the semi-axes
#
# see http://www.ncbi.nlm.nih.gov/pmc/articles/PMC2288654/
# J. R. Minter
# for the center, linear algebra to the rescue
# center is the solution to the pair of equations
# 2ax +  by + d = 0
# bx  + 2cy + e = 0
# or
# | 2a   b |   |x|   |-d|
# |  b  2c | * |y| = |-e|
# or
# A x = b
# or
# x = Ainv b
# or
# x = solve(A) %*% b
A <- matrix(c(2*f, f, f, 2*f), nrow=2, ncol=2, byrow=T )
b <- matrix(c(-f, -f), nrow=2, ncol=1, byrow=T)
soln <- solve(A) %*% b

b2 <- f^2 / 4

center <- c(soln, soln)
names(center) <- c("x", "y")

num  <- 2 * (f * f^2 / 4 + f * f^2 / 4 + f * b2 - f*f*f/4 - f*f*f)
den1 <- (b2 - f*f)
den2 <- sqrt((f - f)^2 + 4*b2)
den3 <- f + f

semi.axes <- sqrt(c( num / (den1 * (den2 - den3)),  num / (den1 * (-den2 - den3)) ))

# calculate the angle of rotation
term <- (f - f) / f
angle <- atan(1 / term) / 2

list(coef=f, center = center, major = max(semi.axes), minor = min(semi.axes), angle = unname(angle))
}

```
Next here is a utility function which takes a fitted ellipse and returns a matrix of vertices for plotting:
```get.ellipse <- function( fit, n=360 )
{
# Calculate points on an ellipse described by
# the fit argument as returned by fit.ellipse
#
# n is the number of points to render

tt <- seq(0, 2*pi, length=n)
sa <- sin(fit\$angle)
ca <- cos(fit\$angle)
ct <- cos(tt)
st <- sin(tt)

x <- fit\$center + fit\$maj * ct * ca - fit\$min * st * sa
y <- fit\$center + fit\$maj * ct * sa + fit\$min * st * ca

cbind(x=x, y=y)
}

```
And finally, some demo code from John:
```create.test.ellipse <- function(Rx=300,         # X-radius
Cx=250,         # X-center
Cy=150,         # Y-center
NoiseLevel=0.5) # Gaussian Noise level
{
set.seed(42)
t <- seq(0, 100, by=1)
x <- Rx * cos(t)
y <- Ry * sin(t)
nx <- x*cos(Rotation)-y*sin(Rotation) + Cx
nx <- nx + rnorm(length(t))*NoiseLevel
ny <- x*sin(Rotation)+y*cos(Rotation) + Cy
ny  <- ny + rnorm(length(t))*NoiseLevel
cbind(x=nx, y=ny)
}

X <- create.test.ellipse()
efit <- fit.ellipse(X)
e <- get.ellipse(efit)
plot(X)
lines(e, col="red")

print(efit)

```
1Halíř 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)