python - Algorithm to find the minimum-area-rectangle for given points in order to compute the major and minor axis length

Question

I have a set of points (black dots in geographic coordinate value) derived from the convex hull (blue) of a polygon (red). see Figure:

[(560023.44957588764,6362057.3904932579), 
 (560023.44957588764,6362060.3904932579), 
 (560024.44957588764,6362063.3904932579), 
 (560026.94957588764,6362068.3904932579), 
 (560028.44957588764,6362069.8904932579), 
 (560034.94957588764,6362071.8904932579), 
 (560036.44957588764,6362071.8904932579), 
 (560037.44957588764,6362070.3904932579), 
 (560037.44957588764,6362064.8904932579), 
 (560036.44957588764,6362063.3904932579), 
 (560034.94957588764,6362061.3904932579), 
 (560026.94957588764,6362057.8904932579), 
 (560025.44957588764,6362057.3904932579), 
 (560023.44957588764,6362057.3904932579)]

I need to calculate the the major and minor axis length following these steps (form this post write in R-project and in Java) or following the this example procedure

1. Compute the convex hull of the cloud.
2. For each edge of the convex hull: 2a. compute the edge orientation, 2b. 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, 2c. Store the orientation corresponding to the minimum area found,
3. Return the rectangle corresponding to the minimum area found.

After that we know the The angle Theta (represented the orientation of the bounding rectangle relative to the y-axis of the image). The minimum and maximum of a and b over all boundary points are found:

• a(xi,yi) = xi*cos Theta + yi sin Theta
• b(xi,yi) = xi*sin Theta + yi cos Theta

The values (a_max - a_min) and (b_max - b_min) defined the length and width, respectively, of the bounding rectangle for a direction Theta.

Answer 1

I just implemented this myself, so 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.

(edit by @heltonbiker) A simple code for plotting:

import matplotlib.pyplot as plt
for n in range(10):
    points = np.random.rand(4,2)
    plt.scatter(points[:,0], points[:,1])
    bbox = minimum_bounding_rectangle(points)
    plt.fill(bbox[:,0], bbox[:,1], alpha=0.2)
    plt.axis('equal')
    plt.show()

(end edit)

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

Link to the same answer over on gis.stackexchange for my own reference.