use crate::Result;
use crate::bulk::DEFAULT_FILL_FACTOR;
const BYTES_PER_CELL: usize = 24;
const NODE_HEADER: usize = 4;
const RND_TOWARDS: f64 = 1.0 - 1.0 / 8_388_608.0;
const RND_AWAY: f64 = 1.0 + 1.0 / 8_388_608.0;
fn coord_down(d: f64) -> f32 {
let f = d as f32;
if f64::from(f) > d {
(d * if d < 0.0 { RND_AWAY } else { RND_TOWARDS }) as f32
} else {
f
}
}
fn coord_up(d: f64) -> f32 {
let f = d as f32;
if f64::from(f) < d {
(d * if d < 0.0 { RND_TOWARDS } else { RND_AWAY }) as f32
} else {
f
}
}
#[derive(Debug, Clone, Copy)]
struct Cell {
id: i64,
bounds: [f32; 4],
}
impl Cell {
fn union(bounds: &[f32; 4], other: &[f32; 4]) -> [f32; 4] {
[
bounds[0].min(other[0]),
bounds[1].max(other[1]),
bounds[2].min(other[2]),
bounds[3].max(other[3]),
]
}
fn centre(&self) -> (f64, f64) {
(
(f64::from(self.bounds[0]) + f64::from(self.bounds[1])) / 2.0,
(f64::from(self.bounds[2]) + f64::from(self.bounds[3])) / 2.0,
)
}
}
pub(crate) trait NodeSink {
fn node(&mut self, nodeno: i64, blob: &[u8]) -> Result<()>;
fn rowid(&mut self, rowid: i64, nodeno: i64) -> Result<()>;
fn parent(&mut self, nodeno: i64, parentnode: i64) -> Result<()>;
}
fn hilbert(x: f64, y: f64, extent: &[f64; 4]) -> u32 {
const SIDE: u32 = 1 << 16;
let norm = |v: f64, lo: f64, hi: f64| -> u32 {
if !v.is_finite() || hi <= lo {
return 0;
}
let t = ((v - lo) / (hi - lo)).clamp(0.0, 1.0);
let scaled = t * f64::from(SIDE - 1);
#[expect(
clippy::cast_sign_loss,
reason = "t is clamped to [0, 1], so scaled lies in [0, SIDE - 1] and is never negative"
)]
let index = scaled as u32;
index.min(SIDE - 1)
};
let mut hx = norm(x, extent[0], extent[1]);
let mut hy = norm(y, extent[2], extent[3]);
let mut d: u32 = 0;
let mut s: u32 = SIDE / 2;
while s > 0 {
let rx = u32::from((hx & s) > 0);
let ry = u32::from((hy & s) > 0);
d = d.wrapping_add(s.wrapping_mul(s).wrapping_mul((3 * rx) ^ ry));
if ry == 0 {
if rx == 1 {
hx = s.wrapping_sub(1).wrapping_sub(hx);
hy = s.wrapping_sub(1).wrapping_sub(hy);
}
std::mem::swap(&mut hx, &mut hy);
}
s /= 2;
}
d
}
fn encode_node_into(blob: &mut Vec<u8>, depth: u16, cells: &[Cell], node_size: usize) {
blob.clear();
blob.reserve(node_size);
blob.extend_from_slice(&depth.to_be_bytes());
let count = u16::try_from(cells.len()).unwrap_or(u16::MAX);
blob.extend_from_slice(&count.to_be_bytes());
for cell in cells {
blob.extend_from_slice(&cell.id.to_be_bytes());
for coord in &cell.bounds {
blob.extend_from_slice(&coord.to_be_bytes());
}
}
blob.resize(node_size, 0);
}
fn level_sizes(entries: usize, per_node: usize) -> Vec<usize> {
let fanout = per_node.max(2);
let mut levels = Vec::new();
let mut count = entries.div_ceil(fanout).max(1);
levels.push(count);
while count > 1 {
count = count.div_ceil(fanout);
levels.push(count);
}
levels
}
pub(crate) fn pack_into<S: NodeSink>(
entries: &[(i64, [f64; 4])],
node_size: usize,
fill_factor: f64,
sink: &mut S,
) -> Result<()> {
let capacity = node_size.saturating_sub(NODE_HEADER) / BYTES_PER_CELL;
debug_assert!(capacity > 0, "node size {node_size} holds no cells");
let per_node = {
let requested = if fill_factor.is_nan() {
DEFAULT_FILL_FACTOR
} else {
fill_factor.clamp(f64::MIN_POSITIVE, 1.0)
};
let scaled = (capacity as f64) * requested;
#[expect(
clippy::cast_sign_loss,
reason = "`requested` is positive and at most 1, so `scaled` lies in [0, capacity] and rounds to a non-negative value"
)]
let rounded = scaled.round() as usize;
rounded.clamp(2.min(capacity), capacity.max(1))
};
let mut blob = Vec::with_capacity(node_size);
if entries.is_empty() {
encode_node_into(&mut blob, 0, &[], node_size);
return sink.node(1, &blob);
}
let mut keyed: Vec<(u32, Cell)> = Vec::with_capacity(entries.len());
let mut extent = [
f64::INFINITY,
f64::NEG_INFINITY,
f64::INFINITY,
f64::NEG_INFINITY,
];
for &(id, [min_x, max_x, min_y, max_y]) in entries {
let cell = Cell {
id,
bounds: [
coord_down(min_x),
coord_up(max_x),
coord_down(min_y),
coord_up(max_y),
],
};
let (cx, cy) = cell.centre();
extent = [
extent[0].min(cx),
extent[1].max(cx),
extent[2].min(cy),
extent[3].max(cy),
];
keyed.push((0, cell));
}
for (key, cell) in &mut keyed {
let (cx, cy) = cell.centre();
*key = hilbert(cx, cy, &extent);
}
keyed.sort_unstable_by_key(|(key, _)| *key);
let levels = level_sizes(keyed.len(), per_node);
let depth = u16::try_from(levels.len().saturating_sub(1)).unwrap_or(u16::MAX);
let mut level_base = Vec::with_capacity(levels.len());
let mut next = 2_i64;
for (index, count) in levels.iter().enumerate() {
if index + 1 == levels.len() {
level_base.push(1_i64);
} else {
level_base.push(next);
next += i64::try_from(*count).unwrap_or(i64::MAX);
}
}
let base = level_base.first().copied().unwrap_or(1);
let mut level: Vec<Cell> = Vec::with_capacity(levels.first().copied().unwrap_or(1));
let mut cells: Vec<Cell> = Vec::with_capacity(per_node);
for (index, chunk) in keyed.chunks(per_node).enumerate() {
let nodeno = base + i64::try_from(index).unwrap_or(0);
cells.clear();
cells.extend(chunk.iter().map(|(_, cell)| *cell));
let mut bounds = cells.first().map_or([0.0; 4], |first| first.bounds);
for cell in &cells {
bounds = Cell::union(&bounds, &cell.bounds);
sink.rowid(cell.id, nodeno)?;
}
encode_node_into(&mut blob, depth, &cells, node_size);
sink.node(nodeno, &blob)?;
level.push(Cell { id: nodeno, bounds });
}
for height in 1..levels.len() {
let base = level_base.get(height).copied().unwrap_or(1);
let mut parents: Vec<Cell> = Vec::with_capacity(levels.get(height).copied().unwrap_or(1));
for (index, chunk) in level.chunks(per_node).enumerate() {
let nodeno = base + i64::try_from(index).unwrap_or(0);
let mut bounds = chunk.first().map_or([0.0; 4], |first| first.bounds);
for cell in chunk {
bounds = Cell::union(&bounds, &cell.bounds);
sink.parent(cell.id, nodeno)?;
}
let node_depth = u16::try_from(height).unwrap_or(0);
encode_node_into(&mut blob, node_depth, chunk, node_size);
sink.node(nodeno, &blob)?;
parents.push(Cell { id: nodeno, bounds });
}
level = parents;
}
Ok(())
}
#[cfg(test)]
mod tests {
use super::*;
#[derive(Default)]
struct Collected {
nodes: Vec<(i64, Vec<u8>)>,
rowid_map: Vec<(i64, i64)>,
parent_map: Vec<(i64, i64)>,
}
impl NodeSink for Collected {
fn node(&mut self, nodeno: i64, blob: &[u8]) -> Result<()> {
self.nodes.push((nodeno, blob.to_vec()));
Ok(())
}
fn rowid(&mut self, rowid: i64, nodeno: i64) -> Result<()> {
self.rowid_map.push((rowid, nodeno));
Ok(())
}
fn parent(&mut self, nodeno: i64, parentnode: i64) -> Result<()> {
self.parent_map.push((nodeno, parentnode));
Ok(())
}
}
fn pack(entries: &[(i64, [f64; 4])], node_size: usize) -> Result<Collected> {
let mut collected = Collected::default();
pack_into(entries, node_size, 1.0, &mut collected)?;
Ok(collected)
}
fn decode_cells(blob: &[u8]) -> Vec<Cell> {
let count = blob
.get(2..4)
.and_then(|b| <[u8; 2]>::try_from(b).ok())
.map_or(0, |b| usize::from(u16::from_be_bytes(b)));
let mut cells = Vec::with_capacity(count);
for i in 0..count {
let start = NODE_HEADER + i * BYTES_PER_CELL;
let Some(raw) = blob.get(start..start + BYTES_PER_CELL) else {
break;
};
let id = <[u8; 8]>::try_from(raw.get(0..8).unwrap_or_default())
.map(i64::from_be_bytes)
.unwrap_or_default();
let mut bounds = [0.0f32; 4];
for (n, slot) in bounds.iter_mut().enumerate() {
let at = 8 + n * 4;
*slot = <[u8; 4]>::try_from(raw.get(at..at + 4).unwrap_or_default())
.map(f32::from_be_bytes)
.unwrap_or_default();
}
cells.push(Cell { id, bounds });
}
cells
}
#[test]
fn all_nodes_have_the_declared_size() {
let entries: Vec<(i64, [f64; 4])> = (0..5000)
.map(|i| {
let f = f64::from(i);
(i64::from(i), [f, f + 1.0, -f, -f + 1.0])
})
.collect();
let packed = pack(&entries, 1228).unwrap();
assert!(packed.nodes.len() > 1, "expected a multi-level tree");
for (_, blob) in &packed.nodes {
assert_eq!(blob.len(), 1228);
}
}
#[test]
fn level_sizes_terminates_for_a_degenerate_fanout() {
for per_node in [0, 1, 2] {
let levels = level_sizes(1000, per_node);
assert_eq!(
levels.last(),
Some(&1),
"level walk did not reach a root for per_node {per_node}"
);
assert!(
levels.len() <= 64,
"level walk produced {} levels for per_node {per_node}",
levels.len()
);
}
}
#[test]
fn fill_factor_controls_node_occupancy() {
let entries: Vec<(i64, [f64; 4])> = (0..1000)
.map(|i| (i64::from(i), [f64::from(i), f64::from(i), 0.0, 0.0]))
.collect();
let full = pack(&entries, 1228).unwrap();
let mut half = Collected::default();
pack_into(&entries, 1228, 0.5, &mut half).unwrap();
let leaf_counts = |c: &Collected| -> Vec<usize> {
let mut counts: Vec<usize> = c
.nodes
.iter()
.map(|(_, blob)| decode_cells(blob).len())
.collect();
counts.sort_unstable();
counts
};
assert_eq!(
leaf_counts(&full).last(),
Some(&51),
"full packs to capacity"
);
assert_eq!(leaf_counts(&half).last(), Some(&26), "half fills half");
assert!(
half.nodes.len() > full.nodes.len(),
"a lower fill factor should need more nodes"
);
for packed in [&full, &half] {
let mut ids: Vec<i64> = packed.rowid_map.iter().map(|&(id, _)| id).collect();
ids.sort_unstable();
ids.dedup();
assert_eq!(ids.len(), entries.len());
}
}
#[test]
fn out_of_range_fill_factors_are_clamped() {
let entries: Vec<(i64, [f64; 4])> = (0..200)
.map(|i| (i64::from(i), [f64::from(i), f64::from(i), 0.0, 0.0]))
.collect();
for factor in [0.0, -1.0, 2.0, f64::NAN] {
let mut collected = Collected::default();
pack_into(&entries, 1228, factor, &mut collected).unwrap();
for (_, blob) in &collected.nodes {
let cells = decode_cells(blob).len();
assert!(cells <= 51, "node over capacity with factor {factor}");
}
let mut ids: Vec<i64> = collected.rowid_map.iter().map(|&(id, _)| id).collect();
ids.sort_unstable();
ids.dedup();
assert_eq!(
ids.len(),
entries.len(),
"entries lost with factor {factor}"
);
}
}
#[test]
fn empty_entry_set_still_has_a_root() {
let packed = pack(&[], 1228).unwrap();
assert_eq!(packed.nodes.len(), 1);
let (nodeno, blob) = packed.nodes.first().unwrap();
assert_eq!(*nodeno, 1);
assert_eq!(blob.len(), 1228);
assert_eq!(u16::from_be_bytes([blob[0], blob[1]]), 0, "depth");
assert_eq!(u16::from_be_bytes([blob[2], blob[3]]), 0, "cell count");
assert!(packed.rowid_map.is_empty());
assert!(packed.parent_map.is_empty());
}
#[test]
fn small_set_is_a_single_root_leaf() {
let entries: Vec<(i64, [f64; 4])> = (1..=10)
.map(|i| (i, [i as f64, i as f64, 0.0, 0.0]))
.collect();
let packed = pack(&entries, 1228).unwrap();
assert_eq!(packed.nodes.len(), 1);
assert_eq!(packed.nodes.first().unwrap().0, 1);
assert!(packed.parent_map.is_empty());
assert_eq!(packed.rowid_map.len(), 10);
assert!(packed.rowid_map.iter().all(|&(_, nodeno)| nodeno == 1));
}
#[test]
fn root_is_node_one() {
let entries: Vec<(i64, [f64; 4])> = (0..20_000)
.map(|i| (i64::from(i), [f64::from(i), f64::from(i), 0.0, 0.0]))
.collect();
let packed = pack(&entries, 1228).unwrap();
assert!(packed.nodes.iter().any(|(nodeno, _)| *nodeno == 1));
assert!(
packed.parent_map.iter().all(|&(child, _)| child != 1),
"the root must not appear as a child"
);
}
#[test]
fn bounds_are_rounded_outward() {
let lo = 0.1_f64;
let hi = 0.300_000_000_000_000_04_f64;
assert!(f64::from(coord_down(lo)) <= lo);
assert!(f64::from(coord_up(hi)) >= hi);
let entries = vec![(1_i64, [lo, hi, -hi, -lo])];
let packed = pack(&entries, 1228).unwrap();
let cells = decode_cells(&packed.nodes.first().unwrap().1);
let cell = cells.first().unwrap();
assert!(f64::from(cell.bounds[0]) <= lo);
assert!(f64::from(cell.bounds[1]) >= hi);
assert!(f64::from(cell.bounds[2]) <= -hi);
assert!(f64::from(cell.bounds[3]) >= -lo);
}
#[test]
fn parent_bounds_contain_child_bounds() {
let entries: Vec<(i64, [f64; 4])> = (0..30_000)
.map(|i| {
let x = f64::from(i % 173) * 1.7;
let y = f64::from(i % 91) * -2.3;
(i64::from(i), [x, x + 0.5, y, y + 0.5])
})
.collect();
let packed = pack(&entries, 1228).unwrap();
let bounds_of: std::collections::HashMap<i64, [f32; 4]> = packed
.nodes
.iter()
.map(|(nodeno, blob)| {
let cells = decode_cells(blob);
let mut bounds = cells.first().map_or([0.0; 4], |c| c.bounds);
for cell in &cells {
bounds = Cell::union(&bounds, &cell.bounds);
}
(*nodeno, bounds)
})
.collect();
for (nodeno, blob) in &packed.nodes {
if packed.parent_map.iter().any(|&(child, _)| child == *nodeno) || *nodeno == 1 {
for cell in decode_cells(blob) {
if let Some(child) = bounds_of.get(&cell.id)
&& packed
.parent_map
.iter()
.any(|&(c, p)| c == cell.id && p == *nodeno)
{
assert!(cell.bounds[0] <= child[0], "min_x not contained");
assert!(cell.bounds[1] >= child[1], "max_x not contained");
assert!(cell.bounds[2] <= child[2], "min_y not contained");
assert!(cell.bounds[3] >= child[3], "max_y not contained");
}
}
}
}
}
#[test]
fn mappings_cover_every_entry_and_node() {
let entries: Vec<(i64, [f64; 4])> = (0..12_345)
.map(|i| (i64::from(i), [f64::from(i), f64::from(i), 0.0, 0.0]))
.collect();
let packed = pack(&entries, 1228).unwrap();
let mut ids: Vec<i64> = packed.rowid_map.iter().map(|&(id, _)| id).collect();
ids.sort_unstable();
ids.dedup();
assert_eq!(ids.len(), entries.len(), "every rowid mapped exactly once");
let mut children: Vec<i64> = packed.parent_map.iter().map(|&(c, _)| c).collect();
children.sort_unstable();
children.dedup();
assert_eq!(
children.len(),
packed.nodes.len() - 1,
"every non-root node has exactly one parent entry"
);
}
}