A captured 3D Gaussian Splatting scene is a few million splats; a lidar scan or a photogrammetry point cloud can be a billion points. In the order they come out of the trainer or the scanner they are scattered: neighbours in space sit megabytes apart in memory. Before you can compress them, stream them, build levels of detail or find neighbours, you put them in spatial order. The workhorse for that is the Morton code (Z-order curve, Morton 1966): quantise each position to an integer grid, interleave the bits of x, y and z, and sort. The sorted array is an octree — every cell at every level is one contiguous run — so a "linear octree" falls out almost for free.
1. Quantise, then interleave
Pick a bounding box and a resolution: with 21 bits per axis the grid has
2²¹ ≈ 2.1 million cells per axis, and three axes fill 63 bits of a uint64_t
(10 bits per axis fit a uint32_t). Quantise each coordinate:
q = clamp(⌊(v − lo) / (hi − lo) · 2²¹⌋, 0, 2²¹ − 1)
The clamp matters: the bounding box of a point cloud is defined by points on
its faces, and for v == hi the formula gives 2²¹ — one past the last cell.
The Morton code takes one bit from each axis in turn, from the least significant
up: bit i of x goes to bit 3i, of y to 3i + 1, of z to 3i + 2.
Two facts make the code useful:
- Prefixes are cells. The top
3Lbits of a point's code are the Morton code of its cell at level L (a grid of 2ᴸ per axis):code >> 3·(21 − L). Dropping the lowest triple halves every coordinate — that's the parent cell. - Sorting by code is sorting by cell, at every level at once. Points in the same cell share a prefix, so they are contiguous after the sort.
2. Spreading bits with magic numbers
The obvious encoder loops over 21 bits per axis. The standard one (see Fabian
Giesen's "Decoding Morton codes") moves all bits at once in ⌈log₂ 21⌉ = 5
steps: each x = (x | x << s) & mask splits the blocks of bits in two and moves
the upper halves up by s; the mask keeps exactly where the blocks should be.
#include <cstdint>
#include <print>
std::uint64_t part1by2(std::uint32_t v) {
std::uint64_t x = v & 0x1fffff; // 64-bit from the start!
x = (x | x << 32) & 0x1f00000000ffff;
x = (x | x << 16) & 0x1f0000ff0000ff;
x = (x | x << 8) & 0x100f00f00f00f00f;
x = (x | x << 4) & 0x10c30c30c30c30c3;
x = (x | x << 2) & 0x1249249249249249;
return x;
}
std::uint32_t compact1by2(std::uint64_t x) { // the same steps, backwards
x &= 0x1249249249249249;
x = (x ^ (x >> 2)) & 0x10c30c30c30c30c3;
x = (x ^ (x >> 4)) & 0x100f00f00f00f00f;
x = (x ^ (x >> 8)) & 0x1f0000ff0000ff;
x = (x ^ (x >> 16)) & 0x1f00000000ffff;
x = (x ^ (x >> 32)) & 0x1fffff;
return static_cast<std::uint32_t>(x);
}
std::uint64_t morton3d(std::uint32_t x, std::uint32_t y, std::uint32_t z) {
return part1by2(x) | part1by2(y) << 1 | part1by2(z) << 2;
}
int main() {
std::println("{} {} {} {}", morton3d(1, 0, 0), morton3d(0, 1, 0), morton3d(0, 0, 1), morton3d(5, 3, 6));
std::println("{:#x}", morton3d(0x1fffff, 0x1fffff, 0x1fffff));
std::uint64_t c = morton3d(123456, 7, 2000000);
std::println("{} {} {}", compact1by2(c), compact1by2(c >> 1), compact1by2(c >> 2));
}Pitfall: the shift width. std::uint32_t x; … x << 32 is undefined
behaviour — the shift count must be less than the width of the promoted left
operand, and x stays 32-bit even if the result is stored in a uint64_t. On
x86 the hardware masks the count to 5 bits, so the bug silently produces wrong
codes; UBSan reports "shift exponent 32 is too large for 32-bit type".
Widen first. The same applies to every loop-based encoder: (v >> i & 1u) << 3*i
is a 32-bit expression.
The BMI2 aside. x86's pdep ("parallel bits deposit") does part1by2 in one
instruction, and pext is compact1by2:
#include <immintrin.h>
__attribute__((target("bmi2")))
std::uint64_t morton3d_bmi2(std::uint32_t x, std::uint32_t y, std::uint32_t z) {
return _pdep_u64(x, 0x1249249249249249) | _pdep_u64(y, 0x2492492492492492) |
_pdep_u64(z, 0x4924924924924924);
}Fast on Intel since Haswell and AMD since Zen 3 (microcoded and slow on Zen 1/2), not portable to ARM — so libraries keep the magic-number version as the default.
3. Sorting: an LSD radix sort
Codes are 64-bit integers, so a radix sort beats comparison sorts: 8 passes, each a stable counting sort on one byte, least significant first. Each pass is O(n) with sequential reads and writes, and stability is essential — it is what keeps the order established by earlier (less significant) passes, and it keeps duplicate points in a deterministic order. Sort an index permutation (or (code, index) pairs) rather than the points themselves; then gather attributes once.
#include <array>
#include <cmath>
#include <cstdint>
#include <print>
#include <vector>
std::uint64_t part1by2(std::uint32_t v) {
std::uint64_t x = v & 0x1fffff;
x = (x | x << 32) & 0x1f00000000ffff;
x = (x | x << 16) & 0x1f0000ff0000ff;
x = (x | x << 8) & 0x100f00f00f00f00f;
x = (x | x << 4) & 0x10c30c30c30c30c3;
x = (x | x << 2) & 0x1249249249249249;
return x;
}
std::vector<std::uint32_t> radix_sort(const std::vector<std::uint64_t>& keys) {
std::vector<std::uint32_t> idx(keys.size()), tmp(keys.size());
for (std::uint32_t i = 0; i < idx.size(); ++i) idx[i] = i;
for (int shift = 0; shift < 64; shift += 8) {
std::array<std::size_t, 257> count{};
for (auto i : idx) ++count[((keys[i] >> shift) & 0xff) + 1];
for (int b = 0; b < 256; ++b) count[b + 1] += count[b];
for (auto i : idx) tmp[count[(keys[i] >> shift) & 0xff]++] = i; // stable scatter
idx.swap(tmp);
}
return idx;
}
struct P { float x, y, z; };
int main() {
// 20000 points on a sphere (like a splat scene's surface), in random order.
std::uint32_t s = 12345;
auto rnd = [&] { s ^= s << 13; s ^= s >> 17; s ^= s << 5; return (s >> 8) * 0x1p-24f; };
std::vector<P> pts;
for (int i = 0; i < 20000; ++i) {
float z = 2 * rnd() - 1, phi = 6.2831853f * rnd(), r = std::sqrt(1 - z * z);
pts.push_back({r * std::cos(phi), r * std::sin(phi), z});
}
std::vector<std::uint64_t> keys;
for (auto p : pts) {
auto q = [](float v) { return std::uint32_t(std::fmin((v + 1) / 2 * 2097152.0, 2097151.0)); };
keys.push_back(part1by2(q(p.x)) | part1by2(q(p.y)) << 1 | part1by2(q(p.z)) << 2);
}
auto order = radix_sort(keys);
auto gap = [&](auto at) {
double sum = 0;
for (std::size_t k = 1; k < pts.size(); ++k) {
P a = pts[at(k - 1)], b = pts[at(k)];
sum += std::sqrt((a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y) + (a.z - b.z) * (a.z - b.z));
}
return sum / (pts.size() - 1);
};
std::println("mean distance between consecutive points: input order {:.4f}, Morton order {:.4f}",
gap([](std::size_t k) { return k; }), gap([&](std::size_t k) { return order[k]; }));
}The Morton-ordered stream steps a tiny fraction of the distance of the input order between consecutive elements — that is the locality caches, delta coders and chunked streaming feed on.
4. The linear octree
With codes sorted, an octree node is just (prefix, level, [begin, end)):
the root covers [0, n); a node's children are the sub-runs where the next
triple, (code >> 3·(21 − L − 1)) & 7, takes each value. Store the nodes in a
flat array, children of a node contiguous and in octant order, plus an 8-bit
child mask saying which octants exist — child k lives at
first_child + popcount(mask & ((1 << k) − 1)). Split only nodes with more
than leaf_size points (and stop at level 21: identical points can never be
separated).
#include <algorithm>
#include <bit>
#include <cstdint>
#include <print>
#include <vector>
struct Node { std::uint64_t prefix; int level; std::uint32_t begin, end, first_child; std::uint8_t mask; };
int main() {
// Already sorted 6-bit codes (a 2-level octree: 4×4×4 grid).
std::vector<std::uint64_t> codes{0b000'001, 0b000'011, 0b000'011, 0b000'110, 0b001'000, 0b001'101,
0b110'010, 0b110'111};
const int depth = 2;
const std::uint32_t leaf_size = 2;
std::vector<Node> nodes{{0, 0, 0, std::uint32_t(codes.size()), 0, 0}};
for (std::size_t i = 0; i < nodes.size(); ++i) { // breadth-first
Node n = nodes[i];
if (n.end - n.begin <= leaf_size || n.level == depth) continue;
int shift = 3 * (depth - (n.level + 1));
std::uint32_t first = std::uint32_t(nodes.size());
std::uint8_t mask = 0;
for (std::uint32_t b = n.begin; b < n.end;) {
std::uint64_t p = codes[b] >> shift;
std::uint32_t e = b;
while (e < n.end && codes[e] >> shift == p) ++e;
mask |= std::uint8_t(1u << (p & 7));
nodes.push_back({p, n.level + 1, b, e, 0, 0});
b = e;
}
nodes[i].first_child = first;
nodes[i].mask = mask;
}
for (std::size_t i = 0; i < nodes.size(); ++i) {
const Node& n = nodes[i];
std::println("node {}: level {} prefix {:0{}b} points [{}, {}) mask {:08b} children {}", i, n.level,
n.prefix, std::max(1, 3 * n.level), n.begin, n.end, n.mask, std::popcount(n.mask));
}
}Karras (2012) goes one step further for GPUs: in a sorted array of n codes, the
internal nodes of a binary radix tree can be built in parallel, one thread per
node, by finding the longest common prefix of neighbouring codes
(countl_zero(a ^ b)). LBVH builders for ray tracing and many GPU octree
builders are this algorithm.
5. Queries
Points in a cell. A point is in level-L cell (cx, cy, cz) iff
code >> 3·(21 − L) == morton3d(cx, cy, cz). That projection is monotone, so
in the sorted array the cell is one run: two binary searches. The classic bug
is comparing a level-L code with full-depth codes — always shift to the same
depth before comparing.
Radius (and k-NN) queries. Walk from the root; a node's cell spans
lo + c·size to lo + (c + 1)·size per axis, with c from decoding its prefix
and size = extent / 2ᴸ. Skip nodes whose box is farther than r from the query
(clamp the query point into the box, measure the distance), and test the points
of the leaves you reach. For k nearest neighbours, keep a max-heap of the best k
and shrink r to the current k-th distance; visit the nearer children first.
Uniform grid / spatial hashing. When the radius is fixed, a hash grid is
simpler: cell size r, key = hash of the integer cell, and a query reads its own
cell plus the 26 around it. Teschner et al. (2003) introduced the classic
(x·73856093) ^ (y·19349663) ^ (z·83492791) mod T hash for this, and
Instant-NGP's multiresolution hash encoding uses the same family of hashes for
its feature grids.
3D Gaussian Splatting initialises each Gaussian's scale from the mean distance
to its 3 nearest neighbours among the SfM points (the reference code calls a
CUDA distCUDA2 built on exactly this: Morton-sort the points, then search
nearby boxes):
#include <cmath>
#include <cstdint>
#include <print>
#include <unordered_map>
#include <vector>
struct P { float x, y, z; };
int main() {
std::uint32_t s = 7;
auto rnd = [&] { s ^= s << 13; s ^= s >> 17; s ^= s << 5; return (s >> 8) * 0x1p-24f; };
std::vector<P> pts(3000);
for (auto& p : pts) p = {rnd(), rnd(), rnd()};
const float cell = 0.15f; // comfortably above the typical 3-NN distance (~0.05 here)
auto key = [&](int x, int y, int z) {
return (std::uint64_t(std::uint32_t(x)) * 73856093u) ^ (std::uint64_t(std::uint32_t(y)) * 19349663u) ^
(std::uint64_t(std::uint32_t(z)) * 83492791u);
};
std::unordered_multimap<std::uint64_t, std::uint32_t> grid;
auto cell_of = [&](float v) { return int(std::floor(v / cell)); };
for (std::uint32_t i = 0; i < pts.size(); ++i)
grid.emplace(key(cell_of(pts[i].x), cell_of(pts[i].y), cell_of(pts[i].z)), i);
double mean_scale = 0;
int mismatches = 0;
for (std::uint32_t i = 0; i < pts.size(); ++i) {
float best[3] = {1e30f, 1e30f, 1e30f};
auto offer = [&](float d2) {
if (d2 < best[2]) { best[2] = d2; if (best[2] < best[1]) std::swap(best[1], best[2]);
if (best[1] < best[0]) std::swap(best[0], best[1]); }
};
int cx = cell_of(pts[i].x), cy = cell_of(pts[i].y), cz = cell_of(pts[i].z);
for (int dx = -1; dx <= 1; ++dx) for (int dy = -1; dy <= 1; ++dy) for (int dz = -1; dz <= 1; ++dz) {
auto [b, e] = grid.equal_range(key(cx + dx, cy + dy, cz + dz));
for (auto it = b; it != e; ++it) {
const P& q = pts[it->second];
if (it->second == i || cell_of(q.x) != cx + dx || cell_of(q.y) != cy + dy || cell_of(q.z) != cz + dz)
continue; // itself, or a hash collision
offer((q.x - pts[i].x) * (q.x - pts[i].x) + (q.y - pts[i].y) * (q.y - pts[i].y) +
(q.z - pts[i].z) * (q.z - pts[i].z));
}
}
if (best[2] > cell * cell) ++mismatches; // 3rd neighbour may lie outside the 27 cells
mean_scale += (std::sqrt(best[0]) + std::sqrt(best[1]) + std::sqrt(best[2])) / 3;
}
std::println("mean 3-NN distance {:.4f}; queries needing a wider search: {}", mean_scale / pts.size(), mismatches);
}The check at the end is the correctness condition (try cell = 0.1f): a grid search is exact only
for neighbours within one cell size; a k-NN search over a grid must widen the
block (or fall back to the tree) when the k-th distance exceeds it.
6. Level of detail and sparse voxel octrees
Because every cell is a run, an LOD pyramid is one sweep per level: for each distinct prefix at level L, merge its points — average positions and colours for a point cloud; for splats, merge Gaussians by moment matching (weights add, means average, covariances add the spread of the means). Streaming renderers (Potree for point clouds, hierarchical 3DGS by Kerbl et al. 2024) load coarse cells first and refine where the screen-space error is large.
A sparse voxel octree (SVO, Laine & Karras 2010) stores only occupied
nodes, each as a child mask plus a pointer to its first child — exactly the
(first_child, child_mask) pair of the build task — so empty space costs
nothing. Ray marching an SVO skips whole empty octants; NeRF accelerators
(PlenOctrees, Yu et al. 2021) bake radiance fields into such trees for
real-time rendering.
In Python (NumPy)
q = np.clip(((p - lo) / (hi - lo) * 2**21).astype(np.int64), 0, 2**21 - 1)
def part1by2(x): # vectorised over an int64 array
x &= 0x1fffff
for s, m in [(32, 0x1f00000000ffff), (16, 0x1f0000ff0000ff), (8, 0x100f00f00f00f00f),
(4, 0x10c30c30c30c30c3), (2, 0x1249249249249249)]:
x = (x | (x << s)) & m
return x
codes = part1by2(q[:, 0]) | part1by2(q[:, 1]) << 1 | part1by2(q[:, 2]) << 2
order = np.argsort(codes, kind="stable") # radix sort for integer keysIn an interview
- "How would you order a billion points for streaming?" — quantise, Morton code, radix sort (out of core: bucket by the top bits first), chunk by prefix; mention Hilbert for better locality.
- "Why 21 bits?" — 3 × 21 = 63 fits 64.
- "Find all neighbours within r" — hash grid with cell r and 27 cells, or an octree/k-d tree for varying radii; know what each costs.
Recap
- Quantise with a clamp; interleave x/y/z into bits 0/1/2 of each triple.
- part1by2 / compact1by2: five shift-mask steps in 64-bit arithmetic (or
pdep/pext). - Sort with a stable LSD radix sort; every octree cell becomes a contiguous run.
- Linear octree = (prefix, level, range, first_child, child_mask); parents are common prefixes.
- Compare codes only at the same depth:
code >> 3·(21 − L). - Radius / k-NN queries prune cells by box distance; hash grids are the fixed-radius shortcut.
- Hilbert curves keep better locality; SVOs and LOD pyramids reuse the same prefix structure.