geometry_algorithm/
minimum_rotated_rect.rs

1//! Minimum-area rotated rectangle around a planar geometry.
2//!
3//! Boost.Geometry has no rotated-envelope entry. This implementation follows
4//! Toussaint's rotating-calipers construction (1983): compute the convex hull,
5//! align an orthogonal frame with every hull edge, and retain the frame with
6//! minimum projected area.
7
8use alloc::vec;
9
10#[cfg(not(feature = "std"))]
11use geometry_coords::math::Float;
12use geometry_model::{Polygon, Ring};
13use geometry_strategy::{ConvexHullStrategy, MonotoneChain};
14use geometry_trait::{Point, PointMut};
15
16use crate::convex_hull::convex_hull;
17
18/// Compute the minimum-area rotated rectangle enclosing `geometry`.
19///
20/// This is Cartesian-only through the convex-hull strategy bound. Empty input
21/// returns an empty polygon; one- and two-point inputs return a closed,
22/// zero-area degenerate rectangle.
23#[inline]
24#[must_use]
25pub fn minimum_rotated_rect<G, P>(geometry: &G) -> Polygon<P>
26where
27    P: Point<Scalar = f64> + PointMut + Default + Copy,
28    MonotoneChain: ConvexHullStrategy<G, Output = Ring<P, true, true>>,
29{
30    let hull = convex_hull(geometry);
31    let mut points = hull.0;
32    while points.len() > 1 && same_xy(points.first(), points.last()) {
33        points.pop();
34    }
35    match points.len() {
36        0 => return Polygon::default(),
37        1 => {
38            let point = points[0];
39            return Polygon::new(Ring::from_vec(vec![point, point, point, point, point]));
40        }
41        2 => {
42            let first = points[0];
43            let second = points[1];
44            return Polygon::new(Ring::from_vec(vec![first, first, second, second, first]));
45        }
46        _ => {}
47    }
48
49    let mut best: Option<Frame> = None;
50    for index in 0..points.len() {
51        let first = points[index];
52        let second = points[(index + 1) % points.len()];
53        let dx = second.get::<0>() - first.get::<0>();
54        let dy = second.get::<1>() - first.get::<1>();
55        let length = dx.hypot(dy);
56        if length <= f64::EPSILON {
57            continue;
58        }
59        let ux = dx / length;
60        let uy = dy / length;
61        let vx = -uy;
62        let vy = ux;
63        let mut min_u = f64::INFINITY;
64        let mut max_u = f64::NEG_INFINITY;
65        let mut min_v = f64::INFINITY;
66        let mut max_v = f64::NEG_INFINITY;
67        for point in &points {
68            let x = point.get::<0>();
69            let y = point.get::<1>();
70            let along = x * ux + y * uy;
71            let across = x * vx + y * vy;
72            min_u = min_u.min(along);
73            max_u = max_u.max(along);
74            min_v = min_v.min(across);
75            max_v = max_v.max(across);
76        }
77        let frame = Frame {
78            ux,
79            uy,
80            vx,
81            vy,
82            min_u,
83            max_u,
84            min_v,
85            max_v,
86        };
87        if best.is_none_or(|current| frame.area() < current.area()) {
88            best = Some(frame);
89        }
90    }
91
92    let Some(frame) = best else {
93        return Polygon::default();
94    };
95    let lower_left = frame.point::<P>(frame.min_u, frame.min_v);
96    let upper_left = frame.point::<P>(frame.min_u, frame.max_v);
97    let upper_right = frame.point::<P>(frame.max_u, frame.max_v);
98    let lower_right = frame.point::<P>(frame.max_u, frame.min_v);
99    Polygon::new(Ring::from_vec(vec![
100        lower_left,
101        upper_left,
102        upper_right,
103        lower_right,
104        lower_left,
105    ]))
106}
107
108#[derive(Clone, Copy)]
109struct Frame {
110    ux: f64,
111    uy: f64,
112    vx: f64,
113    vy: f64,
114    min_u: f64,
115    max_u: f64,
116    min_v: f64,
117    max_v: f64,
118}
119
120impl Frame {
121    fn area(self) -> f64 {
122        (self.max_u - self.min_u) * (self.max_v - self.min_v)
123    }
124
125    fn point<P>(self, along: f64, across: f64) -> P
126    where
127        P: Point<Scalar = f64> + PointMut + Default,
128    {
129        let mut point = P::default();
130        point.set::<0>(along * self.ux + across * self.vx);
131        point.set::<1>(along * self.uy + across * self.vy);
132        point
133    }
134}
135
136#[allow(
137    clippy::float_cmp,
138    reason = "coordinate identity is used only to detect the closing duplicate"
139)]
140fn same_xy<P: Point<Scalar = f64>>(first: Option<&P>, second: Option<&P>) -> bool {
141    first.zip(second).is_some_and(|(first, second)| {
142        first.get::<0>() == second.get::<0>() && first.get::<1>() == second.get::<1>()
143    })
144}
145
146#[cfg(test)]
147mod tests {
148    use geometry_cs::Cartesian;
149    use geometry_model::{MultiPoint, Point2D};
150
151    use super::minimum_rotated_rect;
152    use crate::area::area;
153
154    #[test]
155    fn diamond_has_area_two() {
156        type P = Point2D<f64, Cartesian>;
157        let points = MultiPoint::from_vec(alloc::vec![
158            P::new(0.0, 1.0),
159            P::new(1.0, 0.0),
160            P::new(0.0, -1.0),
161            P::new(-1.0, 0.0),
162        ]);
163        let rectangle = minimum_rotated_rect(&points);
164        assert!((area(&rectangle).abs() - 2.0).abs() < 1e-12);
165    }
166}
geometry_algorithm/minimum_rotated_rect.rs

geometry_algorithm/
minimum_rotated_rect.rs
