Find the smallest rotated rectangle that covers a given polygon

min-rotated-rect.py

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

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

# 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()

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: