[R] Bounding ellipse for any set of points

From: Andy Lyons <ajlyons_at_berkeley.edu>
Date: Thu, 24 Mar 2011 17:23:00 -0700

   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
    R-help_at_r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Received on Fri 25 Mar 2011 - 07:27:45 GMT

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.2.0, at Fri 25 Mar 2011 - 08:40:24 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.

list of date sections of archive