The mvee() function is intended to be released under the BSD license.
Copyright (c) 2009, Nima Moshtagh
Copyright (c) 2011, Andy Lyons
All rights reserved.
http://www.opensource.org/licenses/bsd-license.php
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
Redistributions of source code must retain the above copyright notice,
this list of conditions and the following disclaimer.
Redistributions in binary form must reproduce the above copyright
notice, this list of conditions and the following disclaimer in the
documentation and/or other materials provided with the distribution.
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
POSSIBILITY OF SUCH DAMAGE.
>Date: Thu, 24 Mar 2011 17:23:00 -0700
>From: Andy Lyons <ajlyons at berkeley.edu>
>To: r-help at r-project.org
>Subject: [R] Bounding ellipse for any set of points
>Message-ID: <6.2.1.2.2.20110324051124.117a0560 at
calmail.berkeley.edu>
>Content-Type: text/plain; charset="us-ascii"
>
>
> 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
