WGSL Spatial Algorithms on the GPU

The spatial index a CPU builds with a recursive descent — subdivide a node, sort points into children, recurse — is exactly the shape a GPU hates. Recursion, data-dependent branching, and a growing pointer tree are serial by construction, so a naive port stalls a single thread while thousands sit idle. Yet the algorithms themselves are embarrassingly parallel once you stop expressing them recursively: a quadtree is a sort by interleaved bit code, a viewport cull is one independent bounds test per feature, and building a level-of-detail subset is a scan followed by a scatter. This area is about rewriting the classic spatial toolbox — tree construction, culling, and compaction — as flat, index-free WGSL compute passes whose only sequential dependency is a prefix sum, so a million features are indexed or culled inside a single frame rather than across dozens of them.

These guides sit inside Spatial Compute Shaders & Geometry Pipelines and assume you already treat a compute pass as a first-class GPU program rather than a fragment shader pressed into service — the distinction is drawn in the compute versus render pipeline fundamentals reference. The algorithms here are the structural counterparts to the value-oriented passes elsewhere in this section: where geometry filtering with WGSL compute shaders drops features by predicate and spatial aggregation in GPU memory reduces them into grids, the passes below reorganize features into an index, a survivor list, or a thinned subset. All three lean on the same three primitives — an atomic counter, a parallel prefix sum, and an indirect dispatch — and all three are sensitive to the buffer stride rules in memory alignment for spatial data buffers, because a Morton code or a compaction cursor computed over a misaligned array indexes garbage.

Prerequisites

Before implementing the passes below, you should have:

1. A valid GPUDevice with maxComputeInvocationsPerWorkgroup, maxComputeWorkgroupStorageSize, and maxStorageBufferBindingSize inspected against your feature count — adapter negotiation is covered in initializing WebGPU devices for GIS workloads.
2. A working exclusive prefix sum. Every algorithm here compacts, and compaction needs a scan that maps each survivor’s flag to a unique output slot. A single-workgroup Blelloch scan is enough up to 2 × workgroup_size elements; beyond that you need the block-scan-then-add pattern shown in Step 3.
3. Understanding of WGSL atomic<u32>. The order-independent variant of compaction claims output slots with atomicAdd on a shared cursor. Only atomic<u32> and atomic<i32> exist — there is no atomic float.
4. Structure-of-Arrays inputs. Positions, bounds, and keys live in separate tightly packed arrays. Interleaving them shifts every element off its stride and corrupts the Morton code; the memory alignment for spatial data buffers reference explains the 16-byte rule that bites here.
5. Coordinates already projected into one linear space (Web Mercator metres or normalized tile space). Morton codes and frustum tests both assume a single coordinate system; mixing degrees with metres scrambles the bit interleave.

API and WGSL reference

The types and limits below govern every algorithm in this area. The prefix-sum workgroup budget and the indirect-argument layout are the two constraints that bite first at GIS scale.

These tables scroll horizontally on narrow viewports. For authoritative wording on atomic ordering, indirect arguments, and storage access, consult the W3C WebGPU specification and the WGSL specification.

Implementation walkthrough

Every algorithm here is a variation on one skeleton: a build pass computes a per-feature key, a scan pass turns per-feature keep-flags into output offsets, and a compact pass scatters survivors into a dense buffer. The steps below assemble that skeleton, then specialize it for tree construction and for indirect hand-off.

Step 1 — Compute a per-feature key in parallel

The build pass runs one invocation per feature with no cross-thread communication, so it scales linearly across the whole GPU. For a quadtree it emits a Morton code — the interleaved bits of the quantized x and y — which sorts points into z-order, the traversal that places spatially close points at adjacent indices. The same pass shape emits a single keep-flag for culling or thinning; only the body changes.

@group(0) @binding(0) var<storage, read>       positions: array<vec2<f32>>;
@group(0) @binding(1) var<storage, read_write> codes:     array<u32>;
@group(0) @binding(2) var<uniform>             extent:    vec4<f32>; // min_x,min_y,max_x,max_y

// Spread the low 16 bits of v into even bit positions (Morton part 1).
fn part1by1(x: u32) -> u32 {
    var v = x & 0x0000ffffu;
    v = (v | (v << 8u)) & 0x00ff00ffu;
    v = (v | (v << 4u)) & 0x0f0f0f0fu;
    v = (v | (v << 2u)) & 0x33333333u;
    v = (v | (v << 1u)) & 0x55555555u;
    return v;
}

@compute @workgroup_size(256)
fn build_codes(@builtin(global_invocation_id) gid: vec3<u32>) {
    let i = gid.x;
    if (i >= arrayLength(&positions)) { return; }  // guard the ragged tail workgroup

    let span = vec2<f32>(extent.z - extent.x, extent.w - extent.y);
    let n = (positions[i] - extent.xy) / span;      // normalize to [0,1)
    // Quantize to a 16-bit grid per axis, then interleave: 32-bit Morton key.
    let qx = u32(clamp(n.x, 0.0, 1.0) * 65535.0);
    let qy = u32(clamp(n.y, 0.0, 1.0) * 65535.0);
    codes[i] = part1by1(qx) | (part1by1(qy) << 1u);
}

The spatial-data-specific choice is quantizing to a fixed 16-bit grid per axis before interleaving: it fixes the tree’s maximum depth at 16 levels and makes every code directly comparable, which is what lets the subsequent sort be a plain radix sort over u32 rather than a geometric descent.

Step 2 — Flag survivors with an independent predicate

Culling and level-of-detail thinning both reduce to a boolean per feature: inside the viewport, or selected for this detail level. The build pass writes a 0/1 flag array that the scan consumes. Keeping the flag separate from the payload means the scan touches only 4 bytes per feature, not the whole struct.

@group(0) @binding(0) var<storage, read>       bounds: array<vec4<f32>>; // min_x,min_y,max_x,max_y
@group(0) @binding(1) var<storage, read_write> flags:  array<u32>;
@group(0) @binding(2) var<uniform>             view:   vec4<f32>;        // viewport extent

@compute @workgroup_size(256)
fn flag_visible(@builtin(global_invocation_id) gid: vec3<u32>) {
    let i = gid.x;
    if (i >= arrayLength(&bounds)) { return; }
    let b = bounds[i];
    // Separating-axis overlap test between feature AABB and viewport AABB.
    let outside = b.z < view.x || b.x > view.z || b.w < view.y || b.y > view.w;
    flags[i] = select(1u, 0u, outside);   // 1 = keep, 0 = cull
}

Step 3 — Scan the flags into output offsets

An exclusive prefix sum turns the flag array into per-feature output slots: offset[i] is the count of survivors strictly before i, and the final total is the survivor count. Small arrays scan in one workgroup; large arrays scan per block, then a second pass adds each block’s base. The single-block Blelloch scan below is the reusable core.

const BLOCK: u32 = 256u;
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<atomic<u32>>; // per-block sums

@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 offset = 1u;
    for (var d = BLOCK; d > 0u; d >>= 1u) {        // up-sweep (reduce)
        workgroupBarrier();
        if (lid < d) {
            let ai = offset * (2u * lid + 1u) - 1u;
            let bi = offset * (2u * lid + 2u) - 1u;
            temp[bi] += temp[ai];
        }
        offset <<= 1u;
    }
    if (lid == 0u) {
        atomicStore(&blockTotal[wg.x], temp[BLOCK * 2u - 1u]);
        temp[BLOCK * 2u - 1u] = 0u;                // clear the last element
    }
    for (var d = 1u; d < BLOCK * 2u; d <<= 1u) {   // down-sweep
        offset >>= 1u;
        workgroupBarrier();
        if (lid < d) {
            let ai = offset * (2u * lid + 1u) - 1u;
            let bi = offset * (2u * lid + 2u) - 1u;
            let t = temp[ai];
            temp[ai] = temp[bi];
            temp[bi] += t;
        }
    }
    workgroupBarrier();
    offsets[base + 2u * lid]      = temp[2u * lid];
    offsets[base + 2u * lid + 1u] = temp[2u * lid + 1u];
}

The two workgroupBarrier() phases are load-bearing: the up-sweep must complete before the last element is cleared, and every down-sweep level must publish its writes before the next reads them. Skipping either yields offsets that collide, so two survivors scatter to the same slot and one is silently lost.

Step 4 — Compact survivors into a dense buffer

With offsets in hand, compaction is a pure scatter: each kept feature copies its payload to outBase + offsets[i]. Because every survivor owns a distinct offset, the scatter is race-free without any atomic. The outBase is the running total of survivors in preceding blocks, added from the per-block sums the scan produced.

@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 block sums
@group(0) @binding(3) var<storage, read>       src:      array<vec4<f32>>;
@group(0) @binding(4) var<storage, read_write> dst:      array<vec4<f32>>;

const BLOCK: u32 = 256u;

@compute @workgroup_size(BLOCK * 2u)
fn compact(@builtin(global_invocation_id) gid: vec3<u32>) {
    let i = gid.x;
    if (i >= arrayLength(&flags)) { return; }
    if (flags[i] == 0u) { return; }                       // culled — nothing to write
    let slot = blockBase[i / (BLOCK * 2u)] + offsets[i];  // global output index
    dst[slot] = src[i];                                   // dense, gap-free write
}

Step 5 — Size the next pass from the GPU count with indirect dispatch

The survivor total is only known on the GPU, so writing it into an indirect-argument buffer lets the render pass draw exactly the survivors — no readback, no CPU stall. A tiny pass converts the count into either draw arguments (instanceCount) or dispatch arguments (workgroup count for a following compute pass over the compacted set).

@group(0) @binding(0) var<storage, read>       blockTotal: array<u32>;   // per-block survivor sums
@group(0) @binding(1) var<storage, read_write> drawArgs:   array<u32>;   // 4 x u32 indirect args

const BLOCKS: u32 = 1024u;   // number of scan blocks dispatched

@compute @workgroup_size(1)
fn write_indirect() {
    var total = 0u;
    for (var b = 0u; b < BLOCKS; b++) { total += blockTotal[b]; }
    drawArgs[0] = 6u;        // vertexCount: 6 verts per instanced quad
    drawArgs[1] = total;     // instanceCount: one instance per survivor
    drawArgs[2] = 0u;        // firstVertex
    drawArgs[3] = 0u;        // firstInstance
}

The TypeScript driver records build, scan, compact, and this indirect writer into one command encoder, then the render pass consumes drawArgs with pass.drawIndirect(drawArgs, 0). The buffer dependency chain is honored by WebGPU’s submission model without an explicit barrier, so the four compute passes and the draw run back to back.

Memory and performance implications

The scan is the only pass with a sequential dependency, and it dominates the analysis. A work-efficient Blelloch scan does $O(N)$ work in $O(\log N)$ depth, so over $N$ features the span is proportional to ${\log }_{2}N$ — roughly 20 levels for a million features — against the naive $O(N)$ sequential prefix. The build and compact passes are each one coalesced pass over the feature array with no shared-memory traffic, so their cost is bandwidth, not arithmetic: reading $N$ vec4<f32> bounds is $16N$ bytes, and a compacted output the same size again. Budget two feature-sized buffers (source and compacted destination) plus the small flag, offset, and per-block arrays, which are $4N$ bytes and $4\times \text{blocks}$ bytes respectively.

The indirect hand-off is what keeps the chain off the render-blocking path. Reading the survivor count back to the CPU to size a normal draw would force a queue flush and a PCIe round-trip every frame; drawIndirect reads the count from VRAM in the same submission, so the whole build-scan-compact-draw sequence is a single submit. @workgroup_size(256) is the portable default for the build and compact passes, and the scan block size must be a power of two so the up-sweep and down-sweep index cleanly. Picking the size that actually maximizes occupancy for your feature count and device class is the subject of workgroup occupancy optimization for spatial kernels, and confirming the scan is the bottleneck rather than the scatter is what frame profiling with timestamp queries measures directly. When the same feature set is rebuilt every frame, the buffer-recycling discipline from reducing GPU memory fragmentation during spatial aggregation keeps the source and destination allocations from churning.

Failure modes and diagnostics

• Colliding compaction offsets. A scan with a missing workgroupBarrier() between up-sweep and down-sweep produces offsets that repeat, so two survivors scatter to the same slot and one is overwritten. Detection: the compacted count is correct but the output has duplicates and holes. Fix: keep both barrier phases in scan_block (Step 3) and verify the exclusive scan of [1,0,1,1] yields [0,1,1,2].
• GPUValidationError — indirect buffer missing INDIRECT usage. A buffer written by compute but bound to drawIndirect without GPUBufferUsage.INDIRECT fails validation at draw time. Detection: synchronous error at drawIndirect. Fix: create the argument buffer with STORAGE | INDIRECT | COPY_DST.
• Morton codes clustered at one corner. Feeding un-normalized or mixed-unit coordinates into build_codes collapses the quantized grid, so most points share a code and the tree degenerates to a list. Detection: the code histogram spikes at 0 or 0xffffffff. Fix: project to one linear space before the build pass and confirm the extent uniform matches the data, per memory alignment for spatial data buffers.
• Stale survivor count in the indirect args. Reusing the argument buffer across frames without clearing the accumulator leaves last frame’s total, so the draw renders too many or too few instances. Detection: instance count lags the view by one frame. Fix: zero drawArgs (or recompute total from scratch as in Step 5) before each submission.
• Device lost on an oversized single-block scan. Trying to scan more than 2 × workgroup_size elements in one workgroup silently truncates, or a giant dispatch trips the driver TDR watchdog. Detection: device.lost resolves, or survivors past the first block vanish. Fix: use the block-scan-then-add pattern (Steps 3–4) and cap per-dispatch work; recover through the browser support fallback routing strategies.

Continue in this section

• Generating a Quadtree on the GPU with WGSL — build a linear Morton-code quadtree in parallel, with a radix sort and a range-finding pass instead of a recursive descent.
• On-GPU Viewport Culling for Vector Tiles — test feature bounds against the view extent and compact survivors into an indirect draw with an atomic counter.
• Point-Cloud LOD Compaction with Stream Compaction — thin a point cloud to a level-of-detail subset with a prefix-sum scan and a race-free scatter.

Related

• Geometry Filtering with WGSL Compute Shaders — the predicate-filtering companion whose compaction pattern these algorithms generalize.
• Spatial Aggregation in GPU Memory — the reduction-oriented passes that consume the indexes and survivor lists built here.
• Optimization Flags for Compute Dispatches — dispatch and pipeline tuning that compounds with the workgroup sizing chosen here.
• Workgroup Occupancy Optimization for Spatial Kernels — choosing the scan and scatter block sizes that keep the GPU saturated.
• Memory Alignment for Spatial Data Buffers — the stride rules that keep Morton codes and compaction offsets addressing the right elements.

Up: Spatial Compute Shaders & Geometry Pipelines