Mapping Hilbert values to 3D points

Question

I have a set of Hilbert values (length from the start of the Hilbert curve to the given point).

What is the best way to convert these values to 3D points? Original Hilbert curve was not in 3D, so I guess I have to pick by myself the Hilbert curve rank I need. I do have total curve length though (that is, the maximum value in the set).

Perhaps there is an existing implementation? Some library that would allow me to work with Hilbert curve / values? Language does not matter much.

Answer 1

Not an answer about 3D conversion, but there is a nice algorithm and discussion of Hilbert values here Two-dimensional spatial hashing with space-filling curves

From MIT

4 algorithms for the n-dimensional Hilbert Space-Filling Curve

* A. R. Butz, "Alternative Algorithm for Hilbert's Space-Filling Curve",
  IEEE Trans. Comp., April, 1971, pp 424-426. [Butz 1971]

* S. W. Thomas, "hilbert.c" in the Utah Raster Toolkit circa 1993,
  http://web.mit.edu/afs/athena/contrib/urt/src/urt3.1/urt-3.1b.tar.gz

* D. Moore, Fast Hilbert Curves in C, without Recursion

* J.K.Lawder, Calculation of Mappings Between One and n-dimensional Values Using the Hilbert Space-filling Curve, [JL1_00]

Answer 2

If I get your question right you have some curve-length distance l from start point of Hilbert 3D curve and want to get coordinates corresponding to such point.

if you pre-generate the whole 3D Hilbert curve (covering unit cube) as a Polyline then all the sequenced points are at same distance between previous and next point. So you can compute your point using piecewise linear interpolation.

This is how I generate and render 2D/3D Hilbert curve in C++:

//---------------------------------------------------------------------------
#ifndef _Hilbert_vector_h
#define _Hilbert_vector_h
//---------------------------------------------------------------------------
#include "list.h"
//---------------------------------------------------------------------------
void Hilbert2D(List<double> &pnt,double x,double y,double z,double a,int n)
    {
    int i,j,m;
    double x0,y0,x1,y1,q;
    for (m=4*3,i=1,j=2;j<=n;j++,i+=i+1) m*=4; a/=i; // m = needed size of pnt[]
    pnt.num=0;
    // init generator
          pnt.add(x); pnt.add(y); pnt.add(z);
    y+=a; pnt.add(x); pnt.add(y); pnt.add(z);
    x+=a; pnt.add(x); pnt.add(y); pnt.add(z);
    y-=a; pnt.add(x); pnt.add(y); pnt.add(z);
    x0=x-0.5*a; // center of generator
    y0=y+0.5*a;
    // iterative subdivision
    for (j=2;j<=n;j++)
        {
        // mirror/rotate 2 qudrants
        x1=x0; y1=y0; m=pnt.num;
        for (i=m;i>=3;)
            {
            i--; z=pnt.dat[i]   ;
            i--; y=pnt.dat[i]-y0;
            i--; x=pnt.dat[i]-x0;
            q=x; x=+y; y=-q;    // z+
            pnt.dat[i+0]=(x1+x);
            pnt.dat[i+1]=(y1-y);
            pnt.dat[i+2]=(   z);
            }
        for (y1+=2.0*a,i=m;i>=3;)
            {
            i--; z=pnt.dat[i]   ;
            i--; y=pnt.dat[i]-y0;
            i--; x=pnt.dat[i]-x0;
            q=x; x=-y; y=+q;    // z-
            pnt.add(x1+x);
            pnt.add(y1+y);
            pnt.add(   z);
            }
        // mirror the rest
        x0+=a; y0+=a; m=pnt.num;
        for (i=m;i>=3;)
            {
            i--; z=pnt.dat[i]   ;
            i--; y=pnt.dat[i]-y0;
            i--; x=pnt.dat[i]-x0;
            pnt.add(x0-x);
            pnt.add(y0+y);
            pnt.add(   z);
            }
        a*=2.0;
        }
/*
        // rotations
        q=x; x=+y; y=-q;    // z+
        q=x; x=-y; y=+q;    // z-
*/
    }
//---------------------------------------------------------------------------
void Hilbert3D(List<double> &pnt,double x,double y,double z,double a,int n)
    {
    int i,j,m;
    double x0,y0,z0,x1,y1,z1,q;
    for (m=8*3,i=1,j=2;j<=n;j++,i+=i+1) m*=8; a/=i; // m = needed size of pnt[]
    pnt.num=0;
    // init generator
          pnt.add(x); pnt.add(y); pnt.add(z);
    z-=a; pnt.add(x); pnt.add(y); pnt.add(z);
    x+=a; pnt.add(x); pnt.add(y); pnt.add(z);
    z+=a; pnt.add(x); pnt.add(y); pnt.add(z);
    y+=a; pnt.add(x); pnt.add(y); pnt.add(z);
    z-=a; pnt.add(x); pnt.add(y); pnt.add(z);
    x-=a; pnt.add(x); pnt.add(y); pnt.add(z);
    z+=a; pnt.add(x); pnt.add(y); pnt.add(z);
    x0=x+0.5*a; // center of generator
    y0=y-0.5*a;
    z0=z-0.5*a;
    // iterative subdivision
    for (j=2;j<=n;j++)
        {
        // mirror/rotate qudrants
        x1=x0; y1=y0; z1=z0; m=pnt.num;
        for (i=m;i>=3;)
            {
            i--; z=pnt.dat[i]-z0;
            i--; y=pnt.dat[i]-y0;
            i--; x=pnt.dat[i]-x0;
            q=y; y=-z; z=+q;    // x-
            pnt.dat[i+0]=(x1+x);
            pnt.dat[i+1]=(y1+y);
            pnt.dat[i+2]=(z1-z);
            }
        for (z1-=2.0*a,i=m;i>=3;)
            {
            i--; z=pnt.dat[i]-z0;
            i--; y=pnt.dat[i]-y0;
            i--; x=pnt.dat[i]-x0;
            q=z; z=+x; x=-q;    // y+
            q=y; y=+z; z=-q;    // x+
            pnt.add(x1-x);
            pnt.add(y1+y);
            pnt.add(z1+z);
            }
        for (x1+=2.0*a,i=m;i>=3;)
            {
            i--; z=pnt.dat[i]-z0;
            i--; y=pnt.dat[i]-y0;
            i--; x=pnt.dat[i]-x0;
            q=y; y=+z; z=-q;    // x+
            pnt.add(x1+x);
            pnt.add(y1+y);
            pnt.add(z1+z);
            }
        for (z1+=2.0*a,i=m;i>=3;)
            {
            i--; z=pnt.dat[i]-z0;
            i--; y=pnt.dat[i]-y0;
            i--; x=pnt.dat[i]-x0;
            q=z; z=+x; x=-q;    // y+
            pnt.add(x1-x);
            pnt.add(y1-y);
            pnt.add(z1+z);
            }
        // mirror octants
        x0+=a; y0+=a; z0-=a; m=pnt.num;
        for (i=m;i>=3;)
            {
            i--; z=pnt.dat[i]-z0;
            i--; y=pnt.dat[i]-y0;
            i--; x=pnt.dat[i]-x0;
            pnt.add(x0+x);
            pnt.add(y0-y);
            pnt.add(z0+z);
            }
        a*=2.0;
        }
/*
        // rotations
        q=z; z=+x; x=-q;    // y+
        q=z; z=-x; x=+q;    // y-
        q=y; y=+z; z=-q;    // x+
        q=y; y=-z; z=+q;    // x-
        q=x; x=+y; y=-q;    // z+
        q=x; x=-y; y=+q;    // z-
*/
    }
//---------------------------------------------------------------------------
void pnt_draw2(List<double> &pnt)   // piecewise linear
    {
    int i;
    glBegin(GL_LINE_STRIP);
    for (i=0;i<pnt.num;i+=3) glVertex3dv(pnt.dat+i);
    glEnd();
    }
//---------------------------------------------------------------------------
void pnt_draw4(List<double> &pnt)   // piecewise cubic
    {
    int i,j;
    double  d1,d2,t,tt,ttt,*p0,*p1,*p2,*p3,a0[3],a1[3],a2[3],a3[3],p[3];
    glBegin(GL_LINE_STRIP);
    for (i=-3;i<pnt.num;i+=3)
        {
        j=i-3; if (j>pnt.num-3) j=pnt.num-3; if (j<0) j=0; p0=pnt.dat+j;
        j=i  ; if (j>pnt.num-3) j=pnt.num-3; if (j<0) j=0; p1=pnt.dat+j;
        j=i+3; if (j>pnt.num-3) j=pnt.num-3; if (j<0) j=0; p2=pnt.dat+j;
        j=i+6; if (j>pnt.num-3) j=pnt.num-3; if (j<0) j=0; p3=pnt.dat+j;
        for (j=0;j<3;j++)
            {
            d1=0.5*(p2[j]-p0[j]);
            d2=0.5*(p3[j]-p1[j]);
            a0[j]=p1[j];
            a1[j]=d1;
            a2[j]=(3.0*(p2[j]-p1[j]))-(2.0*d1)-d2;
            a3[j]=d1+d2+(2.0*(-p2[j]+p1[j]));
            }
        for (t=0.0;t<=1.0;t+=0.1)   // single curve patch/segment
            {
            tt=t*t;
            ttt=tt*t;
            for (j=0;j<3;j++) p[j]=a0[j]+(a1[j]*t)+(a2[j]*tt)+(a3[j]*ttt);
            glVertex3dv(p);
            }
        }
    glEnd();
    }
//---------------------------------------------------------------------------
#endif
//---------------------------------------------------------------------------

I used mine dynamic list template so:

List<double> xxx; is the same as double xxx[];
xxx.add(5); adds 5 to end of the list
xxx[7] access array element (safe)
xxx.dat[7] access array element (unsafe but fast direct access)
xxx.num is the actual used size of the array
xxx.reset() clears the array and set xxx.num=0
xxx.allocate(100) preallocate space for 100 items

But you can use dynamic or even static 1D array instead as the number of points of Hilbert curve is easily computable (m at the start of each Hilbert function).

Usage is simple just do something like this:

List<double> pnt;
Hilbert3D(pnt,-0.8,-0.8,+0.8,1.6,n); 

Where n is number of iterations and pnt is linear list of (x,y,z) coordinates for each point (3 numbers per point). the start position and initial size is set to cover cube centered around (0,0,0) with half size 0.8 <-0.8,+0.8>.

Now just compute unit length between points, index of closest Hilbert curve point to the left and parameter (distance to it) from that just linearly interpolate. Here C++ example:

if (pnt.num>=6)
    {
    int i;
    double x,y,z,t,l,dl;
    dl=sqrt(                                                    // base distance between points
            ((pnt[0]-pnt[3])*(pnt[0]-pnt[3]))
           +((pnt[1]-pnt[4])*(pnt[1]-pnt[4]))
           +((pnt[2]-pnt[5])*(pnt[2]-pnt[5]))
            );
    l=double(Form1->sb_t->Position)/double(Form1->sb_t->Max);   // <0,1>
    l*=dl*double((pnt.num/3)-1);                                // <0,Hilbert_curve_lenght>
    i=floor(l/dl); t=(l-(double(i)*dl))/dl; i*=3;               // index in pnt[] and single line segment paramerer
    x=pnt[i+0]+(pnt[i+3]-pnt[i+0])*t;                           // linear interpolation
    y=pnt[i+1]+(pnt[i+4]-pnt[i+1])*t;
    z=pnt[i+2]+(pnt[i+5]-pnt[i+2])*t;
    glColor3f(0.0,1.0,0.0); t=0.05;                             // render of marker
    glBegin(GL_LINES);
    glVertex3d(x-t,y-t,z); glVertex3d(x+t,y+t,z);
    glVertex3d(x+t,y-t,z); glVertex3d(x-t,y+t,z);
    glVertex3d(x,y-t,z-t); glVertex3d(x,y+t,z+t);
    glVertex3d(x,y-t,z+t); glVertex3d(x,y+t,z-t);
    glVertex3d(x-t,y,z-t); glVertex3d(x+t,y,z+t);
    glVertex3d(x+t,y,z-t); glVertex3d(x-t,y,z+t);
    glEnd();
    }

2D preview:

3D preview: