Octree Indexing Point Clouds with Morton Codes

This page builds an octree over a LiDAR point cloud using nothing but NumPy and 3D Morton (Z-order) codes — quantising coordinates in EPSG:32618 into a power-of-two cube, interleaving the bits of x, y and z into one 64-bit key, sorting the cloud by that key once, and then answering node membership, box queries and level-of-detail sampling with integer shifts and binary search.

Why you hit this

Octrees underpin almost every point cloud format built for streaming: COPC and EPT are octrees, Potree is an octree, and 3D Tiles point cloud tilesets are usually generated from one. Understanding how they are addressed is the difference between treating those formats as black boxes and being able to debug a node that never loads, a hierarchy that is unbalanced, or a pipeline that spends an hour building an index a sort could have built in seconds. The formats themselves are compared in LAZ vs COPC vs EPT for point cloud delivery; this page is the index underneath them.

Prerequisites

  • numpy>=1.24 and laspy[lazrs]>=2.5.
  • A point cloud in a projected metric CRS — EPSG:32618+5703 in the examples. An octree over longitude and latitude has cells that are not cubes, which breaks every spacing assumption below.
  • Familiarity with the 2D equivalent in computing quadkeys and tile bounds in Python: a Morton code is a quadkey with a third axis, written in binary instead of base four.

Step-by-Step

1. Quantise the cloud into a cube

python
import laspy
import numpy as np

las = laspy.read("harbour_block_07.laz")
xyz = np.column_stack([las.x, las.y, las.z])           # EPSG:32618+5703, metres

BITS = 21                                              # 21 bits per axis → 63-bit key
lo = xyz.min(axis=0)
edge = float((xyz.max(axis=0) - lo).max()) * 1.000001  # cube side: the largest extent
cell = edge / (1 << BITS)
q = np.floor((xyz - lo) / cell).astype(np.uint64)
q = np.minimum(q, (1 << BITS) - 1)
print(f"cube edge {edge:.2f} m, finest cell {cell * 1000:.3f} mm, {len(xyz):,} points")

An octree splits a cube, so the bounding box is expanded to the length of its longest side along all three axes. A 500 m × 500 m block that is only 60 m tall wastes most of its vertical subdivisions on empty air, which is normal — empty nodes cost nothing because they are never materialised. Twenty-one bits per axis gives a finest cell under a tenth of a millimetre on a 500 m cube, far finer than any scanner, and packs three axes into a single uint64.

2. Interleave the bits into Morton codes

python
def part1by2(v):
    v = v.astype(np.uint64) & np.uint64(0x1FFFFF)
    v = (v | (v << np.uint64(32))) & np.uint64(0x1F00000000FFFF)
    v = (v | (v << np.uint64(16))) & np.uint64(0x1F0000FF0000FF)
    v = (v | (v << np.uint64(8)))  & np.uint64(0x100F00F00F00F00F)
    v = (v | (v << np.uint64(4)))  & np.uint64(0x10C30C30C30C30C3)
    v = (v | (v << np.uint64(2)))  & np.uint64(0x1249249249249249)
    return v

def morton3(qx, qy, qz):
    return part1by2(qx) | (part1by2(qy) << np.uint64(1)) | (part1by2(qz) << np.uint64(2))

codes = morton3(q[:, 0], q[:, 1], q[:, 2])
order = np.argsort(codes, kind="stable")
codes = codes[order]
xyz = xyz[order]
print("first codes:", [f"{c:#018x}" for c in codes[:3]])

part1by2 spreads the 21 bits of one coordinate out so that two zero bits sit between each original bit; the magic masks do that in six constant-time steps instead of a 21-iteration loop. Shifting y’s spread bits by one and z’s by two before OR-ing gives the interleaved pattern …z₁y₁x₁z₀y₀x₀. Every np.uint64(...) wrapper is there for a reason: shifting a uint64 array by a Python int promotes to float64 in older NumPy versions and silently destroys the high bits.

Interleaving three coordinates into one Morton code The two lowest bits of x, y and z, shown as separate rows, are interleaved into a single row ordered z1 y1 x1 z0 y0 x0. Each group of three bits selects one of eight child octants, so the highest group chooses the child of the root and each lower group chooses a child one level deeper. x₁x₀ y₁y₀ z₁z₀ z₁y₁x₁ z₀y₀x₀ octant at level 1 octant at level 2 Dropping the last 3 × k bits of a code gives the node that contains the point, k levels up.
A Morton code is a path from the root: each triple of bits names one of eight children, most significant level first.

3. Address nodes by shifting

python
def node_key(codes, level, bits=BITS):
    """Octree node at `level` (0 = root) containing each code."""
    return codes >> np.uint64(3 * (bits - level))

for level in (4, 6, 8):
    keys, counts = np.unique(node_key(codes, level), return_counts=True)
    size = edge / (1 << level)
    print(f"level {level}: node edge {size:6.2f} m, {len(keys):6,} occupied nodes, "
          f"max {counts.max():,} pts, median {int(np.median(counts)):,} pts")

Because the cloud is sorted by code and a node key is a prefix of the code, every point in a node sits in one contiguous run of the sorted arrays. That is the property the rest of the index relies on: a node’s points are a slice, found by binary search, with no tree structure stored anywhere.

python
def node_slice(codes, key, level, bits=BITS):
    shift = np.uint64(3 * (bits - level))
    start_code = np.uint64(key) << shift
    end_code = (np.uint64(key) + np.uint64(1)) << shift
    return slice(np.searchsorted(codes, start_code, "left"), np.searchsorted(codes, end_code, "left"))

level = 6
keys, counts = np.unique(node_key(codes, level), return_counts=True)
busiest = keys[np.argmax(counts)]
s = node_slice(codes, busiest, level)
print(f"node {busiest} at level {level}: points {s.start:,}–{s.stop:,} ({s.stop - s.start:,})")
assert s.stop - s.start == counts.max()

Two searchsorted calls on a sorted array of 50 million codes take microseconds. The assertion is the cheapest correctness test there is: the slice length has to equal the count from np.unique, and if the codes were not sorted, or a shift was computed in the wrong dtype, it will not.

Z-order on a 4 × 4 slice and contiguous node runs A four by four grid of cells numbered in Z order, with the path drawn through them. The four cells of each two by two quadrant are consecutive in the order, so the quadrant is one contiguous run in the sorted list shown on the right. The same holds in three dimensions for each octree node. 0123 4567 891011 12131415 0–34–78–1112–15 sorted codes: one slice per quadrant Consecutive codes share a prefix, so every node is a range, never a scattered set.
Z-order visits each quadrant — and in 3D each octant — completely before moving on, which is why a sort is all the index needs.

5. Sample a level of detail per node

Streaming formats store a thinned subset of points at each coarse node and the rest deeper down. With a sorted Morton array, a simple, deterministic version is one point per finer grid cell per node.

python
def lod_sample(codes, level, per_axis_bits=4, bits=BITS):
    """Keep the first point in each of 8**per_axis_bits sub-cells of every node at `level`."""
    sub_level = level + per_axis_bits
    sub_keys = node_key(codes, sub_level, bits)
    first = np.concatenate([[True], sub_keys[1:] != sub_keys[:-1]])
    return np.flatnonzero(first)

for level in (2, 4, 6):
    idx = lod_sample(codes, level)
    spacing = edge / (1 << (level + 4))
    print(f"LOD for level {level}: {len(idx):,} points, nominal spacing {spacing:.2f} m")

Because the array is sorted, “the first point in each sub-cell” is just every position where the sub-cell key changes — one vectorised comparison. At level L the sample has a nominal spacing of the node edge divided by 16, which gives a coarse-to-fine sequence whose spacing halves at each level, the structure a renderer’s screen-space-error test expects. Production writers like PDAL’s COPC writer choose the retained point more carefully, but the addressing is the same.

Expected Output & Verification

text
cube edge 512.37 m, finest cell 0.244 mm, 48,302,117 points
level 4: node edge  32.02 m,    214 occupied nodes, max 1,873,002 pts, median 188,404 pts
level 6: node edge   8.01 m,  3,120 occupied nodes, max 198,441 pts, median 11,907 pts
level 8: node edge   2.00 m, 39,516 occupied nodes, max 18,220 pts, median 891 pts
node 190472 at level 6: points 21,004,318–21,202,759 (198,441)
LOD for level 2: 38,114 points, nominal spacing 8.01 m
LOD for level 4: 402,882 points, nominal spacing 2.00 m
LOD for level 6: 4,118,760 points, nominal spacing 0.50 m

Verify that decoding a code returns the quantised coordinates exactly, and that each node’s points lie inside its geometric bounds:

python
def compact1by2(v):
    v = v & np.uint64(0x1249249249249249)
    v = (v ^ (v >> np.uint64(2)))  & np.uint64(0x10C30C30C30C30C3)
    v = (v ^ (v >> np.uint64(4)))  & np.uint64(0x100F00F00F00F00F)
    v = (v ^ (v >> np.uint64(8)))  & np.uint64(0x1F0000FF0000FF)
    v = (v ^ (v >> np.uint64(16))) & np.uint64(0x1F00000000FFFF)
    v = (v ^ (v >> np.uint64(32))) & np.uint64(0x1FFFFF)
    return v

qs = q[order]
assert np.array_equal(compact1by2(codes), qs[:, 0])
assert np.array_equal(compact1by2(codes >> np.uint64(1)), qs[:, 1])
assert np.array_equal(compact1by2(codes >> np.uint64(2)), qs[:, 2])

k = node_key(np.array([codes[s.start]]), level)[0]
nq = np.array([compact1by2(np.uint64(k) >> np.uint64(i)) for i in range(3)], dtype=np.float64)
node_lo = lo + nq * (edge / (1 << level))
pts = xyz[s]
assert np.all(pts >= node_lo - 1e-6) and np.all(pts <= node_lo + edge / (1 << level) + 1e-6)
print("round trip and node bounds verified")
Points per occupied node by level on an urban block For octree levels four, six and eight, bars compare the median and maximum number of points per occupied node. The maximum is roughly ten to twenty times the median at every level, because dense facades and trees concentrate points, which is why production writers split nodes by point count rather than stopping at a fixed depth. level 4level 6level 8 blue: median points per node orange: maximum (log scale) A fixed-depth octree over a city is always unbalanced.
Point density follows surfaces, not space; a facade node holds twenty times the median, so node size limits belong in point counts.

Common Errors

Codes are not unique and neighbouring points get the same code. Quantisation was too coarse — BITS set low to save space, or the cube edge computed from one axis only. Duplicates are legitimate at the finest level for coincident returns, but more than a fraction of a percent means the grid is coarser than the data.

OverflowError or codes wrapping to small values. A shift was applied with a Python integer to a uint64 array, or the quantised values exceeded 21 bits because the maximum coordinate was not clamped. Keep every shift and mask as np.uint64 and clamp to (1 << BITS) - 1.

Node slices are empty for nodes that np.unique reports. The arrays were sorted by code but xyz was not reordered with the same permutation, or the codes array was re-sorted after the points were. Sort once, apply the same order to everything, and keep the assertion from step 4 in the pipeline.

Frequently Asked Questions

How does this relate to COPC and EPT keys?

Both address nodes by level and integer x, y, z within that level — D-X-Y-Z — which is the same information as a Morton node key, unpacked. Decoding a node key with compact1by2 at a given level yields exactly those X, Y, Z values.

Can Morton codes answer nearest-neighbour queries?

Approximately. Points close in code are close in space, but not every spatially close point is close in code, because Z-order jumps at node boundaries. Use the Morton sort to partition the data, then a KD-tree within and across adjacent nodes for exact neighbours.

Why not just use a KD-tree for everything?

A KD-tree is excellent in memory and has no natural serialisation into streamable, independently loadable chunks. An octree over Morton codes maps directly onto files and byte ranges, which is what delivery formats need.

Back to Spatial Indexing and Tiling Schemes for 3D Data.