# How to find the minimum-area-rectangle for given points?

The following is another problem that I am working on it these days. As you see in the figure, (BTW, I tried to make it self-presenter) the question is:

How to find the minimum-area-rectangle (MAR) fitted on the given points?

and a supporting question is:

Is there any analytical solution for the problem?

(A development of the question will be to fit a box (3D) to a cluster of points in a 3D point cloud.)

As a first stage I propose to find the convex-hull for the points which reforms the problem (by removing those points are not involved in the solution) to: fitting a MAR to a polygon. The required method will provide X (center of rectangle), D (two dimensions) and A (angle).

My proposal for solution:

• Find the centroid of the polygon (see this post)
• [S] Fit a simple fitted rectangle i.e., parallel to the axes X and Y
• you may use `minmax` function for X and Y of the given points (e.g., polygon's vertices)
• Store the area of the fitted rectangle
• Rotate the polygon about the centroid by e.g., 1 degree
• Repeat from [S] until a full rotation done
• Report the angle of the minimum area as the result

It looks to me promising, however the following problems exist:

• choose of a good resolution for the angle change could be challenging,
• the computation cost is high,
• the solution is not analytical but experimental.

-
+1 for an interesting problem, and the way you present your question. – stakx Feb 3 '13 at 18:49

Yes, there is an analytical solution for this problem. The algorithm you are looking for is known in polygon generalisation as "smallest surrounding rectangle".

The algorithm you describe is fine but in order to solve the problems you have listed, you can use the fact that the orientation of the MAR is the same as the one of one of the edges of the point cloud convex hull. So you just need to test the orientations of the convex hull edges. You should:

• Compute the convex hull of the cloud.
• For each edge of the convex hull:
• compute the edge orientation (with arctan),
• rotate the convex hull using this orientation in order to compute easily the bounding rectangle area with min/max of x/y of the rotated convex hull,
• Store the orientation corresponding to the minimum area found,
• Return the rectangle corresponding to the minimum area found.

An example of implementation in java is available there.

In 3D, the same applies, except:

• The convex hull will be a volume,
• The orientations tested will be the orientations (in 3D) of the convex hull faces.

Good luck!

-
+1 Very nice answer! I would like to point out that actual rotation of the cloud is unnecessary. First--you probably meant this--only the vertices of the hull have to be considered. Second, instead of rotating, represent the current side as a pair of orthogonal unit vectors. Taking their dot products with the hull vertex coordinates (which could be done as a single matrix operation) gives the rotated coordinates: no trigonometry necessary, fast, and perfectly accurate. – whuber Apr 5 '12 at 15:01
Thanks for the links. Indeed, rotating only for # of edges makes the proposed method very efficient. I could find the paper proves that. Although I marked this as the answer for loyalty to the first good answer (cannot choose two/more great answers :( ) I would like to recommend strongly considering whuber's complete answer below. The efficiency of the given method there (avoiding rotations!) is incredible, and the whole procedure is only a few lines of code. To me it is readily translatable to Python :) – Developer Apr 6 '12 at 0:56
Can you please update java implementation link ? – Myra Aug 23 '12 at 13:49
yes, it is done ! – julien Aug 23 '12 at 15:19
@julien, looks like you might need to update the java implementation link again. Getting a 404. – Fezter Apr 22 '14 at 0:18

To supplement @julien's great solution, here is a working implementation in `R`, which could serve as pseudocode to guide any GIS-specific implementation (or be applied directly in `R`, of course). Input is an array of point coordinates. Output (the value of `mbr`) is an array of the vertices of the minimum bounding rectangle (with the first one repeated to close it). Note the complete absence of any trigonometric calculations.

``````library(alphahull)                                  # Exposes ashape()
MBR <- function(points) {
# Analyze the convex hull edges
a <- ashape(points, alpha=1000)                 # One way to get a convex hull...
e <- a\$edges[, 5:6] - a\$edges[, 3:4]            # Edge directions
norms <- apply(e, 1, function(x) sqrt(x %*% x)) # Edge lengths
v <- diag(1/norms) %*% e                        # Unit edge directions
w <- cbind(-v[,2], v[,1])                       # Normal directions to the edges

# Find the MBR
vertices <- (points) [a\$alpha.extremes, 1:2]    # Convex hull vertices
minmax <- function(x) c(min(x), max(x))         # Computes min and max
x <- apply(vertices %*% t(v), 2, minmax)        # Extremes along edges
y <- apply(vertices %*% t(w), 2, minmax)        # Extremes normal to edges
areas <- (y[1,]-y[2,])*(x[1,]-x[2,])            # Areas
k <- which.min(areas)                           # Index of the best edge (smallest area)

# Form a rectangle from the extremes of the best edge
cbind(x[c(1,2,2,1,1),k], y[c(1,1,2,2,1),k]) %*% rbind(v[k,], w[k,])
}
``````

Here is an example of its use:

``````# Create sample data
set.seed(23)
points <- matrix(rnorm(20*2), ncol=2)                 # Random (normally distributed) points
mbr <- MBR(points)

# Plot the hull, the MBR, and the points
limits <- apply(mbr, 2, function(x) c(min(x),max(x))) # Plotting limits
plot(ashape(points, alpha=1000), col="Gray", pch=20,
xlim=limits[,1], ylim=limits[,2])                # The hull
lines(mbr, col="Blue", lwd=3)                         # The MBR
points(points, pch=19)                                # The points
``````

Timing is limited by the speed of the convex hull algorithm, because the number of vertices in the hull is almost always much less than the total. Most convex hull algorithms are asymptotically O(n*log(n)) for n points: you can compute almost as fast as you can read the coordinates.

-
Indeed, no need to rotate! thanks. – julien Apr 5 '12 at 18:25
+1 What a amazing solution! Such an idea comes only after having long experiences. From now I will be curious to optimize my existing codes being inspired with this great answer. – Developer Apr 6 '12 at 0:59
I wish I could upvote this twice. I am learning R and your answers are a continual source of inspiration. – John Barça Feb 6 '15 at 10:00

I just implemented this myself and posted my answer over on StackOverflow, but I figured I'd drop my version here for others to view:

``````import numpy as np
from scipy.spatial import ConvexHull

def minimum_bounding_rectangle(points):
"""
Find the smallest bounding rectangle for a set of points.
Returns a set of points representing the corners of the bounding box.

:param points: an nx2 matrix of coordinates
:rval: an nx2 matrix of coordinates
"""
from scipy.ndimage.interpolation import rotate
pi2 = np.pi/2.

# get the convex hull for the points
hull_points = points[ConvexHull(points).vertices]

# calculate edge angles
edges = np.zeros((len(hull_points)-1, 2))
edges = hull_points[1:] - hull_points[:-1]

angles = np.zeros((len(edges)))
angles = np.arctan2(edges[:, 1], edges[:, 0])

angles = np.abs(np.mod(angles, pi2))
angles = np.unique(angles)

# find rotation matrices
# XXX both work
rotations = np.vstack([
np.cos(angles),
np.cos(angles-pi2),
np.cos(angles+pi2),
np.cos(angles)]).T
#     rotations = np.vstack([
#         np.cos(angles),
#         -np.sin(angles),
#         np.sin(angles),
#         np.cos(angles)]).T
rotations = rotations.reshape((-1, 2, 2))

# apply rotations to the hull
rot_points = np.dot(rotations, hull_points.T)

# find the bounding points
min_x = np.nanmin(rot_points[:, 0], axis=1)
max_x = np.nanmax(rot_points[:, 0], axis=1)
min_y = np.nanmin(rot_points[:, 1], axis=1)
max_y = np.nanmax(rot_points[:, 1], axis=1)

# find the box with the best area
areas = (max_x - min_x) * (max_y - min_y)
best_idx = np.argmin(areas)

# return the best box
x1 = max_x[best_idx]
x2 = min_x[best_idx]
y1 = max_y[best_idx]
y2 = min_y[best_idx]
r = rotations[best_idx]

rval = np.zeros((4, 2))
rval[0] = np.dot([x1, y2], r)
rval[1] = np.dot([x2, y2], r)
rval[2] = np.dot([x2, y1], r)
rval[3] = np.dot([x1, y1], r)

return rval
``````

Here are four different examples of it in action. For each example, I generated 4 random points and found the bounding box.

It's relatively quick too for these samples on 4 points:

``````>>> %timeit minimum_bounding_rectangle(a)
1000 loops, best of 3: 245 µs per loop
``````
-

There is a tool in Whitebox GAT (http://www.uoguelph.ca/~hydrogeo/Whitebox/) called Minimum Bounding Box for solving this exact problem. There is also a minimum convex hull tool in there too. Several of the tools in the Patch Shape toolbox, e.g. patch orientation and elongation, are based on finding the minimum bounding box.

-

I came across this thread while looking for a Python solution for a minimum-area bounding rectangle.

Here's my implementation, for which the results were verified with Matlab.

Test code is included for simple polygons, and I am using it to find the 2D minimum bounding box and axes directions for a 3D PointCloud.

-
Has your answer been deleted? – Paul Richter Jan 5 at 16:46

Thanks @whuber's answer. It is a great solution, but slow for big point cloud. I found `convhulln` function in R package `geometry` is much faster (138 s vs 0.03 s for 200000 points). I pasted my codes here for anyone is interesting for a faster solution.

``````library(alphahull)                                  # Exposes ashape()
MBR <- function(points) {
# Analyze the convex hull edges
a <- ashape(points, alpha=1000)                 # One way to get a convex hull...
e <- a\$edges[, 5:6] - a\$edges[, 3:4]            # Edge directions
norms <- apply(e, 1, function(x) sqrt(x %*% x)) # Edge lengths
v <- diag(1/norms) %*% e                        # Unit edge directions
w <- cbind(-v[,2], v[,1])                       # Normal directions to the edges

# Find the MBR
vertices <- (points) [a\$alpha.extremes, 1:2]    # Convex hull vertices
minmax <- function(x) c(min(x), max(x))         # Computes min and max
x <- apply(vertices %*% t(v), 2, minmax)        # Extremes along edges
y <- apply(vertices %*% t(w), 2, minmax)        # Extremes normal to edges
areas <- (y[1,]-y[2,])*(x[1,]-x[2,])            # Areas
k <- which.min(areas)                           # Index of the best edge (smallest area)

# Form a rectangle from the extremes of the best edge
cbind(x[c(1,2,2,1,1),k], y[c(1,1,2,2,1),k]) %*% rbind(v[k,], w[k,])
}

MBR2 <- function(points) {
tryCatch({
a2 <- geometry::convhulln(points, options = 'FA')

e <- points[a2\$hull[,2],] - points[a2\$hull[,1],]            # Edge directions
norms <- apply(e, 1, function(x) sqrt(x %*% x)) # Edge lengths

v <- diag(1/norms) %*% as.matrix(e)                        # Unit edge directions

w <- cbind(-v[,2], v[,1])                       # Normal directions to the edges

# Find the MBR
vertices <- as.matrix((points) [a2\$hull, 1:2])    # Convex hull vertices
minmax <- function(x) c(min(x), max(x))         # Computes min and max
x <- apply(vertices %*% t(v), 2, minmax)        # Extremes along edges
y <- apply(vertices %*% t(w), 2, minmax)        # Extremes normal to edges
areas <- (y[1,]-y[2,])*(x[1,]-x[2,])            # Areas
k <- which.min(areas)                           # Index of the best edge (smallest area)

# Form a rectangle from the extremes of the best edge
as.data.frame(cbind(x[c(1,2,2,1,1),k], y[c(1,1,2,2,1),k]) %*% rbind(v[k,], w[k,]))
}, error = function(e) {
assign('points', points, .GlobalEnv)
stop(e)
})
}

# Create sample data
#set.seed(23)
points <- matrix(rnorm(200000*2), ncol=2)                 # Random (normally distributed) points
system.time(mbr <- MBR(points))
system.time(mmbr2 <- MBR2(points))

# Plot the hull, the MBR, and the points
limits <- apply(mbr, 2, function(x) c(min(x),max(x))) # Plotting limits
plot(ashape(points, alpha=1000), col="Gray", pch=20,
xlim=limits[,1], ylim=limits[,2])                # The hull
lines(mbr, col="Blue", lwd=10)                         # The MBR
lines(mbr2, col="red", lwd=3)                         # The MBR2
points(points, pch=19)
``````

Two methods get the same answer (example for 2000 points):

-