Recursive Formulation of Multidimensional Hilbert Space-Filling Curves

	http://people.csail.mit.edu/jaffer/Geometry/RFMDSFF
	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.

Properties

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

space-filling: a continuous function whose range comes arbitrarily close to any given point in the n-dimensional unit hypercube.
multidimensional: exist for all ranks greater than one.
self-similar: can be defined as a nonterminating recurrence.
Lebesgue measure-preserving.: equal length intervals map to equal volumes.
isotropic: 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 R³ 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 2ⁿ to a Hamiltonian cycle on an n-dimensional hypercube. Gray-codes are thus used to generate offsets at each hypercube level.

Rotations

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.

sequence Gray code log2-binary-factors + 2 sequence Gray code log2-binary-factors + 2
0 00000 1
1 00001 2 2 00011 2
3 00010 3 4 00110 3
5 00111 2 6 00101 2
7 00100 4 8 01100 4
9 01101 2 10 01111 2
11 01110 3 12 01010 3
13 01011 2 14 01001 2
15 01000 5 16 11000 5

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 ranknbits (if (null? nbits) (let ((rank^2 ( rank rank))) (* (quotient (+ -1 rank^2 (integer-length scalar)) rank^2) 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)) lst)))) (loop (- rank ranknbits) (ash igry (- rank ranknbits)) 0 0 '()))

(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^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))
                    lst))))
  (loop (- rank rank*nbits) (ash igry (- rank rank*nbits)) 0 0 '()))

Notes

integer->hilbert-coordinates takes an optional argument, nbits.
igry is the Gray-coded representation of input scalar.
rnkmsk is a rank-bit mask: 2^rank-1
rnkhib has just the (high-order) rank bit set: 2^rank-1.
Flipping the high-order bit of a Gray-coded number is the equivalent of taking the ones-complement of the uncoded number.
For 0 <= k < 2*rnkhib:
```
  (= (- (* 2 rnkhib) k 1)
     (gray-code->integer (logxor rnkhib (integer->gray-code k))))
  
```
All but the highest order bdxn get their high-order bit flipped. With flipbit, all bdxn have an even number of flipped bits.
rank*nbits is the number of bits to process from the scalar. If the nbits argument is given, then it is just nbits times rank; otherwise, the smallest integer multiple of rank² greater or equal to the number of bits necessary to represent scalar.
The orientation of nested cubes precesses through all rank axes. To guarantee match of orientation for varied precision inputs, a multiple of rank cubes is processed. There are rank bits in each cube; hence rank² bits are processed.
loop is the recursive function.
If bdxn is positive, then all the digits have been processed; convert the Gray-coded coordinates to regular binary representation and return them.
If bdxn is not positive, then call loop with arguments:
- bdxn is the (negative) number of bits to right shift igry in order to bring the current set of rank bits to lowest order.
- chnk is the rank-bit Gray-coded number currently being processed.
- rotation is the number of bits to rotate the chunk by. Notice that rotation accumulates the earlier rotations.
- flipbit is the bit in the rotated chnk to invert; it is the bit indexed by the calling function's rotation. Note that while rotation is one call delayed from the chnk which generated it, flipbit is two calls delayed.
- The list lst accumulates the results of processing chnk.
The initial chnk is the highest order; it needs no masking.
But if the rank (high) bit of the highest order group is flipped, then zero maps to the other end of the curve from the origin. If the initial flipbit is 1, then zero maps to a point other than the origin.

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.

bxdn chnk rotation flipbit lst
-20 00100 0 00000 ()
-15 00001 4 00001 (00100)
-10 01011 1 10000 (10001 00100)
-5 01111 3 00010 (00110 10001 00100)
0 00100 0 01000 (11001 00110 10001 00100)
5 ~~10000~~ 4 ~~00001~~ (01100 11001 00110 10001 00100)

bxdn	chnk	rotation	flipbit	lst
-20	`00100`	0	`00000`	()
-15	`00001`	4	`00001`	(00100)
-10	`01011`	1	`10000`	(10001 00100)
-5	`01111`	3	`00010`	(00110 10001 00100)
0	`00100`	0	`01000`	(11001 00110 10001 00100)
5	~~`10000`~~	4	~~`00001`~~	(01100 11001 00110 10001 00100)

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

lst reverse delaminate gray-code->integer

01100
11001
00110
10001
00100
==>

00100
10001
00110
11001
01100
==>

01010
00011
10101
00100
01010
==>

01100
00010
11001
00111
01100

(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:

Why is the entire scalar converted to Gray-code, rather than Gray-code converting the individual chnks?
Why do the returned coordinates undergo gray-code->integer; is this related to item 1?
Why is the high bit of each (except the first) chnk flipped; is this related to flipbit?
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?
Why does the rotation calculation using log2-binary-factors ignore changes in the low order bit of the Gray-coded chnk?

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 @ alum.mit.edu	Go Figure!

sequence	Gray code	`log2-binary-factors` + 2	sequence	Gray code	`log2-binary-factors` + 2
			0	`00000`	1
1	`00001`	2	2	`00011`	2
3	`00010`	3	4	`00110`	3
5	`00111`	2	6	`00101`	2
7	`00100`	4	8	`01100`	4
9	`01101`	2	10	`01111`	2
11	`01110`	3	12	`01010`	3
13	`01011`	2	14	`01001`	2
15	`01000`	5	16	`11000`	5