Fitting an ellipse to point data

September 30, 2012
By

(This article was first published on Last Resort Software, and kindly contributed to R-bloggers)

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[1], f[2], f[2], 2*f[3]), nrow=2, ncol=2, byrow=T )
b <- matrix(c(-f[4], -f[5]), nrow=2, ncol=1, byrow=T)
soln <- solve(A) %*% b

b2 <- f[2]^2 / 4

center <- c(soln[1], soln[2])
names(center) <- c("x", "y")

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])
den1 <- (b2 - f[1]*f[3])
den2 <- sqrt((f[1] - f[3])^2 + 4*b2)
den3 <- f[1] + f[3]

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

# calculate the angle of rotation
term <- (f[1] - f[3]) / f[2]
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[1] + fit$maj * ct * ca - fit$min * st * sa
y <- fit$center[2] + 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
Ry=200, # Y-radius
Cx=250, # X-center
Cy=150, # Y-center
Rotation=0.4, # Radians
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)

To leave a comment for the author, please follow the link and comment on his blog: Last Resort Software.

R-bloggers.com offers daily e-mail updates about R news and tutorials on topics such as: visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...



If you got this far, why not subscribe for updates from the site? Choose your flavor: e-mail, twitter, RSS, or facebook...

Comments are closed.