Quadtree and Octree, in depth: PR quadtrees, depth versus point spacing, Morton keys, Barnes-Hut and GPU builds

Point, PR and region quadtrees

There are two families, and introductory descriptions often blur them.

A point quadtree (Finkel and Bentley, 1974) splits at the data points themselves: the first point inserted becomes the root and divides the plane into four quadrants around itself. It is the two-dimensional cousin of a binary search tree, and it shares that tree's weakness. Its shape depends on insertion order, and sorted input produces a long chain. Deleting a point means rebuilding the subtree under it, because the point defines its children's regions.

A point-region (PR) quadtree splits regions at their midpoints, whatever the data. Its cell boundaries are fixed by the domain, so the same set of points always produces the same tree, in any order. Deletion is local. This is the variant almost everyone means today, and the one used below. A region quadtree applies the same midpoint rule to a raster image, splitting until each cell is uniform in colour. An octree is the PR rule applied to a cube, with eight children per node.

Depth depends on spacing, not on n

The common claim that a quadtree has O(log n) depth is false for real data. A PR quadtree keeps splitting until each leaf holds at most its capacity, so two points that are very close force splits until a cell boundary falls between them. With capacity 1, de Berg and colleagues' textbook bound puts the depth at most log2(s / c) + 3/2, where s is the side of the root square and c is the smallest distance between two points. The depth depends on the ratio of the extent to the closest pair, and the number of points does not enter it at all.

Measured with the code below, capacity 8, 100,000 points each:

The node counts are similar, but the clustered tree is almost twice as deep, and a uniform data set would need billions of points to reach that depth. Query cost still tracks the answer size, because the search skips every empty or distant cell. Depth costs you on the descent and in recursion, and it is what blows up when points coincide.

A bucketed PR quadtree in Python

The tree below is a bucketed PR quadtree over a half-open square. Points on a midpoint go to the upper or right child by the >= rule, so every point has exactly one home. The depth guard is not optional, as the failure-modes section shows.

class Node:
    __slots__ = ("x0", "y0", "size", "pts", "kids")
    def __init__(self, x0, y0, size):
        self.x0, self.y0, self.size = x0, y0, size
        self.pts, self.kids = [], None

class PRQuadtree:
    def __init__(self, x0, y0, size, capacity=8, max_depth=24):
        self.root = Node(x0, y0, size)
        self.cap, self.max_depth = capacity, max_depth

    def insert(self, p):
        n, depth = self.root, 0
        while n.kids is not None:
            n, depth = self._child(n, p), depth + 1
        n.pts.append(p)
        if len(n.pts) > self.cap and depth < self.max_depth:
            self._split(n, depth)

    def _child(self, n, p):
        h = n.size / 2
        i = (p[0] >= n.x0 + h) + 2 * (p[1] >= n.y0 + h)
        return n.kids[i]

    def _split(self, n, depth):
        h = n.size / 2
        n.kids = [Node(n.x0, n.y0, h), Node(n.x0 + h, n.y0, h),
                  Node(n.x0, n.y0 + h, h), Node(n.x0 + h, n.y0 + h, h)]
        pts, n.pts = n.pts, []
        for p in pts:
            self._child(n, p).pts.append(p)
        for k in n.kids:            # all points may land in one child: cascade
            if len(k.pts) > self.cap and depth + 1 < self.max_depth:
                self._split(k, depth + 1)

    def range(self, qx0, qy0, qx1, qy1):
        out, stack = [], [self.root]
        while stack:
            n = stack.pop()
            if (n.x0 >= qx1 or n.x0 + n.size <= qx0 or
                    n.y0 >= qy1 or n.y0 + n.size <= qy0):
                continue                     # cell misses the query box
            if n.kids is None:
                out.extend(p for p in n.pts
                           if qx0 <= p[0] < qx1 and qy0 <= p[1] < qy1)
            else:
                stack.extend(n.kids)
        return out

Test it the way the numbers above were checked: compare every range query against a brute-force scan over random boxes, including boxes that touch cell boundaries exactly. Nearest-neighbour search uses the same tree with a priority queue ordered by the distance from the query to each cell, pruning any cell farther away than the best point found so far.

Morton keys and linear quadtrees

Pointer-based trees are awkward on disk and on GPUs. The linear quadtree drops the pointers. Quantise each coordinate to an integer, interleave the bits of x and y into a single Morton key, and sort. Every quadtree cell then corresponds to a contiguous range of keys: the cell at depth d is the set of keys sharing a 2d-bit prefix. A B-tree or a sorted array becomes a quadtree. Many geospatial indexes use exactly this idea under different names.

def part1by1(v):                  # spread the low 16 bits: abcd -> 0a0b0c0d
    v &= 0xFFFF
    v = (v | (v << 8)) & 0x00FF00FF
    v = (v | (v << 4)) & 0x0F0F0F0F
    v = (v | (v << 2)) & 0x33333333
    v = (v | (v << 1)) & 0x55555555
    return v

def morton2(x, y):
    return part1by1(x) | (part1by1(y) << 1)

assert morton2(5, 3) == 27        # x=101, y=011 -> 011011

The inverse, compacting every other bit back out, round-trips over the full 16-bit range in the test harness. For an octree, interleave three coordinates; 21 bits per axis fill a 63-bit key. Two properties make this useful. Nearby points usually get nearby keys, so a sorted array has good cache and disk locality. And a rectangular query becomes a small set of key intervals: the BIGMIN computation of Tropf and Herzog finds the next key that re-enters the box, so a scan can skip ahead instead of filtering. The Z curve has jumps, so a query box can map to many intervals. Hilbert keys have fewer jumps but cost more to compute.

Octrees: Barnes-Hut, point clouds and GPU builds

Octrees appear wherever 3D space is mostly empty. Three uses show the range.

Barnes-Hut N-body simulation. Each internal node stores the total mass and centre of mass of its subtree. To compute the force on a body, walk the tree, and if a node's cell side divided by its distance to the body is below an opening angle theta (0.5 is a common choice), treat the whole node as a single mass. This replaces the all-pairs O(n^2) sum with roughly O(n log n) work per step for reasonably distributed bodies, at a controllable accuracy cost. Fast multipole methods refine the same octree with series expansions.

Point clouds and voxels. Lidar scans and 3D reconstructions are stored in octrees for downsampling (keep one point per leaf at a chosen depth), neighbour search for normal estimation, and change detection between scans. Sparse voxel octrees store only occupied cells, which turns a 1024^3 grid that would hold a billion cells into a structure proportional to the surface area.

GPU construction. Karras (2012) showed how to build a radix tree over sorted Morton codes fully in parallel, then derive BVHs, octrees and k-d trees from it. The pipeline is: compute keys, radix sort, build internal nodes from the longest common prefixes of neighbouring keys, then aggregate bounds or masses bottom up. Sorting dominates, so rebuilding the whole tree every frame is often cheaper than updating it, which is why game physics and simulation codes rebuild rather than rebalance.

Objects with extent

Points have no size; game objects, map features and bounding boxes do. An object that straddles a midpoint belongs to no single child. There are three standard answers. Store it at the smallest node that fully contains it, as the MX-CIF quadtree does. Objects crossing the root's centre lines then pile up at the root and get tested against every query. Duplicate it into every leaf it overlaps, which makes queries fast but forces result deduplication and multiplies update cost. Or use a loose quadtree or octree, where each node's bounds are enlarged (commonly to twice the side) so an object can be filed by its centre and size alone. Above roughly this point, an R-tree, which is built for rectangles, is usually the better tool.

Failure modes

These are the failures that reach production:

• Duplicate points. Five identical points in a capacity-4 tree can never be separated. With the depth guard the tree stops at depth 24 with a chain of 97 nodes. Without it the split recurses until the stack overflows. Points 1e-9 apart behave the same way, since separating them would need about 30 levels. Cap the depth and let the deepest leaf exceed its capacity, or store a count per distinct point.
• Points on the boundary of the domain. With half-open cells a point at exactly x = 1.0 lies outside the unit-square root. Pad the domain or clamp keys, and decide this once, in one function.
• Floating-point midpoints. Once a cell shrinks to the spacing between adjacent doubles (about 52 halvings of a unit square away from the origin), the midpoint stops differing from the corner and children become empty copies of their parent. Quantise to integers and use Morton keys if you need deep trees.
• Objects stuck at the root. See the previous section; profile how many objects live in the top two levels.
• Recursive traversal. A deep clustered tree can exhaust a recursion limit. The range query above uses an explicit stack for that reason.

Operational guidance and trade-offs

Choose the leaf capacity from the cost of a leaf scan versus a node visit; 8 to 64 points per leaf is typical, and larger buckets suit SIMD distance checks. Set the domain and resolution explicitly rather than from the first batch of data. Monitor maximum depth, the number of points per leaf at the depth cap, and query nodes visited per hit; a rising ratio means the data distribution changed. For moving objects, measure rebuild time against incremental update time before choosing; sort-based rebuilds often win.

What to do next

1. Implement the PR quadtree above and check range queries against brute force, including points and boxes on cell boundaries.
2. Measure maximum depth on your real data and compare it to log2(extent / closest pair), not to log n.
3. Add a duplicate-point test and confirm the depth cap holds.
4. Rebuild the index as sorted Morton keys and compare query time and memory with the pointer tree.
5. Compare against a median-split k-d tree and an R-tree on the same queries.
6. Use the tree to speed up neighbour queries in DBSCAN, and read the k-d tree pruning and best-first search article for the median-split alternative.