osm-katana/src/shared/generalise.rs
//! **Generalisation — the ONE place a WKB blob becomes a coarser WKB blob (LAW 5).**
//!
//! `geometry_z6`, `geometry_z9`, `geometry_z12` are Douglas–Peucker run at the rung's
//! epsilon and written into a column. Nothing here is a new simplifier: the DP kernel,
//! its scratch buffers, the per-ring area rule and the epsilon are all
//! [`crate::tiles`]'s, imported rather than re-derived.
//!
//! # Why the columns exist at all — the number, not the intuition
//!
//! MEASURED on oden 2026-08-31, packed `ways.parquet`, central Stockholm
//! `(17.85, 59.22)–(18.30, 59.45)`, warm page cache. The row-group prune is already
//! excellent — 71 of 1 139 groups, 6.23 % — and it costs **18.5 ms / 4.21 MB**. The
//! rest of the 144.0 ms read, **125.5 ms**, is zstd + Arrow decode of 52.38 MB of
//! FULL-PRECISION WKB that a z12 camera cannot resolve. No row-group tuning closes a
//! decode; a narrower column does.
//!
//! # The two rules this module inherits rather than invents
//!
//! 1. **Per RING, never over the flat vertex list.** Concatenated rings put a chord
//! between the outer ring's last vertex and the first hole's first vertex, and DP
//! then measures deviation from an edge the geometry does not have — keeping or
//! dropping real vertices on the strength of it — and the result used to be
//! re-encoded as a `LineString`, so the hole stopped being a hole.
//! [`crate::tiles::simplify_polys`] already paid for that lesson; this module calls
//! it.
//! 2. **Sub-pixel ⇒ a POINT, not a deletion.** A feature whose every ring simplifies
//! below four vertices is genuinely sub-pixel at that rung, and the tiler degrades
//! it to its first vertex rather than dropping the row. Same rule here — a
//! generalisation column must have a cell for every row the table has **that the
//! rung admits**, or a reader projecting it loses features the base column still
//! holds. The admission is a different question, answered by a different column:
//! since DELTA 3 `spatial::pack_batches` writes NULL into `geometry_z<z>` for a row
//! whose `min_zoom > z` — the ladder withholds it at that camera, so there is
//! nothing to degrade. This module never sees those rows; the packer skips them
//! before calling in. MEASURED on the Sweden `ways.parquet` (2026-09-01): the
//! 5 609 616 bare `Point`s the z6 column used to carry were almost all service roads
//! and footpaths no z6 painter draws.
//!
//! # And ONE rule it tightens
//!
//! The tile scan flattens a `MultiLineString`'s runs into one vertex list before DP
//! (`shape.vertices()`), which is defect (1) in its line form: a chord across the seam
//! between two disjoint runs. Here every run is simplified **on its own**, for exactly
//! the reason the ring case gives. `a_ferry_route_keeps_both_of_its_runs` is the guard.
use crate::shared::geometry::{
WkbShape, decode_wkb, encode_wkb_line_into, encode_wkb_multilinestring_into,
encode_wkb_multipolygon_into, push_wkb_point,
};
use crate::tiles::{
DEFAULT_TILE_RASTER_PX, DpScratch, douglas_peucker_into, simplify_polys, simplify_polys_into,
};
/// **The Douglas–Peucker epsilon for generalisation rung `z`, in degrees.**
///
/// This is [`crate::tiles::tile_pixel_tol`] at [`DEFAULT_TILE_RASTER_PX`] and nothing
/// else — the SAME tolerance the pyramid was cut with, so a column read at z12 shows
/// the shape a z12 tile showed. A second constant here would be a second simplification
/// basis, and the two would disagree the first time `DEFAULT_TILE_RASTER_PX` moved (it
/// already moved once, 256 → 2048, and the reason is written up on that constant).
#[must_use]
pub fn generalisation_eps(z: u8) -> f64 {
crate::tiles::tile_pixel_tol(z, DEFAULT_TILE_RASTER_PX)
}
/// Simplify a decoded shape at `eps`, ring-wise and run-wise. `None` only when the
/// shape carries no vertex at all.
///
/// A `Point` comes back unchanged: a point has no vertices to drop, which is also why
/// a POINT TABLE gets no generalisation columns at all (see [`table_generalises`]).
#[must_use]
pub fn generalise_shape(shape: &WkbShape, eps: f64) -> Option<WkbShape> {
match shape {
WkbShape::Point(p) => Some(WkbShape::Point(*p)),
WkbShape::Lines(runs) => {
let mut scratch = DpScratch::default();
let mut thin = Vec::new();
let mut out: Vec<Vec<(f64, f64)>> = Vec::with_capacity(runs.len());
for r in runs {
douglas_peucker_into(r, eps, &mut scratch, &mut thin);
if thin.len() >= 2 {
out.push(std::mem::take(&mut thin));
}
}
if out.is_empty() {
// Every run collapsed to one vertex: the feature IS a point at this
// rung. Degraded, not dropped — the same rule the area case follows.
return shape.vertices().next().map(WkbShape::Point);
}
Some(WkbShape::Lines(out))
}
WkbShape::Polys(ps) => {
let s = simplify_polys(ps, eps);
if s.is_empty() {
return shape.vertices().next().map(WkbShape::Point);
}
Some(WkbShape::Polys(s))
}
}
}
/// **Generalise one WKB blob into a caller-owned buffer** — the hot path (LAW 5; the
/// allocating [`generalise_wkb`] is a wrapper over it).
///
/// `out` is cleared first. Returns `false` — leaving `out` empty — when the blob does
/// not decode or has no vertices, and the caller must then append a NULL rather than an
/// empty cell, so "I could not read this geometry" never wears the same bytes as "this
/// feature has no shape here".
///
/// `scratch` is the DP keep-bitmap and index stack, held per worker for the whole pass:
/// allocating them per row per rung is the `tilemem-2` defect (11.13 GiB of glibc arena
/// free list, measured) in its column form, and there are three rungs per row.
pub fn generalise_wkb_into(
wkb: &[u8],
eps: f64,
scratch: &mut DpScratch,
thin: &mut Vec<(f64, f64)>,
out: &mut Vec<u8>,
) -> bool {
out.clear();
let Some(shape) = decode_wkb(wkb) else {
return false;
};
match &shape {
WkbShape::Point((x, y)) => {
push_wkb_point(out, *x, *y);
true
}
WkbShape::Lines(runs) => {
let mut kept: Vec<Vec<(f64, f64)>> = Vec::with_capacity(runs.len());
for r in runs {
douglas_peucker_into(r, eps, scratch, thin);
if thin.len() >= 2 {
kept.push(std::mem::take(thin));
}
}
match kept.len() {
// Nothing survived as a line ⇒ the feature is a point at this rung.
0 => match shape.vertices().next() {
Some((x, y)) => {
push_wkb_point(out, x, y);
true
}
None => false,
},
// A single run is a LineString, not a one-member MultiLineString —
// the same encoding the base column uses for the same geometry.
1 => {
encode_wkb_line_into(&kept[0], out);
true
}
_ => {
encode_wkb_multilinestring_into(&kept, out);
true
}
}
}
WkbShape::Polys(ps) => {
// Into the caller's pool (row 24). `simplify_polys` — the allocating form —
// was the pack's hot path until this row: a fresh output ring per ring, a
// fresh Douglas–Peucker `Vec` per RETRY, and a clone of the source ring on
// give-up, once per row per rung. The counts land in `scratch.books`, which
// the caller folds once per unit.
simplify_polys_into(ps, eps, scratch);
if scratch.polys.is_empty() {
return match shape.vertices().next() {
Some((x, y)) => {
push_wkb_point(out, x, y);
true
}
None => false,
};
}
encode_wkb_multipolygon_into(scratch.polys.as_slice(), out);
true
}
}
}
/// [`generalise_wkb_into`] with its buffers allocated for you. Tests and one-shot
/// callers only — the pack pass lends buffers.
#[must_use]
pub fn generalise_wkb(wkb: &[u8], eps: f64) -> Option<Vec<u8>> {
let (mut s, mut t, mut o) = (DpScratch::default(), Vec::new(), Vec::new());
generalise_wkb_into(wkb, eps, &mut s, &mut t, &mut o).then_some(o)
}
/// **Does this table's geometry generalise at all?** Read off its `geo` document's
/// declared `geometry_types`.
///
/// `false` for a POINT-only table, and that is not a micro-optimisation. MEASURED on
/// the Sweden cut: `nodes.parquet` holds 102 518 979 points at 655.18 MB of compressed
/// `geometry`; three generalisation columns on it would be **three byte-for-byte copies
/// of that column**, +1.97 GB of disk for a shape no epsilon can change. An absent
/// column is honest (the same argument `ways_schema` makes about `node_refs`); a
/// present duplicate claims a coarser geometry exists when it does not.
///
/// A reader therefore asks the SCHEMA for `geometry_z12` and falls back to `geometry`
/// when it is absent — which is the fallback it needs anyway for a table written before
/// this landed.
#[must_use]
pub fn table_generalises(geometry_types: &[String]) -> bool {
geometry_types
.iter()
.any(|t| !matches!(t.as_str(), "Point" | "MultiPoint"))
}
/// The `geometry_types` a generalisation column may hold.
///
/// Wider than the base column's on purpose, and in both directions: an area that
/// simplifies below a ring becomes a `Point` (the degrade rule above), and every
/// surviving area is re-encoded through the ONE ring encoder, which emits
/// `MultiPolygon`. Declaring only the base column's types here would be a `geo`
/// document that its own data contradicts.
pub const GENERALISED_GEOMETRY_TYPES: &[&str] = &[
"Point",
"LineString",
"MultiLineString",
"Polygon",
"MultiPolygon",
];
/// **The min-area rule: how many rendered pixels an AREA must span to earn a rung.**
///
/// The degrade rule above ("sub-pixel ⇒ a POINT, not a deletion") is right for a row
/// the ladder admits — but the LADDER decides admission from tags alone
/// (`shared::class::min_zoom_of_tags`), and tags cannot see size. MEASURED on the live
/// Sweden cut, 2026-09-01, the Gävle-coast box `(16.8, 60.2)–(17.6, 61.0)` of
/// The owner's country-zoom screenshot: **407 of 454 `natural=water` ways (89.6 %) sat on
/// `min_zoom = 6` and degraded to a bare `Point` in `geometry_z6`** (median extent
/// 0.00195° ≈ 130 m against the rung's 0.00275° epsilon) — and korp paints a classified
/// water point as a 2.5 px dot, so four hundred dots in a ~60 px box render as the
/// solid dark BLOB the archipelago became. The relations leg adds 56 more in the same
/// box.
///
/// Two pixels, because that is the ring-survival threshold in practice: a ring whose
/// extent is under ~2 × epsilon cannot keep 4 vertices through Douglas–Peucker (the
/// probe above: point-degraded p90 = 0.00586° ≈ 2.1 ε, surviving-poly min = 0.00415° ≈
/// 1.5 ε). At or above it, the feature draws as the area it is.
pub const AREA_MIN_EXTENT_PIXELS: f64 = 2.0;
/// **The first rung a feature of `extent_deg` may draw as an AREA** — the size half of
/// the `min_zoom` decision, beside [`crate::shared::class::min_zoom_of_tags`]'s tag
/// half. The packer takes the MAX of the two for polygon rows (never a coarser value:
/// this deepens, it cannot promote), so a sub-pixel islet waits for the rung where it
/// is a shape instead of degrading to country-view confetti.
///
/// Walks the SAME ladder the rung columns are cut on ([`crate::shared::class::RUNGS`])
/// with the SAME epsilon ([`generalisation_eps`]) — a second tolerance here would be a
/// second simplification basis, which is the drift this module's header forbids. An
/// extent too small for even the deepest rung still answers the deepest rung: past the
/// ladder there is full-precision `geometry`, where every feature draws.
#[must_use]
pub fn area_rung_floor(extent_deg: f64) -> u8 {
for &z in &crate::shared::class::RUNGS {
if extent_deg >= AREA_MIN_EXTENT_PIXELS * generalisation_eps(z) {
return z;
}
}
*crate::shared::class::RUNGS.last().expect("the rung ladder is never empty")
}
#[cfg(test)]
mod tests {
use super::*;
use crate::shared::geometry::{encode_wkb_line, encode_wkb_multipolygon};
/// A closed ring, `n` vertices around a circle of radius `r` degrees.
fn ring(cx: f64, cy: f64, r: f64, n: usize) -> Vec<(f64, f64)> {
let mut v: Vec<(f64, f64)> = (0..n)
.map(|i| {
#[allow(clippy::cast_precision_loss)]
let a = std::f64::consts::TAU * i as f64 / n as f64;
(cx + r * a.cos(), cy + r * a.sin())
})
.collect();
v.push(v[0]);
v
}
/// ★ **The min-area ladder: a sub-pixel area waits for the rung where it is a
/// shape** (`waterblob-1`, the Gävle-coast dot-blob).
///
/// Driven with the PROBE's own numbers (2026-09-01, live Sweden cut): the median
/// point-degraded islet (extent 0.00195°) must NOT get rung 6 — it must land on 9,
/// where its extent clears 2 × ε(9) = 0.000687° and survives as a ring — while a
/// real lake (0.02°) keeps rung 6 and a house-sized pond (0.0001°) waits for 12.
/// The floor saturates at the deepest rung: past it sits full-precision
/// `geometry`, where everything draws.
#[test]
fn the_area_floor_walks_the_rung_ladder_by_size() {
assert_eq!(area_rung_floor(0.02), 6, "a real lake earns the country rung");
assert_eq!(area_rung_floor(0.00195), 9, "the median blob islet waits for z9");
assert_eq!(area_rung_floor(0.0001), 12, "a pond waits for z12");
assert_eq!(area_rung_floor(1e-6), 15, "sub-pixel everywhere ⇒ the deepest rung");
// The threshold really is 2 × ε — one epsilon to either side flips the answer.
let eps6 = generalisation_eps(6);
assert_eq!(area_rung_floor(AREA_MIN_EXTENT_PIXELS * eps6), 6);
assert_ne!(area_rung_floor(AREA_MIN_EXTENT_PIXELS * eps6 - eps6 * 0.1), 6);
}
/// ★ **RED-when-broken: the rungs must actually THIN, and thin monotonically.**
///
/// A Douglas–Peucker that drops everything is the fastest DP there is, and one that
/// drops nothing is the safest — both would make the byte gate below read fine for
/// the wrong reason. So this asserts a BAND at every rung: strictly fewer vertices
/// than full precision, strictly more than nothing, and never more vertices at a
/// coarser rung than at a finer one.
///
/// The input is a real Swedish shape: 4 096 vertices along the Bohuslän coast run,
/// with the amplitude of the wiggle chosen so z12 can resolve it and z6 cannot.
#[test]
fn every_rung_thins_and_the_ladder_is_monotone() {
// ~1.4° of coast at 58.5 N, wiggling ±0.002° (~120 m) — visible at z12
// (eps 4.3e-5°), invisible at z6 (eps 2.7e-3°).
let coast: Vec<(f64, f64)> = (0..4096)
.map(|i| {
#[allow(clippy::cast_precision_loss)]
let t = i as f64 / 4096.0;
(
11.2 + t * 1.4,
58.5 + t * 0.9 + 0.002 * (t * 220.0).sin(),
)
})
.collect();
let full = encode_wkb_line(&coast);
let mut prev = coast.len();
let mut seen = 0usize;
for z in crate::shared::schema::GENERALISED_RUNGS.iter().rev() {
let g = generalise_wkb(&full, generalisation_eps(*z)).expect("coast generalises");
let n = match decode_wkb(&g).expect("generalised blob decodes") {
WkbShape::Lines(r) => r.iter().map(Vec::len).sum::<usize>(),
WkbShape::Point(_) => 1,
WkbShape::Polys(_) => panic!("a line became an area"),
};
assert!(
n < coast.len(),
"z{z} kept {n} of {} vertices — the simplifier did not run, and a byte \
gate over this column would be measuring the base column",
coast.len()
);
assert!(
n >= 2,
"z{z} kept {n} vertices — a DP that drops everything is the fastest DP \
there is, and it is not simplification"
);
assert!(
n <= prev,
"z{z} kept {n} vertices but the FINER rung kept {prev} — a coarser rung \
can never keep more"
);
prev = n;
seen += 1;
}
assert_eq!(seen, 3, "the rung ladder is empty; this test asserted nothing");
assert!(
prev < 512,
"the coarsest rung still keeps {prev} of 4 096 vertices — the epsilon is not \
reaching this geometry and the column will not be smaller"
);
}
/// ★ **RED-when-broken: a lake keeps its island, at every rung.**
///
/// The shape of `tiles.rs`'s `a_lake_keeps_its_island_through_the_tile_writer`,
/// pointed at the column writer. A 0.4°-wide lake with a 0.1° island — both far
/// above z6's ~2.7e-3° epsilon, so NEITHER may vanish at any rung.
///
/// **The red arm** is written into the assertion messages: flatten the rings (run DP
/// over `shape.vertices()` and re-encode as a `LineString`, which is what this code
/// replaced) and the result decodes as `Lines`, not `Polys` — the first assert
/// fires. Keep the polygon but drop holes and the inner-ring count goes to 0 — the
/// second fires. Neither can pass by accident.
#[test]
fn a_lake_keeps_its_island_all_the_way_into_the_written_column() {
let lake = vec![(ring(15.0, 58.0, 0.2, 64), vec![ring(15.0, 58.0, 0.05, 32)])];
let wkb = encode_wkb_multipolygon(&lake);
let mut checked = 0usize;
for z in crate::shared::schema::GENERALISED_RUNGS {
let g = generalise_wkb(&wkb, generalisation_eps(z)).expect("the lake generalises");
let shape = decode_wkb(&g).expect("the generalised lake decodes");
let WkbShape::Polys(ps) = &shape else {
panic!(
"z{z}: the lake came back as {shape:?} — not an area at all. A `Lines` \
here is the flat ring-concatenating simplifier: DP took a chord from \
the shore's last vertex to the island's first and the hole stopped \
being a hole."
);
};
assert_eq!(ps.len(), 1, "z{z}: the lake split into {} polygons", ps.len());
assert_eq!(
ps[0].1.len(),
1,
"z{z}: the island is gone — {} inner rings, expected 1. A 0.1° hole is \
{:.0}x this rung's epsilon; nothing about the tolerance can explain \
dropping it.",
ps[0].1.len(),
0.1 / generalisation_eps(z)
);
assert!(
ps[0].1[0].len() >= 4,
"z{z}: the island survived as {} vertices — a ring needs 4 to enclose \
area once closed, so this is a hole that bounds nothing",
ps[0].1[0].len()
);
checked += 1;
}
assert_eq!(checked, 3, "no rung was examined — the loop body never ran");
}
/// ★ **RED-when-broken: two disjoint runs are simplified SEPARATELY.**
///
/// The line-shaped form of the ring lesson, and the one place this module is
/// stricter than `tiles.rs`: the tile scan flattens `shape.vertices()` before DP, so
/// a two-run `MultiLineString` gets a chord across the gap between the runs.
///
/// Two ferry legs with a 0.5° gap between them, each a shallow arc.
///
/// **The damage a flat DP does here is not a vertex count** — that was this test's
/// first control, and it went red on its own arm: flattening these two runs happens
/// to keep exactly as many vertices (34) as simplifying them separately, so a count
/// assertion would have passed for a path that produces the wrong geometry. Kept as
/// a note because the obvious control is the wrong one.
///
/// What flattening actually produces is **one run containing a segment that is not
/// in the input**: the 0.5° seam from the end of leg 1 to the start of leg 2, which
/// draws as a line across open water. So that is what is measured — the longest
/// consecutive-vertex step in the output. Per-run it cannot exceed a leg's own 0.2°
/// span; flat it must be at least the gap.
#[test]
fn a_ferry_route_keeps_both_of_its_runs() {
let leg = |x0: f64| -> Vec<(f64, f64)> {
(0..256)
.map(|i| {
#[allow(clippy::cast_precision_loss)]
let t = i as f64 / 255.0;
(x0 + t * 0.2, 59.0 + 0.01 * (t * std::f64::consts::PI).sin())
})
.collect()
};
let runs = vec![leg(18.0), leg(18.7)];
let mut wkb = Vec::new();
encode_wkb_multilinestring_into(&runs, &mut wkb);
let g = generalise_wkb(&wkb, generalisation_eps(12)).expect("the route generalises");
let WkbShape::Lines(out) = decode_wkb(&g).expect("decodes") else {
panic!("a multi-run line came back as something else");
};
assert_eq!(
out.len(),
2,
"the two legs came back as {} run(s) — they were concatenated before DP and \
the 0.5 deg chord between them is now part of the geometry",
out.len()
);
// The longest step in a vertex run — the seam, if there is one.
let longest = |r: &[(f64, f64)]| -> f64 {
r.windows(2)
.map(|w| (w[1].0 - w[0].0).hypot(w[1].1 - w[0].1))
.fold(0.0_f64, f64::max)
};
let ours = out.iter().map(|r| longest(r)).fold(0.0_f64, f64::max);
assert!(
ours < 0.3,
"the longest step in the generalised route is {ours:.3} deg — longer than a \
leg (0.2 deg), so a seam across the gap is in the output after all"
);
// The control: the FLAT arm, on the same input, through the same DP. Without a
// control that is measurably worse, the assertions above prove nothing.
let flat: Vec<(f64, f64)> = runs.iter().flatten().copied().collect();
let mut s = DpScratch::default();
let mut thin = Vec::new();
douglas_peucker_into(&flat, generalisation_eps(12), &mut s, &mut thin);
let flat_step = longest(&thin);
assert!(
flat_step >= 0.5,
"the flat arm's longest step is {flat_step:.3} deg — it did NOT draw the \
0.5 deg seam, so this test is not measuring the defect it is named for"
);
assert!(
flat_step > ours * 2.0,
"flat {flat_step:.3} deg vs per-run {ours:.3} deg — the two arms agree, and \
a green here would be the geometry, not the code"
);
}
/// ★ **RED-when-broken: a sub-pixel area degrades to a POINT, it does not vanish.**
///
/// A 20 m building is below z6's ~300 m epsilon by two orders of magnitude. The
/// column must still have a cell for it: a generalisation column with fewer non-null
/// rows than the table has rows silently loses features from the coarse view, which
/// is the pyramid defect (87.9 % of z6 areas shipped as `Point`) reappearing as a
/// NULL instead of a point.
#[test]
fn a_subpixel_building_degrades_to_a_point_and_is_not_dropped() {
let b = vec![(ring(18.07, 59.33, 0.0001, 8), vec![])];
let wkb = encode_wkb_multipolygon(&b);
let g = generalise_wkb(&wkb, generalisation_eps(6)).expect("the building has a cell");
match decode_wkb(&g).expect("decodes") {
WkbShape::Point(_) => {}
other => panic!("a 20 m building at z6 came back as {other:?}, not a Point"),
}
// And the control: at z12 the same building is still sub-pixel, but at FULL
// precision it is an area — so the degrade is the epsilon's doing, not the
// encoder's.
assert!(
matches!(
decode_wkb(&generalise_wkb(&wkb, 0.0).unwrap()).unwrap(),
WkbShape::Polys(_)
),
"eps 0 also collapsed the building — the degrade is unconditional, so the \
assertion above proves nothing about the tolerance"
);
}
/// ★ **RED-when-broken: a point table declares no generalisation columns.**
#[test]
fn a_point_only_table_does_not_generalise() {
assert!(
!table_generalises(&[String::from("Point")]),
"a point table would be given three byte-for-byte copies of its geometry \
column — +1.97 GB on the Sweden nodes table for nothing"
);
assert!(
!table_generalises(&[String::from("Point"), String::from("MultiPoint")]),
"MultiPoint has no vertex DP can drop either"
);
assert!(
table_generalises(&[String::from("LineString"), String::from("Polygon")]),
"the WAY table must generalise — without this arm the negative asserts above \
would pass for a function that always returns false"
);
assert!(
!table_generalises(&[]),
"an empty type list is not evidence of extent"
);
}
}