space-filling curve

Recursive Formulation of Multidimensional Hilbert Space-Filling Curves

There is a surfeit of 2-dimensional space-filling curves, but generalizing them to higher rank is not necessarily practical or even possible. [Bongki 2001] attributes generalization of space-filling curves to any rank to [Butz 1969]; who gives algorithms for the Hilbert and the Peano space-filling curves in separate papers.

The rigorous explanations of algorithms for multidimensional space-filling functions tend to be cryptic and difficult to reconcile with their simple geometric interpretations. This article attempts to rectify that.


The properties we desire for space-filling transforms are that they be:

a continuous function whose range comes arbitrarily close to any given point in the n-dimensional unit hypercube.
exist for all ranks greater than one.
can be defined as a nonterminating recurrence.
Lebesgue measure-preserving.
equal length intervals map to equal volumes.
the lengths of axis-aligned projections of the space-filling curve should be equal.
Although there are many shapes which can tile the plane, the only regular polyhedron which tiles R3 is the cube. Isotropy rules out parallelograms and their ilk. The general family of curves will fill squares, cubes, and hypercubes.

Because the transforms are self-similar, the scalar can be considered to be composed of digits in some base, where digits correspond to positions in nested shapes. Because transforms are measure-preserving, the largest copies of the self-similar shape correspond to the highest-order digits of the scalar.

Any real scalar has a most significant digit; because the least-signficant digits could be arbitrarily small, the same can't be said of them. Thus the calculation of transformations between scalars and multidimensional coordinates must start with the highest-order digits; and information must pass from higher-order to lower-order frames.

This suggests a recursive organization for conversion computations. Were the task to sequentially enumerate points on the curve, then recursive computation would have the familiar profile where the number of calls increases exponentially with the depth of calls. But that is only useful for drawing pictures. The other uses convert or compare a few points sparsely scattered in the multidimensional space. The computation will return a single value: a boolean, scalar or vector (coordinates). This is still recursive; but tail-recursive, which is equalent to the interative loops in the extant algorithms.

A Gray code is an ordering of non-negative integers in which exactly one bit differs between each pair of successive elements. There are multiple Gray codings. An n-bit Gray code maps the integers from 0 to 2n to a Hamiltonian cycle on an n-dimensional hypercube. Gray-codes are thus used to generate offsets at each hypercube level.


A cube's rotation is the combination of the rotations of all the cubes it is part of. Thus the rotation is passed to the recursive function and each call augments it.

Although any permutation which cycles all axes through all positions could be used to "rotate" hypercubes in order to align them with the corners they they are replacing, cyclic permuations are sufficient and used here. The rotate-bit-field function is used to rotate fields of rank bits.

But how much rotation should be applied? The coordinates of a vertex alone do not tell which direction the next edge is to be drawn. That direction is controlled by the Gray-code sequence. The log2-binary-factors machinations appear to compute the additional rotation based on the next bit to be changed in the Gray-code sequence.

log2-binary-factors is a simple but non-obvious function. I encountered it in Doug Moore's "hilbert.c" ([Moore 1999]). Here is the scheme version:

(define (log2-binary-factors n)
  (+ -1 (integer-length (logand n (- n)))))

The log2-binary-factors function returns the number of factors-of-2 of chnk, which is also the index of the lowest-order 1 bit in chnk. But more importantly for this application, log2-binary-factors is related to the next-bit-changed in incrementing Gray-codes.

For Gray-codes with odd numbers of 1 bits (which correspond to odd integers), log2-binary-factors returns one less than the index of the bit different in the next Gray-code in the sequence.

For Gray-codes with even numbers of bits (which correspond to even integers), log2-binary-factors returns one less than the index of the bit different in the previous Gray-code in the sequence.

sequenceGray codelog2-binary-factors + 2sequenceGray codelog2-binary-factors + 2
1000012 2000112
3000103 4001103
5001112 6001012
7001004 8011004
9011012 10011112
11011103 12010103
13010112 14010012
15010005 16110005

The other points in the sequence all are changes of the low-order bit.

The Algorithm

Recursive Formulation of integer->hilbert-coordinates
(define (integer->hilbert-coordinates scalar rank . nbits)
  (define igry (integer->gray-code scalar))
  (define rnkmsk (lognot (ash -1 rank)))
  (define rnkhib (ash 1 (+ -1 rank)))
  (define rank*nbits
    (if (null? nbits)
        (let ((rank^2 (* rank rank)))
          (* (quotient (+ -1 rank^2 (integer-length scalar)) rank^2)
        (* rank (car nbits))))
  (define (loop bdxn chnk rotation flipbit lst)
    (if (positive? bdxn)
        (map gray-code->integer (delaminate-list rank (reverse lst)))
        (loop (+ rank bdxn)
              (logxor rnkhib (logand (ash igry (+ rank bdxn))) rnkmsk)
              (modulo (+ rotation (log2-binary-factors chnk) 2) rank)
              (ash 1 rotation)
              (cons (logxor flipbit (rotate-bit-field chnk rotation 0 rank))
  (loop (- rank rank*nbits) (ash igry (- rank rank*nbits)) 0 0 '()))


Example Calculation

Convert the scalar #b0011100001011010101011000 to rank 5 Hilbert coordinates.

rank*nbits is 25.

(integer->gray-code #b0011100001011010101011000)
                ==> #b0010010001110111111110100
    xor rnkhib  ==> #b0010000001010110111100100
The high bit of each rank-bit group except the highest order group are flipped (by rnkhib) and are underlined. The leftmost number in lst is derived from the chnk rotated by rotation, then exclusive-ored with the flipbit of the preceeding row.

-1001011110000(10001 00100)
-501111300010(00110 10001 00100)
000100001000(11001 00110 10001 00100)
510000400001(01100 11001 00110 10001 00100)

From left to right, lst gets reversed (top to bottom), delaminated (matrix transpose), and Gray-decoded.


(integer->hilbert-coordinates #b0011100001011010101011000 5)
                ==> (#b01100 #b00010 #b11001 #b00111 #b01100)

Further Work

Although the SLIB code works and replicates Butz's example calculations, my attempt to fully understand the algorithm has failed. Questions remain:
  1. Why is the entire scalar converted to Gray-code, rather than Gray-code converting the individual chnks?

  2. Why do the returned coordinates undergo gray-code->integer; is this related to item 1?

  3. Why is the high bit of each (except the first) chnk flipped; is this related to flipbit?

  4. The chnk is already rotated by rotation when the bit indexed by rotation is flipped (by flipbit). What relation does this bit position have to this or the previous chnk?

  5. Why does the rotation calculation using log2-binary-factors ignore changes in the low order bit of the Gray-coded chnk?

Copyright © 2005 Aubrey Jaffer

I am a guest and not a member of the MIT Computer Science and Artificial Intelligence Laboratory.  My actions and comments do not reflect in any way on MIT.
Space-Filling Curves
agj @
Go Figure!