osm-katana/src/polygon_mask.rs
//! **A raster acceleration mask for point-in-polygon** — O(1) per node instead of
//! O(vertices), with the polygon left exactly as it was.
//!
//! ## Why
//! A region clip is `nodes × vertices`. The planet is 10.58 billion nodes; Sweden's
//! real border (relation 52822) is 14 602 vertices where the downloaded outline it
//! replaced was 146. The bbox pre-filter only helps for nodes outside the box, so
//! every node that lands inside Sweden's bounding box pays a full ray-cast — and
//! the honest reading of that is a 100× cost increase for the fidelity gained.
//!
//! There are three ways to attack `nodes × vertices`. Simplifying the ring attacks
//! `vertices` — and deleting vertices from a coastline is exactly how the old
//! 48-vertex ring lost Gotland, Öland and 4.52 % of Sweden's land. Fanning the
//! test across cores attacks `nodes` — and it is **already done**: `clip_keeps`
//! runs inside the gatling pass-1 workers, measured at 29.9 of 32 cores busy.
//! This attacks the third factor, the cost of one test, and it is the only one of
//! the three that changes no geometry at all.
//!
//! ## How
//! Lay a grid over the polygon's bbox. Every cell is one of three things:
//!
//! * **Outside** — no edge crosses it and its interior is outside the polygon.
//! * **Inside** — no edge crosses it and its interior is inside.
//! * **Boundary** — at least one polygon edge passes through it.
//!
//! A query maps `(lon, lat)` to its cell in constant time. Outside and Inside
//! answer immediately; only a Boundary cell falls through to the real ray-cast.
//! Since no edge crosses a non-boundary cell, every point in it has the same
//! verdict, so answering from the cell IS the ray-cast's answer — the mask cannot
//! disagree with [`crate::Polygon::contains_exact`], and a test asserts that over
//! millions of real planet nodes rather than trusting the argument.
//!
//! ## Building it is a scanline, not a raster of ray-casts
//! Classifying each clean cell by ray-casting its centre would be `cells ×
//! vertices` — 1024² cells over 14 602 vertices is 1.5 × 10¹⁰ operations. Instead
//! each grid ROW is done as one even-odd scanline: intersect the row's centre
//! latitude with every edge once, sort the crossings, and every cell centre in
//! that row is classified by a binary search. That is `rows × vertices` — 1.5 × 10⁷
//! for the same grid, and it is the same even-odd rule `contains` uses, so the two
//! agree by construction rather than by luck.
//!
//! Rows are independent, so the build is a `gatling` fork-join (ROOT LAW #0 —
//! never rayon), mirroring the 192×192 coverage grid in korp-collectors'
//! `osm_tiles_verify` (which measures 733 ms over a 10 000-vertex Sweden polygon).
//! That one answers a different question — "does the tile pyramid cover the clip"
//! — and lives in a leaf crate that nothing may depend on, so the technique is
//! read from it, not imported. The shared primitive belongs here, next to the
//! `Polygon` it accelerates. The scanline build costs **2 ms** at 1024² over
//! Sweden's 14 602-vertex border, against that 733 ms — the difference is
//! per-cell ray-casts versus one sweep per row.
//!
//! ## Measured (oden, 2026-08-03, `benches/clip_containment`)
//! 52 931 795 REAL node coordinates sampled from the planet (stride 200), against
//! `korp/assets/sweden.poly` — relation 52822, 2 rings, 14 602 vertices. 2.273 %
//! of planet nodes land inside Sweden's bbox and therefore pay for a ray-cast.
//!
//! | | ns/node | boundary cells | vs exact | nodes moved |
//! |---|---|---|---|---|
//! | exact ray-cast | 264.6 | — | 1.0× | 0 |
//! | mask 256² (64 KiB) | 7.4 | 2.35 % | 35.6× | **0** |
//! | mask 512² (256 KiB) | 5.9 | 1.18 % | 44.7× | **0** |
//! | mask 1024² (1 MiB) | 5.2 | 0.60 % | 51.3× | **0** |
//! | mask 2048² (4 MiB) | 4.8 | 0.30 % | 54.7× | **0** |
//!
//! Extrapolated over the planet's 10.58 billion nodes across 31 gatling workers,
//! that is 91 s of clip wall time today against **1.8 s** at 1024² — about what
//! the 146-vertex DOWNLOADED outline this border replaced would have cost, at the
//! full border's fidelity. 1024² is the default because 2048² buys 0.4 ns for 4×
//! the memory.
//!
//! Simplification, for comparison, buys its speed by MOVING nodes: Douglas–Peucker
//! at `eps = 0.0001` is 2.5× faster and already gives 41 of those 52.9 M nodes a
//! different verdict; at `eps = 0.05` it is 57× faster, drops Treriksröset out of
//! Sweden, and moves 1 883. The mask is faster than any of them and moves none.
use crate::Bounds;
/// A cell that no polygon edge passes through, wholly outside the polygon.
const CELL_OUTSIDE: u8 = 0;
/// A cell that no polygon edge passes through, wholly inside the polygon.
const CELL_INSIDE: u8 = 1;
/// A cell at least one polygon edge passes through — the only cells that pay for
/// a real ray-cast.
const CELL_BOUNDARY: u8 = 2;
/// Raster acceleration mask over one [`crate::Polygon`]. See the module docs.
#[derive(Clone, Debug)]
pub struct PolygonMask {
/// `nx * ny` cell states, row-major, `ny` rows of `nx`.
cells: Vec<u8>,
nx: usize,
ny: usize,
/// The polygon's bbox — the mask's extent.
bbox: Bounds,
/// Cell size in degrees.
dlon: f64,
dlat: f64,
}
impl PolygonMask {
/// Build a `steps × steps` mask over `rings`' bbox.
///
/// `steps` is clamped to at least 1. A degenerate bbox (zero width or height)
/// yields a mask whose cells are all BOUNDARY, i.e. every query falls through
/// to the exact test — slower, never wrong.
pub fn build(rings: &[Vec<(f64, f64)>], bbox: Bounds, steps: usize) -> Self {
let nx = steps.max(1);
let ny = steps.max(1);
let width = bbox.max_lon - bbox.min_lon;
let height = bbox.max_lat - bbox.min_lat;
if !(width > 0.0) || !(height > 0.0) {
return Self {
cells: vec![CELL_BOUNDARY; nx * ny],
nx,
ny,
bbox,
dlon: f64::INFINITY,
dlat: f64::INFINITY,
};
}
let dlon = width / nx as f64;
let dlat = height / ny as f64;
// ── 1. Mark every cell any edge passes through. ───────────────────────
// Conservative by construction: a segment is walked cell by cell (a
// supercover DDA), so a cell an edge merely clips a corner of is still
// BOUNDARY. Over-marking costs a ray-cast; under-marking would be a wrong
// answer, so the walk always rounds outward.
let mut cells = vec![CELL_OUTSIDE; nx * ny];
for ring in rings {
for pair in ring.windows(2) {
let (x0, y0) = pair[0];
let (x1, y1) = pair[1];
mark_segment(&mut cells, nx, ny, &bbox, dlon, dlat, x0, y0, x1, y1);
}
// A ring given without its repeated closing vertex still closes.
if let (Some(&first), Some(&last)) = (ring.first(), ring.last()) {
if first != last {
mark_segment(
&mut cells, nx, ny, &bbox, dlon, dlat, last.0, last.1, first.0, first.1,
);
}
}
}
// ── 2. Classify the clean cells, one even-odd scanline per row. ───────
// Row-parallel: rows share nothing (each writes only its own slice), so
// this is a gatling fork-join, never rayon.
let row_states: Vec<Vec<u8>> = crate::par::par_map(ny, |row| {
let lat = bbox.min_lat + (row as f64 + 0.5) * dlat;
// Every edge crossing the row's centre latitude, as x positions —
// the same test `Polygon::contains` applies, hoisted out of the loop.
let mut xs: Vec<f64> = Vec::new();
for ring in rings {
let n = ring.len();
if n < 3 {
continue;
}
let mut j = n - 1;
for i in 0..n {
let (xi, yi) = ring[i];
let (xj, yj) = ring[j];
if (yi > lat) != (yj > lat) {
xs.push((xj - xi) * (lat - yi) / (yj - yi) + xi);
}
j = i;
}
}
xs.sort_unstable_by(|a, b| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal));
let base = row * nx;
let mut out = vec![CELL_OUTSIDE; nx];
for col in 0..nx {
if cells[base + col] == CELL_BOUNDARY {
out[col] = CELL_BOUNDARY;
continue;
}
let lon = bbox.min_lon + (col as f64 + 0.5) * dlon;
// Even-odd: an odd number of crossings strictly left of the point
// means inside — `Polygon::contains`'s rule, counted rather than
// toggled.
let crossings = xs.partition_point(|&x| x <= lon);
out[col] = if crossings % 2 == 1 {
CELL_INSIDE
} else {
CELL_OUTSIDE
};
}
out
});
for (row, states) in row_states.into_iter().enumerate() {
cells[row * nx..(row + 1) * nx].copy_from_slice(&states);
}
Self {
cells,
nx,
ny,
bbox,
dlon,
dlat,
}
}
/// Cells per axis.
pub fn steps(&self) -> usize {
self.nx
}
/// Mask size in bytes — one byte per cell.
pub fn bytes(&self) -> usize {
self.cells.len()
}
/// How many cells are BOUNDARY, INSIDE, OUTSIDE — the mask's whole story: the
/// boundary share is the fraction of the bbox where a query still pays for a
/// ray-cast.
pub fn cell_counts(&self) -> (usize, usize, usize) {
let mut b = 0;
let mut i = 0;
let mut o = 0;
for &c in &self.cells {
match c {
CELL_BOUNDARY => b += 1,
CELL_INSIDE => i += 1,
_ => o += 1,
}
}
(b, i, o)
}
/// The BOUNDARY cells' row-major indices (`row * steps + col`), ascending.
///
/// The mask's whole content in sparse form: at 1024² over Sweden's border
/// these are 0.60 % of the cells, and everything else is decided by sign. It
/// is what a sparse membership index — a sorted-key structure such as
/// `znippy_zoomies::stree` — would have to hold, so this is the honest input
/// to "could a search tree be the mask" rather than an estimate of it.
pub fn boundary_cells(&self) -> Vec<i64> {
self.cells
.iter()
.enumerate()
.filter(|&(_, &c)| c == CELL_BOUNDARY)
.map(|(i, _)| i as i64)
.collect()
}
/// The row-major cell index of `(lon, lat)`, or `None` outside the bbox — the
/// key a sparse index would be queried with.
#[inline]
pub fn cell_key(&self, lon: f64, lat: f64) -> Option<i64> {
if !self.bbox.contains(lon, lat) {
return None;
}
let cx = (((lon - self.bbox.min_lon) / self.dlon) as usize).min(self.nx.saturating_sub(1));
let cy = (((lat - self.bbox.min_lat) / self.dlat) as usize).min(self.ny.saturating_sub(1));
Some((cy * self.nx + cx) as i64)
}
/// Answer for a precomputed [`PolygonMask::cell_key`] — the same lookup with
/// the coordinate arithmetic hoisted out, so a dense array index and a sparse
/// key-search can be compared on identical inputs.
#[inline]
pub fn lookup_key(&self, key: i64) -> Option<bool> {
match self.cells.get(usize::try_from(key).ok()?) {
Some(&CELL_INSIDE) => Some(true),
Some(&CELL_OUTSIDE) => Some(false),
_ => None,
}
}
/// Answer for `(lon, lat)`, or `None` when the point lands in a BOUNDARY cell
/// and only the exact ray-cast can decide.
#[inline]
pub fn lookup(&self, lon: f64, lat: f64) -> Option<bool> {
if !self.bbox.contains(lon, lat) {
return Some(false);
}
let cx = ((lon - self.bbox.min_lon) / self.dlon) as usize;
let cy = ((lat - self.bbox.min_lat) / self.dlat) as usize;
// A point exactly on the max edge floors to `n`; clamp into the last cell,
// which is the cell that contains it.
let cx = cx.min(self.nx.saturating_sub(1));
let cy = cy.min(self.ny.saturating_sub(1));
match self.cells.get(cy * self.nx + cx) {
Some(&CELL_INSIDE) => Some(true),
Some(&CELL_OUTSIDE) => Some(false),
_ => None, // BOUNDARY (or out of range) → ray-cast
}
}
}
/// Mark every cell the segment `(x0,y0)–(x1,y1)` passes through as BOUNDARY.
///
/// A cell-walking (supercover) DDA: step to whichever of the next vertical or
/// horizontal cell boundary comes first, marking as it goes, so a segment that
/// only clips a cell's corner still marks it. The alternative — marking a
/// segment's bounding box — is also conservative but marks a whole rectangle for
/// a long diagonal maritime edge, which is most of Sweden's outer ring.
#[allow(clippy::too_many_arguments)]
fn mark_segment(
cells: &mut [u8],
nx: usize,
ny: usize,
bbox: &Bounds,
dlon: f64,
dlat: f64,
x0: f64,
y0: f64,
x1: f64,
y1: f64,
) {
let cell_of = |v: f64, min: f64, d: f64, n: usize| -> isize {
let i = ((v - min) / d).floor();
if !i.is_finite() {
return 0;
}
i.clamp(0.0, n.saturating_sub(1) as f64) as isize
};
let mut cx = cell_of(x0, bbox.min_lon, dlon, nx);
let mut cy = cell_of(y0, bbox.min_lat, dlat, ny);
let ex = cell_of(x1, bbox.min_lon, dlon, nx);
let ey = cell_of(y1, bbox.min_lat, dlat, ny);
let mark = |cells: &mut [u8], cx: isize, cy: isize| {
if cx >= 0 && cy >= 0 && (cx as usize) < nx && (cy as usize) < ny {
cells[cy as usize * nx + cx as usize] = CELL_BOUNDARY;
}
};
mark(cells, cx, cy);
if cx == ex && cy == ey {
return;
}
let dx = x1 - x0;
let dy = y1 - y0;
let step_x: isize = if dx > 0.0 {
1
} else if dx < 0.0 {
-1
} else {
0
};
let step_y: isize = if dy > 0.0 {
1
} else if dy < 0.0 {
-1
} else {
0
};
// Parametric distance (t in 0..1 along the segment) to the next cell border.
let next_border = |c: isize, step: isize, min: f64, d: f64| -> f64 {
if step > 0 {
min + (c + 1) as f64 * d
} else {
min + c as f64 * d
}
};
let mut t_max_x = if step_x == 0 {
f64::INFINITY
} else {
(next_border(cx, step_x, bbox.min_lon, dlon) - x0) / dx
};
let mut t_max_y = if step_y == 0 {
f64::INFINITY
} else {
(next_border(cy, step_y, bbox.min_lat, dlat) - y0) / dy
};
let t_delta_x = if step_x == 0 {
f64::INFINITY
} else {
(dlon / dx).abs()
};
let t_delta_y = if step_y == 0 {
f64::INFINITY
} else {
(dlat / dy).abs()
};
// Bounded: a segment crosses at most nx + ny cell borders, and the guard
// means a NaN or a pathological coordinate ends the walk instead of spinning.
let max_steps = nx.saturating_add(ny).saturating_add(2);
for _ in 0..max_steps {
if cx == ex && cy == ey {
break;
}
if t_max_x < t_max_y {
if !t_max_x.is_finite() {
break;
}
cx += step_x;
t_max_x += t_delta_x;
} else {
if !t_max_y.is_finite() {
break;
}
cy += step_y;
t_max_y += t_delta_y;
}
mark(cells, cx, cy);
}
mark(cells, ex, ey);
}
#[cfg(test)]
mod tests {
use super::*;
use crate::Polygon;
/// A ring with a deep inlet and a detached island — the shape a mask must not
/// smooth over (it is the Gotland case in miniature).
fn coastline() -> Polygon {
let mainland = vec![
(0.0, 0.0),
(10.0, 0.0),
(10.0, 10.0),
(6.0, 10.0),
(6.0, 2.0),
(4.0, 2.0),
(4.0, 10.0), // the inlet
(0.0, 10.0),
(0.0, 0.0),
];
let island = vec![
(12.0, 4.0),
(14.0, 4.0),
(14.0, 6.0),
(12.0, 6.0),
(12.0, 4.0),
];
Polygon::new(vec![mainland, island]).expect("rings")
}
/// **The mask must agree with the exact ray-cast everywhere.**
///
/// Swept over a dense grid of the bbox — including the inlet, the island and
/// the gap between them — so a mask that filled the inlet or swallowed the
/// island is red. This is the property the planet-scale test in
/// `tests/clip_mask_equivalence.rs` asserts over real node coordinates.
#[test]
fn mask_agrees_with_the_exact_raycast_over_the_whole_bbox() {
let poly = coastline();
for steps in [8usize, 32, 256] {
let mask = PolygonMask::build(poly.rings(), poly.bbox(), steps);
let mut answered = 0usize;
let mut fell_through = 0usize;
for iy in 0..=400 {
for ix in 0..=400 {
let lon = -1.0 + ix as f64 * 16.0 / 400.0;
let lat = -1.0 + iy as f64 * 12.0 / 400.0;
let exact = poly.contains_exact(lon, lat);
match mask.lookup(lon, lat) {
Some(v) => {
answered += 1;
assert_eq!(
v, exact,
"mask({steps}) says {v} at ({lon:.4}, {lat:.4}), exact says {exact}",
);
}
None => fell_through += 1,
}
}
}
assert!(
answered > 0,
"the mask must answer something at {steps} steps"
);
assert!(
fell_through > 0,
"a mask over a real coastline must have boundary cells at {steps} steps; \
if none fall through, the edge rasterisation did nothing",
);
}
}
/// The mask is only worth building if it ANSWERS most queries. A finer grid
/// must shrink the boundary share — that is the whole speed-up.
#[test]
fn finer_grids_shrink_the_boundary_share() {
let poly = coastline();
let share = |steps: usize| {
let m = PolygonMask::build(poly.rings(), poly.bbox(), steps);
let (b, _, _) = m.cell_counts();
b as f64 / (steps * steps) as f64
};
let coarse = share(16);
let fine = share(256);
assert!(
fine < coarse * 0.5,
"a 16× finer grid should more than halve the boundary share, got {coarse:.4} → {fine:.4}",
);
assert!(
fine < 0.10,
"at 256 steps most cells must be decidable, got {fine:.4}"
);
}
/// A degenerate (zero-area) bbox must not produce a mask that answers — every
/// query has to fall through to the exact test.
#[test]
fn degenerate_extent_falls_through_instead_of_guessing() {
let flat = vec![(1.0, 5.0), (3.0, 5.0), (2.0, 5.0), (1.0, 5.0)];
let bbox = Bounds::new(1.0, 5.0, 3.0, 5.0);
let mask = PolygonMask::build(&[flat], bbox, 16);
assert_eq!(mask.lookup(2.0, 5.0), None, "no cell may claim to know");
}
}