Eight children at a time: a sparse octree over a point cloud

A voxel grid over a point cloud is the obvious spatial index and the obvious waste: a scan of a room at 2 cm resolution needs a grid a few hundred cells on a side, and almost all of those cells are empty air. An octree fixes the memory problem by subdividing only where there’s something to subdivide, but the textbook version pays for it on the way in, every insertion walks from the root, one comparison per level, for however many levels deep the tree goes. In the middle of my masters I wrote a C++ octree (LibTorch for its resizable typed arrays, pybind11 to call it from Python) that does neither: it finds a point’s node with a handful of bit operations, no walk, no comparisons, and it never allocates a node that nothing has visited yet. This post rebuilds that structure in the browser, over the Stanford Bunny, so you can drag a depth slider and watch millions of “addressable” cells collapse to the few thousand that actually hold a point.

With JavaScript on, this becomes a live WebGL octree explorer: drag the depth slider and watch the count of occupied leaves against the 8^depth addressable cells at that depth, toggle the leaves as wireframe cubes and the traversal order as a curve through their centres, click a cube to reveal its own children, or drop in your own small PLY or CSV point cloud. Without JavaScript, here is the shape of the argument in numbers: at depth 6 over the bundled, downsampled bunny, 6,498 of 262,144 addressable cells are occupied (2.5%); at depth 8, 7,269 of 16,777,216 (0.04%).

Peeling three bits at a time

Voxelise a point by dividing its position by the voxel size and flooring, same as any voxel grid. The trick is in how the octree finds that voxel’s node. At depth d, splitting space in half along each axis d times, a voxel’s integer coordinate is exactly d bits wide per axis. The bit at position level of that coordinate says which half of the parent cell the voxel falls in at that level, so reading the coordinate’s bits from the coarsest level down to the finest is the descent, there’s no separate step where you compute which child to visit; the coordinate already encodes the whole path from the root.

That’s what the actual descent loop does. Trimmed for the parts that matter (comments and variable names are the original’s):

// worldmap.cpp:124-172 (WorldMap::locateOnLevelCPU, function signature trimmed)
for (; outIdxPtr < endPtr; outIdxPtr++){
    uint8_t sfcOffset = 0;
    uint8_t chunkIdx = *chunkIdxPtr;
    for (auto dim = SPATIAL_DIMS-1; dim >= 0; dim--){
        // Floor(vox_idx/2**level)
        sfcOffset |= ((*voxIdxPtr >> level) & 0b1) << dim;
        voxIdxPtr++;
    }
    long idx = sfcOffset + ((*outIdxPtr)<<SPATIAL_DIMS);
    long nextIdx = octreeLevelIdx[chunkIdx][idx];
    if (nextIdx < 0) {
        levelMaxIdx[chunkIdx][0] += 1;
        nextIdx = levelMaxIdx[chunkIdx][0];
        octreeLevelIdx[chunkIdx][idx] = nextIdx;
        if (childParentIdx[chunkIdx] != nullptr){
            childParentIdx[chunkIdx][nextIdx] = idx;
        }
    }
    *outIdxPtr = nextIdx;
    chunkIdxPtr++;
}

sfcOffset is 3 bits, one per axis, built by pulling bit level out of each of x, y and z (SPATIAL_DIMS = 3, SPATIAL_COMB = 1<<SPATIAL_DIMS = 8, worldmap.h:8-9). idx then shifts the previous level’s index left by 3 and ORs in this level’s 3 bits. Do that once per level from the coarsest down to the finest and the final idx is the bits of x, y and z interleaved, x in the high position of every 3-bit group, which is a Morton (Z-order) code, computed with no encoding step at all, it falls out of reading the coordinate’s own bits in the right order. That’s also why nearby octree indices land near each other in memory: neighbouring cells share most of their high bits, so their Morton codes are close together too. The widget’s traversal-order toggle draws exactly this, a line through every occupied leaf’s centre in that order, and it visibly clusters into local runs rather than jumping around the volume.

The other half of the trick is nextIdx < 0: a level’s array starts full of -1, and a node is only ever created the first time some point’s descent actually reaches it (levelMaxIdx[chunkIdx][0] += 1, a running high-water mark, one per level per chunk). An 8-level tree over an empty volume allocates nothing; it grows exactly as many nodes as points touch, one bit-shift and one array write per level per point, no root-to-leaf comparison chain.

Reading a leaf back to a position

Descending gives you an index, not a coordinate, insertion never needed one. Getting a coordinate back out (which the widget needs, to draw a cube) means walking that same path in reverse, using the parent pointers locateOnLevelCPU wrote down on the way in:

// worldmap.cpp:94-118 (WorldMap::getVoxelData, trimmed)
for (auto lvl_it = level; lvl_it < m_tree_depth; lvl_it++){
    for (long idx_chk = 0; idx_chk < octree_idx_length; idx_chk++){
        for (auto dim = 0; dim < SPATIAL_DIMS; dim ++){
            out_acc[out_idx][dim] += ((*parent_it>>(SPATIAL_DIMS-1-dim))&0b1) << lvl_it;
        }
        if (lvl_it < m_tree_depth -1) {
            *parent_it = lvl_parent_idx[*parent_it>>SPATIAL_DIMS];
        }
    }
}

Same bit-peeling, run backwards: each level’s own 3 bits get OR’d into the output coordinate at the position that level owns, then the parent pointer moves up one level. My widget doesn’t reproduce this reverse walk, and that’s a deliberate simplification worth being upfront about, not a hidden shortcut. The lazy per-level arrays are a memory optimisation, not a different mapping: whichever way you get there, the voxel a point belongs to at depth d is always floor(position / voxel_size(d)). So the widget computes occupancy directly with a hash grouping over that formula (src/widgets/sparse-octree/octree.ts, occupiedCells), and only reproduces the order faithfully, via the same bit-interleaving mortonCode shown above, because that order is the actual subject of the post. I checked this equivalence against the repo’s own ground truth before trusting it: test/test_octree.py’s four single-point cases (test_simple_block, test_x_correct, test_y_correct, test_z_correct) each locate one point in a depth-1 octree and assert the exact occupied voxel coordinate; the widget runs the same four cases against its own occupiedCells at load time and warns to the console if they ever disagree.

Chunks: a hash lookup instead of one giant tree

A single octree over an unbounded scene has to pick a depth for the whole thing up front, which either wastes memory far from where the sensor’s been or runs out of precision once you’ve walked far enough. My octree sidesteps this by paying for it at a coarser level: space is partitioned into chunk_size-metre cubes, each its own independent, fixed-depth octree, keyed by chunk coordinate in a hash map. Finding the right chunk is one hash lookup instead of a tree walk, and a chunk that’s never been visited simply doesn’t exist in the map.

The wrinkle is what happens at a chunk’s edge. A batch of points (a frame of a depth image, say) can straddle a real chunk boundary, and re-hashing per point would be wasteful when almost every point in the batch lands in the same one or two chunks. So loadChunks loads 8 chunks at once, the 2x2x2 block around the batch’s own minimum corner, and each point is then routed to 1-of-8 by comparing its coordinate against one fixed midpoint per axis:

// worldmap.cpp:239-249 (WorldMap::getChunkIdx)
void WorldMap::getChunkIdx(const long* voxIdxPtr, uint8_t* outChunkIdxPtr, int64_t numPts) {
    const uint8_t* const endChnkPtr = &(outChunkIdxPtr[numPts]);
    for (; outChunkIdxPtr < endChnkPtr; outChunkIdxPtr++){
        uint8_t chunkIdxVal = 0;
        for (auto dim = 0; dim<SPATIAL_DIMS; dim++){
            chunkIdxVal |= (*voxIdxPtr >= m_loadedChunksMidpoint[dim]) * (FIRST_COORDINATE_MULTIPLIER >> dim);
            voxIdxPtr++;
        }
        *outChunkIdxPtr = chunkIdxVal;
    }
}

That’s the whole routing decision, three comparisons and three multiplies, reused for every point in the batch without a second hash lookup. locate (worldmap.cpp:331-356) guards the assumption this depends on directly: TORCH_CHECK(maxRange<=m_chunk_size, ...), the batch has to fit within one chunk’s width of its own minimum corner, or the 8 loaded chunks can’t possibly cover it.

The widget above has a small second canvas under the main view for exactly this: drag the square (the batch’s own minimum corner) and the circle (a query point, kept within one chunk of the square) to see the routing decision live, with a toggle to reveal or hide the actual chunk boundary underneath it.

What the widget actually shows

The main view loads the bundled bunny (or a small PLY/CSV you drop in) and rebuilds the octree’s occupancy from scratch on every depth change, drawing the occupied leaves as wireframe cubes over the points. The depth slider is the whole point of the post made draggable: at depth 1 every one of the 8 cells is occupied, at depth 6 it’s 6,498 of 262,144 (2.5%), and by depth 8 it’s 7,269 of 16,777,216, 0.04%, a dense grid at that depth would need one byte per cell just to say “empty” and still cost more memory than the actual point cloud. Click any occupied cube and its own up-to-8 children at the next depth light up, real occupancy, not an illustration, since the widget queries the same point data the parent cube did. The traversal-order toggle draws the space-filling curve described above through the occupied leaves’ centres, in the order the real descent would visit them.

Everything here is plain JavaScript, typed arrays and one hand-rolled WebGLRenderingContext, one vertex shader for points, one for lines, a handful of gl.POINTS/gl.LINES draw calls, no three.js. Bit-shifting and hashing a few thousand points at interactive rates is squarely the “plain JS is faster to ship” case rather than a WASM one; there’s no per-pixel image work or heavy linear algebra here to justify the round trip.

Where this goes next

An octree like this only answers “is this cell occupied”: a bit, nothing more. The natural next step, which I built into the same structure, is to keep a running mean and covariance of the points inside each occupied voxel instead of just a bit, updated one point at a time with no need to ever revisit old data. That’s the piece this one is missing, and it’s the subject of its own post.