Clarification about ellipsoidhull in R

Question

I wonder if someone can help me with some questions about the function ellipsoidhull in the R package cluster. I am using it to find a minimum ellipse containing a series of 2d points. e.g.

library(cluster)
d <- matrix(c(1,2,3,1,3,2),ncol=2)
e <- ellipsoidhull(d)

The function computes the ellipse OK, returning a structure that includes the center of the ellipse and a covariance matrix.

summary(e)
## 'ellipsoid' in 2 dimensions:
## center = ( 2 2 ); squared ave.radius d^2 =  2 
## and shape matrix =
##        [,1]    [,2]
##   [1,] 0.66667 0.33333
##   [2,] 0.33333 0.66667
##     hence, area  =  3.6276 

Questions

a) How can I use this data to check if a given point belongs to the ellipse?

b) How can I use this data to compute the distance from a given point to the ellipse?

Answer 1

We can try the following:

library(cluster)
d <- matrix(c(1,2,3,1,3,2),ncol=2)
e <- ellipsoidhull(d)
eg <- eigen(e$cov)
axes <- sqrt(eg$values)
angle <- atan(eg$vectors[1,1]/eg$vectors[2,1]) # angle of major axis with x axis

# check if the point (xp, yp) belongs to the ellipse with parameters a,b,... with tolerance eps
belongs.to <- function (xp, yp, a, b, x0, y0, alpha, eps=1e-3) {
  return(abs((cos(alpha)*(xp-x0)+sin(alpha)*(yp-y0))^2/a^2+(sin(alpha)*(xp-x0)-cos(alpha)*(yp-y0))^2/b^2 - 1) <= eps)
} 

# check if the point (xp, yp) is inside the ellipse with parameters a,b,...
is.inside <- function (xp, yp, a, b, x0, y0, alpha) {
  return((cos(alpha)*(xp-x0)+sin(alpha)*(yp-y0))^2/a^2+(sin(alpha)*(xp-x0)-cos(alpha)*(yp-y0))^2/b^2 <= 1)
}

# plot ellipse
plot(e$loc, xlim=c(0,4), ylim=c(0,4), main = "ellipsoidhull", xlab='x', ylab='y')
lines(predict(e), col="blue")
points(rbind(e$loc), col = "red", cex = 3, pch = 13)

x0 <- e$loc[1] # centroid locations
y0 <- e$loc[2]  
a <- sqrt(e$d2) * axes[1]  # major axis length
b <- sqrt(e$d2) * axes[2]  # minor axis length

alpha <- angle
xp <- 3
yp <- 2.9
is.inside(xp, yp, a, b, x0, y0, alpha)
# [1] TRUE
points(xp, yp, pch=19, col='green')
xp <- 3
yp <- 3.1
is.inside(xp, yp, a, b, x0, y0, alpha)
# [1] FALSE
points(xp, yp, pch=19, col='blue')
xp <- 3
yp <- 3
belongs.to(xp, yp, a, b, x0, y0, alpha)
# [1] TRUE
points(xp, yp, pch=19, col='pink')


# distance of a point from the center of the ellipse
sqrt((xp-x0)^2+(yp-y0)^2)