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

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

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

  • 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.

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