axiolid_reference/
convex_hull.rs

1//! Deterministic 2D convex hulls and oriented bounding rectangles.
2use crate::orient2d;
3use axiolid_contracts::{GeomError, GeomResult, Sign};
4use axiolid_core::{Point2, Rectangle2, Scalar};
5
6/// The oriented rectangle this module used to define itself.
7///
8/// Kept as an alias rather than a distinct type: it stored four corners plus a
9/// cached area and side pair, which is the same rectangle `Rectangle2` models
10/// as an origin and two edge vectors, and having two spellings of one concept
11/// forced every caller to convert. The edge-vector form is also the one that
12/// cannot drift -- four independently stored corners can be edited into a
13/// non-parallelogram, and the cached area can disagree with them.
14pub type OrientedRectangle2 = Rectangle2;
15
16/// Side lengths of an oriented rectangle, shortest edge first.
17///
18/// A free function rather than an inherent method because `Rectangle2` lives
19/// in `axiolid-core`, which holds data and no algorithms. Ordering is fixed so
20/// the result does not depend on which hull edge the caliper happened to stop
21/// on.
22#[must_use]
23pub fn side_lengths(rectangle: &Rectangle2) -> [Scalar; 2] {
24    let mut sides = [rectangle.x.length(), rectangle.y.length()];
25    sides.sort_by(|left, right| left.total_cmp(right));
26    sides
27}
28
29pub fn strict_convex_hull(points: &[Point2]) -> GeomResult<Vec<usize>> {
30    for (index, point) in points.iter().enumerate() {
31        if !point.is_finite() {
32            return Err(GeomError::InvalidInput(format!(
33                "point {index} is not finite"
34            )));
35        }
36    }
37    let mut ordered: Vec<usize> = (0..points.len()).collect();
38    ordered.sort_by(|&a, &b| {
39        points[a]
40            .x
41            .total_cmp(&points[b].x)
42            .then_with(|| points[a].y.total_cmp(&points[b].y))
43            .then(a.cmp(&b))
44    });
45    ordered.dedup_by(|a, b| points[*a] == points[*b]);
46    if ordered.len() < 3 {
47        return Err(GeomError::Degenerate("need three distinct points".into()));
48    }
49    let mut lower = Vec::new();
50    for &index in &ordered {
51        push_strict(&mut lower, index, points);
52    }
53    let mut upper = Vec::new();
54    for &index in ordered.iter().rev() {
55        push_strict(&mut upper, index, points);
56    }
57    lower.pop();
58    upper.pop();
59    lower.extend(upper);
60    if lower.len() < 3 {
61        return Err(GeomError::Degenerate("points are collinear".into()));
62    }
63    Ok(lower)
64}
65
66pub fn minimum_area_rectangle(points: &[Point2]) -> GeomResult<OrientedRectangle2> {
67    let hull = strict_convex_hull(points)?;
68    let mut best: Option<Rectangle2> = None;
69    // Tracked beside the rectangle because `Rectangle2` deliberately stores no
70    // cached area: recomputing it per candidate is a multiply, and a cached
71    // field is one more thing that can disagree with the geometry.
72    let mut best_area: Option<Scalar> = None;
73    for i in 0..hull.len() {
74        let a = points[hull[i]];
75        let b = points[hull[(i + 1) % hull.len()]];
76        let edge = b - a;
77        let width = edge.length();
78        let u = edge / width;
79        let v = Point2::new(-u.y, u.x);
80        let (mut ulo, mut uhi, mut vlo, mut vhi) = (
81            Scalar::INFINITY,
82            Scalar::NEG_INFINITY,
83            Scalar::INFINITY,
84            Scalar::NEG_INFINITY,
85        );
86        for &index in &hull {
87            let point = points[index];
88            let pu = point.dot(u);
89            let pv = point.dot(v);
90            ulo = ulo.min(pu);
91            uhi = uhi.max(pu);
92            vlo = vlo.min(pv);
93            vhi = vhi.max(pv);
94        }
95        let sides = [uhi - ulo, vhi - vlo];
96        let area = sides[0] * sides[1];
97        // Origin plus two edge vectors, rather than four corners: the
98        // parallelogram property then holds by construction instead of being
99        // an invariant four separately stored points could violate.
100        let rectangle = Rectangle2::new(u * ulo + v * vlo, u * sides[0], v * sides[1]);
101        if best_area.is_none_or(|current| area < current) {
102            best_area = Some(area);
103            best = Some(rectangle);
104        }
105    }
106    Ok(best.expect("non-empty strict hull has an edge"))
107}
108
109fn push_strict(hull: &mut Vec<usize>, index: usize, points: &[Point2]) {
110    while hull.len() >= 2 {
111        let n = hull.len();
112        if sign(orient2d(
113            points[hull[n - 2]],
114            points[hull[n - 1]],
115            points[index],
116        )) == Sign::Positive
117        {
118            break;
119        }
120        hull.pop();
121    }
122    hull.push(index);
123}
124fn sign(value: axiolid_contracts::Certified) -> Sign {
125    value.sign().expect("orient2d is total")
126}
axiolid_reference/convex_hull.rs

axiolid_reference/
convex_hull.rs
