A compact Hilbert curve routine

January 30, 2017

Here is a relatively compact Mathematica routine for generating the b-th iterate of a Hilbert curve in [0,1]^n. The algorithm is due to Skilling:

HilbertCurve =
     Compile[{{x, _Integer}, {b, _Integer}, {n, _Integer}},
             Module[{t = BitXor[x, Quotient[x, 2]], p = 2, k, q, r, xx},
                    xx = Total[Table[BitAnd[Quotient[t, 2^(r (n - 1) + k)], 2^r],
                                     {r, b - 1, 0, -1}, {k, n - 1, 0, -1}]];
                    Do[q = p - 1;
                       Do[t = xx[[k]];
                          If[BitAnd[t, p] != 0,
                             xx[[1]] = xx[[1]] ~BitXor~ q,
                             t = BitAnd[t ~BitXor~ xx[[1]], q];
                             xx[[{1, k}]] = BitXor[xx[[{1, k}]], t]],
                          {k, n, 2, -1}];
                       If[BitAnd[xx[[1]], p] != 0,
                          xx[[1]] = xx[[1]] ~BitXor~ q];
                       p *= 2, {r, b - 1}];
                    2 xx/(p - 1) - 1], RuntimeAttributes -> {Listable}]

I’ve seen a variety of routines for generating the Hilbert curve before, but none were as compact or as general as this one. If you want a version that generates exact rational coordinates, it is easy to take out the requisite parts from the Compile[] function.

Here are the 2D and 3D versions of the curve:

With[{p = 6}, 
     Graphics[{ColorData[97, 1],
               Line[HilbertCurve[Range[0, 2^(2 p) - 1], p, 2]]}, Frame -> True]]

h2d

With[{p = 3}, 
     Graphics3D[Tube[HilbertCurve[Range[0, 2^(3 p) - 1], p, 3]], Axes -> True]]

h3d

Here’s a perspective projection of a four-dimensional Hilbert curve:

With[{p = 3, k = 2, f = 4, d = 3/2, r = 1/8}, 
     Graphics3D[{Directive[EdgeForm[], CapForm[None], GrayLevel[1/5], 
                           Glow[ColorData["Legacy", "SteelBlue"]],
                           Specularity[ColorData["Legacy", "Gold"], 16]], 
                 Tube[(f Delete[#, k])/(d - Extract[#, k]) & /@ 
                      HilbertCurve[Range[0, 2^(4 p) - 1], p, 4], r]}, 
                Boxed -> False, Lighting -> "Neutral"]]

h4d

(This visualization is admittedly a bit messy-looking; if you have suggestions for making a nice display of the four-dimensional version, tell me about it in the comments!)