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}