Morton Code Generation for Spatial Sorting in WGSL

A Morton code — a Z-order curve index — interleaves the bits of two quantised coordinates into a single integer with one useful property: numerically close keys are spatially close. That single property is what lets a plain linear radix sort produce spatial locality, what makes a quadtree node a contiguous range of a sorted array, and what turns a range query into a pair of binary searches. Generating one is a branch-free ladder of shifts and masks costing a handful of ALU operations. This page is that ladder, the quantisation it depends on, the 64-bit variant, and the two mistakes that silently destroy the locality it exists to create. It is one stage of GPU sorting and prefix sums for spatial data.

Two 8-bit coordinates becoming one 16-bit key A single strip of sixteen bit positions showing the interleave. Reading from the most significant bit down, the positions alternate between the y coordinate and the x coordinate: y7, x7, y6, x6, y5, x5 and so on to y0, x0. Because the bits alternate, the top two bits of the key identify a quadrant of the extent, the top four identify a quadrant of that quadrant, and in general the top two times k bits are the quadtree path to depth k. BIT POSITION 0 2 4 6 8 10 12 14 16 y7 x7 y6 x6 y5 x5 y4 x4 y3 x3 y2 x2 y1 x1 y0 x0 16-bit key from 8-bit x, y the top 2k bits are exactly the quadtree path to depth k A 32-bit key takes 16 bits per axis; a 64-bit key takes 32 and doubles the sort passes.
The alternation is the whole design. It is what makes a numeric comparison between two keys equivalent to a comparison of their positions on a Z-order curve, and what makes a subtree a contiguous range.

Runnable reference implementation

The interleave works by spreading each coordinate’s bits apart — inserting a zero between every pair — and then shifting one of them left by one and combining. The spreading is the classic magic-number ladder, which is branch-free and constant-time.

wgsl
/// Spread the low 16 bits of v so each occupies an even position:
/// abcd efgh ijkl mnop -> 0a0b 0c0d 0e0f 0g0h 0i0j 0k0l 0m0n 0o0p
fn spread_bits_16(value : u32) -> u32 {
  var v = value & 0x0000FFFFu;
  v = (v | (v << 8u))  & 0x00FF00FFu;
  v = (v | (v << 4u))  & 0x0F0F0F0Fu;
  v = (v | (v << 2u))  & 0x33333333u;
  v = (v | (v << 1u))  & 0x55555555u;
  return v;
}

/// 32-bit Morton code from two 16-bit quantised coordinates.
/// y occupies the odd bit positions, x the even ones.
fn morton_2d(qx : u32, qy : u32) -> u32 {
  return (spread_bits_16(qy) << 1u) | spread_bits_16(qx);
}

The quantisation is the part that carries the spatial meaning, and it has to use the same extent for every point in the sort or the ordering is meaningless.

wgsl
struct Extent {
  min_xy : vec2<f32>,
  inv_span : vec2<f32>,     // 1 / (max - min), precomputed on the host
};
@group(0) @binding(0) var<uniform> extent : Extent;
@group(0) @binding(1) var<storage, read>       points : array<vec2<f32>>;
@group(0) @binding(2) var<storage, read_write> keys   : array<u32>;

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

  // Normalise into [0, 1] against the SHARED extent, then quantise to 16 bits.
  let n = clamp((points[i] - extent.min_xy) * extent.inv_span,
                vec2<f32>(0.0), vec2<f32>(1.0));
  let q = vec2<u32>(n * 65535.0);
  keys[i] = morton_2d(q.x, q.y);
}

The clamp matters. A point marginally outside the extent — a coordinate that arrived from a neighbouring tile, or a rounding artefact at a boundary — produces a quantised value above 65535, whose high bits corrupt the interleave and land the point somewhere arbitrary in the ordering. Clamping puts it on the edge of the grid, which is both correct and boring.

Parameter reference

Value Setting here Guidance
Bits per axis 16 Gives a 32-bit key and a four-pass radix sort. At a 40 km extent that is 0.6 m of resolution.
Key width u32 64-bit keys need eight sort passes; use them only when 16 bits per axis is genuinely too coarse.
Extent shared, per sort Every point in one sort must be quantised against the same extent, or the ordering means nothing.
Extent source the tile key Derive it from the tile, never from the data bounds — data-derived extents change when the data does.
Clamp required An out-of-extent point otherwise corrupts the interleave rather than landing at an edge.
The order a Z curve visits an 8×8 grid An eight by eight grid of cells with sixteen cells shaded to show the first quarter of the Z-order traversal. The curve fills the bottom-left four by four block completely before moving to the bottom-right block, then the top-left, then the top-right, and within each block it recurses in the same order. The shaded region is therefore one contiguous run of keys and one quadtree node, which is the property the sort exploits. Z-ORDER · FIRST QUADRANT SHADED a quadrant is a key range recursion is bit prefixes a range query is a few searches candidates still need a bbox test keys 0–15: one quadtree node later quadrants The Z curve jumps between quadrants, so a rectangle maps to several ranges, not one.
Every quadrant is a contiguous run of keys, at every depth. That is what turns a tree traversal into a range lookup and a pointer chase into a binary search.

Why the top bits are a quadtree path

The interleave has a structural consequence that is easy to miss and is the reason Morton codes are used rather than any other space-filling curve that is cheaper to compute.

Take the top two bits of the key. One came from the top bit of y and one from the top bit of x, so together they say which quadrant of the extent the point falls in: 00 is the bottom-left, 01 the bottom-right, 10 the top-left, 11 the top-right. The next two bits say which quadrant of that quadrant — and so on down.

The top 2k bits of a Morton key are therefore exactly the path from the root of a quadtree to the node containing that point at depth k. Which means that after a sort, every node at every depth is a contiguous range of the array, and finding it is a matter of comparing prefixes rather than following pointers. That is the whole of GPU quadtree construction: sort, then scan for prefix changes.

It also gives range queries a cheap first approximation. A rectangle in space maps to a set of key ranges, and while that set is conservative — the Z curve jumps, so a rectangle is not one contiguous range — a small number of ranges covers it, and a binary search finds each one. The candidates then need a real bounding-box test, because the key range is a superset, never an exact match.

Choosing how many bits per axis

The bit budget is the one real tuning decision, and it trades sort cost against spatial resolution.

Sixteen bits per axis gives a 32-bit key, a four-pass radix sort, and 65 536 cells along each axis of the extent. Over a 40-kilometre extent that is 0.6 metres per cell, which is finer than most vector tile sources carry and comfortably finer than a viewport-scale query needs. It is the right default.

Thirty-two bits per axis gives a 64-bit key, an eight-pass sort — twice the most expensive primitive in the pipeline — and resolution nobody has data for at map scale. It earns its place only in survey or point-cloud work where the extent is small and the source data genuinely resolves millimetres, and even then it is worth checking whether shrinking the extent would have been the cheaper answer.

Below sixteen bits, keys start colliding. Two points in the same cell get the same key, the sort puts them in arbitrary relative order, and any structure built on the ordering treats them as interchangeable. For clustering and level-of-detail work that is often fine — points in the same cell are exactly the ones the algorithm wants to merge — but for a picking index it is not, and the failure looks like a hit test returning the wrong feature from a pair.

The useful reframing is that bits per axis is a resolution decision expressed in a strange unit. Pick the ground resolution the application needs, divide the extent by it, take the base-two logarithm, and round up.

Failure modes

  • The sort is correct but neighbouring entries are far apart. x and y were spread and then combined without the shift, so both occupy the same bit positions and the key is a bitwise OR rather than an interleave. Detection: keys collide far more often than they should.
  • Points from different tiles interleave wrongly. Two extents were used in one sort. The extent must be shared across everything being sorted together.
  • A few points land at arbitrary positions. No clamp, and a coordinate outside the extent overflowed 16 bits into the neighbouring axis’s bit positions.
  • The locality is right but the quadtree nodes overlap. The boundary scan compared the wrong number of top bits: it must be 2 × depth, not depth.
  • Everything degrades as the dataset grows. 16 bits per axis is 65 536 cells; past a few million points in a small extent, many points share a key and the ordering within a cell is arbitrary. Fix: more bits, and accept the extra sort passes.
Morton key mistakes and what each destroys A table of four Morton key mistakes with their symptom and fix. Omitting the shift when combining the two spread values makes the key a bitwise OR rather than an interleave, so keys collide and locality is lost. Using two different quantisation extents in one sort makes the ordering meaningless across the boundary. Omitting the clamp lets an out-of-extent coordinate overflow into the neighbouring axis bits. Comparing depth bits rather than two times depth bits when scanning for node boundaries produces overlapping nodes. KEY MISTAKE · SYMPTOM · FIX Symptom Fix no shift on combine keys collide shift y left by 1 two extents in one sort order is meaningless share the extent no clamp stray points scattered clamp to [0,1] wrong prefix width nodes overlap use 2 × depth bits A quick check: sample 100 sorted pairs and confirm the median distance between neighbours is small.
None of these raises an error — every one produces a key that sorts fine and means the wrong thing. Checking that sorted neighbours are spatially close is the assertion worth writing.

Backend / Python interop note

Keys can be generated server-side, and there is a real argument for doing so: a tile whose points arrive already sorted by Morton code needs no sort on the client at all, and the ordering compresses better than the original because neighbouring rows now share high bits.

python
import numpy as np

def spread_bits_16(v: np.ndarray) -> np.ndarray:
    v = v.astype(np.uint32) & 0x0000FFFF
    v = (v | (v << 8)) & 0x00FF00FF
    v = (v | (v << 4)) & 0x0F0F0F0F
    v = (v | (v << 2)) & 0x33333333
    v = (v | (v << 1)) & 0x55555555
    return v

def morton_2d(qx: np.ndarray, qy: np.ndarray) -> np.ndarray:
    return (spread_bits_16(qy) << 1) | spread_bits_16(qx)

Two cautions. The extent used for quantisation has to travel with the data — in the record batch metadata, alongside the origin from relative-to-eye encoding — because a client that re-derives it from the data it received will get a different extent from the one the server used and every key will be wrong. And numpy’s default integer type is platform-dependent; the explicit uint32 casts above are load-bearing, and omitting them produces 64-bit intermediates whose shifts overflow differently.