diff --git a/geometry/rotating_calipers.py b/geometry/rotating_calipers.py new file mode 100644 index 000000000000..2209faf14af1 --- /dev/null +++ b/geometry/rotating_calipers.py @@ -0,0 +1,207 @@ +""" +Rotating Calipers Algorithm for Convex Polygon Diameter. + +References: +- https://en.wikipedia.org/wiki/Rotating_calipers +- https://cp-algorithms.com/geometry/convex-hull-kernel.html +- Toussaint, G. T. (1983). "Solving geometric problems with the rotating calipers". + Proceedings of IEEE MELECON '83, Athens, Greece. + +The rotating calipers paradigm allows computing the diameter (the maximum Euclidean +distance between any pair of points) of a set of 2D points in O(n log n) time +(O(n log n) for the convex hull and O(n) for the calipers sweep). +""" + +from __future__ import annotations + +import math +from typing import NamedTuple + + +class Point(NamedTuple): + """ + A 2D point with real-valued coordinates. + + >>> Point(0.0, 0.0) + Point(x=0.0, y=0.0) + >>> Point(1.5, -2.0) + Point(x=1.5, y=-2.0) + """ + + x: float + y: float + + +def cross_product(origin: Point, point_a: Point, point_b: Point) -> float: + """ + Compute the 2D cross product of vectors (origin -> point_a) and (origin -> point_b). + + The return value represents twice the signed area of triangle + (origin, point_a, point_b): + > 0 : Counter-clockwise turn (left turn) + < 0 : Clockwise turn (right turn) + = 0 : Collinear points + + >>> cross_product(Point(0.0, 0.0), Point(1.0, 0.0), Point(1.0, 1.0)) + 1.0 + >>> cross_product(Point(0.0, 0.0), Point(1.0, 1.0), Point(1.0, 0.0)) + -1.0 + >>> cross_product(Point(0.0, 0.0), Point(1.0, 1.0), Point(2.0, 2.0)) + 0.0 + """ + return (point_a.x - origin.x) * (point_b.y - origin.y) - (point_a.y - origin.y) * ( + point_b.x - origin.x + ) + + +def distance_squared(point_a: Point, point_b: Point) -> float: + """ + Compute the squared Euclidean distance between point_a and point_b. + + >>> distance_squared(Point(0.0, 0.0), Point(3.0, 4.0)) + 25.0 + >>> distance_squared(Point(1.0, 1.0), Point(1.0, 1.0)) + 0.0 + >>> distance_squared(Point(-1.0, -1.0), Point(2.0, 3.0)) + 25.0 + """ + return (point_a.x - point_b.x) ** 2 + (point_a.y - point_b.y) ** 2 + + +def convex_hull(points: list[Point]) -> list[Point]: + """ + Compute the convex hull of a set of 2D points in counter-clockwise order + using Andrew's monotone chain algorithm. + + Time Complexity: O(n log n) where n is the number of points. + Space Complexity: O(n) + + >>> convex_hull([Point(0.0, 0.0), Point(1.0, 1.0)]) + [Point(x=0.0, y=0.0), Point(x=1.0, y=1.0)] + >>> convex_hull([ + ... Point(0.0, 0.0), + ... Point(3.0, 0.0), + ... Point(3.0, 3.0), + ... Point(0.0, 3.0), + ... Point(1.0, 1.0), + ... ]) + [Point(x=0.0, y=0.0), Point(x=3.0, y=0.0), Point(x=3.0, y=3.0), Point(x=0.0, y=3.0)] + >>> convex_hull([Point(0.0, 0.0), Point(1.0, 1.0), Point(2.0, 2.0)]) + [Point(x=0.0, y=0.0), Point(x=2.0, y=2.0)] + >>> convex_hull([Point(1.0, 1.0)]) + [Point(x=1.0, y=1.0)] + """ + unique_points = sorted(set(points)) + if len(unique_points) <= 1: + return unique_points + + lower_hull: list[Point] = [] + for candidate_point in unique_points: + while ( + len(lower_hull) >= 2 + and cross_product(lower_hull[-2], lower_hull[-1], candidate_point) <= 0.0 + ): + lower_hull.pop() + lower_hull.append(candidate_point) + + upper_hull: list[Point] = [] + for candidate_point in reversed(unique_points): + while ( + len(upper_hull) >= 2 + and cross_product(upper_hull[-2], upper_hull[-1], candidate_point) <= 0.0 + ): + upper_hull.pop() + upper_hull.append(candidate_point) + + return lower_hull[:-1] + upper_hull[:-1] + + +def rotating_calipers(points: list[Point]) -> tuple[float, tuple[Point, Point]]: + """ + Find the maximum Euclidean distance (polygon diameter) and an antipodal pair + of points for a given set of 2D points using the rotating calipers algorithm. + + Time Complexity: O(n log n) for convex hull construction + + O(n) for the calipers sweep. + Space Complexity: O(n) for the convex hull. + + Raises: + ValueError: If fewer than 2 points are provided. + + >>> points = [ + ... Point(0.0, 0.0), + ... Point(3.0, 0.0), + ... Point(3.0, 4.0), + ... Point(0.0, 4.0), + ... ] + >>> max_dist, pair = rotating_calipers(points) + >>> max_dist + 5.0 + >>> pair in [ + ... (Point(0.0, 0.0), Point(3.0, 4.0)), + ... (Point(3.0, 4.0), Point(0.0, 0.0)), + ... (Point(3.0, 0.0), Point(0.0, 4.0)), + ... (Point(0.0, 4.0), Point(3.0, 0.0)), + ... ] + True + >>> rotating_calipers([Point(0.0, 0.0), Point(0.0, 5.0)]) + (5.0, (Point(x=0.0, y=0.0), Point(x=0.0, y=5.0))) + >>> rotating_calipers([Point(1.0, 1.0), Point(1.0, 1.0)]) + (0.0, (Point(x=1.0, y=1.0), Point(x=1.0, y=1.0))) + >>> rotating_calipers([ + ... Point(0.0, 0.0), + ... Point(1.0, 1.0), + ... Point(2.0, 2.0), + ... Point(3.0, 3.0), + ... ])[0] + 4.242640687119285 + >>> rotating_calipers([Point(1.0, 1.0)]) + Traceback (most recent call last): + ... + ValueError: At least 2 points are required to compute polygon diameter. + """ + if len(points) < 2: + raise ValueError("At least 2 points are required to compute polygon diameter.") + + hull = convex_hull(points) + hull_size = len(hull) + + if hull_size == 1: + return 0.0, (hull[0], hull[0]) + if hull_size == 2: + return math.hypot(hull[0].x - hull[1].x, hull[0].y - hull[1].y), ( + hull[0], + hull[1], + ) + + max_dist_squared = 0.0 + best_pair = (hull[0], hull[1]) + + # Find initial antipodal point furthest from edge hull[0]-hull[1] + antipodal_idx = 1 + while cross_product( + hull[0], hull[1], hull[(antipodal_idx + 1) % hull_size] + ) > cross_product(hull[0], hull[1], hull[antipodal_idx]): + antipodal_idx = (antipodal_idx + 1) % hull_size + + for current_idx in range(hull_size): + next_idx = (current_idx + 1) % hull_size + while cross_product( + hull[current_idx], hull[next_idx], hull[(antipodal_idx + 1) % hull_size] + ) > cross_product(hull[current_idx], hull[next_idx], hull[antipodal_idx]): + antipodal_idx = (antipodal_idx + 1) % hull_size + + for p in (hull[current_idx], hull[next_idx]): + for candidate_idx in (antipodal_idx, (antipodal_idx + 1) % hull_size): + dist_sq = distance_squared(p, hull[candidate_idx]) + if dist_sq > max_dist_squared: + max_dist_squared = dist_sq + best_pair = (p, hull[candidate_idx]) + + return math.sqrt(max_dist_squared), best_pair + + +if __name__ == "__main__": + import doctest + + doctest.testmod()