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 points distributed over cells. Under a skewed distribution a single hot cell can absorb a fraction of all points, so a global-atomic scheme serializes increments on one memory location. That serialized run is the floor on pass time. Privatizing into workgroups, each with its own local histogram, caps the contention on any local bin at the workgroup size , and the subsequent global flush contends at most writers per bin. The expected serialization on the hottest global bin falls from

and because (256) and (points over ) are both far below 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 four phases of a privatized binning kernel Four phases run left to right inside one dispatch. The zero phase clears the workgroup-private bin array, with each lane clearing its own stride of bins, then a barrier. The accumulate phase has every lane read its point, compute a bin index from the coordinate, and atomically increment the private bin. A second barrier ensures every lane has finished. The merge phase has each lane atomically add one private bin into the corresponding global bin, so a workgroup contributes one global atomic per bin instead of one per point. ONE DISPATCH · FOUR PHASES Zero clear var<workgroup> bins strided across lanes Barrier workgroupBarrier() all clears visible Accumulate bin index from coordinate atomicAdd on private bins Merge one atomicAdd per bin into the global grid A second barrier before the merge is equally required, or a lane merges bins still being written.
The two barriers are not optional and they are not free — they are the price of the privatization. Skipping the first one is the classic bug: a lane starts accumulating into bins another lane has not finished clearing, and the heatmap comes out with random hot cells.

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.

wgsl
@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.

typescript
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

How many bins one workgroup can privatize A table of bin counts against workgroup storage cost and verdict, assuming a 16 kibibyte workgroup storage limit and four bytes per atomic u32 bin. 256 bins cost one kibibyte and fit comfortably. 1024 bins cost four kibibytes and still fit. 4096 bins cost sixteen kibibytes, exactly at the limit, and will fail on any device that reserves storage for anything else. 16384 bins cost 64 kibibytes and cannot be privatized at all, so that grid must be split into several passes over bin ranges. PRIVATE BINS · 4 B EACH · 16 KiB LIMIT Workgroup storage Verdict 256 bins 1 KiB comfortable 1024 bins 4 KiB fits 4096 bins 16 KiB at the limit 16384 bins 64 KiB impossible — split the pass Check maxComputeWorkgroupStorageSize on the adapter — 16 KiB is the guaranteed floor, not the value.
The bin count is a hard architectural limit, not a tuning knob. A heatmap grid finer than about 64 by 64 cells cannot be privatized in one pass, and the honest answer there is several passes over bin ranges rather than falling back to global atomics.
Parameter Typical value Guidance
TILE 32 Local histogram is TILE² × 4 bytes; 32 gives 4 KiB, safely under the 16 KiB maxComputeWorkgroupStorageSize floor. Raise only if occupancy allows.
@workgroup_size (WG) 256 Threads per workgroup; a multiple of 32/64 warp width. More threads zero and flush the local bins faster but raise register pressure.
WEIGHT_SCALE 256.0 Fixed-point factor encoding an f32 weight as u32. Lower it if a hot cell’s summed count risks overflowing u32.
GRID_W × GRID_H 1024 × 1024 Output resolution. VRAM scales with resolution, not point count — a 1024² u32 grid is 4 MiB.
Dispatch shape tilesX × tilesY One workgroup per TILE-sized region of the grid. Must cover the full grid, so round each axis up.
Local-bin count TILE² Bins per workgroup. Governs both shared-memory footprint and the number of global atomics in the flush.

Before tuning anything here, confirm the totals: sum the global bins on the CPU once and compare the result against the number of points dispatched. If those two numbers disagree, no amount of performance work matters yet, because the kernel is dropping data — and the size of the gap usually points straight at which of the failure modes below is responsible.

Failure modes

Reading a wrong heatmap from its pattern A table linking three heatmap defects to their cause and fix. Random hot cells scattered through the grid mean a missing barrier after the clear phase, and the fix is to add workgroupBarrier before accumulation. A grid whose totals are consistently low means bin indices are being clamped or dropped at the edges, and the fix is to clamp the coordinate rather than the index. A grid that is correct but blocky in a repeating pattern means the bin index is computed from a truncated coordinate, and the fix is to compute the index in double precision on the host or subtract a tile origin first. HEATMAP DEFECT · CAUSE · FIX Cause Fix random hot cells missing barrier workgroupBarrier after clear totals run low edge points dropped clamp the coordinate repeating blockiness f32 coordinate truncation subtract a tile origin Sum the global bins and compare against the input count — a mismatch localises the bug in seconds.
All three look like data problems and none of them are. The third in particular gets blamed on the source data for a long time, because a blocky heatmap is exactly what coarse input would produce.
  • 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.

Up: Spatial Aggregation in GPU Memory