R-tree

Deep dive

How it works

In: boxes with ids (insert), a query window or a point and k · Out: ids whose box meets the window, or the k nearest · State: the tree of nested boxes · Cost: search O(log_M N + visited); insert O(M) per level

Entry = (box, pointer): leaf pointer = object id, inner pointer = child node. All leaves at one depth; every non-root node holds m..M entries (m about 40 % of M).

1. Search(q): from the root, visit entries whose box intersects q; report matches at leaves.
2. Insert(r): descend to the child needing least enlargement to cover r (ties: smaller area); add to the leaf; refresh parent boxes on the way up.
3. Split an overflowing node (Guttman's quadratic split): seeds = the pair wasting most area together; then repeatedly give the entry with the strongest preference to the group it enlarges less, forcing the rest into a group that would fall below m. The parent gains an entry and may split; a root split adds a level.
4. Nearest(p, k): best-first (Hjaltason and Samet). A priority queue of nodes and items keyed by minimum distance from p to their box; an item popped is the next neighbour, a node is expanded. No node farther than the k-th answer is touched.
5. Bulk load (STR): sort by x centre, cut into ceil(sqrt(N/M)) slices, sort each by y, pack runs of M into leaves; repeat on the leaf boxes until one root remains. Nodes are about 100 % full.
6. R*-tree refines 2-3: least overlap enlargement at the leaf level, split by smallest margin sum then least overlap, and forced reinsertion of the farthest entries on first overflow per level.

         root: [A][B]
        /            \
 A: [a1][a2][a3]   B: [b1][b2]     A and B may overlap:
     |    |    |       |    |      a window touching both
   leaf leaf leaf    leaf leaf     descends into both

Code

2-D tree in an arena: quadratic-split insertion, window search, best-first k-nearest. M = 16; tested against brute force on 20,000 random points.

use std::cmp::Reverse;
use std::collections::BinaryHeap;

#[derive(Clone, Copy, Debug)]
pub struct Rect { pub lo: [f64; 2], pub hi: [f64; 2] }

impl Rect {
    pub fn point(x: f64, y: f64) -> Rect { Rect { lo: [x, y], hi: [x, y] } }
    fn area(&self) -> f64 { (self.hi[0] - self.lo[0]) * (self.hi[1] - self.lo[1]) }
    fn union(&self, o: &Rect) -> Rect {
        Rect { lo: [self.lo[0].min(o.lo[0]), self.lo[1].min(o.lo[1])],
               hi: [self.hi[0].max(o.hi[0]), self.hi[1].max(o.hi[1])] }
    }
    fn enlargement(&self, o: &Rect) -> f64 { self.union(o).area() - self.area() }
    fn intersects(&self, o: &Rect) -> bool {
        self.lo[0] <= o.hi[0] && o.lo[0] <= self.hi[0] && self.lo[1] <= o.hi[1] && o.lo[1] <= self.hi[1]
    }
    /// Squared distance from a point to the rectangle (0 inside).
    fn mindist2(&self, p: [f64; 2]) -> f64 {
        (0..2).map(|d| { let e = (self.lo[d] - p[d]).max(p[d] - self.hi[d]).max(0.0); e * e }).sum()
    }
}

type Entry = (Rect, usize); // leaf: (box, item id); internal: (child MBR, child node index)
struct Node { leaf: bool, kids: Vec<Entry> }

pub struct RTree { nodes: Vec<Node>, root: usize }
const M: usize = 16;       // max entries per node
const MIN: usize = M * 2 / 5; // min fill 40 %

fn mbr(es: &[Entry]) -> Rect { es.iter().skip(1).fold(es[0].0, |a, e| a.union(&e.0)) }

/// Guttman's quadratic split: returns the second group; `es` keeps the first.
fn split(es: &mut Vec<Entry>) -> Vec<Entry> {
    let n = es.len();
    let (mut si, mut sj, mut worst) = (0, 1, f64::MIN);
    for i in 0..n { for j in i + 1..n {
        let w = es[i].0.union(&es[j].0).area() - es[i].0.area() - es[j].0.area();
        if w > worst { worst = w; si = i; sj = j; }
    }}
    let mut rest = std::mem::take(es);
    let b0 = rest.swap_remove(sj); // sj > si, so si stays valid
    let a0 = rest.swap_remove(si);
    let (mut ga, mut gb) = (vec![a0], vec![b0]);
    let (mut ra, mut rb) = (a0.0, b0.0);
    while !rest.is_empty() {
        let k = if ga.len() + rest.len() == MIN { Some(true) } else if gb.len() + rest.len() == MIN { Some(false) } else { None };
        let pick = (0..rest.len()).max_by(|&x, &y| {
            let d = |i: usize| (ra.enlargement(&rest[i].0) - rb.enlargement(&rest[i].0)).abs();
            d(x).total_cmp(&d(y)) }).unwrap();
        let e = rest.swap_remove(pick);
        let (ea, eb) = (ra.enlargement(&e.0), rb.enlargement(&e.0));
        let to_a = k.unwrap_or(ea < eb || (ea == eb && (ra.area() < rb.area() || (ra.area() == rb.area() && ga.len() <= gb.len()))));
        if to_a { ra = ra.union(&e.0); ga.push(e); } else { rb = rb.union(&e.0); gb.push(e); }
    }
    *es = ga;
    gb
}

impl RTree {
    pub fn new() -> RTree { RTree { nodes: vec![Node { leaf: true, kids: vec![] }], root: 0 } }

    pub fn insert(&mut self, r: Rect, id: usize) {
        if let Some(sib) = self.insert_rec(self.root, r, id) {
            let old = self.root;
            let kids = vec![(mbr(&self.nodes[old].kids), old), (mbr(&self.nodes[sib].kids), sib)];
            self.nodes.push(Node { leaf: false, kids });
            self.root = self.nodes.len() - 1;
        }
    }
    fn insert_rec(&mut self, n: usize, r: Rect, id: usize) -> Option<usize> {
        if self.nodes[n].leaf { self.nodes[n].kids.push((r, id)); } else {
            // least enlargement, ties by smaller area
            let i = (0..self.nodes[n].kids.len()).min_by(|&a, &b| {
                let (ka, kb) = (&self.nodes[n].kids[a].0, &self.nodes[n].kids[b].0);
                (ka.enlargement(&r), ka.area()).partial_cmp(&(kb.enlargement(&r), kb.area())).unwrap() }).unwrap();
            let c = self.nodes[n].kids[i].1;
            let sib = self.insert_rec(c, r, id);
            self.nodes[n].kids[i].0 = mbr(&self.nodes[c].kids);
            if let Some(s) = sib { let e = (mbr(&self.nodes[s].kids), s); self.nodes[n].kids.push(e); }
        }
        if self.nodes[n].kids.len() <= M { return None; }
        let mut kids = std::mem::take(&mut self.nodes[n].kids);
        let other = split(&mut kids);
        let leaf = self.nodes[n].leaf;
        self.nodes[n].kids = kids;
        self.nodes.push(Node { leaf, kids: other });
        Some(self.nodes.len() - 1)
    }

    pub fn search(&self, q: &Rect, out: &mut Vec<usize>) {
        let mut stack = vec![self.root];
        while let Some(n) = stack.pop() {
            for (r, c) in &self.nodes[n].kids {
                if r.intersects(q) { if self.nodes[n].leaf { out.push(*c) } else { stack.push(*c) } }
            }
        }
    }

    /// Best-first k-nearest (Hjaltason-Samet). Non-negative f64 bits order like the values.
    pub fn nearest(&self, p: [f64; 2], k: usize) -> Vec<usize> {
        let mut heap = BinaryHeap::new(); // (dist bits, is_item, node-or-id)
        heap.push(Reverse((0f64.to_bits(), false, self.root)));
        let mut out = vec![];
        while let Some(Reverse((_, is_item, x))) = heap.pop() {
            if is_item { out.push(x); if out.len() == k { break; } continue; }
            for (r, c) in &self.nodes[x].kids {
                heap.push(Reverse((r.mindist2(p).to_bits(), self.nodes[x].leaf, *c)));
            }
        }
        out
    }
}

STR (step 5) replaces insert for static data: sort by centre-x, chunks, sort by centre-y, chunks(M), push a Node per chunk, repeat on the chunk boxes. I ran this variant through the same tests: for 20,000 random points it built 1,335 nodes versus 2,021 for one-by-one insertion above (measured, M = 16).

Complexity and constants

• Search: O(log_M N + visited); visited grows with overlap, and a window returning T items ideally touches order of T/M leaves.
• Insert: O(M) choose-subtree per level, plus an O(M^2) quadratic split amortised over about M inserts. Linear split is O(M) but gives worse trees.
• Fill: packed about 100 %; dynamic lower (above: 1.5x the nodes of STR). Cost driver: one cache miss or page read per visited node.
• Bytes: f64 2-D box 32 B + 8 B pointer = 40 B per entry; M = 16 is 640 B = 10 cache lines. f32 boxes halve the box part.

Efficiency

• Cost: search O(log_M N + visited); one cache miss or page read per visited node, so node visits, not comparisons, are the bound. An f64 2-D entry is 40 B (32 B box + 8 B pointer), so M = 16 is 640 B, 10 cache lines. A dynamic tree needs about 1.5x the nodes of an STR-packed one (measured above).
• Levers:
  • bulk load (STR or Hilbert sort) — full nodes, little overlap; rbush's README states bulk insertion is usually 2-3x faster than one by one — static data (see bulk loading, space-filling curves)
  • flat buffer with implicit child offsets (flatbush, geo-index) — no pointers or allocation, mmap-able; geo-index documents it as smaller than the dynamic rstar — static data
  • compact boxes: f32 (SQLite stores f32 pairs), rounded outward so no hit is lost — more entries per page, fewer levels — large indexes
  • best-first k-NN with LIMIT k — visits a minimal node set, never range-scan and sort — nearest queries
• Naive → standard → state of the art:
  • linear scan of all boxes · O(N)
  • dynamic R*-tree, insertion with forced reinsertion · O(log_M N + visited)
  • STR or Hilbert packed flat tree · full nodes, lowest overlap; for points only, a cell index or a Hilbert-sorted array costs less

Implementation notes

• M = disk page (4-16 KiB) on disk; 8-32 entries in memory. rbush defaults to 9 entries per node, rstar to 6 max / 3 min.
• Scan one node's boxes linearly (SIMD-friendly with structure-of-arrays boxes); keep child boxes in the parent so a failed test costs no child read.
• Squared distances as to_bits keys: non-negative IEEE floats order like their bits, so the heap needs no NaN handling.
• Delete: shrink boxes upward and reinsert entries of underfull nodes (Guttman), or leave them until empty like rbush.

Pitfalls

• Huge or long thin boxes, many duplicates: overlap explodes; keep big objects in a separate list.
• Sorted insertion order with plain quadratic split gives thin nodes; bulk load or use R*.
• Rounding f32 boxes to nearest instead of outward drops hits; geographic lon/lat boxes break at the antimeridian.
• The box test only filters; run the exact predicate on candidates.

Variants and alternatives

In production

• Databases: the SQLite R*Tree module indexes 1 to 5 dimensions and its docs show start and end times of events as an example (docs); PostGIS uses an R-tree implemented on top of GiST (docs).
• GIS and computational geometry: JTS STRtree; Boost.Geometry offers index::rtree, its only spatial index (docs); libspatialindex behind Python rtree.
• Web maps: rbush and flatbush in browser mapping code; rstar in the Rust georust stack.

Implementations

Sources

• Antonin Guttman, R-trees: A Dynamic Index Structure for Spatial Searching, SIGMOD 1984. paper
• Norbert Beckmann, Hans-Peter Kriegel, Ralf Schneider, Bernhard Seeger, The R*-tree: An Efficient and Robust Access Method for Points and Rectangles, SIGMOD 1990. paper
• Scott T. Leutenegger, Mario A. Lopez, Jeffrey Edgington, STR: A Simple and Efficient Algorithm for R-tree Packing, ICDE 1997. paper
• Gisli R. Hjaltason, Hanan Samet, Distance Browsing in Spatial Databases, ACM TODS 1999. paper
• Lars Arge, Mark de Berg, Herman Haverkort, Ke Yi, The Priority R-tree: A Practically Efficient and Worst-Case Optimal R-tree, SIGMOD 2004. paper
• SQLite documentation, The SQLite RTree Module*. docs