Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Find the smallest rotated rectangle that covers a given polygon
import numpy as np
# Gift wrapping algorithm, from https://en.wikipedia.org/wiki/Gift_wrapping_algorithm
def find_convex_hull(poly):
"""
Returns the convex hull of the given polygon
Can be replaced using OpenCV:
```
hull = cv2.convexHull(poly, clockwise=True, returnPoints=False)
return poly[hull[:, 0]]
```
"""
n = len(poly)
p_i = p_0 = np.argmin(poly[:, 0])
hull = []
while p_i != p_0 or len(hull) == 0:
hull.append(p_i)
p_e = 0
# For each vertex, test if p_j is left of line between p_i <-> p_e
for p_j in range(n):
if p_e == p_i or np.cross(poly[p_e], poly[p_j] - poly[p_i]) > np.cross(poly[p_i], poly[p_j]):
p_e = p_j
p_i = p_e
return np.array([poly[k] for k in hull])
def find_min_rotated_rect(poly):
"""Returns the four points of the smallest rotated rect that covers the given poly"""
ch = find_convex_hull(poly)
edges = ch - np.roll(ch, 1, 0)
edges = edges / np.linalg.norm(edges, 2, -1)[:, None]
normals = np.stack((-edges[:, 1], edges[:, 0]), 1)
basis = np.stack((edges, normals), -1)
ps = np.matmul(ch, basis)
ps0 = np.amin(ps, 1)
ps1 = np.amax(ps, 1)
areas = np.prod(ps1 - ps0, -1)
k = np.argmin(areas)
rotated_rect = np.array(
[[ps0[k, 0], ps0[k, 1]],
[ps0[k, 0], ps1[k, 1]],
[ps1[k, 0], ps1[k, 1]],
[ps1[k, 0], ps0[k, 1]]])
return rotated_rect @ basis[k].T
@janhuenermann

This comment has been minimized.

Copy link
Owner Author

@janhuenermann janhuenermann commented May 19, 2021

Test cases

from matplotlib import pyplot as plt

# Only for visualisation
def plot_poly(poly):
    poly = np.concatenate([poly, poly[0:1]], 0)
    plt.plot(poly[:, 0], poly[:, 1])
# Test case 1
poly = np.array([[0., 1.], [0., 2.], [2., 2.], [4., 6.], [6., 3.], [4., 0.], [2., 0.]])
poly = np.roll(poly, 1, 0)
rect = find_min_rotated_rect(poly)
rect = norm_rect_with_ll_point(rect, poly)

plot_poly(poly)
plot_poly(rect)

plt.scatter(poly[0, 0], poly[0, 1])
plt.scatter(rect[0, 0], rect[0, 1])

plt.axis('scaled')
plt.grid('on')
plt.show()

out-1

# Test case 2
poly = np.array([[0.5, 0.5], [-0.5, -0.2], [0.5, -0.5], [1.5, 0.2], [0.5, 0.5]])
rect = find_min_rotated_rect(poly)
rect = norm_rect_with_ll_point(rect, poly)

plot_poly(poly)
plot_poly(rect)

plt.scatter(poly[0, 0], poly[0, 1])
plt.scatter(rect[0, 0], rect[0, 1])

plt.axis('scaled')
plt.grid('on')
plt.show()

out-2

@janhuenermann

This comment has been minimized.

Copy link
Owner Author

@janhuenermann janhuenermann commented May 19, 2021

Here's another test case

poly = np.zeros((n, 2))
for i in range(1, n):
    poly[i] = poly[i-1] + np.random.normal(size=(2,)) * 0.33
    
poly[:, 0] = 1. * (poly[:, 0] - np.min(poly[:, 0])) / (np.max(poly[:, 0]) - np.min(poly[:, 0])) - .5
poly[:, 1] = 1. * (poly[:, 1] - np.min(poly[:, 1])) / (np.max(poly[:, 1]) - np.min(poly[:, 1])) - .5

rect = find_min_rotated_rect(poly)
rect = norm_rect_with_ll_point(rect, poly)

plot_poly(poly)
plot_poly(rect)

plt.xlim(-1., 1.)
plt.ylim(-1., 1.)

plt.scatter(poly[0, 0], poly[0, 1])
plt.scatter(rect[0, 0], rect[0, 1])

plt.grid('on')
plt.show()

And an animated version of the outputs:

out

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment