**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 Flusser^{1} 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)

^{1}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)

**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...