# Convex hull, (minimum) bounding box, minimum enclosing ellipse, and minimum enclosing circle

## Convex hull

chull() returns the ordered indices of the matrix rows corresponding to corners of the convex hull.

set.seed(123)
xy      <- matrix(rnorm(24, 100, 15), ncol=2)
hullIdx <- chull(xy)
hullIdx
[1] 3 6 9 8 4
plot(xy, xlab="x", ylab="y", asp=1, type="n")
polygon(xy[hullIdx, ], border="blue", lwd=2)
points(xy, pch=16, cex=1.5)

## Bounding box

getBoundingBox <- function(xy) {
stopifnot(is.matrix(xy), is.numeric(xy), nrow(xy) >= 2, ncol(xy) == 2)
x   <- range(xy[ , 1])
y   <- range(xy[ , 2])
pts <- c(xleft=x[1], ybottom=y[1], xright=x[2], ytop=y[2])
return(list(pts=pts, width=abs(diff(x)), height=abs(diff(y))))
}
bb <- getBoundingBox(xy)
plot(xy, xlab="x", ylab="y", asp=1, type="n")
rect(bb$pts[1], bb$pts[2], bb$pts[3], bb$pts[4], border="blue", lwd="2")
points(xy, pch=16, cex=1.5)
bb$width * bb$height
[1] 2516.849

## Minimum bounding box

### Rotating calipers algorithm

getMinBBox <- function(xy) {
stopifnot(is.matrix(xy), is.numeric(xy), nrow(xy) >= 2, ncol(xy) == 2)

## rotating calipers algorithm using the convex hull
H    <- chull(xy)                    # hull indices, vertices ordered clockwise
n    <- length(H)                    # number of hull vertices
hull <- xy[H, ]                      # hull vertices

## unit basis vectors for all subspaces spanned by the hull edges
hDir  <- diff(rbind(hull, hull[1,])) # account for circular hull vertices
hLens <- sqrt(rowSums(hDir^2))       # length of basis vectors
huDir <- diag(1/hLens) %*% hDir      # scaled to unit length

## unit basis vectors for the orthogonal subspaces
## rotation by 90 deg -> y' = x, x' = -y
ouDir <- cbind(-huDir[ , 2], huDir[ , 1])

## project hull vertices on the subspaces spanned by the hull edges, and on
## the subspaces spanned by their orthogonal complements - in subspace coords
projMat <- rbind(huDir, ouDir) %*% t(hull)

## range of projections and corresponding width/height of bounding rectangle
rangeH  <- matrix(numeric(n*2), ncol=2)   # hull edge
rangeO  <- matrix(numeric(n*2), ncol=2)   # orth subspace
widths  <- numeric(n)
heights <- numeric(n)
for(i in seq(along=H)) {
rangeH[i, ] <- range(projMat[  i, ])
rangeO[i, ] <- range(projMat[n+i, ])  # orth subspace is in 2nd half
widths[i]   <- abs(diff(rangeH[i, ]))
heights[i]  <- abs(diff(rangeO[i, ]))
}

## extreme projections for min-area rect in subspace coordinates
eMin  <- which.min(widths*heights)   # hull edge leading to minimum-area
hProj <- rbind(   rangeH[eMin, ], 0)
oProj <- rbind(0, rangeO[eMin, ])

## move projections to rectangle corners
hPts <- sweep(hProj, 1, oProj[ , 1], "+")
oPts <- sweep(hProj, 1, oProj[ , 2], "+")

## corners in standard coordinates, rows = x,y, columns = corners
## in combined (4x2)-matrix: reverse point order to be usable in polygon()
basis <- cbind(huDir[eMin, ], ouDir[eMin, ])  # basis formed by hull edge and orth
hCorn <- basis %*% hPts
oCorn <- basis %*% oPts
pts   <- t(cbind(hCorn, oCorn[ , c(2, 1)]))

return(list(pts=pts, width=widths[eMin], height=heights[eMin]))
}

### Draw the minimum bounding box

mbb <- getMinBBox(xy)       ## minimum bounding box
H   <- chull(xy)            ## convex hull

# plot original points, convex hull, and minimum bounding box
plot(xy, xlab="x", ylab="y", asp=1, type="n",
xlim=range(c(xy[ , 1], mbb$pts[ , 1])), ylim=range(c(xy[ , 2], mbb$pts[ , 2])))
polygon(xy[H, ], col=NA)    ## show convex hull
polygon(mbb$pts, border="blue", lwd=2) points(xy, pch=16, cex=1.5) mbb$width * mbb$height ## box area [1] 2127.255 ## Minimum enclosing ellipse ### Kachiyan’s algorithm Adapted from Matlab code by Jacob on Stack Overflow. getMinEllipse <- function(xy, tol=0.001, max_iter=1000) { stopifnot(is.matrix(xy), is.numeric(xy), nrow(xy) >= 2, ncol(xy) == 2) H <- chull(xy) # convex hull indices (vertices ordered clockwise) hPts <- t(xy[H, ]) # points that make up the convex hull p <- nrow(hPts) # dimension N <- ncol(hPts) # number of points on hull Q <- rbind(hPts, rep(1, N)) iter <- 1 # iteration err <- 1 # remaining error u <- rep(1/N, N) # Khachiyan's algorithm while((err > tol) && (iter < max_iter)) { X <- Q %*% diag(u) %*% t(Q) m <- diag(t(Q) %*% solve(X) %*% Q) maximum <- max(m) j <- which.max(m) step_size <- (maximum - p-1) / ((p+1)*(maximum-1)) new_u <- (1 - step_size)*u new_u[j] <- new_u[j] + step_size err <- sqrt(sum((new_u-u)^2)) iter <- iter + 1 u <- new_u } if(iter >= max_iter) { warning(paste("Maximum number of iterations reached. Error is still:", err)) } ctr <- hPts %*% u # ellipse center E <- (1/p) * solve((hPts %*% diag(u) %*% t(hPts)) - tcrossprod(ctr)) shape <- solve(E) area0 <- pi # area unit circle area <- area0 * sqrt(det(shape)) # area ellipse rad <- sqrt(eigen(shape)$values) # length ellipse semi axes

}

### Draw the minimum enclosing ellipse

me     <- getMinEllipse(xy)
CF     <- chol(me$shape, pivot=TRUE) # Cholesky-factor CFord <- order(attr(CF, "pivot")) angles <- seq(0, 2*pi, length.out=100) # angles in radians ell <- 1 * cbind(cos(angles), sin(angles)) %*% CF[ , CFord] # ellipse ellCtr <- sweep(ell, 2, me$ctr, "+")      # move ellipse to center

# determine axis limits so that the ellipse will be visible
xLims <- me$ctr[1] + c(-me$rad[1], me$rad[1]) yLims <- me$ctr[2] + c(-me$rad[1], me$rad[1])

# draw points and ellipse
plot(xy, xlab="x", ylab="y", xlim=xLims, ylim=yLims, asp=1, type="n")
points(xy, pch=16)
points(me$ctr[1], me$ctr[2], pch=4, cex=1.5, lwd=2, col="blue")  # center
polygon(ellCtr, border="blue")         # ellipse

## Minimum enclosing circle

Skyum algorithm based on the convex hull

### Identify collinear points on convex hull

## assumes that points in xy are already oredered
idCollinear <- function(xy) {
if(!is.matrix(xy))  { stop("xy must be a matrix") }
if(!is.numeric(xy)) { stop("xy must be numeric") }
if(nrow(xy) < 3L)   { stop("xy must have at least 3 points") }
if(ncol(xy) != 2L)  { stop("xy must have 2 columns") }

n    <- nrow(xy)
idx  <- seq_len(n)
post <- (idx %% n) + 1              # next point in S
prev <- idx[order(post)]            # previous point in S

del <- integer(0)
## check all sets of 3 consecutive points if they lie in 1D sub-space
for(i in idx) {
pts <- rbind(xy[prev[i], ],
xy[i, ],
xy[post[i], ])
pts_rank <- qr(scale(pts, center=TRUE, scale=FALSE))$rank if(pts_rank < 2L) { del[length(del) + 1] <- i } } del } ### Circle defined by three points getCircleFrom3 <- function(xy) { stopifnot(is.matrix(xy), is.numeric(xy), nrow(xy) == 3, ncol(xy) == 2) aa <- xy[1, ] bb <- xy[2, ] cc <- xy[3, ] y <- xy[ , 2] xDeltaA <- bb[1] - aa[1] yDeltaA <- bb[2] - aa[2] xDeltaB <- cc[1] - bb[1] yDeltaB <- cc[2] - bb[2] xDeltaC <- cc[1] - aa[1] yDeltaC <- cc[2] - aa[2] ## check if the points are collinear: qr(xy)$rank == 1, or:
## determinant of difference matrix = 0, no need to use det()
dMat <- rbind(c(xDeltaA, yDeltaA), c(xDeltaB, yDeltaB))
if(isTRUE(all.equal(dMat[1,1]*dMat[2,2] - dMat[1,2]*dMat[2,1], 0, check.attributes=FALSE))) {
## define the circle as the one that's centered between the points
rangeX <- range(c(aa[1], bb[1], cc[1]))
rangeY <- range(c(aa[2], bb[2], cc[2]))
ctr    <- c(rangeX[1] + 0.5*diff(rangeX), rangeY[1] + 0.5*diff(rangeY))
} else {
v1  <- rowSums(xy^2)                    # first vector in the numerator
v2x <- c( xDeltaB, -xDeltaC,  xDeltaA)  # 2nd vector numerator for Mx
v2y <- c(-yDeltaB,  yDeltaC, -yDeltaA)  # 2nd vector numerator for My
ctr <- c(t(v1) %*% v2y, t(v1) %*% v2x) / (2 * (t(y) %*% v2x))  # center
}

}

### Vertex that produces the circle with the maximum radius

Used later in getMinCircle()

getMaxRad <- function(xy, S) {
stopifnot(is.matrix(xy), is.numeric(xy), nrow(xy) >= 2, ncol(xy) == 2)
stopifnot(is.numeric(S), length(S) >= 2, length(S) <= nrow(xy))

n    <- length(S)                    # number of points
Sidx <- seq(along=S)                 # index for points
post <- (Sidx %% n) + 1              # next point in S
prev <- Sidx[order(post)]            # previous point in S
for(i in Sidx) {
pts     <- rbind(xy[S[prev[i]], ], xy[S[i], ], xy[S[post[i]], ])
rad    <- maxPD$d / 2 # half the distance -> radius rangeX <- c(hPts[idx[1], 1], hPts[idx[2], 1]) rangeY <- c(hPts[idx[1], 2], hPts[idx[2], 2]) ctr <- c(rangeX[1] + 0.5*diff(rangeX), rangeY[1] + 0.5*diff(rangeY)) ## check if circle centered between hPts[pt1Idx, ] and hPts[pt2Idx, ] ## contains all points (all distances <= rad) dst2ctr <- dist(rbind(ctr, hPts[-idx, ])) # distances to center if(all(as.matrix(dst2ctr)[-1, 1] <= rad)) { # if all <= rad, we're done tri <- rbind(hPts[idx, ], ctr) return(getCircleFrom3(tri)) } ## min circle touches hull in three points - Skyum algorithm S <- H # copy of hull indices that will be changed while(length(S) >= 2) { n <- length(S) # number of remaining hull vertices Sidx <- seq(along=S) # index for vertices post <- (Sidx %% n) + 1 # next vertex in S prev <- Sidx[order(post)] # previous vertex in S mIdx <- getMaxRad(xy, S) # idx for maximum radius ## triangle where mIdx is vertex B in ABC Smax <- rbind(xy[S[prev[mIdx]], ], xy[S[mIdx], ], xy[S[post[mIdx]], ]) ## if there's only two hull vertices, we're done if(n <= 2) { break } ## check if angle(ABC) is > 90 ## if so, eliminate B - if not, we're done if(isBiggerThan90(Smax)) { S <- S[-mIdx] } else { break } } return(getCircleFrom3(Smax)) } ### Draw the minimal enclosing circle mc <- getMinCircle(xy) angles <- seq(0, 2*pi, length.out=200) circ <- cbind(mc$ctr[1] + mc$rad*cos(angles), mc$ctr[2] + mc$rad*sin(angles)) # determine axis limits so that the circle will be visible xLims <- mc$ctr[1] + c(-mc$rad, mc$rad)
yLims <- mc$ctr[2] + c(-mc$rad, mc\$rad)
plot(xy, xlab="x", ylab="y", xlim=xLims, ylim=yLims, asp=1, type="n")
lines(circ, col="blue", lwd=2)
points(xy, pch=16, cex=1.5)

## Notes

Note that the presented implementation of the minimum bounding box and of the minimum enclosing circle were written ad-hoc “to get the job done” with a limited amount of mostly random point data. The implementations are not rigorously tested for numerical stability, and are not optimized. Advanced implementations (with generalizations to higher dimensions) may be found in solid computer graphics libraries.