GPU Sorting and Prefix Sums for Spatial Data

Almost every non-trivial spatial algorithm on a GPU decomposes into two primitives that have nothing spatial about them. A prefix sum turns a per-element pass/fail mask into a per-element output offset, which is what makes stream compaction, bucket construction and variable-length output possible at all. A sort imposes an ordering, and on spatial data that ordering is a space-filling curve, which turns scattered points into contiguous runs that share cache lines and bound the same subtree. Quadtree construction, tile binning, level-of-detail compaction and nearest-neighbour search are all compositions of these two. This page covers what each primitive does, the WGSL shapes they take, how they compose, and where their cost lands. It sits under spatial compute shaders and geometry pipelines.

Prerequisites

1. A compute pipeline and a storage buffer. Everything here is a compute pass writing storage buffers; the descriptor surface is covered under compute vs render pipeline fundamentals.
2. Workgroup storage. Both primitives use var<workgroup> arrays and barriers, which are subject to maxComputeWorkgroupStorageSize — 16 KiB is the guaranteed floor.
3. A tail guard habit. Both primitives run over arrays whose length is rarely a multiple of the workgroup size, and both corrupt silently without one.
4. Morton keys, for the spatial half. Interleaving two coordinates into one integer is the encoding that makes a linear sort spatially meaningful.

The prefix sum, and what it is for

An exclusive prefix sum over [1, 0, 1, 1, 0, 1] is [0, 1, 1, 2, 3, 3]: each element holds the total of everything before it. Applied to a pass/fail mask, that total is the output index — element three passed, and two elements before it passed, so it belongs at output slot two. One pass produces every output offset at once, with no contention and no atomics.

That is the difference between a prefix sum and the atomic counter used in simpler compaction. An atomic gives each survivor a unique slot; a prefix sum gives each survivor its slot, deterministically, in input order. For point clouds the order rarely matters and the atomic is simpler. For anything where the output has to be reproducible run to run — a test fixture, a picking index, a bucket structure another pass will binary-search — the prefix sum is required.

The workgroup-level form is a Hillis–Steele scan: log₂(n) rounds of “add the value from n lanes back”, separated by barriers.

const WG : u32 = 256u;
var<workgroup> scratch : array<u32, WG>;

@group(0) @binding(0) var<storage, read>       flags : array<u32>;
@group(0) @binding(1) var<storage, read_write> offsets : array<u32>;

@compute @workgroup_size(256)
fn scan(@builtin(global_invocation_id) gid : vec3<u32>,
        @builtin(local_invocation_id)  lid : vec3<u32>) {
  let i = gid.x;
  let t = lid.x;
  scratch[t] = select(0u, flags[i], i < arrayLength(&flags));
  workgroupBarrier();

  // Hillis-Steele: log2(WG) rounds, each reading a value the previous round wrote.
  for (var stride : u32 = 1u; stride < WG; stride = stride << 1u) {
    var add : u32 = 0u;
    if (t >= stride) { add = scratch[t - stride]; }
    workgroupBarrier();                 // everyone has read before anyone writes
    scratch[t] = scratch[t] + add;
    workgroupBarrier();                 // everyone has written before the next read
  }

  // Convert the inclusive result to exclusive: subtract this element's own flag.
  if (i < arrayLength(&offsets)) {
    offsets[i] = scratch[t] - select(0u, flags[i], i < arrayLength(&flags));
  }
}

The two barriers inside the loop are both required and are the most common thing to get wrong. Without the first, a lane may write scratch[t] while another is still reading it; without the second, a lane may read the next round’s value before it exists. A single-barrier version passes on some hardware and produces wrong offsets on others, which is the worst possible failure profile.

Scanning an array larger than one workgroup

A workgroup scan handles 256 elements. Scanning millions takes three passes, and the structure is worth knowing because it recurs in every GPU library that implements it.

The first pass scans each block of 256 independently and writes both the per-element offsets within the block and the block’s total. The second pass scans the array of block totals — which is n / 256 elements, small enough for a single workgroup at moderate sizes, or recursively scanned when it is not. The third pass adds each block’s scanned total back into every element of that block, producing global offsets.

The cost is three passes over the data rather than one, which sounds worse than it is: passes two and three are bandwidth-bound and the middle one touches a four-hundredth of the array. In practice a full scan over four million elements costs a small multiple of a single pass.

The sort, and why it is a Morton sort

Sorting on a GPU means a radix sort: several passes, each stable, each sorting on a few bits of the key using — unsurprisingly — a prefix sum to compute output positions. Four passes of eight bits sorts a 32-bit key. It is the expensive primitive on this page, and everything that uses it is designed to sort once and reuse the ordering.

What makes it spatial is the key. Interleaving the bits of a quantised x and y produces a Morton code, and Morton codes have the property that numerically adjacent values are spatially adjacent. Sorting by that key therefore produces an array where neighbouring entries are neighbouring points — which is exactly what makes a quadtree node a contiguous range, a tile’s features a contiguous range, and a cache line worth fetching.

The last row is the reason both belong in the same page: nearly every spatial structure is a sort followed by a scan, and knowing that is more useful than knowing either algorithm in isolation.

Where the cost lands

Both primitives are bandwidth-bound rather than arithmetic-bound, which has three practical consequences.

The first is that key width matters more than element count. A 64-bit Morton key needs eight radix passes rather than four, so halving the key — by quantising coordinates to 16 bits per axis rather than 32 — halves the sort. For viewport-scale work, 16 bits per axis is about 0.3 metres at zoom 14, which is finer than most vector tile sources carry anyway.

The second is that the scan’s workgroup size should sit in the middle of the occupancy band, because it is memory-bound and its var<workgroup> array grows with the size. 256 is the usual answer and 512 is worth measuring.

The third is that sorting once and reusing beats re-sorting. A point cloud whose ordering is a property of the data, not of the camera, is sorted at load and never again; a per-frame sort by depth for transparency is a different and much more expensive proposition, and is usually better replaced by an order-independent technique.

Composing the two into a spatial structure

The reason both primitives live on one page is that they are almost always used together, and the composition follows the same shape in every algorithm built on them.

Quadtree construction is the canonical case. Quantise each point’s coordinates and interleave the bits into a Morton key; radix sort by that key; scan a flag array marking every position where the top 2k bits change, which is where one node at depth k ends and the next begins; the scan’s output is the node table. Three primitives, no tree, no pointer, no allocation — and the resulting structure is a sorted array plus an index, which is far friendlier to a GPU than a pointer graph would be.

Tile binning is the same shape with a different key. Replace the Morton code with a tile index, sort, and scan for boundaries; the result is, for each tile, a contiguous range of the sorted array. A pass that then processes tiles in parallel reads a range per workgroup with perfect coalescing.

Deterministic compaction is the scan alone. A predicate produces a mask, the scan turns the mask into offsets, and a scatter writes survivors to their offsets. This is where the scan earns its cost over an atomic counter: the output order matches the input order exactly, every run.

Nearest-neighbour search uses both indirectly — it wants the sorted order for locality and the node table from the boundary scan to bound its search — which is why it is usually built on top of the quadtree rather than from the primitives directly.

The pattern worth extracting is that neither primitive is spatial. The spatial content is entirely in the key, and choosing the key is where the design work is. Everything after it is the same two passes regardless of whether the data is points, polygon centroids, or tile addresses.

Buffer budgeting for the two primitives

Both primitives need scratch space, and sizing it is a decision that has to be made at allocation time rather than discovered at dispatch time.

A multi-block scan needs one buffer of block sums — the element count divided by the workgroup size, so four million elements at 256 per block is 15 625 u32 values, about 62 kilobytes. That is small enough to allocate once at start-up for the largest array the pipeline will ever scan and forget about. If the scan recurses, each level needs its own, but the sizes fall by a factor of the workgroup size each time, so the whole chain is a rounding error on the first level.

A radix sort needs considerably more: a second key buffer and a second payload buffer of the same size as the originals, because each pass reads from one and writes to the other. Sorting four million points with a 32-bit key and a 32-bit payload therefore needs 32 megabytes of working space on top of the 32 megabytes the data occupies. Budgeting for one buffer and discovering the need for two mid-implementation is a common and avoidable surprise.

Both scratch allocations should live in the same pool as everything else and be counted in the same byte budget. The temptation is to treat them as transient — allocated per sort, freed after — and that is exactly the churn that fragments an allocator over a long session. Allocate once for the worst case, reuse, and never free until the pipeline is torn down.

Failure modes and diagnostics

• Offsets are wrong on one machine and right on another. A missing barrier in the scan loop. Both barriers are required; hardware that happens to execute a workgroup in lockstep hides the bug.
• The last partial workgroup corrupts the output. No tail guard. Out-of-bounds reads return zero, which silently contributes a zero to the scan and shifts every later offset.
• The sort is correct but the locality is not. The key interleaved the bits in the wrong order, or quantised x and y over different extents. Detection: sorted neighbours are spatially far apart.
• A sort over floats produces nonsense. Radix sort operates on unsigned integers; negative floats compare backwards under a bit-pattern sort. Fix: flip the sign bit for positives and invert entirely for negatives before sorting.
• Memory use spikes during the sort. A radix sort needs a second buffer of the same size to ping-pong between. Budget for two, not one.

When not to reach for either

Both primitives are general enough to be over-applied, and two cases are worth naming.

The first is a small array. A scan over four hundred elements is one workgroup and a handful of microseconds, but so is the dispatch overhead around it, and a CPU loop over four hundred elements finishes before the command buffer is submitted. The threshold moves with hardware; below a few thousand elements it is worth measuring rather than assuming.

The second is a case where an atomic counter is genuinely sufficient. Compaction whose output order does not matter — point clouds, particle survivors, anything the renderer draws as an unordered set — is exactly what an atomic is for, and it is one pass rather than three. The scan earns its extra cost only when the output has to be reproducible or when the offsets are needed as data by a later pass, and paying it by reflex is a real waste on a hot path.

There is also a third case that looks like a job for these primitives and is not: sorting for transparency. Depth-sorting translucent geometry per frame is a full radix sort every frame, which is the most expensive thing on this page applied at the highest possible frequency. Order-independent transparency techniques exist precisely to avoid it, and on a map with translucent overlays they are almost always the better answer.

Verifying a scan or a sort

Both primitives are easy to get subtly wrong and easy to verify exactly, which is an unusually good combination and worth exploiting.

A scan has a single global invariant: the last output plus the last input equals the sum of all inputs. Reading back three numbers — the final offset, the final flag and a separately computed total — is enough to catch every off-by-one and every missing barrier that affects the result. It costs one small readback in a development build and can run every frame without being noticed.

A sort has an equally simple one: the output is a permutation of the input, and it is non-decreasing. Checking non-decreasing is a single pass with a comparison per element; checking permutation is a checksum — sum the keys before and after, and compare. Together they catch the two failure modes that matter, a comparison that sorts backwards and a scatter that drops or duplicates elements.

Both checks belong behind a flag rather than in production, and both are worth running against a fixture in continuous integration with a CPU reference implementation for comparison. A few thousand elements is enough; the bugs in these primitives are structural rather than scale-dependent, and a small fixture finds them just as reliably as a large one while running in milliseconds.

Continue in this section

• Morton code generation for spatial sorting in WGSL — the bit interleave, its quantisation, and the extent it depends on.
• Implementing a workgroup prefix sum for stream compaction — the scan in full, including the multi-block form.
• Radix sorting point clouds on the GPU — the four-pass sort, the ping-pong buffers, and float key handling.
• Binning features into tiles with atomic counters — the cheaper alternative when order does not matter.

Related

• Spatial compute shaders and geometry pipelines — the section this topic belongs to.
• WGSL spatial algorithms on the GPU — the algorithms these primitives compose into.
• Generating a quadtree on the GPU with WGSL — a sort followed by a scan, in one page.
• Point cloud LOD compaction with stream compaction — the compaction the scan makes deterministic.
• Workgroup occupancy optimization for spatial kernels — how to size the scan’s workgroup.