Skip to main content

map_tile_toolkit/
geom.rs

1//! Exact integer geometry predicates shared by clipping, ring orientation, and point-in-ring tests.
2//!
3//! Every predicate is exact across the full `i32` coordinate range and never panics (no floats, no
4//! division). A 2-D cross product of full-`i32` points can reach `2^65`, so exactness needs 128-bit
5//! arithmetic *somewhere*; following [`i_overlay`]/`i_float` — which widen `i32 → i64` and
6//! cross-multiply the edge *difference* vectors within a bounded coordinate domain — [`orient`]
7//! evaluates the cross in `i64` on the difference vectors and only widens to `i128` when a difference
8//! does not fit `i32` (the rare near-full-`i32`-span case). Ordinary tile-local geometry therefore
9//! never pays the 128-bit width, while the crate keeps its full-`i32` contract.
10//!
11//! [`i_overlay`]: https://crates.io/crates/i_overlay
12
13use core::cmp::Ordering;
14
15use geo_types::Coord;
16
17use crate::vertex::Vertex;
18
19/// Orientation of the turn `a → b → c`: the sign of the cross product `(b − a) × (c − a)`.
20///
21/// [`Ordering::Greater`] is a left turn (counter-clockwise), [`Ordering::Less`] a right turn
22/// (clockwise), and [`Ordering::Equal`] means the three points are collinear.
23#[inline]
24pub(crate) fn orient(a: Coord<i32>, b: Coord<i32>, c: Coord<i32>) -> Ordering {
25    // Edge vectors as `i64` — an `i32 − i32` difference always fits `i64`.
26    let ex1 = i64::from(b.x) - i64::from(a.x);
27    let ey1 = i64::from(b.y) - i64::from(a.y);
28    let ex2 = i64::from(c.x) - i64::from(a.x);
29    let ey2 = i64::from(c.y) - i64::from(a.y);
30    // Fast path: when every component fits `i32`, each product is `≤ (2^31 − 1)^2 < 2^62`, so the two
31    // `i64` products compare exactly with no risk of overflow. Otherwise widen the two products to
32    // `i128`. Comparing the products (rather than subtracting) avoids a possible `i64` overflow in the
33    // difference and needs no wider type on the fast path.
34    if fits_i32(ex1) && fits_i32(ey1) && fits_i32(ex2) && fits_i32(ey2) {
35        (ex1 * ey2).cmp(&(ey1 * ex2))
36    } else {
37        (i128::from(ex1) * i128::from(ey2)).cmp(&(i128::from(ey1) * i128::from(ex2)))
38    }
39}
40
41/// Whether an `i64` value fits back into an `i32`.
42#[inline]
43fn fits_i32(v: i64) -> bool {
44    i32::try_from(v).is_ok()
45}
46
47/// The smallest integer `≥ x`, where `x` is the point at which segment `a`–`b` crosses the horizontal
48/// line at height `y` (`a.y ≠ b.y`, `y` between them). Exact: `x = a.x + dx·t/dy` is evaluated as a
49/// rounded-up integer quotient in `i64` when the factors fit `i32`, widening to `i128` otherwise —
50/// the same fast path / fallback split as [`orient`]. The result lies between `a.x` and `b.x`.
51#[expect(
52    clippy::cast_possible_truncation,
53    reason = "the offset is dx·t/dy with |t| ≤ |dy|, so it fits i64"
54)]
55pub(crate) fn crossing_x_ceil(a: Coord<i32>, b: Coord<i32>, y: i32) -> i64 {
56    let dx = i64::from(b.x) - i64::from(a.x);
57    let mut dy = i64::from(b.y) - i64::from(a.y);
58    let mut t = i64::from(y) - i64::from(a.y);
59    // A positive divisor lets `div_euclid` floor; ceil(n/d) = −floor(−n/d).
60    if dy < 0 {
61        dy = -dy;
62        t = -t;
63    }
64    let offset = if fits_i32(dx) && fits_i32(t) {
65        -(-(dx * t)).div_euclid(dy)
66    } else {
67        (-(-(i128::from(dx) * i128::from(t))).div_euclid(i128::from(dy))) as i64
68    };
69    i64::from(a.x) + offset
70}
71
72/// Twice the signed area of `ring` (the shoelace sum `Σ xᵢ·yᵢ₊₁ − xᵢ₊₁·yᵢ`), exact for any `i32`
73/// ring. Works whether or not the ring repeats its first vertex, and is `0` for a degenerate ring.
74///
75/// The sign is the ring's winding: **positive** turns counter-clockwise with `y` up — which is
76/// *clockwise on screen* with `y` down, i.e. an MVT exterior ring in tile coordinates — and
77/// **negative** the other way. The slicer preserves each input ring's winding (synthetic corners
78/// follow it), so a caller that needs a fixed convention (MVT: exterior positive, holes negative)
79/// normalizes once per feature before slicing, reversing any ring whose sign disagrees with its role.
80///
81/// Accumulated in `i128`: each term is a cross product of two `i32`-difference vectors (`< 2^65`), so
82/// the sum cannot overflow below `2^62` vertices.
83#[must_use]
84pub fn signed_area_2x<V: Vertex>(ring: &[V]) -> i128 {
85    let Some(v0) = ring.first().map(Vertex::position) else {
86        return 0;
87    };
88    // Measured about the first vertex, so the closing edge `v_{n-1} → v0` and the opening edge
89    // `v0 → v1` both contribute zero: iterating `windows(2)` gives the full doubled area whether or not
90    // the ring repeats its first vertex, and smaller magnitudes keep the terms well inside `i128`.
91    let mut area2: i128 = 0;
92    for w in ring.windows(2) {
93        let (p, q) = (w[0].position(), w[1].position());
94        let px = i128::from(p.x) - i128::from(v0.x);
95        let py = i128::from(p.y) - i128::from(v0.y);
96        let qx = i128::from(q.x) - i128::from(v0.x);
97        let qy = i128::from(q.y) - i128::from(v0.y);
98        area2 += px * qy - py * qx;
99    }
100    area2
101}
102
103/// The winding of a closed `ring` from the sign of [`signed_area_2x`]: [`Ordering::Greater`] for
104/// counter-clockwise (y up), [`Ordering::Less`] for clockwise, [`Ordering::Equal`] for a degenerate
105/// ring (zero area).
106pub(crate) fn ring_orientation<V: Vertex>(ring: &[V]) -> Ordering {
107    signed_area_2x(ring).cmp(&0)
108}
109
110/// Whether point `p` lies inside `ring` (a closed polygon ring), counting the boundary as **inside**.
111///
112/// Even-odd (crossing-number) ray cast along `+x`, evaluated with [`orient`] so there is no division
113/// and it is exact for all `i32` inputs. A half-open rule on each edge's `y` span (`[y_lo, y_hi)`)
114/// counts every crossing once, so a ray grazing a shared vertex is handled consistently. Points on an
115/// edge or vertex return `true`. The ring is treated cyclically; a repeated closing vertex is fine.
116pub(crate) fn point_in_ring(p: Coord<i32>, ring: &[Coord<i32>]) -> bool {
117    let n = ring.len();
118    if n < 3 {
119        return false;
120    }
121    // Drop a repeated closing vertex so the cyclic walk visits each edge once.
122    let m = if ring[n - 1] == ring[0] { n - 1 } else { n };
123    if m < 3 {
124        return false;
125    }
126    let mut inside = false;
127    let mut prev = m - 1;
128    for i in 0..m {
129        let tail = ring[prev];
130        let head = ring[i];
131        // On-boundary check: collinear and within the edge's bounding box → treat as inside.
132        if orient(tail, head, p) == Ordering::Equal
133            && p.x >= tail.x.min(head.x)
134            && p.x <= tail.x.max(head.x)
135            && p.y >= tail.y.min(head.y)
136            && p.y <= tail.y.max(head.y)
137        {
138            return true;
139        }
140        // Does the edge straddle the horizontal ray at `p.y`? Half-open: exactly one endpoint is
141        // strictly above `p.y`.
142        if (tail.y > p.y) != (head.y > p.y) {
143            // The crossing is to the right of `p` iff `p` is on the correct side of the directed edge.
144            // For an upward edge a left turn (`Greater`) puts `p` left of it → the ray crosses to
145            // `p`'s right; a downward edge flips the sense.
146            let left = orient(tail, head, p) == Ordering::Greater;
147            if left == (head.y > tail.y) {
148                inside = !inside;
149            }
150        }
151        prev = i;
152    }
153    inside
154}
155
156#[cfg(test)]
157mod tests {
158    use geo_types::Coord;
159
160    use super::*;
161
162    fn c(x: i32, y: i32) -> Coord<i32> {
163        Coord { x, y }
164    }
165
166    #[test]
167    fn degenerate_ring_contains_nothing() {
168        // Fewer than three vertices, or only two distinct ones once the closing copy is dropped.
169        assert!(!point_in_ring(c(0, 0), &[c(0, 0), c(5, 5)]));
170        assert!(!point_in_ring(c(0, 0), &[c(0, 0), c(5, 5), c(0, 0)]));
171    }
172
173    #[test]
174    fn orient_basic_and_extremes() {
175        assert_eq!(orient(c(0, 0), c(1, 0), c(0, 1)), Ordering::Greater); // CCW / left
176        assert_eq!(orient(c(0, 0), c(1, 0), c(0, -1)), Ordering::Less); // CW / right
177        assert_eq!(orient(c(0, 0), c(2, 2), c(1, 1)), Ordering::Equal); // collinear
178        // Full-`i32`-span vectors take the `i128` fallback and must stay exact.
179        assert_eq!(
180            orient(
181                c(i32::MIN, i32::MIN),
182                c(i32::MAX, i32::MIN),
183                c(i32::MIN, i32::MAX)
184            ),
185            Ordering::Greater
186        );
187        assert_eq!(
188            orient(
189                c(i32::MIN, i32::MIN),
190                c(i32::MAX, i32::MAX),
191                c(i32::MIN + 1, i32::MIN + 1)
192            ),
193            Ordering::Equal
194        );
195    }
196
197    #[test]
198    fn crossing_x_ceil_is_exact() {
199        // x = 0 + 10·(1/3) = 3.33… → 4; reversed direction gives the same point.
200        assert_eq!(crossing_x_ceil(c(0, 0), c(10, 3), 1), 4);
201        assert_eq!(crossing_x_ceil(c(10, 3), c(0, 0), 1), 4);
202        // Exact integer crossings stay put, negative offsets round toward +∞.
203        assert_eq!(crossing_x_ceil(c(0, 0), c(10, 2), 1), 5);
204        assert_eq!(crossing_x_ceil(c(0, 0), c(-10, 3), 1), -3);
205        // Full-`i32` spans take the `i128` path.
206        let x = crossing_x_ceil(c(i32::MIN, i32::MIN), c(i32::MAX, i32::MAX), 0);
207        assert_eq!(x, 0);
208        // The anti-diagonal `x + y = −1` crosses `y = 1` exactly at `x = −2`; one unit to the right
209        // of an exact crossing must not round further.
210        let x = crossing_x_ceil(c(i32::MIN, i32::MAX), c(i32::MAX, i32::MIN), 1);
211        assert_eq!(x, -2);
212        let x = crossing_x_ceil(c(i32::MIN, i32::MAX), c(i32::MAX - 1, i32::MIN), 1);
213        assert_eq!(x, -2); // the slightly steeper line crosses just left of −2
214    }
215
216    #[test]
217    fn ring_orientation_ccw_cw() {
218        let ccw = [c(0, 0), c(4, 0), c(4, 4), c(0, 4), c(0, 0)];
219        let cw = [c(0, 0), c(0, 4), c(4, 4), c(4, 0), c(0, 0)];
220        assert_eq!(ring_orientation(&ccw), Ordering::Greater);
221        assert_eq!(ring_orientation(&cw), Ordering::Less);
222        assert_eq!(ring_orientation(&ccw[..2]), Ordering::Equal); // degenerate
223        // Open and closed rings agree, and the magnitude is twice the area.
224        assert_eq!(signed_area_2x(&ccw), 32);
225        assert_eq!(signed_area_2x(&ccw[..4]), 32);
226        assert_eq!(signed_area_2x(&cw), -32);
227        assert_eq!(signed_area_2x::<Coord<i32>>(&[]), 0);
228        // Full-`i32` extremes stay exact.
229        let big = [
230            c(i32::MIN, i32::MIN),
231            c(i32::MAX, i32::MIN),
232            c(i32::MAX, i32::MAX),
233        ];
234        let side = i128::from(i32::MAX) - i128::from(i32::MIN);
235        assert_eq!(signed_area_2x(&big), side * side);
236    }
237
238    #[test]
239    fn point_in_ring_square() {
240        let sq = [c(0, 0), c(10, 0), c(10, 10), c(0, 10), c(0, 0)];
241        assert!(point_in_ring(c(5, 5), &sq)); // interior
242        assert!(!point_in_ring(c(15, 5), &sq)); // exterior
243        assert!(point_in_ring(c(0, 5), &sq)); // on the left edge → boundary counts as inside
244        assert!(point_in_ring(c(10, 10), &sq)); // corner vertex
245        assert!(!point_in_ring(c(-1, 5), &sq));
246    }
247
248    #[test]
249    fn point_in_ring_concave() {
250        // A C-shape (concave): the notch on the right is outside.
251        let c_shape = [
252            c(0, 0),
253            c(10, 0),
254            c(10, 3),
255            c(3, 3),
256            c(3, 7),
257            c(10, 7),
258            c(10, 10),
259            c(0, 10),
260            c(0, 0),
261        ];
262        assert!(point_in_ring(c(1, 5), &c_shape)); // in the spine
263        assert!(!point_in_ring(c(6, 5), &c_shape)); // in the notch → outside
264        assert!(point_in_ring(c(6, 1), &c_shape)); // in the bottom arm
265    }
266}