GLSL Spatial Hashing for 3D Particle Grids

Spatial hashing in GLSL is an optimization technique that maps continuous 3D particle positions into discrete cells within a uniform grid, collapsing three-dimensional space into a one-dimensional index using a hash function. This approach allows GPU shaders to rapidly partition space, reducing pairwise particle interaction checks from an \(O(N^2)\) brute-force complexity to an efficient \(O(N)\) local neighborhood search. This article explains the core mathematical formulation, the standard GLSL implementation, and the process for querying neighbor cells.

Coordinate Discretization

The first step converts a floating-point world-space position into integer grid cell coordinates. Dividing the continuous position vector by the uniform cell spacing size and taking the floor of each component yields the discrete 3D index:

ivec3 getCellCoords(vec3 position, float cellSize) {
    return ivec3(floor(position / cellSize));
}

The cellSize parameter is typically set to match or slightly exceed the interaction radius (or smoothing length) of the particles to ensure that interacting particles reside either in the same cell or in immediately adjacent cells.

The 3D Hash Function

Once the discrete integer coordinates (x, y, z) are established, they must be converted into a single integer hash value. When working with unbounded domains or fixed-size hash tables (such as Shader Storage Buffer Objects), large prime numbers are used to distribute the 3D coordinates evenly across the table while minimizing collisions:

uint hashCell(ivec3 cell, uint tableSize) {
    // Large prime numbers for bitwise mixing
    const uint p1 = 73856093u;
    const uint p2 = 19349663u;
    const uint p3 = 83492791u;

    // Convert negative coordinates safely to unsigned integers
    uvec3 uCell = uvec3(cell);

    // Compute hash via XOR and prime multiplication
    uint hash = (uCell.x * p1) ^ (uCell.y * p2) ^ (uCell.z * p3);

    return hash % tableSize;
}

If the domain is strictly bounded within a fixed resolution (for example, a grid of size \(W \times H \times D\)), a linear flat-index can be calculated directly without hashing collisions:

uint getLinearIndex(ivec3 cell, ivec3 gridDim) {
    return uint(cell.x + cell.y * gridDim.x + cell.z * (gridDim.x * gridDim.y));
}

Neighborhood Querying

To detect collisions or compute inter-particle forces (such as in Smoothed Particle Hydrodynamics), a thread must inspect the current cell and all 26 adjacent neighboring cells. This is achieved using nested loops covering the offset range from -1 to +1 along each axis:

void findNeighbors(vec3 position, float cellSize, uint tableSize) {
    ivec3 baseCell = getCellCoords(position, cellSize);

    for (int z = -1; z <= 1; ++z) {
        for (int y = -1; y <= 1; ++y) {
            for (int x = -1; x <= 1; ++x) {
                ivec3 neighborCell = baseCell + ivec3(x, y, z);
                uint cellHash = hashCell(neighborCell, tableSize);

                // Retrieve particle list or head pointer associated with cellHash
                // and evaluate distance checks
            }
        }
    }
}

Memory Layout and Collision Management

Because multiple cells can map to the same hash value when using a modulo hash function, the GPU implementation must handle collisions. Standard GPU spatial binning architectures resolve this using two passes:

1. Sorting/Counting Pass: Each particle computes its cell hash. The particle indices are then sorted by their hash value using a GPU-friendly parallel sort (such as Bitonic or Radix sort).
2. Cell Offset Buffer: A secondary buffer stores the starting index and particle count for each unique hash bucket, enabling direct lookups during physics or rendering passes.