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 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.

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

 ```(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: 2rank-1

• rnkhib has just the (high-order) rank bit set: 2rank-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 rank2 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 rank2 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.

bxdnchnkrotationflipbitlst
-2000100000000()
-1500001400001(00100)
-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.

lstreversedelaminategray-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:
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?