Hilbert Curve

Why a continuous space-filling curve matters

The fundamental question of space-filling curves is: can you visit every point in a 2D grid by following a path that never jumps, and that produces a 1D ordering whose nearby indices correspond to nearby cells? Morton's Z-order is one answer, but it has discontinuities — consecutive indices can leap across the grid. The Hilbert curve is the better answer. Its path is continuous: every consecutive pair of cells shares an edge.

David Hilbert constructed the curve in 1891 to demonstrate that the unit interval and the unit square have the same cardinality — a counterintuitive result for the 19th century. The construction is recursive and produces a curve with infinite length but zero area, packed into a finite square. For computer science, we don't take the limit; we use the finite-order approximations.

Recursive construction

The order-1 Hilbert curve in a 2×2 grid visits four cells in a U shape:

order 1 (4 cells)        order 2 (16 cells)

  1  2                     5  6   9 10
  |  |                     |  |   |  |
  0  3                     4  7---8 11
                            \         /
                             3       12
                             |        |
                             2  1---14 13
                                |   |
                                0  15

To go from order n to order n+1, divide the square into four quadrants. Inside each quadrant, place a smaller copy of the order-n curve, rotated or reflected so that the four sub-curves connect end-to-end into a single continuous path. The standard convention:

• Bottom-left quadrant: rotated 90° clockwise.
• Top-left quadrant: unrotated.
• Top-right quadrant: unrotated.
• Bottom-right quadrant: rotated 90° counter-clockwise.

The exit of each sub-curve aligns with the entry of the next. Result: an order n+1 curve through 4ⁿ⁺¹ cells, still continuous at every step.

Encoding: coordinate → Hilbert index

The textbook algorithm walks bit-by-bit from MSB to LSB. At each level, you look at one bit of x and one bit of y, consult a small state-transition table, and emit two bits of the Hilbert index. The state tracks the rotation/reflection of the current quadrant.

// Encode (x, y) on a 2^order x 2^order grid into Hilbert index.
uint64_t hilbert_encode(uint32_t x, uint32_t y, int order) {
    uint64_t d = 0;
    for (int s = (1u << (order - 1)); s > 0; s >>= 1) {
        uint32_t rx = (x & s) ? 1 : 0;
        uint32_t ry = (y & s) ? 1 : 0;
        d += (uint64_t)s * (uint64_t)s * ((3 * rx) ^ ry);
        // rotate the quadrant
        if (ry == 0) {
            if (rx == 1) {
                x = s - 1 - x;
                y = s - 1 - y;
            }
            uint32_t t = x; x = y; y = t;
        }
    }
    return d;
}

It's about 50 lines including the inverse direction. Runtime per encode is roughly 20 nanoseconds on a modern x86 — slower than Morton's 3 ns because the state machine prevents the single-instruction PDEP trick. For workloads dominated by range scans rather than encoding throughput, the locality advantage usually wins.

Decoding: Hilbert index → coordinate

Symmetric: walk the bits of the index, two at a time, applying the inverse rotation. Each pair of index bits identifies the quadrant; the rotation state advances to the sub-quadrant's frame.

void hilbert_decode(uint64_t d, int order, uint32_t* x, uint32_t* y) {
    uint32_t cx = 0, cy = 0;
    for (uint32_t s = 1; s < (1u << order); s <<= 1) {
        uint32_t rx = 1 & (d / 2);
        uint32_t ry = 1 & (d ^ rx);
        if (ry == 0) {
            if (rx == 1) {
                cx = s - 1 - cx;
                cy = s - 1 - cy;
            }
            uint32_t t = cx; cx = cy; cy = t;
        }
        cx += s * rx;
        cy += s * ry;
        d /= 4;
    }
    *x = cx; *y = cy;
}

Application: R-tree indexing

An R-tree groups spatial objects into nested bounding rectangles for fast range and nearest-neighbor queries. The quality of the index depends on rectangle clustering — overlapping rectangles cost time during queries because the search must descend into both. Hilbert R-trees, introduced by Kamel and Faloutsos in 1994, insert objects in Hilbert order so spatially adjacent objects end up in the same leaf. Bounding rectangles stay small and roughly square. Empirically, Hilbert R-trees outperform Morton-ordered R-trees by 10-30% on range queries and 5-15% on k-NN.

Production databases use this. PostGIS's GiST indexes can be built with Hilbert-ordered packing. SQL Server's spatial indexes use a related "tessellation" scheme that's effectively Hilbert order.

Application: image compression

JPEG-LS, the "near-lossless" JPEG variant, optionally orders pixels along a Hilbert curve before entropy coding. Because adjacent pixels on the curve are spatially adjacent, the local correlation is preserved and the entropy coder gets a more compressible sequence. Compression ratios improve by 2-5% over row-major scanning for natural images.

The same idea appears in academic work on neural-image codecs and in some game-asset compression tools. The cost is the Hilbert encode/decode per pixel; the benefit is fewer bits per coefficient at the entropy stage.

Application: geospatial database keys

If you need to put a billion latitude/longitude points into a B-tree, the question is how to define a 1D sort key that makes range queries fast. Hilbert encoding is the gold-standard answer: pack (lat, lon) into a single 64-bit key by walking the recursive construction. Range queries on the key are mostly contiguous; out-of-range gaps are bounded and small.

Uber's H3 and Google's S2 are sphere-aware variants that solve the same problem on a sphere (Hilbert curves on the icosahedron and on the cube faces respectively). The fundamental trick is the same: a locality-preserving 1D order over 2D coordinates lets you reuse B-tree machinery for spatial queries.

Hilbert vs Morton vs row-major

	Row-major	Morton (Z-order)	Hilbert
Adjacent in 1D ⇒ adjacent in 2D	Within rows only	Usually, but with jumps	Always
Max 2D distance per index step	width	up to 2ⁿ⁻¹	1 cell
Encode time (modern CPU)	1 ns	3 ns (BMI2)	~20 ns
Code length	1 line	~10 lines	~50 lines
Range scan locality (cache hit %)	~50%	~92%	~95%
R-tree pack quality	poor	good	best
Generalizes to n-D	No	Trivially	Yes (Butz 1971)
Best for	Sequential row scans	GPU swizzle, octrees	R-trees, DB indexing

Python — encode, decode, render

def hilbert_encode(x, y, order):
    d = 0
    s = 1 << (order - 1)
    while s > 0:
        rx = 1 if (x & s) else 0
        ry = 1 if (y & s) else 0
        d += s * s * ((3 * rx) ^ ry)
        # rotate
        if ry == 0:
            if rx == 1:
                x = s - 1 - x
                y = s - 1 - y
            x, y = y, x
        s //= 2
    return d

def hilbert_decode(d, order):
    x, y = 0, 0
    s = 1
    while s < (1 << order):
        rx = 1 & (d >> 1)
        ry = 1 & (d ^ rx)
        if ry == 0:
            if rx == 1:
                x = s - 1 - x
                y = s - 1 - y
            x, y = y, x
        x += s * rx
        y += s * ry
        d //= 4
        s *= 2
    return x, y

# Order 2 -> 16 cells. Walk them in Hilbert order:
for i in range(16):
    print(i, hilbert_decode(i, 2))

Variants and curiosities

Peano curve

The Peano curve (1890) is the original space-filling curve, predating Hilbert by one year. It uses a 3×3 base shape instead of 2×2 — recursion produces 9ⁿ cells. Same locality property as Hilbert, but the implementation is less common because powers of 3 are awkward for digital representations.

Moore curve

A closed variant: like Hilbert but the curve returns to its starting cell at the end. Useful when you want a cyclic 1D ordering, e.g., for traveling-salesman heuristics on a grid.

Sierpinski curve

Triangular-cell space-filling curve. Used in some triangular-mesh applications and as a base for finite-element ordering.

Hilbert curve in higher dimensions

Butz 1971 gives the n-D construction. Implementations exist for d ≤ 6 in public repositories. Beyond that, the number of orientations grows combinatorially and pre-computed tables become impractical.

Common pitfalls

• Assuming Hilbert is faster than row-major. Hilbert encoding takes 20+ nanoseconds; raw indexing is a single multiply-add. Hilbert wins only when the application benefits from locality (range scans, R-trees) — not for sequential memory traversal.
• Mixing up orientations. The recursive construction depends on which quadrant rotation convention you use. Two implementations using different conventions will produce different Hilbert indices for the same coordinate — they're both valid Hilbert curves, just different ones.
• Believing Hilbert preserves distance. It preserves adjacency, not Euclidean distance. Two cells far apart in 1D can still be near in 2D (one of Hilbert's strengths). Two cells near in 1D are always near in 2D (its other strength). But the index difference is not a metric.
• Generalizing badly to non-power-of-two grids. Like Morton, Hilbert encoding assumes 2ⁿ × 2ⁿ. For 1080×1920, you must pad to 2048×2048 and reject out-of-range indices.
• Using Hilbert for GPU texture swizzling. Morton is still the right choice there — encoding cost dominates and the texture cache locality is already 90%+ with Morton.
• Reinventing the encoder. The state-machine encoding is a famous source of off-by-one bugs. Use a tested public library (libhilbert, Skilling's algorithm, or scikit-image's _hilbert module) instead of writing from scratch.