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);