jeudi 13 août 2015

N dimension hilbert index in Delphi

I am trying calculate a 3d coordinate into a Hilbert index and back. I found this c code and I am trying to to convert it to Delphi. The AxestoTranspose function in the code below takes an N dimension and is supposed to output the hilbert index. There is a similar C# conversion of it here on this site.

I think having the 3 axes as individual integers instead of encoding the axes in the format in the code might more intuitive since the code is to use with 3 dimensions only.

//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
// Functions: TransposetoAxes
//            AxestoTranspose
//
// Purpose:   Transform between Hilbert transpose and geometrical axes
//
// Example:   b=5 bits for each of n=3 coordinates
//            Hilbert transpose
//             X[0] = A D b e 3                  X[1]|  
//             X[1] = B E c 1 4    <------->         |  /X[2]
//             X[2] = C a d 2 5                axes  | /
//                    high  low                      |/______
//                                                         X[0]
//            Axes are stored conventially as b-bit integers.
//         
// History:   John Skilling  20 Apr 2001, 3 Sep 2003, 14 Oct 2003
//-----------------------------------------------------------------------------
// 
static void TransposetoAxes(
coord_t* X,            // I O  position   [n]
int      b,            // I    # bits
int      n)            // I    dimension
{
    coord_t  M, P, Q, t;
    int      i;

// Gray decode by  H ^ (H/2)
    t = X[n-1] >> 1;
    for( i = n-1; i; i-- )
        X[i] ^= X[i-1];
    X[0] ^= t;

// Undo excess work
    M = 2 << (b - 1);
    for( Q = 2; Q != M; Q <<= 1 )
    {
        P = Q - 1;
        for( i = n-1; i; i-- )
            if( X[i] & Q ) X[0] ^= P;                              // invert
            else{ t = (X[0] ^ X[i]) & P;  X[0] ^= t;  X[i] ^= t; } // exchange
        if( X[0] & Q ) X[0] ^= P;                                  // invert
    }
} 
static void AxestoTranspose(
coord_t* X,            // I O  position   [n]
int      b,            // I    # bits
int      n)            // I    dimension
{
    coord_t  P, Q, t;
    int      i;

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

// Gray encode (inverse of decode)
    for( i = 1; i < n; i++ )
        X[i] ^= X[i-1];
    t = X[n-1];
    for( i = 1; i < b; i <<= 1 )
        X[n-1] ^= X[n-1] >> i;
    t ^= X[n-1];
    for( i = n-2; i >= 0; i-- )
        X[i] ^= t;
}



via Chebli Mohamed

Aucun commentaire:

Enregistrer un commentaire