Point-Cloud LOD Compaction with Stream Compaction

Rendering a hundred-million-point LiDAR cloud at every zoom level wastes bandwidth: most points project to the same pixel and cost fill rate for nothing. The fix is a level-of-detail subset — keep one point per screen-space cell, or every -th point along a Morton curve — but that subset has to be dense in a buffer for the draw to be efficient, and it must be rebuilt whenever the view changes. Stream compaction is the primitive that turns a sparse per-point keep-flag into a gap-free survivor buffer in two GPU passes: an exclusive prefix sum assigns each survivor a unique output slot, then a scatter writes it there. Unlike an atomic-counter cull, the prefix-sum scatter is race-free and order-preserving, so the thinned cloud keeps its spatial ordering and stays coherent for a following pass. This page gives the complete WGSL and TypeScript, the order-preserving counterpart to the atomic compaction in on-GPU viewport culling for vector tiles and a concrete case of the compaction skeleton in WGSL spatial algorithms on the GPU.

The core is the exclusive prefix sum. Given a flag array where marks whether point survives, the output slot for a survivor is the count of survivors strictly before it:

and the total survivor count is . Because each survivor’s offset is unique, the scatter dst[offset_i] = src[i] never collides — no atomic is needed once the scan is correct.

Runnable reference implementation

A sparse survivor mask collapsing into a dense prefix Two strips of sixteen slots each. The upper strip is the survivor mask after the level-of-detail predicate: slots 0, 2, 3, 6, 9, 12 and 13 are marked keep and the other nine are marked drop. The lower strip is the compacted output: the seven surviving indices have been written into slots 0 through 6 in their original order, and slots 7 through 15 are unused capacity that was allocated for the worst case. The draw call reads only the dense prefix, so nine ninths of the dropped work never reaches the rasterizer. SLOT 0 2 4 6 8 10 12 14 16 0 · 2 3 · · 6 · · 9 · · 12 13 · · Survivor mask 7 of 16 kept the predicate is evaluated for every point — only the write is conditional 0 2 3 6 9 12 13 unused capacity — allocated for the worst case Compacted dense prefix the draw call reads the prefix only; the tail is never touched The atomic counter that hands out prefix slots is also the draw call’s vertex count — write it to the indirect buffer.
The buffer is allocated for the worst case and used for the actual case, which is what lets the whole operation run without a CPU round trip: nothing has to know the survivor count until the indirect draw reads it.

Three WGSL entry points. flag_lod marks survivors: it keeps the first point that lands in each screen-space cell, which decimates dense regions while preserving sparse detail. scan_block is the work-efficient Blelloch prefix sum over a block, emitting per-point offsets and a per-block total. scatter_lod then adds each block’s base (itself the prefix over block totals) and writes survivors densely.

wgsl
// ---- Pass 1: flag one survivor per screen-space cell ----
@group(0) @binding(0) var<storage, read>       points:   array<vec4<f32>>;   // xyz + intensity
@group(0) @binding(1) var<storage, read_write> flags:    array<u32>;
@group(0) @binding(2) var<storage, read_write> occupied: array<atomic<u32>>; // one flag per cell
@group(0) @binding(3) var<uniform>             lod:      vec4<f32>;           // min_x,min_y,span,cells

@compute @workgroup_size(256)
fn flag_lod(@builtin(global_invocation_id) gid: vec3<u32>) {
    let i = gid.x;
    if (i >= arrayLength(&points)) { return; }            // guard the ragged tail
    let p = points[i].xy;
    let cells = u32(lod.w);
    let n = (p - lod.xy) / lod.z;                         // normalize to [0,1)
    if (n.x < 0.0 || n.x >= 1.0 || n.y < 0.0 || n.y >= 1.0) { flags[i] = 0u; return; }
    let cx = u32(n.x * f32(cells));
    let cy = u32(n.y * f32(cells));
    let cell = cy * cells + cx;
    // First writer of this cell wins; the rest are decimated. exchange returns the old value.
    let prev = atomicExchange(&occupied[cell], 1u);
    flags[i] = select(1u, 0u, prev != 0u);               // 1 only if this thread claimed the cell
}

// ---- Pass 2: exclusive Blelloch scan over a 512-element block ----
const BLOCK: u32 = 256u;                                  // 2*BLOCK = 512 elements per block
var<workgroup> temp: array<u32, BLOCK * 2u>;

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

@compute @workgroup_size(BLOCK)
fn scan_block(@builtin(local_invocation_index) lid: u32,
              @builtin(workgroup_id) wg: vec3<u32>) {
    let base = wg.x * BLOCK * 2u;
    temp[2u * lid]      = flags[base + 2u * lid];
    temp[2u * lid + 1u] = flags[base + 2u * lid + 1u];

    var off = 1u;
    for (var d = BLOCK; d > 0u; d >>= 1u) {               // up-sweep
        workgroupBarrier();
        if (lid < d) {
            let ai = off * (2u * lid + 1u) - 1u;
            let bi = off * (2u * lid + 2u) - 1u;
            temp[bi] += temp[ai];
        }
        off <<= 1u;
    }
    if (lid == 0u) {
        blockTotal[wg.x] = temp[BLOCK * 2u - 1u];         // this block's survivor count
        temp[BLOCK * 2u - 1u] = 0u;                       // seed the down-sweep
    }
    for (var d = 1u; d < BLOCK * 2u; d <<= 1u) {          // down-sweep
        off >>= 1u;
        workgroupBarrier();
        if (lid < d) {
            let ai = off * (2u * lid + 1u) - 1u;
            let bi = off * (2u * lid + 2u) - 1u;
            let t = temp[ai];
            temp[ai] = temp[bi];
            temp[bi] += t;
        }
    }
    workgroupBarrier();
    offsets[base + 2u * lid]      = temp[2u * lid];       // exclusive per-point offset
    offsets[base + 2u * lid + 1u] = temp[2u * lid + 1u];
}

// ---- Pass 3: scatter survivors into a dense, order-preserving buffer ----
@group(0) @binding(0) var<storage, read>       flags:     array<u32>;
@group(0) @binding(1) var<storage, read>       offsets:   array<u32>;
@group(0) @binding(2) var<storage, read>       blockBase: array<u32>;        // prefix over blockTotal
@group(0) @binding(3) var<storage, read>       points:    array<vec4<f32>>;
@group(0) @binding(4) var<storage, read_write> lodPoints: array<vec4<f32>>;

const BLOCK: u32 = 256u;

@compute @workgroup_size(BLOCK * 2u)
fn scatter_lod(@builtin(global_invocation_id) gid: vec3<u32>) {
    let i = gid.x;
    if (i >= arrayLength(&flags)) { return; }
    if (flags[i] == 0u) { return; }                       // decimated point
    let slot = blockBase[i / (BLOCK * 2u)] + offsets[i];  // unique, collision-free slot
    lodPoints[slot] = points[i];                          // preserves input order
}

The TypeScript driver clears the per-cell occupancy grid each frame, runs flag, scan, a short prefix over the block totals, then scatter. The blockBase prefix is a single small scan over blockTotal — for up to a few thousand blocks a one-workgroup scan covers it.

typescript
function buildLod(
  device: GPUDevice,
  pipelines: {
    flag: GPUComputePipeline; scan: GPUComputePipeline;
    blockScan: GPUComputePipeline; scatter: GPUComputePipeline;
  },
  binds: {
    flag: GPUBindGroup; scan: GPUBindGroup;
    blockScan: GPUBindGroup; scatter: GPUBindGroup;
  },
  occupied: GPUBuffer,
  lod: Float32Array,           // [min_x, min_y, span, cellsPerAxis]
  lodUbo: GPUBuffer,
  pointCount: number,
): void {
  device.queue.writeBuffer(lodUbo, 0, lod);
  const groups = Math.ceil(pointCount / 512);            // 512 elements per scan block

  const enc = device.createCommandEncoder();
  enc.clearBuffer(occupied);                             // reset cell claims; stale claims drop points

  const p1 = enc.beginComputePass();
  p1.setPipeline(pipelines.flag); p1.setBindGroup(0, binds.flag);
  p1.dispatchWorkgroups(Math.ceil(pointCount / 256)); p1.end();

  const p2 = enc.beginComputePass();
  p2.setPipeline(pipelines.scan); p2.setBindGroup(0, binds.scan);
  p2.dispatchWorkgroups(groups); p2.end();               // per-block offsets + block totals

  const p3 = enc.beginComputePass();
  p3.setPipeline(pipelines.blockScan); p3.setBindGroup(0, binds.blockScan);
  p3.dispatchWorkgroups(1); p3.end();                    // prefix over block totals → blockBase

  const p4 = enc.beginComputePass();
  p4.setPipeline(pipelines.scatter); p4.setBindGroup(0, binds.scatter);
  p4.dispatchWorkgroups(groups); p4.end();               // dense, ordered survivors

  device.queue.submit([enc.finish()]);
  // lodPoints is now bindable as vertex data — no readback of the survivor count.
}

Parameter reference

Points surviving each level of detail at four zoom levels A bar chart in millions of points showing how many of a twelve million point cloud survive the level-of-detail predicate at four zoom levels. At zoom 8 the predicate keeps about 0.4 million points, at zoom 11 about 1.9 million, at zoom 14 about 5.6 million, and at zoom 17, where the viewport covers a small area at full density, about 11.2 million of the resident points survive but the viewport itself holds far fewer. SURVIVORS BY ZOOM · MILLIONS OF POINTS Zoom 8 0.4 M Zoom 11 1.9 M Zoom 14 5.6 M Zoom 17 11.2 M 0 2 4 6 8 10 12 M comfortable watch the budget needs viewport culling too Chain the two passes: cull to the viewport, then compact by level of detail on the survivors.
Compaction alone stops helping at high zoom, because the predicate keeps nearly everything. That is the point at which the viewport cull has to run first, so the compaction is deciding detail among visible points rather than among all of them.

Every tunable value in the implementation above. These tables scroll horizontally on narrow viewports.

Parameter Typical value Guidance
cellsPerAxis 256–2048 Screen-space decimation grid. More cells keep more points; match it to viewport pixels over point radius.
BLOCK 256 (512 elements) Scan block half-width. Must be a power of two; the shared temp array is 2 × BLOCK u32 = 2 KiB.
occupied size cellsPerAxis² u32 One atomic flag per cell. A 1024² grid is 4 MiB and independent of point count.
@workgroup_size (scan) 256 Fixed by BLOCK; the scatter uses 2 × BLOCK = 512 to cover a full block per workgroup.
point stride 16 bytes vec4<f32> (xyz + intensity). Keep it 16-byte aligned; a vec3 payload still pads to 16.
blockBase scan 1 workgroup Valid while blocks ≤ 2 × workgroup_size; beyond that, recurse the block scan.
lodPoints capacity ≈ cellsPerAxis² At most one survivor per cell, so the output never exceeds the cell count.

Failure modes

Compaction bugs and the artefact each one produces A table of three compaction defects with cause and fix. Points disappearing at high density mean the atomic counter saturated past the buffer capacity and the out-of-bounds writes were silently discarded, and the fix is to allocate for the worst case or guard the write. A visibly shuffled point order means the scatter used the counter value without accounting for the fact that atomic order is not deterministic, and the fix is to accept unordered output or use a prefix sum instead. A draw call rendering stale geometry means the indirect buffer was written after the draw was recorded, and the fix is to write it in the same encoder before the render pass. COMPACTION DEFECT · CAUSE · FIX Cause Fix points vanish at density counter past capacity allocate worst case point order shuffles atomic order is not stable use a prefix sum stale geometry drawn indirect written too late write before the pass Read the counter back once per second in development and compare it against the buffer capacity.
The middle row is worth accepting rather than fixing in most cases. Point clouds rarely care about ordering, and a prefix sum costs a pass that the atomic counter did not.
  • Colliding offsets, lost points. A scan_block missing a workgroupBarrier() between up-sweep and down-sweep emits repeated offsets, so two survivors scatter to the same slot and one is overwritten. Detection: lodPoints has duplicates and stale entries even though the survivor count looks right. Fix: keep both barrier phases; verify the exclusive scan of [1,0,1,1] gives [0,1,1,2].
  • Dropped survivors past the first block. Skipping the blockBase prefix, or scanning more than 2 × workgroup_size in one workgroup, leaves every block writing from offset 0, so blocks overwrite each other. Detection: only the first ~512 survivors appear. Fix: run the block-totals prefix (Pass 3) and add blockBase[block] in the scatter.
  • Stale decimation across frames. Not clearing occupied leaves last frame’s cell claims, so every cell reads as already taken and flag_lod keeps nothing. Detection: the thinned cloud is empty after the first frame. Fix: clearBuffer(occupied) before the flag pass.
  • Device lost on a huge scan. A single dispatch over a hundred million points can exceed maxComputeWorkgroupsPerDimension or trip the TDR watchdog. Detection: device.lost resolves, or the dispatch is silently clamped. Fix: cap points per dispatch and recover via the browser support fallback routing strategies; size dispatches against the limits negotiated during WebGPU device initialization for GIS workloads.

Worth stating plainly before the failure modes: the output buffer here is sized for the worst case and used for the actual case, so its allocated length and its meaningful length are different numbers that live in different places. The length is in the atomic counter; the capacity is in the buffer descriptor. Nearly every bug on this page comes from code that read one and meant the other.

Backend / Python interop note

The point payload must arrive as a contiguous vec4<f32> array so the scatter writes and the downstream draw read the same 16-byte stride; a packed xyz-only vec3 layout drifts off alignment and corrupts the compacted output, the failure the memory alignment for spatial data buffers reference details. Pre-sort the cloud into Morton order server-side so that even before compaction, points that land in the same screen cell are contiguous in the buffer — the flag_lod first-writer rule then keeps a spatially representative point per cell rather than an arbitrary one, and the resulting subset is coherent for streaming. Cap each streamed batch to the reserved staging size using the ring-buffer discipline in reducing GPU memory fragmentation during spatial aggregation.

Up: WGSL Spatial Algorithms on the GPU