Spatial Hashing for O(1) Neighbour Lookup

Naïve neighbour search checks every pair of particles: O(N²). A spatial hash grid maps 3D positions to a flat hash table using only integer arithmetic, giving O(1) average lookups per query and enabling smooth 100k-particle simulations in real time.

1The core idea: cells and hash function

The world is conceptually divided into uniform cells of size h (the interaction radius). Each 3D coordinate maps to an integer cell index, which is then hashed to a table slot:

          const h = 0.5; // cell size = interaction radius // 3D position →
          integer cell coordinates function cellOf(x, y, z) { return
          [Math.floor(x / h), Math.floor(y / h), Math.floor(z / h)]; } //
          Integer cell → flat table index (large primes minimise collision
          clustering) const TABLE = 1 << 20; // 1M slots — must be power of 2
          for cheap modulo const P1 = 73856093, P2 = 19349663, P3 = 83492791;
          function hashCell(cx, cy, cz) { // XOR of scaled coordinates — fast
          and reasonably uniform return ((cx * P1) ^ (cy * P2) ^ (cz * P3)) &
          (TABLE - 1); }
        

Cell size should equal the interaction radius h. Querying the 3×3×3 (27 cells in 3D) neighborhood then guarantees every potential neighbor is found. In 2D it's 3×3 = 9 cells.

2Hash table layout — flat arrays

Instead of JavaScript objects (slow allocation), use two integer-typed arrays — a count array (how many particles in each slot) and a particle list array (flat list of particle indices per slot).

          const N = 100_000; // max particles // Flat sorted-list approach
          (fastest rebuild each frame): const cellStart = new Int32Array(TABLE +
          1); // prefix-sum start index const cellCount = new Int32Array(TABLE);
          // particles per cell const sortedIds = new Int32Array(N); // particle
          indices sorted by cell // Alternatively: chained buckets (simpler to
          update, slightly slower iterate) const buckets = new
          Array(TABLE).fill(null).map(() => []);
        

3Build phase: assign particles to cells

          // particles: Float32Array of [x,y,z, ...] with stride 3 function
          build(particles, n) { // 1. Count particles per cell
          cellCount.fill(0); for (let i = 0; i < n; i++) { const [cx, cy, cz]
          = cellOf(particles[i*3], particles[i*3+1], particles[i*3+2]);
          cellCount[hashCell(cx, cy, cz)]++; } // 2. Prefix sum → start index
          per cell cellStart[0] = 0; for (let b = 0; b < TABLE; b++) {
          cellStart[b + 1] = cellStart[b] + cellCount[b]; } // 3. Fill sortedIds
          (stable assignment) const tempCount = new Int32Array(TABLE); //
          working counter for (let i = 0; i < n; i++) { const [cx, cy, cz] =
          cellOf(particles[i*3], particles[i*3+1], particles[i*3+2]); const slot
          = hashCell(cx, cy, cz); sortedIds[cellStart[slot] + tempCount[slot]] =
          i; tempCount[slot]++; } } // Rebuild every frame for moving particles
          — O(N) time
        

4Query phase: iterate over neighbour cells

          // Iterate all particles within radius h of point (qx, qy, qz)
          function query(particles, qx, qy, qz, callback) { const [bcx, bcy,
          bcz] = cellOf(qx, qy, qz); const h2 = h * h; // Check the 3×3×3
          neighbourhood (27 cells) for (let dz = -1; dz <= 1; dz++) { for
          (let dy = -1; dy <= 1; dy++) { for (let dx = -1; dx <= 1; dx++)
          { const slot = hashCell(bcx + dx, bcy + dy, bcz + dz); const start =
          cellStart[slot]; const end = start + cellCount[slot]; for (let k =
          start; k < end; k++) { const j = sortedIds[k]; const rx =
          particles[j*3] - qx; const ry = particles[j*3+1] - qy; const rz =
          particles[j*3+2] - qz; const r2 = rx*rx + ry*ry + rz*rz; if (r2 <
          h2) callback(j, Math.sqrt(r2), rx, ry, rz); } } } } } // Usage: for
          each particle, find its neighbours for (let i = 0; i < N; i++) {
          query(particles, particles[i*3], particles[i*3+1], particles[i*3+2],
          (j, dist) => { // j is a neighbour index, dist is the distance }); }
        

Hash collisions (two cells mapping to the same slot) cause false-positive candidates — particles that hash to the same slot but are in different cells. The exact distance check r2 < h2 filters these out, so correctness is preserved at the cost of a few extra candidates.

5SPH fluid density with spatial hashing

          // SPH kernel: poly6 function poly6(r, h) { if (r >= h) return 0;
          const diff = h*h - r*r; return (315 / (64 * Math.PI * h**9)) *
          diff**3; } const mass = 0.02; const densities = new Float32Array(N);
          function computeDensities(particles) { build(particles, N); for (let i
          = 0; i < N; i++) { densities[i] = 0; query(particles,
          particles[i*3], particles[i*3+1], particles[i*3+2], (j, dist) => {
          densities[i] += mass * poly6(dist, h); } ); } } // Without spatial
          hashing: O(N²) = 10^10 operations for 100k particles // With spatial
          hashing: O(N·k) where k ≈ 27 × avg_per_cell ≈ 50..100
        

6Extensions and limits

Non-uniform densities — for wildly varying particle densities, consider a BVH or k-d tree instead (spatial hashing is best for roughly-uniform distributions).
Multiple radii — build separate hash grids for small/large radii, or always use the largest with distance filtering.
GPU spatial hashing — on the GPU, use a counting sort (prefix sum) in compute shaders/transform feedback. GPUSorting libraries exist for WebGPU.
Compact hashing — for sparse distributions, use a JavaScript Map<number, number[]> — lower memory for sparse grids, but slower iteration than typed arrays.

          // Tuning advice: // TABLE should be 2-4× N for good load factor (low
          collision rate) // Cell size h should match the interaction radius
          exactly // For 2D (particles confined to a plane): only query 3×3 = 9
          cells const TABLE_2D = 1 << 18; // 256k slots for 2D sim with
          50k particles function hashCell2D(cx, cy) { return ((cx * P1) ^ (cy *
          P2)) & (TABLE_2D - 1); } // Query 9 cells in 2D for (let dy = -1;
          dy <= 1; dy++) { for (let dx = -1; dx <= 1; dx++) { const slot =
          hashCell2D(bcx + dx, bcy + dy); // ... iterate slot } }
        