As you see in the figure, 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 Finding center of geometry of object?)
  • [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.

enter image description here


8 Answers 8


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 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!

  • 13
    +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
    Commented Apr 5, 2012 at 15:01
  • 4
    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
    Commented Apr 6, 2012 at 0:56
  • 2
    Note that the extension into 3D is a bit more complicated than that. Each face of the 3D convex hull defines a possible orientation of one face of the bounding box, but not the orientation of faces perpendicular to it. The problem of how to rotate the box in that plane becomes the 2D minimum-bounding-rectangle problem in the plane of that face. For each edge of the convex hull of the cloud projected onto a given plane you can draw a bounding box which will give you a different volume in 3D.
    – Will
    Commented Dec 5, 2017 at 11:12
  • 1
    @julien Do you know of a way to solve the "minimum bounding square" problem? In that case it seems the "the orientation of the MAR is the same as the one of one of the edges" assumption no longer holds. Commented Jul 12, 2022 at 19:58
  • 1
    @ThomasAhle Indeed ! Already for a triangle, the problem does not seem obvious. Maybe something to think about this summer !
    – julien
    Commented Jul 14, 2022 at 7:17

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.

MBR <- function(p) {
  # Analyze the convex hull edges     
  a <- chull(p)                                   # Indexes of extremal points
  a <- c(a, a[1])                                 # Close the loop
  e <- p[a[-1],] - p[a[-length(a)], ]             # Edge directions
  norms <- sqrt(rowSums(e^2))                     # Edge lengths
  v <- e / norms                                  # Unit edge directions
  w <- cbind(-v[,2], v[,1])                       # Normal directions to the edges

  # Find the MBR
  vertices <- p[a, ]                              # Convex hull vertices
  x <- apply(vertices %*% t(v), 2, range)         # Extremes along edges
  y <- apply(vertices %*% t(w), 2, range)         # 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
p <- matrix(rnorm(20*2), ncol=2)                 # Random (normally distributed) points
mbr <- MBR(p)

# Plot the hull, the MBR, and the points
limits <- apply(mbr, 2, range) # Plotting limits
plot(p[(function(x) c(x, x[1]))(chull(p)), ], 
     type="l", asp=1, bty="n", xaxt="n", yaxt="n",
     col="Gray", pch=20, 
     xlab="", ylab="",
     xlim=limits[,1], ylim=limits[,2])                # The hull
lines(mbr, col="Blue", lwd=3)                         # The MBR
points(p, 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.

  • 1
    +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
    Commented Apr 6, 2012 at 0:59
  • I wish I could upvote this twice. I am learning R and your answers are a continual source of inspiration. Commented Feb 6, 2015 at 10:00
  • 1
    @retrovius The bounding rectangle of a set of (rotated) points is determined by four numbers: the smallest x coordinate, the largest x coordinate, the smallest y coordinate, and the largest y coordinate. That's what the "extremes along edges" refers to.
    – whuber
    Commented May 29, 2019 at 11:13
  • 1
    @retrovius The origin plays no role in these calculations, because everything is based on differences of coordinates except at the end, where the best rectangle as computed in rotated coordinates is simply rotated back. Although it's a smart idea to use a coordinate system in which the origin is close to the points (to minimize loss of floating point precision), the origin otherwise is irrelevant.
    – whuber
    Commented May 29, 2019 at 21:56
  • 1
    @Retrovius You can interpret this in terms of a property of rotations: namely, the matrix of a rotation is orthogonal. Thus, one kind of resource would be a study of linear algebra (generally) or analytic Euclidean geometry (specifically). However, I have found that the easiest way to deal with rotations (and translations and rescalings) in the plane is to view the points as complex numbers: rotations are simply carried out by multiplying values by unit-length numbers.
    – whuber
    Commented May 30, 2019 at 13:06

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([
#     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.

enter image description here

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
  • 1
    Hi JesseBuesking, are you able to generate rectangles with 90 deg corners? Your code is working great for getting parallelograms but 90 deg corners is required in my specific use case. Could you recommend how your code can be modified to reach that? Thanks! Commented Nov 7, 2017 at 5:38
  • @NaderAlexan If you're asking whether it can handle squares, then yes it certainly can! I just tried it on a unit square points = np.array([[0, 0], [0, 1], [1, 0], [1, 1]]), and the output is array([[1.00000000e+00, 6.12323400e-17], [0.00000000e+00, 0.00000000e+00], [6.12323400e-17, 1.00000000e+00], [1.00000000e+00, 1.00000000e+00]]) which is the unit square itself (including some floating point rounding errors). Note: a square is just a rectangle with equal sides, so I'm assuming if it can handle a square it generalizes to all rectangles. Commented Nov 8, 2017 at 15:32
  • 1
    thank you for your answer. Yes, it is working great but I am attempting to force it to always produce a rectangle (4 sides with 90 deg angles for each side) over any other 4-sided polygon, although in certain cases it does produce a rectangle it does not seem to be a constant constraint, do you know of how to modify the code to add this constraint? Thanks! Commented Nov 9, 2017 at 5:05
  • Maybe gis.stackexchange.com/a/22934/48041 might guide you toward a solution, given their answer appears to have this constraint? Once you do find a solution, you should contribute it as I'm sure others will find it useful. Good luck! Commented Nov 10, 2017 at 13:42
  • 3
    @NaderAlexan I would assume that just the plots were created with different axis scaling for x and y, thus producing rectangles that look like parallelograms.
    – moooeeeep
    Commented Nov 12, 2019 at 10:51

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.

enter image description here


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.


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) {
        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)

# Create sample data
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):

enter image description here

  • 1
    Is it possible to extend this implementation to 3d space (i.e. find a minimum volume box which includes all given points in 3d space)?
    – Sasha
    Commented Dec 20, 2016 at 17:37

I simply recommend the OpenCV's build-in function minAreaRect, which finds a rotated rectangle of the minimum area enclosing the input 2D point set. To see how to use this function, one may refer to this tutorial.


There's another cool solution for this: https://stackoverflow.com/questions/32892932/create-the-oriented-bounding-box-obb-with-python-and-numpy

import numpy as np
from shapely.geometry.base import BaseGeometry
from shapely.geometry import Point, Polygon, MultiPolygon, GeometryCollection

def pca_eigenvectors(pts: np.ndarray) -> np.ndarray:
    Returns the principal axes of a set of points.
    Method is essentially running a PCA on the points.

    pts : array_like
    ca = np.cov(pts, y=None, rowvar=False, bias=True)
    val, vect = np.linalg.eig(ca)

    return np.transpose(vect)

def oriented_bounding_box(pts: np.ndarray) -> np.ndarray:
    Returns the oriented bounding box width set of points.

    Based on [Create the Oriented Bounding-box (OBB) with Python and NumPy](https://stackoverflow.com/questions/32892932/create-the-oriented-bounding-box-obb-with-python-and-numpy).

    pts : array_like
    tvect = pca_eigenvectors(pts)
    rot_matrix = np.linalg.inv(tvect)

    rot_arr = np.dot(pts, rot_matrix)

    mina = np.min(rot_arr, axis=0)
    maxa = np.max(rot_arr, axis=0)
    diff = (maxa - mina) * 0.5

    center = mina + diff

    half_w, half_h = diff
    corners = np.array([
        center + [-half_w, -half_h],
        center + [half_w, -half_h],
        center + [half_w, half_h],
        center + [-half_w, half_h],

    return np.dot(corners, tvect)

def polygon_from_obb(obb: np.ndarray) -> Polygon:
    Returns the oriented bounding box width set of points.

    obb : array_like
    obb = np.vstack((obb, obb[0]))
    return Polygon(obb)

This essentially performs a mini Principal Component Analysis to pull out the two perpendicular axes best describing the points, which can then be used to draw. It's all implemented with numpy for the most part so performant for most use cases.

I've pulled that example's implementation and added some utilities for working with shapely and geopandas here: https://github.com/raphaellaude/geo-obb/blob/main/geoobb/obb.py

See example usage: https://github.com/raphaellaude/geo-obb/blob/main/examples/parcel_obbs.ipynb

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.