{{ message }}

Instantly share code, notes, and snippets.

# janhuenermann/min-rotated-rect.py

Last active May 21, 2021
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 commented May 19, 2021 • edited

 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()``` ### janhuenermann commented May 19, 2021 • edited

 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: 