Parallel Histogram Binning for GPU Heatmaps

The bottleneck this page targets is atomic contention on a spatial histogram. A heatmap bins millions of points into a 2D grid by hashing each point to a cell and incrementing that cell’s counter — and when thousands of points land in the same hot cell, a dense city centre against a sparse rural extent, every one of those increments serializes on a single global atomicAdd. The arithmetic is trivial; the pass time is dominated entirely by how long the hardware queues those contending atomics against one memory location. On a skewed distribution a naive global-atomic histogram can run an order of magnitude slower than the memory bandwidth would allow, because the hottest cell alone dictates the schedule. The fix is privatization: give each workgroup its own copy of the histogram (or a tile of it) in fast var<workgroup> memory, accumulate there where contention is bounded by the workgroup size rather than the whole dispatch, then flush each workgroup’s partial histogram to the global grid with one atomic per non-empty bin. Contention on any global bin drops from the total point count to the number of workgroups, which is what makes the pass bandwidth-bound instead of contention-bound.

This is the binning core of spatial aggregation in GPU memory, specialized for the heatmap output: where that reference covers the full scatter-normalize-smooth chain, this page zooms into making the scatter itself scale under skew. It assumes packed Structure-of-Arrays inputs from the memory alignment for spatial data buffers rules, and the privatized flush is a direct application of the workgroup occupancy optimization for spatial kernels trade-off between shared-memory footprint and thread count.

Why privatization cuts contention

Model the scatter as $N$ points distributed over $C$ cells. Under a skewed distribution a single hot cell can absorb a fraction $\rho$ of all points, so a global-atomic scheme serializes $\rho N$ increments on one memory location. That serialized run is the floor on pass time. Privatizing into $G$ workgroups, each with its own local histogram, caps the contention on any local bin at the workgroup size $W$, and the subsequent global flush contends at most $G$ writers per bin. The expected serialization on the hottest global bin falls from

and because $W$ (256) and $G$ (points over $W$) are both far below $\rho N$ at GIS scale, the hot-cell penalty effectively disappears. The cost paid for it is the var<workgroup> storage for the local histogram plus two workgroupBarrier synchronizations — cheap against the serialization it removes.

Runnable reference implementation

The kernel tiles the histogram: each workgroup owns a TILE × TILE local sub-histogram in workgroup memory. Points are pre-binned by their owning workgroup (partitioned by spatial tile on the host side), so every point a workgroup processes falls inside that workgroup’s tile. Threads cooperatively zero the local bins, accumulate with fast workgroup-local atomics, barrier, then flush non-empty bins to the global grid.

@group(0) @binding(0) var<storage, read>       positions: array<vec2<f32>>;
@group(0) @binding(1) var<storage, read>       weights:   array<f32>;
@group(0) @binding(2) var<storage, read_write> grid:      array<atomic<u32>>; // global histogram
@group(0) @binding(3) var<uniform>             extent:    vec4<f32>;          // min_x,min_y,max_x,max_y

const GRID_W: u32 = 1024u;
const GRID_H: u32 = 1024u;
const TILE: u32 = 32u;                 // 32*32 u32 local bins = 4 KiB, well within budget
const WG: u32 = 256u;
const WEIGHT_SCALE: f32 = 256.0;       // fixed-point: no atomic<f32> exists

// One private histogram per workgroup, in fast shared memory.
var<workgroup> local_bins: array<atomic<u32>, TILE * TILE>;

@compute @workgroup_size(WG)
fn histogram(@builtin(global_invocation_id) gid: vec3<u32>,
             @builtin(workgroup_id)         wid: vec3<u32>,
             @builtin(local_invocation_index) lid: u32) {
    // 1. Cooperatively zero the private histogram before any accumulation.
    for (var b = lid; b < TILE * TILE; b += WG) {
        atomicStore(&local_bins[b], 0u);
    }
    workgroupBarrier();                 // all local bins zeroed before scatter

    // 2. Hash this point to a global cell, then to a tile-local bin.
    let i = gid.x;
    if (i < arrayLength(&positions)) {
        let p = positions[i];
        let span = vec2<f32>(extent.z - extent.x, extent.w - extent.y);
        let n = (p - extent.xy) / span;
        if (n.x >= 0.0 && n.x < 1.0 && n.y >= 0.0 && n.y < 1.0) {
            let gx = u32(n.x * f32(GRID_W));
            let gy = u32(n.y * f32(GRID_H));
            // Tile origin for this workgroup; points are partitioned so gx,gy fall inside it.
            let ox = wid.x * TILE;
            let oy = wid.y * TILE;
            let lx = gx - ox;
            let ly = gy - oy;
            if (lx < TILE && ly < TILE) {
                let amount = u32(weights[i] * WEIGHT_SCALE);
                atomicAdd(&local_bins[ly * TILE + lx], amount); // cheap, low-contention
            }
        }
    }
    workgroupBarrier();                 // every local add visible before the flush

    // 3. Flush non-empty local bins to the global grid: one global atomic each.
    for (var b = lid; b < TILE * TILE; b += WG) {
        let v = atomicLoad(&local_bins[b]);
        if (v != 0u) {
            let lx = b % TILE;
            let ly = b / TILE;
            let gx = wid.x * TILE + lx;
            let gy = wid.y * TILE + ly;
            if (gx < GRID_W && gy < GRID_H) {
                atomicAdd(&grid[gy * GRID_W + gx], v); // contention capped at workgroup count
            }
        }
    }
}

The host side dispatches one workgroup per tile of the grid, partitions points into per-tile batches so each workgroup only sees points inside its tile, and clears the global grid before the pass. Read the resulting grid as fixed-point counts and divide by WEIGHT_SCALE in the downstream normalize pass.

const GRID_W = 1024, GRID_H = 1024, TILE = 32;
const tilesX = Math.ceil(GRID_W / TILE);   // 32 workgroups across
const tilesY = Math.ceil(GRID_H / TILE);   // 32 workgroups down

// One u32 accumulator per global cell; cleared each frame, never reallocated.
const grid = device.createBuffer({
  size: GRID_W * GRID_H * 4,
  usage: GPUBufferUsage.STORAGE | GPUBufferUsage.COPY_SRC | GPUBufferUsage.COPY_DST,
});

function binHistogram(
  encoder: GPUCommandEncoder,
  pipeline: GPUComputePipeline,
  bindGroup: GPUBindGroup,
): void {
  encoder.clearBuffer(grid);                 // stale counts brighten the map forever
  const pass = encoder.beginComputePass();
  pass.setPipeline(pipeline);
  pass.setBindGroup(0, bindGroup);
  // 2D dispatch: one workgroup per grid tile, so wid.x/wid.y index the tile.
  pass.dispatchWorkgroups(tilesX, tilesY);
  pass.end();
}

Recording this into the same command encoder as the downstream normalize and smooth passes lets the driver honor the buffer dependency without an explicit barrier. Sequencing the whole chain off the render-blocking path is the job of pipelining compute dispatches with double buffering.

Parameter reference

Failure modes

• Points binned into the wrong tile. If the host partition does not match the wid-derived tile origin in the shader, a point whose global cell lies outside the workgroup’s tile is silently dropped by the lx < TILE && ly < TILE guard. Detection: total binned weight is lower than the input weight sum, and tile seams show as faint dark grid lines. Fix: partition points on the host by the identical floor(cell / TILE) rule the shader uses, and assert per-tile point counts sum to the input.
• Atomic counter overflow in a hot cell. A dense cell’s fixed-point sum can exceed u32 range and wrap toward zero, so the hottest region reads as cold. Detection: the densest area shows a dark hole surrounded by a bright ring. Fix: lower WEIGHT_SCALE, or accumulate the flush with a saturating check, or narrow TILE so each workgroup contributes a smaller partial sum.
• Missing barrier causing undercounts. Dropping either workgroupBarrier — before scatter or before flush — lets some threads read local bins the others have not finished writing, or flush before zeroing completes. Detection: cell totals differ run to run for identical input. Fix: keep both barriers exactly as shown; they are the correctness contract of privatization.
• Local histogram exceeds workgroup storage. Setting TILE too large pushes TILE² × 4 past maxComputeWorkgroupStorageSize, so pipeline creation fails. Detection: a GPUValidationError at createComputePipeline naming workgroup storage. Fix: keep TILE² × 4 ≤ 16384, or query device.limits.maxComputeWorkgroupStorageSize and size TILE against it.

Related

• Spatial Aggregation in GPU Memory — the full scatter, normalize, and smooth chain this binning feeds.
• Reducing GPU Memory Fragmentation During Spatial Aggregation — keep the grid and staging buffers churn-free across zoom levels.
• Workgroup Occupancy Optimization for Spatial Kernels — balance the shared-memory footprint of the local histogram against occupancy.
• Memory Alignment for Spatial Data Buffers — the packed position and weight layout the kernel reads.

Up: Spatial Aggregation in GPU Memory