Octree

How an octree works

Imagine a cube of empty 3D space and you want to drop in fifty floating points. You could store them in a flat list — but answering "what's near this query point?" costs O(n). The octree, introduced by Donald Meagher in 1980, partitions space recursively: every time a cube contains more than a small capacity of points, it splits into eight equal child cubes, and the points get reassigned.

The split is the same one a child draws when asked to cut a cake into eight pieces with three slices: one cut perpendicular to each axis through the center. Each child is named by its sign relative to the parent center: +X+Y+Z, +X+Y-Z, +X-Y+Z, +X-Y-Z, -X+Y+Z, etc. Each child can split further independently. Empty regions stop subdividing immediately, so dense areas get fine-grained cells and empty space stays a single cube.

Two operations dominate:

1. Insert(point). Descend from the root, picking the child whose region contains the point. If a leaf exceeds capacity, subdivide and redistribute. O(log n) on balanced data.
2. Range query (box). Recursively visit any cell whose region intersects the query box. Skip everything else. O(log n + k) where k is the number of points returned.

When to use an octree

• 3D collision broad-phase. Insert each entity into the leaf containing its center; per-frame collision checks only run between entities sharing a leaf or adjacent leaves. Cuts an O(n²) physics tick to roughly O(n log n).
• LIDAR and point clouds. A city scan can have a billion points; brute-force nearest-neighbor is hopeless. Octrees built once at load time give k-NN, range, and surface queries in milliseconds.
• Voxel terrain (Minecraft, voxel engines). Empty space — most of a Minecraft world is air — collapses to a single empty node. Sparse voxel octrees compress voxel data by 10³–10⁴×.
• Ray-tracing acceleration. A ray walks the octree, testing only triangles or voxels in cells it actually passes through. NVIDIA's RTX hardware traces rays through BVHs, but voxel global-illumination (VXGI, sparse voxel octree GI) still uses octrees.
• Surface reconstruction. The marching cubes and dual contouring algorithms typically run on octree-organized voxel data so that resolution scales with surface detail.
• 3D adaptive numerical meshes. Finite-element and fluid simulations refine the octree where the solution changes quickly; coarse where the field is smooth.

Octree vs BVH vs uniform grid vs k-d tree

	Octree	BVH (3D)	Uniform grid	k-d tree (3D)
Partitions	Space (regular)	Objects (irregular)	Space (fixed)	Space (data-driven)
Best for	Voxels, sparse, static	Triangles, dynamic	Uniform density, real-time	Static points
Insert	O(log n)	O(log n) + refit	O(1)	O(log n) at build
3D range query	O(log n + k)	O(log n + k)	O(cells × density)	O(n^(2/3) + k)
Dynamic updates	Re-insert + cascading splits	Refit bounding boxes	O(1) bucket swap	Rebuild required
Memory (sparse)	Excellent — empty = 1 node	O(n)	O(W·H·D)	O(n)
Memory (dense)	O(n)	O(n)	O(volume)	O(n)
Used by	Minecraft, VXGI, PCL	Embree, Unreal, OptiX	NPC AI, FX particles	scikit-learn kNN

Rule of thumb: octree for static voxel-aligned worlds and sparse 3D point data. BVH for dynamic mesh-based scenes (every modern ray tracer). Uniform grid for dense uniform fields. k-d tree for static point queries when you don't need axis-aligned cells.

Pseudo-code

// Octree with point payload, recursive insert, range query.

class OctreeNode:
    bounds   // 3D AABB { min, max }
    capacity // e.g. 8
    depth, maxDepth
    points   // null in internal nodes
    children // [8] or null

insert(node, p):
    if not contains(node.bounds, p): return false
    if node.children:
        for child in node.children:
            if insert(child, p): return true
        return false
    node.points.append(p)
    if len(node.points) > capacity and node.depth < maxDepth:
        subdivide(node)
        for q in node.points: insert(node, q)
        node.points = null
    return true

subdivide(node):
    cx = (node.bounds.min + node.bounds.max) / 2
    node.children = [octant 0..7]  // each is the AABB of one quadrant

rangeQuery(node, box, found):
    if not intersects(node.bounds, box): return
    if node.children:
        for child in node.children: rangeQuery(child, box, found)
    else:
        for p in node.points:
            if contains(box, p): found.append(p)

JavaScript implementation

class Octree {
  constructor(bounds, capacity = 8, depth = 0, maxDepth = 8) {
    // bounds = { min: {x,y,z}, max: {x,y,z} }
    this.bounds = bounds;
    this.capacity = capacity;
    this.depth = depth;
    this.maxDepth = maxDepth;
    this.points = [];
    this.children = null;
  }

  contains(p) {
    const { min, max } = this.bounds;
    return p.x >= min.x && p.x < max.x &&
           p.y >= min.y && p.y < max.y &&
           p.z >= min.z && p.z < max.z;
  }

  intersects(box) {
    const a = this.bounds, b = box;
    return !(b.min.x > a.max.x || b.max.x < a.min.x ||
             b.min.y > a.max.y || b.max.y < a.min.y ||
             b.min.z > a.max.z || b.max.z < a.min.z);
  }

  subdivide() {
    const { min, max } = this.bounds;
    const cx = (min.x + max.x) / 2,
          cy = (min.y + max.y) / 2,
          cz = (min.z + max.z) / 2;
    const oct = (xlo, ylo, zlo, xhi, yhi, zhi) => new Octree(
      { min: { x: xlo, y: ylo, z: zlo }, max: { x: xhi, y: yhi, z: zhi } },
      this.capacity, this.depth + 1, this.maxDepth
    );
    this.children = [
      oct(min.x, min.y, min.z, cx,    cy,    cz   ),  // ---
      oct(cx,    min.y, min.z, max.x, cy,    cz   ),  // +--
      oct(min.x, cy,    min.z, cx,    max.y, cz   ),  // -+-
      oct(cx,    cy,    min.z, max.x, max.y, cz   ),  // ++-
      oct(min.x, min.y, cz,    cx,    cy,    max.z),  // --+
      oct(cx,    min.y, cz,    max.x, cy,    max.z),  // +-+
      oct(min.x, cy,    cz,    cx,    max.y, max.z),  // -++
      oct(cx,    cy,    cz,    max.x, max.y, max.z),  // +++
    ];
    for (const p of this.points) {
      for (const c of this.children) if (c.contains(p)) { c.insert(p); break; }
    }
    this.points = null;
  }

  insert(p) {
    if (!this.contains(p)) return false;
    if (this.children) {
      for (const c of this.children) if (c.insert(p)) return true;
      return false;
    }
    if (this.points.length < this.capacity || this.depth >= this.maxDepth) {
      this.points.push(p);
      return true;
    }
    this.subdivide();
    return this.insert(p);
  }

  rangeQuery(box, found = []) {
    if (!this.intersects(box)) return found;
    if (this.children) {
      for (const c of this.children) c.rangeQuery(box, found);
    } else {
      for (const p of this.points) {
        if (p.x >= box.min.x && p.x < box.max.x &&
            p.y >= box.min.y && p.y < box.max.y &&
            p.z >= box.min.z && p.z < box.max.z) found.push(p);
      }
    }
    return found;
  }
}

// Usage
const oct = new Octree({ min: { x: 0, y: 0, z: 0 }, max: { x: 100, y: 100, z: 100 } });
for (let i = 0; i < 10000; i++) oct.insert({ x: Math.random()*100, y: Math.random()*100, z: Math.random()*100 });
const near = oct.rangeQuery({ min: { x: 40, y: 40, z: 40 }, max: { x: 60, y: 60, z: 60 } });

Python implementation

from dataclasses import dataclass, field
from typing import List, Optional, Tuple

@dataclass
class AABB:
    min: Tuple[float, float, float]
    max: Tuple[float, float, float]

    def contains(self, p):
        return (self.min[0] <= p[0] < self.max[0] and
                self.min[1] <= p[1] < self.max[1] and
                self.min[2] <= p[2] < self.max[2])

    def intersects(self, other):
        return not (other.min[0] > self.max[0] or other.max[0] < self.min[0] or
                    other.min[1] > self.max[1] or other.max[1] < self.min[1] or
                    other.min[2] > self.max[2] or other.max[2] < self.min[2])


class Octree:
    def __init__(self, bounds: AABB, capacity=8, depth=0, max_depth=8):
        self.bounds = bounds
        self.capacity = capacity
        self.depth = depth
        self.max_depth = max_depth
        self.points = []
        self.children: Optional[List['Octree']] = None

    def subdivide(self):
        mn, mx = self.bounds.min, self.bounds.max
        c = ((mn[0]+mx[0])/2, (mn[1]+mx[1])/2, (mn[2]+mx[2])/2)
        boxes = []
        for ix in range(2):
            for iy in range(2):
                for iz in range(2):
                    lo = (mn[0] if ix==0 else c[0],
                          mn[1] if iy==0 else c[1],
                          mn[2] if iz==0 else c[2])
                    hi = (c[0] if ix==0 else mx[0],
                          c[1] if iy==0 else mx[1],
                          c[2] if iz==0 else mx[2])
                    boxes.append(AABB(lo, hi))
        self.children = [Octree(b, self.capacity, self.depth + 1, self.max_depth) for b in boxes]
        for p in self.points:
            for ch in self.children:
                if ch.bounds.contains(p):
                    ch.insert(p); break
        self.points = None

    def insert(self, p):
        if not self.bounds.contains(p):
            return False
        if self.children is not None:
            return any(ch.insert(p) for ch in self.children)
        if len(self.points) < self.capacity or self.depth >= self.max_depth:
            self.points.append(p); return True
        self.subdivide()
        return self.insert(p)

    def range_query(self, box: AABB, out=None):
        if out is None:
            out = []
        if not self.bounds.intersects(box):
            return out
        if self.children is not None:
            for ch in self.children:
                ch.range_query(box, out)
        else:
            for p in self.points:
                if box.contains(p):
                    out.append(p)
        return out

Production-grade Python octree work uses the PCL Python bindings or open3d.geometry.Octree — both written in C++ with vectorized queries that beat any pure-Python implementation by 100× on dense point clouds.

Variants

• Point-region (PR) octree. Points stored in leaves with capacity threshold. The everyday flavor described above.
• Region octree. Cells store occupancy or color. Used for 3D image compression and voxel rendering. Empty regions collapse to single nodes.
• Sparse voxel octree (SVO). Region octree with explicit pointer-only storage and no node for empty space. Compresses a 4096³ voxel volume from 64 GB to a few hundred MB. Used in NVIDIA VXGI and id Tech Megavoxel research.
• Loose octree. Each child's region extends beyond its parent's strict octant by a slack factor. Eliminates the "object straddles a boundary" problem — every entity lives in exactly one cell.
• Linear octree (Morton-encoded). Encode each leaf by a Z-order key and store the leaves in a sorted array. Disk-friendly; the basis of LIDAR file formats like LAS and Entwine.
• Dual octree. Each cell stores both occupancy and the dual structure linking adjacent cells. Used by dual contouring and isosurface algorithms.

Performance — costed claims

• Sparse compression. A typical Minecraft world is 99% air. Encoding it as a sparse octree gives ~10³–10⁴× memory reduction versus a flat 3D array — the difference between 1 GB and 100 KB per chunk.
• Range query. With capacity 8 and depth 8 (16M leaves theoretical, far fewer in practice), a typical k-NN query over a million points takes ~10–100 microseconds in C++ (PCL or Open3D benchmarks).
• Construction. Bulk-building a Morton-ordered octree from a sorted point list is O(n) once the points are sorted (O(n log n) total). Open3D builds a 10M-point octree in under one second on modern desktop hardware.
• Ray traversal. A ray walks an octree in O(log n) hops per traversal; modern CPUs handle 50–200 million ray-cell tests per second per core for VXGI-style sparse octrees.
• Memory per node. Roughly 80–120 bytes (8 child pointers + bounds + metadata). At 100k nodes, that's 10 MB of node overhead.

Common bugs and edge cases

• Boundary tiebreaks. A point on the exact center of a cell can plausibly land in any of the eight children. Pick one rule (< on low edges, <= on high) and use it everywhere; otherwise points silently disappear.
• Pathological clustering. A thousand coincident points force infinite subdivision; recursion only stops at maxDepth, after which the leaf becomes a long list. Always cap depth.
• Floating-point drift. Subdividing a cube repeatedly can produce children whose computed bounds don't quite tile the parent. Use double precision or recompute bounds from the parent's min/max each time rather than from the previous child.
• Forgetting to clear points on subdivide. Internal nodes that retain their old point lists double-count during queries.
• Bounding-box objects in a point octree. A 3D sprite straddling multiple cells will be assigned to one cell only, missing collisions across boundaries. Use a loose octree or insert into all overlapping cells.
• Frame-by-frame rebuild on dynamic data. Cheap when n is small, ruinous at n = 100k. For dynamic scenes use a refitted BVH or a uniform grid; reserve octrees for static or low-churn data.