Mapping N-dimensional value to a point on Hilbert curve
Algorithm for mapping from n->1 and 1->n given here "Calculation of Mappings Between One and n-dimensional Values Using the Hilbert Space-filling Curve" J K Lawder
If you Google for "SFC module and Kademlia overlay" youl find a group who claim to use it as part of their system. If you view the source you can probably extract the relevant function.
I spent a little time translating Paul Chernoch's code to Java and cleaning it up. It's possible that there is a bug in my code, especially because I don't have access to the paper that it's originally from. However, it passes what unit tests I was able to write. It is below.
Note that I have evaluated both Z-Order and Hilbert curves for spatial indexing on largish data sets. I have to say that Z-Order provides far better quality. But feel free to try for yourself.
/**
* Convert the Hilbert index into an N-dimensional point expressed as a vector of uints.
*
* Note: In Skilling's paper, this function is named TransposetoAxes.
* @param transposedIndex The Hilbert index stored in transposed form.
* @param bits Number of bits per coordinate.
* @return Point in N-space.
*/
static long[] HilbertAxes(final long[] transposedIndex, final int bits) {
final long[] result = transposedIndex.clone();
final int dims = result.length;
grayDecode(result, dims);
undoExcessWork(result, dims, bits);
return result;
}
static void grayDecode(final long[] result, final int dims) {
final long swap = result[dims - 1] >>> 1;
// Corrected error in Skilling's paper on the following line. The appendix had i >= 0 leading to negative array index.
for (int i = dims - 1; i > 0; i--)
result[i] ^= result[i - 1];
result[0] ^= swap;
}
static void undoExcessWork(final long[] result, final int dims, final int bits) {
for (long bit = 2, n = 1; n != bits; bit <<= 1, ++n) {
final long mask = bit - 1;
for (int i = dims - 1; i >= 0; i--)
if ((result[i] & bit) != 0)
result[0] ^= mask; // invert
else
swapBits(result, mask, i);
}
}
/**
* Given the axes (coordinates) of a point in N-Dimensional space, find the distance to that point along the Hilbert curve.
* That distance will be transposed; broken into pieces and distributed into an array.
*
* The number of dimensions is the length of the hilbertAxes array.
*
* Note: In Skilling's paper, this function is called AxestoTranspose.
* @param hilbertAxes Point in N-space.
* @param bits Depth of the Hilbert curve. If bits is one, this is the top-level Hilbert curve.
* @return The Hilbert distance (or index) as a transposed Hilbert index.
*/
static long[] HilbertIndexTransposed(final long[] hilbertAxes, final int bits) {
final long[] result = hilbertAxes.clone();
final int dims = hilbertAxes.length;
final long maxBit = 1L << (bits - 1);
inverseUndo(result, dims, maxBit);
grayEncode(result, dims, maxBit);
return result;
}
static void inverseUndo(final long[] result, final int dims, final long maxBit) {
for (long bit = maxBit; bit != 0; bit >>>= 1) {
final long mask = bit - 1;
for (int i = 0; i < dims; i++)
if ((result[i] & bit) != 0)
result[0] ^= mask; // invert
else
swapBits(result, mask, i);
} // exchange
}
static void grayEncode(final long[] result, final int dims, final long maxBit) {
for (int i = 1; i < dims; i++)
result[i] ^= result[i - 1];
long mask = 0;
for (long bit = maxBit; bit != 0; bit >>>= 1)
if ((result[dims - 1] & bit) != 0)
mask ^= bit - 1;
for (int i = 0; i < dims; i++)
result[i] ^= mask;
}
static void swapBits(final long[] array, final long mask, final int index) {
final long swap = (array[0] ^ array[index]) & mask;
array[0] ^= swap;
array[index] ^= swap;
}
It's not clear to me how this will do what you want. Consider this trival 3D case:
001 ------ 101
|\ |\
| \ | \
| 011 ------ 111
| | | |
| | | |
000 -|---- 100 |
\ | \ |
\ | \ |
010 ------ 110
which can be "Hilbertized" by the following path:
001 -----> 101
\ \
\ \
011 111
^ |
| |
000 | 100 |
\ | \ |
\ | \ V
010 110
into the 1D order:
000 -> 010 -> 011 -> 001 -> 101 -> 111 -> 110 -> 100
Here's the nasty bit. Consider the list of pairs and 1D distances below:
000 : 100 -> 7
010 : 110 -> 5
011 : 111 -> 3
001 : 101 -> 1
In all cases, the left- and right-hand values are the same 3D distance from each other (+/- 1 in the first position), which seem to imply similar "spatial locality". But linearizing by any choice of dimensional ordering (y, then z, then z, in the above example) breaks that locality.
Another way of saying this is that taking a starting point and ordering the remaining points by their distance from that starting point will provide significantly different results. Taking 000 as the start, for example:
1D ordering : distance 3D ordering : distance
---------------------- ----------------------
010 : 1 001,010,100 : 1
011,101,110 : sqrt(2)
111 : sqrt(3)
011 : 2
001 : 3
101 : 4
111 : 5
110 : 6
100 : 7
This effect grows exponentially with the number of dimensions (assuming that each dimension has the same "size").
I finally broke down and shelled out some money. The AIP (American Institute of Physics) has a nice, short article with source code in C. "Programming the Hilbert curve" by John Skilling (from the AIP Conf. Proc. 707, 381 (2004)) has an appendix with code for mappings in both directions. It works for any number of dimensions > 1, is not recursive, does not use state-transition lookup tables that gobble up huge amounts of memory, and mostly uses bit operations. Thus it is reasonably fast and has a good memory footprint.
If you choose to purchase the article, I discovered an error in the source code.
The following line of code (found in the function TransposetoAxes) has the error:
for( i = n-1; i >= 0; i-- ) X[i] ^= X[i-1];
The correction is to change the greater than or equal (>=) to a greater than (>). Without this correction, the X array is accessed using a negative index when the variable “i” becomes zero, causing the program to fail.
I recommend reading the article (which is seven pages long, including the code), as it explains how the algorithm works, which is far from obvious.
I translated his code into C# for my own use. The code follows. Skilling performs the transformation in place, overwriting the vector that you pass in. I chose to make a clone of the input vector and return a new copy. Also, I implemented the methods as extension methods.
Skilling's code represents the Hilbert index as a transpose, stored as an array. I find it more convenient to interleave the bits and form a single BigInteger (more useful in Dictionaries, easier to iterate over in loops, etc), but I optimized that operation and its inverse with magic numbers, bit operations and the like, and the code is lengthy, so I have omitted it.
namespace HilbertExtensions
{
/// <summary>
/// Convert between Hilbert index and N-dimensional points.
///
/// The Hilbert index is expressed as an array of transposed bits.
///
/// Example: 5 bits for each of n=3 coordinates.
/// 15-bit Hilbert integer = A B C D E F G H I J K L M N O is stored
/// as its Transpose ^
/// X[0] = A D G J M X[2]| 7
/// X[1] = B E H K N <-------> | /X[1]
/// X[2] = C F I L O axes |/
/// high low 0------> X[0]
///
/// NOTE: This algorithm is derived from work done by John Skilling and published in "Programming the Hilbert curve".
/// (c) 2004 American Institute of Physics.
///
/// </summary>
public static class HilbertCurveTransform
{
/// <summary>
/// Convert the Hilbert index into an N-dimensional point expressed as a vector of uints.
///
/// Note: In Skilling's paper, this function is named TransposetoAxes.
/// </summary>
/// <param name="transposedIndex">The Hilbert index stored in transposed form.</param>
/// <param name="bits">Number of bits per coordinate.</param>
/// <returns>Coordinate vector.</returns>
public static uint[] HilbertAxes(this uint[] transposedIndex, int bits)
{
var X = (uint[])transposedIndex.Clone();
int n = X.Length; // n: Number of dimensions
uint N = 2U << (bits - 1), P, Q, t;
int i;
// Gray decode by H ^ (H/2)
t = X[n - 1] >> 1;
// Corrected error in Skilling's paper on the following line. The appendix had i >= 0 leading to negative array index.
for (i = n - 1; i > 0; i--)
X[i] ^= X[i - 1];
X[0] ^= t;
// Undo excess work
for (Q = 2; Q != N; Q <<= 1)
{
P = Q - 1;
for (i = n - 1; i >= 0; i--)
if ((X[i] & Q) != 0U)
X[0] ^= P; // invert
else
{
t = (X[0] ^ X[i]) & P;
X[0] ^= t;
X[i] ^= t;
}
} // exchange
return X;
}
/// <summary>
/// Given the axes (coordinates) of a point in N-Dimensional space, find the distance to that point along the Hilbert curve.
/// That distance will be transposed; broken into pieces and distributed into an array.
///
/// The number of dimensions is the length of the hilbertAxes array.
///
/// Note: In Skilling's paper, this function is called AxestoTranspose.
/// </summary>
/// <param name="hilbertAxes">Point in N-space.</param>
/// <param name="bits">Depth of the Hilbert curve. If bits is one, this is the top-level Hilbert curve.</param>
/// <returns>The Hilbert distance (or index) as a transposed Hilbert index.</returns>
public static uint[] HilbertIndexTransposed(this uint[] hilbertAxes, int bits)
{
var X = (uint[])hilbertAxes.Clone();
var n = hilbertAxes.Length; // n: Number of dimensions
uint M = 1U << (bits - 1), P, Q, t;
int i;
// Inverse undo
for (Q = M; Q > 1; Q >>= 1)
{
P = Q - 1;
for (i = 0; i < n; i++)
if ((X[i] & Q) != 0)
X[0] ^= P; // invert
else
{
t = (X[0] ^ X[i]) & P;
X[0] ^= t;
X[i] ^= t;
}
} // exchange
// Gray encode
for (i = 1; i < n; i++)
X[i] ^= X[i - 1];
t = 0;
for (Q = M; Q > 1; Q >>= 1)
if ((X[n - 1] & Q)!=0)
t ^= Q - 1;
for (i = 0; i < n; i++)
X[i] ^= t;
return X;
}
}
}
I have posted working code in C# to github.
See https://github.com/paulchernoch/HilbertTransformation
UPDATE: I just published (Fall 2019) a Rust crate on crates.io called "hilbert". It also uses Skilling's algorithm. See https://crates.io/crates/hilbert