Morton ordering: decompose a grid window into code ranges

Last updated: 3 Oct 20269 min read

tutorial

A Morton code interleaves bits from x and y to order points in a fixed integer grid. It turns each aligned square cell into one contiguous code interval, but an arbitrary rectangle generally needs many intervals. This static index sorts point records by code. A range query recursively visits aligned square cells; it skips cells outside the half-open rectangle and uses binary search to read the interval for every fully covered cell. Descending partial cells continues until they are inside, outside, or one grid unit wide. Multiple IDs at the same coordinate share a code and are all retained. The method reports exact rectangle membership here because it only scans fully covered cells, not because nearby planar points always receive adjacent codes.

Operational case

The same four integer depot points used in the quadtree example are sorted by their interleaved eight-bit coordinates. Querying x from 15 through values below 60 and y from 45 through values below 65 returns pump-47, valve-19, and sensor-61. The grid-83 point is outside. A thin diagonal or narrow rectangle may require many small cell intervals even if few points are returned, so this method should not be advertised as one binary search for every rectangle. It is especially useful as an explanation of how a multidimensional window becomes several one-dimensional runs.

Working Python program

from bisect import bisect_left


def morton_code(x, y, bits):
    code = 0
    for position in range(bits):
        code |= ((x >> position) & 1) << (2 * position)
        code |= ((y >> position) & 1) << (2 * position + 1)
    return code


class MortonDepotIndex:
    def __init__(self, points, bits=8):
        if bits < 1 or bits > 16:
            raise ValueError("bits must be between one and sixteen")
        self.bits = bits
        self.width = 1 << bits
        if len({depot_id for depot_id, _, _ in points}) != len(points):
            raise ValueError("duplicate depot ID")
        if any(not (0 <= x < self.width and 0 <= y < self.width)
               for _, x, y in points):
            raise ValueError("point outside integer grid")
        self.records = sorted((morton_code(x, y, bits), depot_id, x, y)
                              for depot_id, x, y in points)
        self.codes = [record[0] for record in self.records]

    def range_report(self, left, bottom, right, top):
        if not (0 <= left <= right <= self.width and
                0 <= bottom <= top <= self.width):
            raise ValueError("invalid query rectangle")
        found = []

        def visit(x0, y0, width, prefix, remaining):
            x1, y1 = x0 + width, y0 + width
            if left >= x1 or right <= x0 or bottom >= y1 or top <= y0:
                return
            if (left <= x0 and x1 <= right and bottom <= y0 and y1 <= top):
                start = bisect_left(self.codes, prefix)
                end = bisect_left(self.codes, prefix + (1 << (2 * remaining)))
                found.extend(record[1] for record in self.records[start:end])
                return
            if remaining == 0:
                return
            half = width // 2
            for quadrant in range(4):
                visit(x0 + (quadrant & 1) * half,
                      y0 + (quadrant >> 1) * half, half,
                      prefix | (quadrant << (2 * (remaining - 1))), remaining - 1)

        visit(0, 0, self.width, 0, self.bits)
        return sorted(found)


index = MortonDepotIndex([
    ("pump-47", 47, 52), ("valve-19", 19, 61),
    ("sensor-61", 52, 58), ("grid-83", 180, 70),
])
print("window=", index.range_report(15, 45, 60, 65), sep="")

Output

window=['pump-47', 'sensor-61', 'valve-19']

Time, space, and tradeoff

Building N records with B bits per coordinate takes O(NB + N log N) time and O(N) stored records and codes. A query costs O(C log N + R log R), where C fully covered cells in the decomposition and R reported points; recursive boundary visits add work proportional to the visited grid cells. C can grow with grid resolution and rectangle shape. The fixed domain is zero through 2^B minus one per coordinate; inserts after construction require rebuilding this immutable list. This implementation supports only exact integer-grid windows and no proximity or geographic-distance guarantee. Morton order preserves some locality but does not preserve all two-dimensional neighbor relations.

Common Mistakes

Do not assume an arbitrary rectangle maps to one Morton interval.
Do not discard points that share one code.
Do not treat code proximity as proof of geometric proximity.
Do not insert outside the configured integer domain without rebuilding the index.

Connected lessons

Compare its geometry and update costs with Spatial hash grids: move points between occupied cells, Point-region quadtrees: subdivide crowded cells, Bounding-volume hierarchies: prune box overlap searches, then run the spatial audit and contract quiz.

data structures

range-query-structures