Algorithm to find the minimum-area-rectangle for given points in order to compute the major and minor axis length
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.
Given a clockwise-ordered list of n points in the convex hull of a set of points, it is an O(n) operation to find the minimum-area enclosing rectangle. (For convex-hull finding, in O(n log n) time, see activestate.com recipe 66527 or see the quite compact Graham scan code at tixxit.net.)
The following python program uses techniques similar to those of the usual O(n) algorithm for computing maximum diameter of a convex polygon. That is, it maintains three indexes (iL, iP, iR) to the leftmost, opposite, and rightmost points relative to a given baseline. Each index advances through at most n points. Sample output from the program is shown next (with an added header):
i iL iP iR Area
0 6 8 0 203.000
1 6 8 0 211.875
2 6 8 0 205.800
3 6 10 0 206.250
4 7 12 0 190.362
5 8 0 1 203.000
6 10 0 4 201.385
7 0 1 6 203.000
8 0 3 6 205.827
9 0 3 6 205.640
10 0 4 7 187.451
11 0 4 7 189.750
12 1 6 8 203.000
For example, the i=10 entry indicates that relative to the baseline from point 10 to 11, point 0 is leftmost, point 4 is opposite, and point 7 is rightmost, yielding an area of 187.451 units.
Note that the code uses mostfar() to advance each index. The mx, my parameters to mostfar() tell it what extreme to test for; as an example, with mx,my = -1,0, mostfar() will try to maximize -rx (where rx is the rotated x of a point), thus finding the leftmost point. Note that an epsilon allowance probably should be used when if mx*rx + my*ry >= best is done in inexact arithmetic: when a hull has numerous points, rounding error may be a problem and cause the method to incorrectly not advance an index.
Code is shown below. The hull data is taken from the question above, with irrelevant large offsets and identical decimal places elided.
#!/usr/bin/python
import math
hull = [(23.45, 57.39), (23.45, 60.39), (24.45, 63.39),
(26.95, 68.39), (28.45, 69.89), (34.95, 71.89),
(36.45, 71.89), (37.45, 70.39), (37.45, 64.89),
(36.45, 63.39), (34.95, 61.39), (26.95, 57.89),
(25.45, 57.39), (23.45, 57.39)]
def mostfar(j, n, s, c, mx, my): # advance j to extreme point
xn, yn = hull[j][0], hull[j][1]
rx, ry = xn*c - yn*s, xn*s + yn*c
best = mx*rx + my*ry
while True:
x, y = rx, ry
xn, yn = hull[(j+1)%n][0], hull[(j+1)%n][1]
rx, ry = xn*c - yn*s, xn*s + yn*c
if mx*rx + my*ry >= best:
j = (j+1)%n
best = mx*rx + my*ry
else:
return (x, y, j)
n = len(hull)
iL = iR = iP = 1 # indexes left, right, opposite
pi = 4*math.atan(1)
for i in range(n-1):
dx = hull[i+1][0] - hull[i][0]
dy = hull[i+1][1] - hull[i][1]
theta = pi-math.atan2(dy, dx)
s, c = math.sin(theta), math.cos(theta)
yC = hull[i][0]*s + hull[i][1]*c
xP, yP, iP = mostfar(iP, n, s, c, 0, 1)
if i==0: iR = iP
xR, yR, iR = mostfar(iR, n, s, c, 1, 0)
xL, yL, iL = mostfar(iL, n, s, c, -1, 0)
area = (yP-yC)*(xR-xL)
print ' {:2d} {:2d} {:2d} {:2d} {:9.3f}'.format(i, iL, iP, iR, area)
Note: To get the length and width of the minimal-area enclosing rectangle, modify the above code as shown below. This will produce an output line like
Min rectangle: 187.451 18.037 10.393 10 0 4 7
in which the second and third numbers indicate the length and width of the rectangle, and the four integers give index numbers of points that lie upon sides of it.
# add after pi = ... line:
minRect = (1e33, 0, 0, 0, 0, 0, 0) # area, dx, dy, i, iL, iP, iR
# add after area = ... line:
if area < minRect[0]:
minRect = (area, xR-xL, yP-yC, i, iL, iP, iR)
# add after print ... line:
print 'Min rectangle:', minRect
# or instead of that print, add:
print 'Min rectangle: ',
for x in ['{:3d} '.format(x) if isinstance(x, int) else '{:7.3f} '.format(x) for x in minRect]:
print x,
print
There is a module that does this already on github. https://github.com/BebeSparkelSparkel/MinimumBoundingBox
All you need to do is insert your point cloud into it.
from MinimumBoundingBox import minimum_bounding_box
points = ( (1,2), (5,4), (-1,-3) )
bounding_box = minimum_bounding_box(points) # returns namedtuple
You can get the major and minor axis lengths by:
minor = min(bounding_box.length_parallel, bounding_box.length_orthogonal)
major = max(bounding_box.length_parallel, bounding_box.length_orthogonal)
It also returns area, rectangle center, rectangle angle, and corner points.