Radix Sorting Point Clouds on the GPU

A radix sort is the only practical way to order millions of points on a GPU, and it is built entirely from the primitive on the previous page: each pass histograms one digit of the key, scans the histogram to get output offsets, and scatters every element to its offset. Four passes of eight bits sorts a 32-bit Morton key, each pass is stable, and the result is a point cloud whose memory order matches its spatial order. This page is the pass structure, the buffers it needs, the float-key transformation, and the places where a sort silently stops being stable. It is one stage of GPU sorting and prefix sums for spatial data.

Runnable reference implementation

Each pass is a histogram, a scan, and a scatter. The histogram is a privatised count per workgroup, exactly as in parallel histogram binning, because 256 buckets is small enough to hold in workgroup memory.

const RADIX : u32 = 256u;                  // one 8-bit digit per pass
const WG    : u32 = 256u;

struct PassInfo { shift : u32 };           // 0, 8, 16, 24 across the four passes
@group(0) @binding(0) var<uniform> info : PassInfo;
@group(0) @binding(1) var<storage, read>       keys_in  : array<u32>;
@group(0) @binding(2) var<storage, read_write> counts   : array<atomic<u32>>;

var<workgroup> local_counts : array<atomic<u32>, RADIX>;

@compute @workgroup_size(256)
fn histogram(@builtin(global_invocation_id) gid : vec3<u32>,
             @builtin(local_invocation_id)  lid : vec3<u32>) {
  atomicStore(&local_counts[lid.x], 0u);
  workgroupBarrier();

  let i = gid.x;
  if (i < arrayLength(&keys_in)) {
    let digit = (keys_in[i] >> info.shift) & (RADIX - 1u);
    atomicAdd(&local_counts[digit], 1u);   // contention stays inside the workgroup
  }
  workgroupBarrier();

  // One global atomic per bucket per workgroup, not one per element.
  atomicAdd(&counts[lid.x], atomicLoad(&local_counts[lid.x]));
}

The scatter reads the scanned offsets and writes each element to its destination. The key detail is that the destination buffer is the other half of a ping-pong pair, because a sort cannot be done in place.

@group(0) @binding(0) var<uniform> info : PassInfo;
@group(0) @binding(1) var<storage, read>       keys_in     : array<u32>;
@group(0) @binding(2) var<storage, read>       payload_in  : array<u32>;
@group(0) @binding(3) var<storage, read>       offsets     : array<u32>;   // scanned
@group(0) @binding(4) var<storage, read_write> keys_out    : array<u32>;
@group(0) @binding(5) var<storage, read_write> payload_out : array<u32>;
@group(0) @binding(6) var<storage, read_write> cursor      : array<atomic<u32>>;

@compute @workgroup_size(256)
fn scatter(@builtin(global_invocation_id) gid : vec3<u32>) {
  let i = gid.x;
  if (i >= arrayLength(&keys_in)) { return; }

  let key   = keys_in[i];
  let digit = (key >> info.shift) & (RADIX - 1u);
  // Claim the next slot within this digit's range. Sequential within a digit,
  // which is what preserves stability across passes.
  let dst   = offsets[digit] + atomicAdd(&cursor[digit], 1u);
  keys_out[dst]    = key;
  payload_out[dst] = payload_in[i];
}

After each pass the two buffers swap roles. Four passes means four swaps, so the sorted result ends up back in the original buffer — a small convenience that disappears if the key width ever changes to an odd number of passes.

Parameter reference

Sorting keys that are not unsigned integers

A radix sort compares bit patterns, which is only equivalent to comparing values for unsigned integers. Morton codes already are unsigned, which is convenient and is one more reason to sort by them rather than by a raw coordinate. Where a float key is unavoidable — sorting by distance from the camera, say — the bit pattern has to be transformed first.

IEEE 754 floats have the useful property that positive values compare correctly as unsigned integers, and the awkward property that negative values compare backwards, because the sign bit is set and the magnitude bits run the wrong way. The standard fix is one branchless expression:

/// Map a float's bit pattern to a u32 that sorts in the same order.
fn float_key(value : f32) -> u32 {
  let bits = bitcast<u32>(value);
  // If the sign bit is set, invert everything; otherwise set the sign bit.
  let mask = select(0x80000000u, 0xFFFFFFFFu, (bits & 0x80000000u) != 0u);
  return bits ^ mask;
}

The inverse, needed if the key has to be read back as a float, applies the opposite mask. In practice it rarely does — the payload carries the index, and the index is what later passes want.

A distance key also raises a design question worth answering before implementing it. Sorting by distance changes every frame the camera moves, which means a per-frame radix sort over the whole cloud: the most expensive primitive in the pipeline at the highest possible frequency. For transparency, order-independent techniques avoid it entirely. For level-of-detail selection, a Morton sort done once at load plus a per-frame filter achieves the same visual result at a fraction of the cost.

Where the four passes actually go

Instrumenting a radix sort with per-pass timestamps is worth doing once, because the distribution is not what most people expect.

The histogram passes are cheap: one read per element and an atomic into workgroup memory, which is close to a pure streaming read. The scans between them are cheaper still, because they operate on 256 counts rather than on the data. The scatter passes dominate, and they dominate for a reason that is structural rather than fixable — a scatter writes each element to an address determined by its data, so consecutive lanes write to addresses that are not consecutive, and the write coalescing that makes the histogram fast is unavailable.

That gives a rough split of one part histogram, four parts scatter, per pass. It has two practical consequences. Reducing the number of passes helps far more than optimising any individual one, which is the argument for a 32-bit key over a 64-bit one restated in measurements rather than in principle. And carrying a smaller payload helps proportionally, because the payload is scattered alongside the key: sorting an index rather than a full record, and gathering the records once at the end, moves four bytes per element per pass instead of thirty-two.

The second of those is the optimisation worth doing by default. Sort (key, index) pairs, then perform one gather pass at the end to reorder the actual data — five passes instead of four, each moving a fraction of the bytes.

Failure modes

• The sort is correct for one pass and wrong after four. A pass is not stable — usually because the scatter used a global atomic across all digits rather than a per-digit cursor, so elements with equal digits reorder and earlier passes’ work is undone.
• The output has holes. The cursor buffer was not reset between passes, so the second pass started counting from where the first finished and wrote past its range.
• Memory use doubles unexpectedly. The ping-pong buffers. They are not optional; budget for two of everything the sort touches.
• Negative distances sort to the wrong end. A float key used without the sign transformation.
• The sort is slower than expected on one vendor. The histogram’s global atomics are contending. Check that the privatised per-workgroup counts are actually being used — a version that skips them does one global atomic per element rather than one per bucket per workgroup.

Backend / Python interop note

Sorting on the server is worth considering seriously, because a point cloud’s spatial ordering is a property of the data rather than of the view, and a tile that arrives already Morton-ordered needs no sort on the client at all.

import numpy as np

def morton_sort(xs: np.ndarray, ys: np.ndarray, payload: np.ndarray):
    """Return the arrays reordered by Morton key, plus the keys themselves."""
    keys = morton_2d(quantise(xs), quantise(ys))     # see the Morton page
    order = np.argsort(keys, kind="stable")
    return keys[order], xs[order], ys[order], payload[order]

The wire-format benefit is real: a Morton-ordered coordinate column has neighbouring rows that share high bits, which any general-purpose compressor exploits, and delta encoding within the column becomes genuinely effective. A tile can shrink by a third or more purely from being sorted, before any other encoding decision.

The cost is that the ordering has to be preserved end to end. A pipeline that sorts server-side and then filters client-side must keep the survivors in order — which is exactly what the deterministic scan on the previous page provides and what an atomic counter does not. That is the clearest case where the scan’s extra pass pays for itself.

Related

• GPU sorting and prefix sums for spatial data — the topic this page belongs to.
• Morton code generation for spatial sorting in WGSL — the keys this sort orders.
• Implementing a workgroup prefix sum for stream compaction — the primitive each pass is built from.
• Generating a quadtree on the GPU with WGSL — what a sorted array becomes.
• VRAM budget management across tile zoom levels — where the ping-pong buffers have to be budgeted.