Spatial hashing

Deep dive

How it works

In: n object positions and a cell size h ≥ the largest object diameter · Out: candidate neighbour pairs, or a neighbour list per object · State: none; the grid (bucket starts plus sorted ids) is rebuilt each step · Cost: build O(n), query O(n + k)

1. Choose cell size h ≥ 2 · r_max (largest bounding-sphere diameter, or largest AABB extent).
2. For each object: cell = floor(center / h) per axis; bucket = hash(cell) & (table_size - 1).
3. Counting sort by bucket: count per bucket, exclusive prefix sum gives cell_start[], scatter ids into sorted[]. Objects of one bucket are now contiguous.
4. For each object, visit the 27 buckets around its cell. For each id j > i in those buckets, run the exact test (sphere or AABB overlap). The j > i check reports each pair once.

cell (x,y,z) ──hash──► bucket b ──► sorted[cell_start[b] .. cell_start[b+1]]
                                       │ ids of all objects whose cell hashes to b
query: 3×3×3 neighbour cells ──► 27 buckets ──► exact test

Variants of step 4: scan only 13 "forward" neighbour cells plus the own cell to halve the work (same pair coverage, no j > i filter needed across cells), or, if h ≥ 4 r_max, choose the 8 cells of the octant the center lies in.

Code

/// Broadphase over spheres with radius <= r_max. Cell size = 2 * r_max, so each sphere
/// is stored once (by its center) and partners can only be in the 27 surrounding cells.
pub struct SpatialHash {
    inv_cell: f32,
    mask: usize,             // table size is a power of two
    cell_start: Vec<u32>,    // bucket b holds sorted[cell_start[b]..cell_start[b+1]]
    sorted: Vec<u32>,        // object ids grouped by bucket
    keys: Vec<u32>,          // bucket per object (scratch)
}

#[inline]
fn hash_cell(x: i32, y: i32, z: i32) -> u32 {
    // Teschner et al. 2003 primes; wrapping keeps it defined for negative coordinates.
    (x.wrapping_mul(73856093) ^ y.wrapping_mul(19349663) ^ z.wrapping_mul(83492791)) as u32
}

impl SpatialHash {
    pub fn new(r_max: f32, n_objects: usize) -> Self {
        let size = (2 * n_objects).next_power_of_two().max(16); // load factor <= 0.5
        SpatialHash { inv_cell: 1.0 / (2.0 * r_max), mask: size - 1,
            cell_start: vec![0; size + 1], sorted: vec![0; n_objects], keys: vec![0; n_objects] }
    }

    #[inline]
    fn cell(&self, p: [f32; 3]) -> [i32; 3] {
        [(p[0] * self.inv_cell).floor() as i32, (p[1] * self.inv_cell).floor() as i32,
         (p[2] * self.inv_cell).floor() as i32]
    }

    /// Counting sort of objects by bucket: 2 passes, no allocation after `new`.
    pub fn build(&mut self, pos: &[[f32; 3]]) {
        self.cell_start.iter_mut().for_each(|c| *c = 0);
        for (i, &p) in pos.iter().enumerate() {
            let c = self.cell(p);
            let b = hash_cell(c[0], c[1], c[2]) as usize & self.mask;
            self.keys[i] = b as u32;
            self.cell_start[b + 1] += 1;
        }
        for b in 0..=self.mask { self.cell_start[b + 1] += self.cell_start[b]; }
        let mut fill = self.cell_start.clone(); // write cursor per bucket
        for i in 0..pos.len() {
            let b = self.keys[i] as usize;
            self.sorted[fill[b] as usize] = i as u32;
            fill[b] += 1;
        }
    }

    /// Calls `f(i, j)` once per overlapping pair (i < j).
    pub fn pairs(&self, pos: &[[f32; 3]], rad: &[f32], mut f: impl FnMut(u32, u32)) {
        let mut seen = [usize::MAX; 27]; // buckets already scanned for this object
        for i in 0..pos.len() {
            let c = self.cell(pos[i]);
            let mut n_seen = 0;
            for dz in -1..=1 { for dy in -1..=1 { for dx in -1..=1 {
                let b = hash_cell(c[0] + dx, c[1] + dy, c[2] + dz) as usize & self.mask;
                // two neighbour cells may hash to the same bucket: scan it once
                if seen[..n_seen].contains(&b) { continue; }
                seen[n_seen] = b; n_seen += 1;
                let (s, e) = (self.cell_start[b] as usize, self.cell_start[b + 1] as usize);
                for &j in &self.sorted[s..e] {
                    if j as usize <= i { continue; }           // each pair once
                    let pj = pos[j as usize];
                    let d = [pj[0] - pos[i][0], pj[1] - pos[i][1], pj[2] - pos[i][2]];
                    let r = rad[i] + rad[j as usize];
                    // exact test also rejects hash collisions from far-away cells
                    if d[0] * d[0] + d[1] * d[1] + d[2] * d[2] <= r * r { f(i as u32, j); }
                }
            }}}
        }
    }
}

Complexity and constants

• Build: O(n + table_size). Two streaming passes plus one scattered write per object.
• Query: O(n · 27 · m) candidate tests, m = objects per occupied cell. With h = 2 r_max and uniform density, m is constant and the cost is O(n + k).
• Memory: 4 bytes (sorted) + 4 bytes (keys) per object + 4 bytes per bucket; with a 2× table, order of 16-24 bytes per object.
• Cache behaviour: the 27 bucket lookups are random accesses. Reordering the object arrays themselves by bucket (or by Morton code of the cell) each frame makes neighbour data contiguous; GPU particle solvers do exactly this.

Efficiency

• Cost: build is two streaming passes plus one scattered write per object; a query makes 27 random bucket reads per object (9 in 2D). Memory is order of 16-24 bytes per object with a 2x table (this page). The bound is memory latency in the scatter and the 27-cell scan.
• Levers:
  • counting sort instead of linked lists — no pointer chase, parallel histograms (this page) — all builds
  • reorder the objects by cell or Morton code each step — neighbour data becomes contiguous (Green 2010) — GPU particle solvers
  • scan 13 forward cells plus own cell — half the lookups — symmetric pair reporting
  • dense grid instead of hash — no collisions, adjacent cells adjacent in memory — bounded worlds
  • keep the neighbour list per particle — one build serves all solver iterations (Ihmsen 2011) — SPH
• Naive → standard → state of the art:
  • compare all pairs · O(n²)
  • hashed grid with counting sort and 27-cell scan · O(n + k), order of 16-24 B per object
  • cell-sorted particles on GPU with radix sort and reused neighbour lists · linear, coalesced reads (Green 2010; Ihmsen 2011)

Implementation notes

• Counting sort, not linked lists. Per-cell linked lists cost a pointer chase per object; the prefix-sum layout is two linear passes and parallelises (per-thread histograms, then scan). On GPU: compute (bucket, id), radix sort by bucket, find bucket starts by comparing neighbours.
• Table size: power of two, order of 1-2× the object count. Smaller tables raise collisions; collisions are harmless for correctness, they only add candidates.
• Dense grid instead of a hash when the world is bounded and small: index = x + nx * (y + ny * z) has no collisions and neighbouring cells are adjacent in memory.
• Morton (Z-order) cell keys instead of a hash keep spatially close cells close in the sorted array, improving locality of the 27-cell scan.
• Objects larger than a cell: insert into every overlapped cell and deduplicate pairs (e.g. report only from the cell containing the pair's min-corner), or keep big objects in a separate list tested against everything.
• Fixed-radius neighbour search (SPH): same structure with h = kernel support radius; keep the neighbour list per particle for the solver iterations.

Pitfalls

• Truncating with as i32 instead of floor: cells -1 and 0 merge around the origin and the 27-cell guarantee breaks for negative coordinates.
• Scanning a bucket twice when two of the 27 cells hash to the same bucket: duplicate pairs. Deduplicate buckets per query (as in the code) or deduplicate pairs.
• Cell size chosen from average instead of maximum size: large objects miss partners two cells away.
• One giant object in the grid: it lands in many cells or forces a huge h, and the whole broadphase turns quadratic.
• Integer overflow in the hash with plain * in debug builds panics in Rust; use wrapping_mul.

Variants and alternatives

In production

Molecular dynamics: LAMMPS builds neighbour lists by binning, which scales linearly with atoms per processor and is almost always faster than the all-pairs style (docs). Fluids and graphics: GPU particle and SPH solvers sort particles by cell hash every step (NVIDIA's CUDA "particles" sample, Green 2010); the SPlisHSPlasH fluid framework uses the CompactNSearch library; NVIDIA Warp offers a HashGrid for GPU point queries; Teschner et al. used it for deformable-object collision. Games and robot simulators: Chipmunk2D can switch its broadphase to a spatial hash (cpSpaceUseSpatialHash); ODE has hash spaces (dHashSpaceCreate, user guide).

Implementations

Sources

• Matthias Teschner, Bruno Heidelberger, Matthias Müller, Danat Pomeranets, Markus Gross, Optimized Spatial Hashing for Collision Detection of Deformable Objects, VMV 2003.
• Christer Ericson, Real-Time Collision Detection, Morgan Kaufmann, 2005 (ch. 7, grids and hierarchical grids).
• Simon Green, Particle Simulation using CUDA, NVIDIA whitepaper, 2010.
• Markus Ihmsen, Nadir Akinci, Markus Becker, Matthias Teschner, A Parallel SPH Implementation on Multi-Core CPUs, Computer Graphics Forum, 2011. paper
• M. P. Allen, D. J. Tildesley, Computer Simulation of Liquids, 2nd ed., Oxford University Press, 2017. doi
• LAMMPS manual, neighbor command. docs