facett-core/src/engine/cull.rs
//! **The cull index** — Ragnar's S+tree, not a Hilbert R-tree.
//!
//! Gemini's spec proposed a packed Hilbert R-tree for the CPU-side slicing.
//! Rejected: the constellation already owns a **measured** better structure. The
//! S(+)-tree in `znippy_zoomies::stree32` is Ragnar Groot Koerkamp's cache-optimal
//! implicit tree — one 64-byte cache line per node, a branchless AVX2
//! `count_lt`, and a software-pipelined prefetch batch traversal. Its published
//! headline is **40× faster than binary search**, and it is already a dependency
//! of `facett-core` (the `gatling` edge), so adopting it costs zero new deps.
//!
//! ## How a *range* query becomes a *key* query
//! The S+tree answers "where does key `q` sit in this sorted `u32` array?"
//! ([`lower_bound`](znippy_zoomies::stree32::STree32::lower_bound)). A viewport
//! query is turned into that shape the way every tile pyramid already does it:
//!
//! 1. Quantise each feature's centre into a `2^k × 2^k` lattice over the dataset
//! bbox and interleave the bits — a **Morton key**.
//! 2. Sort features by Morton key once, at build time. Now spatial locality is
//! array locality.
//! 3. A viewport covers `rows × cols` lattice cells. For **each row**, the cells
//! `x0..=x1` lie between `morton(x0, y)` and `morton(x1, y)` — a contiguous
//! key interval. Two `lower_bound`s per row give a contiguous slice.
//! 4. Those slices are a **superset** (Morton rows interleave), so each survivor
//! gets one exact AABB test. `candidates` records exactly how many — the
//! anti-lie field, the same discipline as
//! [`NearestHit::candidates`](crate::legibility::NearestHit::candidates).
//!
//! The lattice level `k` is chosen **per query** so the viewport spans at most
//! [`MAX_ROWS`] rows. That is not a tuning constant dressed up as a bound: it is
//! the same "pick the pyramid level that matches the zoom" decision the
//! [`Hierarchy`](super::source::Hierarchy) makes, which is the point — *drill-down
//! and LOD are the same operation*.
//!
//! ## What this does NOT yet do
//! Features are indexed by their **centre**. A feature far wider than a lattice
//! cell (a motorway, a country boundary) can therefore be missed by a viewport
//! that clips it without containing its centre. [`CullIndex::build`] detects that
//! case and records it in [`CullIndex::oversize`]; the engine currently draws
//! oversize features unconditionally, which is correct but not yet cheap. The
//! fix is a per-level index (each feature entered at the level where it fits) —
//! the same `Bands` structure, applied to the culler.
use super::source::{WorldAabb, WorldPos};
use znippy_zoomies::stree32::STree32;
/// The most lattice rows one viewport query will decompose into. Above this the
/// lattice coarsens (fewer, larger cells) so the query count stays bounded and
/// the superset grows instead — the classic index trade, made explicit.
pub const MAX_ROWS: u32 = 64;
/// What a [`CullIndex::query`] found, and **how much work it did**.
#[derive(Clone, Debug, Default, PartialEq)]
pub struct CullReport {
/// Features returned.
pub visible: usize,
/// Features that survived the Morton slice and got an exact AABB test. A
/// linear scan reports `candidates == len`; an honest index reports far less.
pub candidates: usize,
/// Contiguous key intervals the viewport decomposed into (2 `lower_bound`s each).
pub ranges: usize,
/// Oversize features drawn without a spatial test (see the module note).
pub oversize_forced: usize,
/// The lattice level the query chose — the LOD/zoom band, made observable.
pub level: u32,
}
impl CullReport {
/// `state_json` fragment.
#[must_use]
pub fn state_json(&self) -> serde_json::Value {
serde_json::json!({
"visible": self.visible,
"candidates": self.candidates,
"ranges": self.ranges,
"oversize_forced": self.oversize_forced,
"level": self.level,
})
}
}
/// A static, Morton-ordered, S+tree-backed spatial index over world-space AABBs.
pub struct CullIndex {
/// The dataset bbox the lattice is normalised against.
bbox: WorldAabb,
/// Morton keys, ascending. Index `i` here is `order[i]` in caller id space.
keys: Vec<u32>,
/// Caller feature ids, in Morton order.
order: Vec<u32>,
/// AABBs in Morton order, for the exact re-test.
boxes: Vec<WorldAabb>,
/// `None` for an empty dataset — `STree32::new` asserts on an empty slice,
/// and an empty layer is a normal frame (a pane before its data arrives),
/// not a crash.
tree: Option<STree32>,
/// Features wider than one lattice cell at the finest level — see module doc.
oversize: Vec<u32>,
/// Bits per axis in the stored keys (`LEVEL_BITS`).
bits: u32,
}
// `STree32` is a packed SIMD array with no `Debug`; the index reports its shape
// rather than its bytes (which are never useful in a failure message anyway).
impl std::fmt::Debug for CullIndex {
fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
f.debug_struct("CullIndex")
.field("len", &self.order.len())
.field("oversize", &self.oversize.len())
.field("bits", &self.bits)
.finish()
}
}
/// Bits per axis. `16 + 16 = 32` — exactly a `u32` Morton key, exactly the
/// S+tree's native width.
const LEVEL_BITS: u32 = 16;
#[inline]
fn part1by1(mut v: u32) -> u32 {
v &= 0x0000_ffff;
v = (v ^ (v << 8)) & 0x00ff_00ff;
v = (v ^ (v << 4)) & 0x0f0f_0f0f;
v = (v ^ (v << 2)) & 0x3333_3333;
v = (v ^ (v << 1)) & 0x5555_5555;
v
}
/// Interleave two 16-bit lattice coordinates into a 32-bit Morton key.
#[inline]
#[must_use]
pub fn morton(x: u32, y: u32) -> u32 {
part1by1(x) | (part1by1(y) << 1)
}
impl CullIndex {
/// Build over `(id, aabb)` pairs. `O(n log n)` once; queries are then
/// `O(rows · log n + candidates)`.
#[must_use]
pub fn build(items: &[(u32, WorldAabb)]) -> Self {
let mut bbox = WorldAabb::empty();
for (_, b) in items {
if b.min.x.is_finite() && b.min.y.is_finite() {
bbox = bbox.union(*b);
}
}
let n = 1u32 << LEVEL_BITS;
let (sx, sy) = Self::scales(&bbox, n);
let mut rows: Vec<(u32, u32, WorldAabb)> = Vec::with_capacity(items.len());
let mut oversize = Vec::new();
// A cell edge at the FINEST level, in world units.
let cell_w = 1.0 / sx.max(f64::MIN_POSITIVE);
let cell_h = 1.0 / sy.max(f64::MIN_POSITIVE);
for &(id, b) in items {
if !b.min.x.is_finite() || !b.min.y.is_finite() {
continue;
}
let c = b.centre();
let gx = (((c.x - bbox.min.x) * sx) as i64).clamp(0, (n - 1) as i64) as u32;
let gy = (((c.y - bbox.min.y) * sy) as i64).clamp(0, (n - 1) as i64) as u32;
rows.push((morton(gx, gy), id, b));
if (b.max.x - b.min.x) > cell_w * (MAX_ROWS as f64) || (b.max.y - b.min.y) > cell_h * (MAX_ROWS as f64) {
oversize.push(id);
}
}
rows.sort_unstable_by_key(|r| r.0);
let keys: Vec<u32> = rows.iter().map(|r| r.0).collect();
let order: Vec<u32> = rows.iter().map(|r| r.1).collect();
let boxes: Vec<WorldAabb> = rows.iter().map(|r| r.2).collect();
let tree = (!keys.is_empty()).then(|| STree32::new(&keys));
Self { bbox, keys, order, boxes, tree, oversize, bits: LEVEL_BITS }
}
fn scales(bbox: &WorldAabb, n: u32) -> (f64, f64) {
let span_x = (bbox.max.x - bbox.min.x).max(f64::MIN_POSITIVE);
let span_y = (bbox.max.y - bbox.min.y).max(f64::MIN_POSITIVE);
(n as f64 / span_x, n as f64 / span_y)
}
#[must_use]
pub fn len(&self) -> usize {
self.order.len()
}
#[must_use]
pub fn is_empty(&self) -> bool {
self.order.is_empty()
}
/// Features too wide for the lattice — drawn unconditionally (see module doc).
#[must_use]
pub fn oversize(&self) -> &[u32] {
&self.oversize
}
/// The dataset bbox the lattice normalises against.
#[must_use]
pub fn bbox(&self) -> WorldAabb {
self.bbox
}
/// **Every feature whose AABB overlaps `view`**, ascending by Morton order,
/// plus the work report. Exact: no false positives, no false negatives (for
/// non-oversize features).
pub fn query(&self, view: WorldAabb, out: &mut Vec<u32>) -> CullReport {
out.clear();
let mut rep = CullReport { level: self.bits, ..Default::default() };
if self.order.is_empty() || view.is_empty() {
return rep;
}
let full = 1u32 << self.bits;
let (sx, sy) = Self::scales(&self.bbox, full);
// Fine lattice cell range covered by the viewport.
let clampi = |v: i64| v.clamp(0, (full - 1) as i64) as u32;
let fx0 = clampi(((view.min.x - self.bbox.min.x) * sx).floor() as i64);
let fx1 = clampi(((view.max.x - self.bbox.min.x) * sx).ceil() as i64);
let fy0 = clampi(((view.min.y - self.bbox.min.y) * sy).floor() as i64);
let fy1 = clampi(((view.max.y - self.bbox.min.y) * sy).ceil() as i64);
// Coarsen until the row count fits the budget. Shifting BOTH axes by
// `shift` keeps the Morton interleave valid: a coarse cell (gx>>s, gy>>s)
// is exactly the contiguous key block `[m<<2s, ((m+1)<<2s)-1]`.
let mut shift = 0u32;
while ((fy1 >> shift) - (fy0 >> shift) + 1) > MAX_ROWS && shift < self.bits {
shift += 1;
}
rep.level = self.bits - shift;
let (cx0, cx1) = (fx0 >> shift, fx1 >> shift);
let (cy0, cy1) = (fy0 >> shift, fy1 >> shift);
let block = 2 * shift; // key bits consumed by one coarse cell
let Some(tree) = self.tree.as_ref() else { return rep };
// Per-row key intervals. Morton is NOT row-major, so consecutive rows'
// intervals OVERLAP (row 0 of a 4-wide span is keys [0..5], row 1 is
// [2..7]). Scanning them unmerged emits the same feature twice — a
// duplicated draw, a duplicated pick, and a `visible` count that lies.
// Merge first; the merge also shrinks the candidate set.
let mut spans: Vec<(usize, usize)> = Vec::with_capacity((cy1 - cy0 + 1) as usize);
for cy in cy0..=cy1 {
let lo = (morton(cx0, cy) as u64) << block;
let hi = (((morton(cx1, cy) as u64) + 1) << block) - 1;
let a = tree.lower_bound(lo.min(u32::MAX as u64) as u32);
// `hi + 1` as a u64 so the top-of-range case does not wrap.
let hi_next = hi + 1;
let b = if hi_next > u32::MAX as u64 {
self.keys.len()
} else {
tree.lower_bound(hi_next as u32)
};
let b = b.min(self.keys.len());
if a < b {
spans.push((a, b));
}
}
spans.sort_unstable();
let mut merged: Vec<(usize, usize)> = Vec::with_capacity(spans.len());
for (a, b) in spans {
match merged.last_mut() {
Some(last) if a <= last.1 => last.1 = last.1.max(b),
_ => merged.push((a, b)),
}
}
rep.ranges = merged.len();
let mut candidates = 0usize;
for (a, b) in merged {
for i in a..b {
candidates += 1;
if self.boxes[i].overlaps_2d(&view) {
out.push(self.order[i]);
}
}
}
rep.candidates = candidates;
rep.visible = out.len();
rep
}
/// The reference answer — a full linear scan. Public so the S+tree path can
/// be **red-proven** against it on random data rather than trusted.
#[must_use]
pub fn query_linear(&self, view: WorldAabb) -> Vec<u32> {
let mut v: Vec<u32> =
self.boxes.iter().enumerate().filter(|(_, b)| b.overlaps_2d(&view)).map(|(i, _)| self.order[i]).collect();
v.sort_unstable();
v
}
}
/// The viewport, in world space, as an AABB — what the engine hands the culler.
#[must_use]
pub fn view_aabb(centre: WorldPos, half_w: f64, half_h: f64) -> WorldAabb {
WorldAabb {
min: WorldPos::flat(centre.x - half_w, centre.y - half_h),
max: WorldPos::flat(centre.x + half_w, centre.y + half_h),
}
}
#[cfg(test)]
mod tests {
use super::*;
fn pt(x: f64, y: f64) -> WorldAabb {
WorldAabb::point(WorldPos::flat(x, y))
}
/// Deterministic LCG — no dev-dependency for a random test.
fn lcg(seed: &mut u64) -> f64 {
*seed = seed.wrapping_mul(6_364_136_223_846_793_005).wrapping_add(1_442_695_040_888_963_407);
((*seed >> 11) as f64) / ((1u64 << 53) as f64)
}
#[test]
fn stree_query_equals_linear_scan_on_random_points() {
let mut seed = 0x5eed_1234_u64;
let items: Vec<(u32, WorldAabb)> =
(0..20_000u32).map(|i| (i, pt(lcg(&mut seed) * 1000.0, lcg(&mut seed) * 1000.0))).collect();
let idx = CullIndex::build(&items);
let mut out = Vec::new();
let mut total_cand = 0usize;
let mut total_vis = 0usize;
for _ in 0..200 {
let (x, y) = (lcg(&mut seed) * 1000.0, lcg(&mut seed) * 1000.0);
let (w, h) = (lcg(&mut seed) * 60.0 + 1.0, lcg(&mut seed) * 60.0 + 1.0);
let view = view_aabb(WorldPos::flat(x, y), w, h);
let rep = idx.query(view, &mut out);
let mut got = out.clone();
got.sort_unstable();
assert_eq!(got, idx.query_linear(view), "index disagreed with linear scan at {x},{y}");
assert_eq!(rep.visible, got.len());
total_cand += rep.candidates;
total_vis += rep.visible;
}
// The anti-lie assertion: the index must actually be an index. A linear
// scan would report 200 * 20_000 = 4_000_000 candidates.
assert!(total_cand < 400_000, "index tested {total_cand} candidates — that is a linear scan in disguise");
assert!(total_vis > 0, "the probe never hit anything — the test proves nothing");
}
#[test]
fn empty_and_degenerate_inputs_do_not_panic() {
let idx = CullIndex::build(&[]);
let mut out = Vec::new();
let rep = idx.query(view_aabb(WorldPos::flat(0.0, 0.0), 1.0, 1.0), &mut out);
assert_eq!(rep.visible, 0);
let same: Vec<(u32, WorldAabb)> = (0..100u32).map(|i| (i, pt(7.0, 7.0))).collect();
let idx = CullIndex::build(&same);
let rep = idx.query(view_aabb(WorldPos::flat(7.0, 7.0), 1.0, 1.0), &mut out);
assert_eq!(rep.visible, 100, "every coincident point must survive a viewport containing them");
}
#[test]
fn row_budget_bounds_the_range_count() {
let mut seed = 99u64;
let items: Vec<(u32, WorldAabb)> =
(0..5_000u32).map(|i| (i, pt(lcg(&mut seed) * 1000.0, lcg(&mut seed) * 1000.0))).collect();
let idx = CullIndex::build(&items);
let mut out = Vec::new();
// The whole world — the worst case for row decomposition.
let rep = idx.query(view_aabb(WorldPos::flat(500.0, 500.0), 600.0, 600.0), &mut out);
assert!(rep.ranges <= MAX_ROWS as usize + 1, "ranges={} exceeded the budget", rep.ranges);
assert_eq!(rep.visible, 5_000, "a world-sized viewport must return everything");
}
#[test]
fn oversize_features_are_flagged_not_silently_dropped() {
let mut items: Vec<(u32, WorldAabb)> = (0..100u32).map(|i| (i, pt(i as f64, i as f64))).collect();
items.push((
999,
WorldAabb { min: WorldPos::flat(0.0, 0.0), max: WorldPos::flat(100.0, 100.0) },
));
let idx = CullIndex::build(&items);
assert!(idx.oversize().contains(&999), "a dataset-spanning feature must be flagged oversize");
}
}