Mapping N-dimensional value to a point on Hilbert curve

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

Another possibility would be to build a kd-tree on your data, and then to an in-order traversal of the tree to get the ordering. Constructing the kd-tree only requires you to have a good median-finding algorithm, of which there are many.

I don't see how you can use a Hilbert curve in one dimension.

If you are interested in mapping points to a lower dimension while preserving distances (with minimum error) then you can look into "Multidimensional Scaling" algorithms.

Simulated annealing is one approach.

Edit: Thanks for the comment. I see what you meant by the Hilbert Curve approach now. However, this is a hard problem, and given N=100 and 10 million data points I don't think any approach will preserve locality well and run in a reasonable amount of time. I don't think kd-trees will work here.

If finding a total ordering is not important to you, then you can look into locality-based hashing and other approximate nearest neighbor schemes. Hierarchical multidimensional scaling with buckets of points to reduce the input size might give you a good ordering, but again it's doubtful in such a high dimension.

Here is John Skilling's original C code for encode/decode of Hilbert coordinates in arbitrary dimensions. This is from the paper cited by Paul Chernoch above, Programming the Hilbert curve" by John Skilling (from the AIP Conf. Proc. 707, 381 (2004)).

This code has the bug fix applied.

I also extended main() to show both encoding and decoding. I also added functions interleaveBits() and uninterleaveBits() that convert the Hilbert-transposed coordinates into a single code and back, which is what most people are interested in.

Skilling's code works for arbitrary dimensions, but my interleaveBits() functions are specific to three dimensions. Easy to extend, though.

//+++++++++++++++++++++++++++ PUBLIC-DOMAIN SOFTWARE ++++++++++++++++++++++++++
// Functions: TransposetoAxes AxestoTranspose
// Purpose: Transform in-place between Hilbert transpose and geometrical axes
// Example: b=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]|
//        X[1] = B E H K N <------->       | /X[1]
//        X[2] = C F I L O            axes |/
//               high  low                 0------ X[0]
// Axes are stored conventially as b-bit integers.
// Author: John Skilling 20 Apr 2001 to 11 Oct 2003
//-----------------------------------------------------------------------------

#include <cstdio>

typedef unsigned int coord_t; // char,short,int for up to 8,16,32 bits per word

void TransposetoAxes(coord_t* X, int b, int n) // Position, #bits, dimension
{
    coord_t N = 2 << (b - 1), P, Q, t;

    // 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 (int 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 (int i = n - 1; i >= 0; i--)
            if (X[i] & Q) // Invert
                X[0] ^= P;
            else { // Exchange
                t = (X[0] ^ X[i]) & P;
                X[0] ^= t;
                X[i] ^= t;
            }
    }
}

void AxestoTranspose(coord_t* X, int b, int n) // Position, #bits, dimension
{
    coord_t M = 1 << (b - 1), P, Q, t;

    // Inverse undo
    for (Q = M; Q > 1; Q >>= 1) {
        P = Q - 1;
        for (int i = 0; i < n; i++)
            if (X[i] & Q) // Invert
                X[0] ^= P;
            else { // Exchange
                t = (X[0] ^ X[i]) & P;
                X[0] ^= t;
                X[i] ^= t;
            }
    }

    // Gray encode
    for (int i = 1; i < n; i++) X[i] ^= X[i - 1];
    t = 0;
    for (Q = M; Q > 1; Q >>= 1)
        if (X[n - 1] & Q) t ^= Q - 1;
    for (int i = 0; i < n; i++) X[i] ^= t;
}

int interleaveBits(coord_t* X, int b, int n) // Position, #bits, dimension
{
    unsigned int codex = 0, codey = 0, codez = 0;

    const int nbits2 = 2 * b;

    for (int i = 0, andbit = 1; i < nbits2; i += 2, andbit <<= 1) {
        codex |= (unsigned int)(X[0] & andbit) << i;
        codey |= (unsigned int)(X[1] & andbit) << i;
        codez |= (unsigned int)(X[2] & andbit) << i;
    }

    return (codex << 2) | (codey << 1) | codez;
}

// From https://github.com/Forceflow/libmorton/blob/main/include/libmorton/morton3D.h
void uninterleaveBits(coord_t* X, int b, int n, unsigned int code) // Position, #bits, dimension
{
    X[0] = X[1] = X[2] = 0;

    for (unsigned int i = 0; i <= b; ++i) {
        unsigned int selector = 1;
        unsigned int shift_selector = 3 * i;
        unsigned int shiftback = 2 * i;
        X[2] |= (code & (selector << shift_selector)) >> (shiftback);
        X[1] |= (code & (selector << (shift_selector + 1))) >> (shiftback + 1);
        X[0] |= (code & (selector << (shift_selector + 2))) >> (shiftback + 2);
    }
}

int main()
{
    coord_t X[3] = {5, 10, 20}; // Any position in 32x32x32 cube
    printf("Input coords = %d,%d,%d\n", X[0], X[1], X[2]);

    AxestoTranspose(X, 5, 3); // Hilbert transpose for 5 bits and 3 dimensions
    printf("Hilbert coords = %d,%d,%d\n", X[0], X[1], X[2]);

    unsigned int code = interleaveBits(X, 5, 3);

    printf("Hilbert integer = %d = %d%d%d %d%d%d %d%d%d %d%d%d %d%d%d = 7865 check\n", code,
           X[0] >> 4 & 1, X[1] >> 4 & 1, X[2] >> 4 & 1, 
           X[0] >> 3 & 1, X[1] >> 3 & 1, X[2] >> 3 & 1, 
           X[0] >> 2 & 1, X[1] >> 2 & 1, X[2] >> 2 & 1, 
           X[0] >> 1 & 1, X[1] >> 1 & 1, X[2] >> 1 & 1, 
           X[0] >> 0 & 1, X[1] >> 0 & 1, X[2] >> 0 & 1);

    uninterleaveBits(X, 5, 3, code);
    printf("Reconstructed Hilbert coords = %d,%d,%d\n", X[0], X[1], X[2]);

    TransposetoAxes(X, 5, 3);
    printf("Orig coords = %d,%d,%d\n", X[0], X[1], X[2]);

    return 0;
}

EDIT: I took this code and paired it with similar code for other space-filling curves (Morton, Raster, Boustrophedonic, and Tiled) and made them all available on GitHub. I include both the forward and inverse transforms for all curves, and some code that measures them against each other for various quality metrics. See https://github.com/davemc0/DMcTools/blob/main/Math/SpaceFillCurve.h and for the testing code https://github.com/davemc0/DMcTools/blob/main/Test/SpaceFillCurveTest.cpp .