[R] Bounding ellipse for any set of points

[Andy Lyons]

   After a lot of effort I developed the following function to compute the
   bounding ellipse (also known a as minimum volume enclosing ellipsoid) for
   any set of points. This script is limited to two dimensions, but I believe
   with minor modification the algorithm should work for 3 or more dimensions.
   It seems to work well (although I don't know if can be optimized to run
   faster) and hope it may be useful to someone else. Andy
   ######################################################################
   ## mvee()
   ## Uses the Khachiyan Algorithm to find the the minimum volume enclosing
   ## ellipsoid (MVEE) of a set of points. In two dimensions, this is just
   ## the bounding ellipse (this function is limited to two dimensions).
   ## Adapted by Andy Lyons from Matlab code by Nima Moshtagh.
   ## Copyright (c) 2009, Nima Moshtagh
   ##         [1]http://www.mathworks.com/matlabcentral/fileexchange/9542
   ##         [2]http://www.mathworks.com/matlabcentral/fileexchange/13844
   ##         [3]http://stackoverflow.com/questions/1768197/bounding-ellipse
   ##
   ## Parameters
   ## xy          a two-column data frame containing x and y coordinates
   ##              if  NULL then a random sample set of 10 points will be
   generated
   ## tolerance   a tolerance value (default = 0.005)
   ## plotme      FALSE/TRUE. Plots the points and ellipse. Default TRUE.
   ## max.iter    The maximum number of iterations. If the script tries this
   ##             number of iterations but still can't get to the tolerance
   ##             value, it displays an error message and returns NULL
   ## shiftxy     TRUE/FALSE. If True, will apply a shift to the coordinates to
   make them
   ##             smaller and speed up the matrix calculations, then reverse
   the shift
   ##             to the center point of the resulting ellipoid. Default TRUE
   ## Output:     A list containing the "center form" matrix equation of the
   ##              ellipse.  i.e.  a  2x2 matrix "A" and a 2x1 vector "C"
   representing
   ##             the center of the ellipse such that:
   ##             (x - C)' A (x - C) <= 1
   ##              Also in the list is a 2x1 vector elps.axes.lngth whose
   elements
   ##              are  one-half  the lengths of the major and minor axes
   (variables
   ##             a and b
   ##             Also in list is alpha, the angle of rotation
   ######################################################################
   mvee <- function(xy=NULL, tolerance = 0.005, plotme = TRUE, max.iter = 500,
   shiftxy = TRUE) {
       if (is.null(xy)) {
           xy <- data.frame(x=runif(10,100,200), y=runif(10,100,200))
       } else if (ncol(xy) != 2) {
           warning("xy must be a two-column data frame")
           return(NULL)
       }
       ## Number of points
       n = nrow(xy)
       ## Dimension of the points (2)
       d = ncol(xy)
       if (n <= d) return(NULL)
        ##  Apply  a  uniform shift to the x&y coordinates to make matrix
   calculations computationally
       ## simpler (if x and y are very large, for example UTM coordinates, this
   may be necessary
       ## to prevent a 'compuationally singular matrix' error
       if (shiftxy) {
           xy.min <- sapply(xy, FUN = "min")
       } else {
           xy.min <- c(0,0)
       }
       xy.use <- xy - rep(xy.min, each = n)
       ## Add a column of 1s to the (n x 2) matrix xy - so it is now (n x 3)
       Q <- t(cbind(xy.use, rep(1,n)))
       ## Initialize
       count <- 1
       err <- 1
       u <- rep(1/n, n)
       ## Khachiyan Algorithm
       while (err > tolerance)
       {
           ## see
   [4]http://stackoverflow.com/questions/1768197/bounding-ellipse
           ## for commented code
           X <- Q %*% diag(u) %*% t(Q)
           M <- diag(t(Q) %*% solve(X) %*% Q)
           maximum <- max(M)
           j <- which(M == maximum)
           step_size = (maximum - d -1) / ((d+1)*(maximum-1))
           new_u <- (1 - step_size) * u
           new_u[j] <- new_u[j] + step_size
           err <- sqrt(sum((new_u - u)^2))
           count <- count + 1
           if (count > max.iter) {
               warning(paste("Iterated", max.iter, "times and still can't find
   the bounding ellipse. \n", sep=""))
               warning("Either increase the tolerance or the maximum number of
   iterations")
               return(NULL)
           }
           u <- new_u
       }
       ## Put the elements of the vector u into the diagonal of a matrix
       U <- diag(u)
       ## Take the transpose of xy
       P <- t(xy.use)
       ## Compute the center, adding back the shifted values
       c <- as.vector((P %*% u) + xy.min)
       ## Compute the A-matrix
       A <- (1/d) * solve(P %*% U %*% t(P) - (P %*% u) %*% t(P %*% u))
       ## Find the Eigenvalues of matrix A which will be used to get the major
   and minor axes
       A.eigen <- eigen(A)
       ## Calculate the length of the semi-major and semi-minor axes
       ## from the Eigenvalues of A.
       semi.axes <- sort(1 / sqrt(A.eigen$values), decreasing=TRUE)
       ## Calculate the rotation angle from the first Eigenvector
       alpha <- atan2(A.eigen$vectors[2,1], A.eigen$vectors[1,1]) - pi/2
       if(plotme) {
           ## Plotting commands adapted from code by Alberto Monteiro
           ## [5]https://stat.ethz.ch/pipermail/r-help/2006-October/114652.html
           ## Create the points for the ellipse
           theta <- seq(0, 2 * pi, length = 72)
           a <- semi.axes[1]
           b <- semi.axes[2]
           elp.plot.xs <- c[1] + a * cos(theta) * cos(alpha) - b * sin(theta) *
   sin(alpha)
           elp.plot.ys <- c[2] + a * cos(theta) * sin(alpha) + b * sin(theta) *
   cos(alpha)
           ## Plot the ellipse with the same scale on each axis
           plot(elp.plot.xs, elp.plot.ys, type = "l", lty="dotted", col="blue",
   asp=1,
                main="minimum volume enclosing ellipsoid", xlab=names(xy)[1],
   ylab=names(xy)[2])
           ## Plot the original points
           points(xy[,1], xy[,2], type="p", pch=16)
           ## add the center of the ellipse using a triangle symbol
           points(c[1], c[2], pch=2, col="blue")
       }
       ellipse.params <- list("A" = A, "c" = c, "ab" = semi.axes, alpha=alpha)
   }

References

   1. http://www.mathworks.com/matlabcentral/fileexchange/9542
   2. http://www.mathworks.com/matlabcentral/fileexchange/13844
   3. http://stackoverflow.com/questions/1768197/bounding-ellipse
   4. http://stackoverflow.com/questions/1768197/bounding-ellipse
   5. https://stat.ethz.ch/pipermail/r-help/2006-October/114652.html