Algorithm for generating a 3D Hilbert space-filling curve in Python

Question

I'd like to map points in a RGB color cube to a one-dimensional list in Python, in a way that makes the list of colors look nice and continuous.

I believe using a 3D Hilbert space-filling curve would be a good way to do this, but I've searched and haven't found very helpful resources for this problem. Wikipedia in particular only provides example code for generating 2D curves.

Answer 1

This paper seems to have quite a discussion: An inventory of three-dimensional Hilbert space-filling curves.

Quoting from the abstract:

Hilbert's two-dimensional space-filling curve is appreciated for its good locality properties for many applications. However, it is not clear what is the best way to generalize this curve to filling higher-dimensional spaces. We argue that the properties that make Hilbert's curve unique in two dimensions, are shared by 10694807 structurally different space-filling curves in three dimensions.

Answer 2

I came across your question while trying to do the same thing in javascript. I figured it out on my own. Here is a recursive function that breaks a cube in 8 parts and rotates each part so that it traverses a hilbert curve in order. The arguments represent the size:s, location:xyz, and 3 vectors for the rotated axes of the cube. The example call uses a 256^3 cube and assumes red,green,blue arrays have length 256^3.

It should be easy to adapt this code to python or other procedural languages.

Adapted from pictures here: http://www.math.uwaterloo.ca/~wgilbert/Research/HilbertCurve/HilbertCurve.html

    function hilbertC(s, x, y, z, dx, dy, dz, dx2, dy2, dz2, dx3, dy3, dz3)
    {
        if(s==1)
        {
            red[m] = x;
            green[m] = y;
            blue[m] = z;
            m++;
        }
        else
        {
            s/=2;
            if(dx<0) x-=s*dx;
            if(dy<0) y-=s*dy;
            if(dz<0) z-=s*dz;
            if(dx2<0) x-=s*dx2;
            if(dy2<0) y-=s*dy2;
            if(dz2<0) z-=s*dz2;
            if(dx3<0) x-=s*dx3;
            if(dy3<0) y-=s*dy3;
            if(dz3<0) z-=s*dz3;
            hilbertC(s, x, y, z, dx2, dy2, dz2, dx3, dy3, dz3, dx, dy, dz);
            hilbertC(s, x+s*dx, y+s*dy, z+s*dz, dx3, dy3, dz3, dx, dy, dz, dx2, dy2, dz2);
            hilbertC(s, x+s*dx+s*dx2, y+s*dy+s*dy2, z+s*dz+s*dz2, dx3, dy3, dz3, dx, dy, dz, dx2, dy2, dz2);
            hilbertC(s, x+s*dx2, y+s*dy2, z+s*dz2, -dx, -dy, -dz, -dx2, -dy2, -dz2, dx3, dy3, dz3);
            hilbertC(s, x+s*dx2+s*dx3, y+s*dy2+s*dy3, z+s*dz2+s*dz3, -dx, -dy, -dz, -dx2, -dy2, -dz2, dx3, dy3, dz3);
            hilbertC(s, x+s*dx+s*dx2+s*dx3, y+s*dy+s*dy2+s*dy3, z+s*dz+s*dz2+s*dz3, -dx3, -dy3, -dz3, dx, dy, dz, -dx2, -dy2, -dz2);
            hilbertC(s, x+s*dx+s*dx3, y+s*dy+s*dy3, z+s*dz+s*dz3, -dx3, -dy3, -dz3, dx, dy, dz, -dx2, -dy2, -dz2);
            hilbertC(s, x+s*dx3, y+s*dy3, z+s*dz3, dx2, dy2, dz2, -dx3, -dy3, -dz3, -dx, -dy, -dz);
        }
    }
    m=0;
    hilbertC(256,0,0,0,1,0,0,0,1,0,0,0,1);

Answer 3

What's frequently used in engineering practice is not strictly Hilbert (Peano) curves - it's Morton code.

https://en.wikipedia.org/wiki/Z-order_curve

Much easier to compute.