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.
[1] 3 6 9 8 4plot(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)
[1] 2516.849Minimum 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)
[1] 2127.255Minimum 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
list(ctr=ctr, shape=shape, area=area, rad=rad)
}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))
rad <- sqrt((0.5*diff(rangeX))^2 + (0.5*diff(rangeY))^2)
} else {
rad <- prod(dist(xy)) / (2 * abs(det(cbind(xy, 1)))) # circle radius
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
}
return(list(ctr=ctr, rad=rad))
}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
rads <- numeric(n) # radii for all circles
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]], ])
rads[i] <- getCircleFrom3(pts)$rad # circle radius
}
return(which.max(rads))
}Check if the angle at B in triangle ABC exceeds 90 degrees
isBiggerThan90 <- function(xy) {
stopifnot(is.matrix(xy), is.numeric(xy), nrow(xy) == 3, ncol(xy) == 2)
d <- dist(xy)
dAB <- d[1]
dAC <- d[2]
dBC <- d[3]
return((dAB^2 + dBC^2 - dAC^2) < 0)
}Maximum pairwise distance between two 2D-points
getMaxPairDist <- function(xy) {
stopifnot(is.matrix(xy), is.numeric(xy), ncol(xy) == 2, nrow(xy) >= 2)
# 2D -> only convex hull is relevant
H <- chull(xy) # convex hull indices (vertices ordered clockwise)
pts <- xy[H, ] # points that make up the convex hull
N <- nrow(pts) # number of points on hull
dMat <- dist(pts, method="euclidean") # distance matrix
idx <- which.max(as.matrix(dMat)) # maximum distance
i <- (idx-1) %/% N+1 # column -> point 1
j <- (idx-1) %% N+1 # row -> point 2
mPts <- H[c(i, j)] # rows with max distance
dst <- max(dMat) # max distance
return(list(d=dst, idx=mPts))
}Minimal enclosing circle
getMinCircle <- function(xy) {
stopifnot(is.matrix(xy), is.numeric(xy), nrow(xy) >= 2, ncol(xy) == 2)
H <- chull(xy) # convex hull indices (vertices ordered clockwise)
hPts <- xy[H, ] # points that make up the convex hull
## remove collinear points on convex hull, if any
del <- idCollinear(hPts)
if(length(del) >= 1L) {
H <- H[-del]
hPts <- hPts[-del, ]
if(length(H) < 2L) {
stop("less than 2 points left after removing collinear points on convex hull")
}
}
## min circle may touch convex hull in only two points
## if so, it is centered between the hull points with max distance
maxPD <- getMaxPairDist(hPts)
idx <- maxPD$idx # index of points with max distance
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.
- Computational Geometry Algorithms Library. 2021. CGAL Chapter 65: Bounding Volumes
- Fischer, K.; Gaertner, B.; Kutz, M. 2003. Fast smallest-enclosing-ball computation in high dimensions. In: Proceedings of the 11th European Symposium on Algorithms (ESA), 630-641
- Gaertner, B. 2021. Miniball: Smallest Enclosing Balls of Points
- Skyum, S. 1991. A simple algorithm for computing the smallest enclosing circle. Information Processing Letters 37(3), 121-125.
- Toussaint, G. T. 1983. Solving geometric problems with the rotating calipers. In: Proceedings of the 1983 IEEE MELECON. Athens, Greece: IEEE Computer Society.
- Welzl, E. 1991. Smallest enclosing disks (balls and ellipsoids). In: Maurer H. (eds), New Results and New Trends in Computer Science 555, 359-370. doi 10.1007/BFb0038202.
- Todd MJ and Yildrim EA. On Khachiyan’s Algorithm for the Computation of Minimum Volume Enclosing Ellipsoids. Discrete Applied Mathematics 2007: 155: 1731-1744. PDF
Get the article source from GitHub
R markdown - markdown - R code - all posts
top