<< Computer AlgorithmsUse  Ctrl +  or    to enlarge or reduce text size.

The Hilbert Curve in Three Dimensions

The Hilbert curve in two dimensions is the limit of a certain sequence of paths inside a square. As the sequence progresses the path passes closer and closer to every point of the square. In the limit it passes through every point of it, and it is called a “space filling” curve. The limit path is rather academic and we will leave it to those interested in such things. Here we focus instead on the sequence of paths leading up to it.

Consider any one of the paths in the sequence. It passes through each point of a lattice of points inside the square (see the diagram below). That by itself is not remarkable. A raster scan, passing through one horizontal row of points after the other, does the same. What makes the Hilbert path interesting is that two points near one another as measured along the meandering path are on average even nearer measured “as the crow flies,” that is, along the shortest line between them in the plane. (The converse isn’t true. Some points are near one another in the plane but far from one another along the curve between them.) Whereas with the raster scan path, two points near one another measured along the path are usually the same distance measured in the plane.

The sequence of Hilbert paths is defined inductively, that is, a path is determined by the one before, as follows.

Imagine the square divided into four quadrants numbered 1 to 4 starting at the lower left. As you go from 1 to 2 to 3 to 4, only one direction (either up-down or left-right) changes at a time.
                                                           
Level  0Level  1Level  2Level  3

Consider a sequence of such squares numbered  0, 1, 2 ... , one for each curve.  Level 0 starts things off; it is merely a point in the center of the square. Now suppose we have the  kth curve in the  kth square. We will say the curve is at level  k.  We can draw the curve at level  k + 1  as follows.  Make four half-sized copies of level  k  and place them in the four quadrants of the square for level  k + 1,  oriented a certain way.  In detail: 
      Start at the lower left, quadrant 1. Flip the copy about the ╱ diagonal.
      Go up, to quadrant 2. Leave the copy as is.
            Go right, to quadrant 3. Leave the copy as is.
            Go down, to quadrant 4. Flip the copy about the ╲ diagonal.
And as you do this connect the end of one copy to the beginning of the next. The purpose of changing the orientation is to make this connection as short as possible. These connections are always parallel to the steps in level 1.

So much for the Hilbert curve in two dimensions. We took care to describe the inductive procedure in a way that can be generalized to higher dimensions. Now we construct a similar sequence of paths in three dimensions, that is do for the cube what we just did for the square. As before we start with one point and proceed inductively.

Imagine the cube divided into eight octants numbered 1 to 8 starting at the lower left forward corner. As you step from 1 to 2 to ... to 8, only one coordinate (left-right X, forward-rearward Y, or up-down Z) changes each time.

Consider a sequence of such cubes numbered  0, 1, 2 ... , one for each curve.  As we said, level 0 starts things off; a point in the center of the cube. Now suppose we have the kth curve in the kth cube. We can determine the curve in the level  k + 1 cube as follows.  Make eight half-sized copies of level k and place them in the eight quadrants of the cube for level  k + 1, oriented to allow their easy connection. This last, finding suitable orientations, is the only difficult part. (The solution is not unique.)

We can specify the orientations by a rotation about an axis. We shall need several axes. To avoid describing them repeatedly use the following nomenclature: 
X  means the ray  0 0 0 →   1 0 0,  that is, from the center to the right.
D0 means the ray  0 0 0 →   1 1 1,  that is, from the center to right rear up.
D4 means the ray  0 0 0 → –1 1 1,  that is, from the center to left rear up.

Here, then, is how to proceed from level  k  to level  k + 1  after making eight half-sized copies of level k:
0Start at octant 1. Rotate the copy  –120° about D0.
0Go up, to octant 2. Rotate the copy   120° about D0.
Go rearward, to octant 3. Rotate ditto.
Go down, to octant 4. Rotate the copy 180° about X.
Go right, to octant 5. Rotate ditto.
0Go up, to octant 6. Rotate  –120° about D4.
      Go forward, to octant 7. Rotate ditto.
0Go down, to octant 8. Rotate   120° about D4.

As you do this connect the end of one copy to the beginning of the next. The purpose of the rotations is to make this connection as short as possible. These connections are always parallel to the steps in level 1.
                                       
Level  0Level  1Level  2

Not drawn is the enclosing cube analogous to the enclosing square drawn for the two dimensional case. At each level beyond zero the edge of that cube is a step length longer than the edge of the smallest cube the curve fits in. In other words, each of its sides is half a step length out from the apparent smallest enclosure. In the diagram for level 2 above, the shaded cubical areas are half the size of an octant of the enclosing cube the curve, in the limit as  level → ∞,  fills.

The program below is written for PBCC. The parts that compute the coordinates of the curve can be used almost as is for Visual Basic, FORTRAN and any other BASIC-like compiler; the graphical display would require extra work.

This particular program does not use recursion in the programming sense, that is, a subroutine calling itself. Computing the next level requires all the points of the previous level. It proceeds “bottom up” rather than “top down.” Recursion can be used to compute the Hilbert curve in two dimensions but I don’t know about three.

Hilbert-3D.exe
'---------------------------------------------------------
' Only the subroutine PBMain is used to
' compute the 3D coordinates of the curve.
' The other subroutines help display it.
'---------------------------------------------------------

Global x(), y(), z() As Single      'dimensioned in PBMain
Global Xbuffer(), Ybuffer(), Zbuffer() As Single

Global t11, t12      As Single      'tilt matrix used by DrawMain (t13 = 0)
Global t21, t22, t23 As Single
Global t31, t32, t33 As Single

Global xtilt, ztilt As Single       'angles for tiltmatrix in radians
Global xtiltin, ztiltin As Long     'in degrees
Global radianf As Single            'conversion factor to go from degrees to radians

Global level As Long                'iteration level
Global n As Long                    'number of points of the curve in Main
Global stepsize As Single           'step size
Global pq As Single                 'pi / 4

%stepsize0 = 1000                   'initial stepsize
%maxlevel = 5                       'maximum level
'---------------------------------------------------------

'Furnishes t11, ... t33 used by DrawMain.  Depends on the angles xtilt, ztilt.
Sub GetTiltMatrix
 Local cxt, sxt As Single
 Local sr, v As Single
 Local e11, e12, e13, e21, e22, e23, e31, e32, e33 As Single

 cxt = Cos(xtilt) :sxt = Sin(xtilt)
 sr  = Sin(ztilt)

 e11 = Cos(ztilt)  :e12 = -cxt*sr          :e13 = -sxt*sr
 v = 1 - e11
 e21 = -e12        :e22 = e11 + sxt*sxt*v  :e23 = -sxt*cxt*v
 e31 = -e13        :e32 = e23              :e33 = e11 + cxt*cxt*v

 t11 =  e11
 t12 =  e12*cxt + e13*sxt
 't13 = 0, not used
 t21 =  e21
 t22 =  e22*cxt + e23*sxt
 t23 = -e22*sxt + e23*cxt
 t31 =  e31
 t32 =  e32*cxt + e33*sxt
 t33 = -e32*sxt + e33*cxt
End Sub
'---------------------------------------------------------

Sub NormalizeAngle(angle As Long)    'put angle back in the range -180 to 180
 If     angle >   180 Then
   angle = angle - 360
 ElseIf angle <= -180 Then
   angle = angle + 360
 End If
End Sub

Sub XtiltInc(ByVal delta As Long)
 xtiltin = xtiltin + delta
 NormalizeAngle xtiltin
 xtilt = xtiltin * radianf
 GetTiltMatrix
End Sub

Sub ZtiltInc(ByVal delta As Long)
 ztiltin = ztiltin + delta
 NormalizeAngle ztiltin
 ztilt = ztiltin * radianf
 GetTiltMatrix
End Sub
'---------------------------------------------------------

%viewerdistance = %stepsize0*5         'viewer distance from screen for perspective projection
%colour = %Red
%delay = 200                   'millisecond wait between drawing steps

'Display  x(i),y(i),z(i)  i = 1 to n  moved by the tilt matrix t11, t12, etc.
'Besides the tilt matrix it uses the globals: level, n.
Sub DrawMain
 Local i As Long
 Local f As Single
 Local xt, yt, zt As Single
 Local xx, yy, zz As Single
 Local s As Single

 s = 2^(level - 2)

 'first point
 xt = x(1)  :yt = y(1)  :zt = z(1)
 xx = t11*xt + t12*yt
 yy = t21*xt + t22*yt + t23*zt
 zz = t31*xt + t32*yt + t33*zt
 f = s * %viewerdistance / (%viewerdistance + yy)
 Graphic Ellipse (f*xx-5, f*zz-5) - (f*xx+6, f*zz+6), %colour, %colour
 Graphic Set Pos (f*xx, f*zz)

 'succeeding points
 For i = 2 To n
   xt = x(i)  :yt = y(i)  :zt = z(i)
   xx = t11*xt + t12*yt
   yy = t21*xt + t22*yt + t23*zt
   zz = t31*xt + t32*yt + t33*zt
   f = s * %viewerdistance / (%viewerdistance + yy)
   Sleep 1
   Graphic Line -(f*xx, f*zz), %colour
   Sleep %delay / level
   If InStat Then WaitKey$ :Exit Sub
   Graphic ReDraw
 Next
 Graphic Ellipse (f*xx-5, f*zz-5) - (f*xx+6, f*zz+6), %colour, %colour
 Graphic ReDraw
End Sub
'---------------------------------------------------------

%leftkey  = 75     'extended key codes for arrow keys
%upkey    = 72
%rightkey = 77
%downkey  = 80

' This breaks the curve up into line segments, sorts them on depth,
' then draws the furthest first so that what's in front covers what's
' in back.  Besides the tilt matrix it uses the globals: level, n.
Function DrawMainFinal As Long                 'return 1 if user presses Esc to quit
 Local i, k As Long
 Local f As Single
 Local xt, yt, zt As Single
 Local s As Single
 Dim xx1(n) As Long, zz1(n) As Long
 Dim xx2(n) As Long, zz2(n) As Long
 Dim yy(n) As Long
 Local yyA, yyB As Single
 Dim index(n) As Long
 Local x1, z1, x2, z2 As Long
 Local vx, wx, vz, wz As Single
 Local vlen As Single
 Local xdot1, zdot1, xdotn, zdotn As Single    'beginning and end dots
 Local wk As String
 Local halocolor As Long

 s = 2^(level - 2)

 Do

   For i = 1 To n                                 'for tag sort
    index(i) = i
   Next

   xt = x(1)  :yt = y(1)  :zt = z(1)
   yyA = t21*xt + t22*yt + t23*zt
   f = s * %viewerdistance / (%viewerdistance + yyA)
   xdot1 = f * (t11*xt + t12*yt)           :xx1(1) = xdot1
   zdot1 = f * (t31*xt + t32*yt + t33*zt)  :zz1(1) = zdot1

   For i = 1 To n - 1
     xt = x(i)  :yt = y(i)  :zt = z(i)
     yyB = t21*xt + t22*yt + t23*zt
     f = s * %viewerdistance / (%viewerdistance + yyB)
     xx2(i) = f * (t11*xt + t12*yt)
     zz2(i) = f * (t31*xt + t32*yt + t33*zt)
     xx1(i+1) = xx2(i)
     zz1(i+1) = zz2(i)
     yy(i) = yyA + yyB
     yyB = yyA
     If InStat Then
      wk = WaitKey$
      If wk = $Cr Or wk = $Esc Then Exit For
     End If
   Next
   xt = x(i)  :yt = y(i)  :zt = z(i)              'here i = n
   yyB = t21*xt + t22*yt + t23*zt
   yy(i) = yyA + yyB
   f = s * %viewerdistance / (%viewerdistance + yyB)
   xdotn = f * (t11*xt + t12*yt)           :xx2(i) = xdotn
   zdotn = f * (t31*xt + t32*yt + t33*zt)  :zz2(i) = zdotn

   Array Sort yy(1), TagArray index(), Descend    'sort back to front

   halocolor = RGB(212,208,200)                   'window background color

   For i = 1 To n
     k = index(i)

     x1 = xx1(k)  :z1 = zz1(k)
     x2 = xx2(k)  :z2 = zz2(k)

     wx = x1 - x2  :wz = z1 - z2                  'backlight
     vlen = Sqr(wx*wx + wz*wz)
     If vlen Then
       vx = -wz*3/vlen                            'perpendicular to line
       vz =  wx*3/vlen
       wx = -vz*2.35                              'along line
       wz =  vx*2.35
       Graphic Line(x1+vx+wx,z1+vz+wz)-(x2+vx-wx,z2+vz-wz),halocolor
       Graphic Line(x1-vx+wx,z1-vz+wz)-(x2-vx-wx,z2-vz-wz),halocolor

     End If

     Graphic Line (x1, z1)-(x2, z2), %colour
   Next

   Graphic Ellipse (xdot1-5, zdot1-5) - (xdot1+6, zdot1+6), %colour, %colour
   Graphic Ellipse (xdotn-5, zdotn-5) - (xdotn+6, zdotn+6), %colour, %colour
   Graphic ReDraw

   Do                                                  'wait for user to press a key
     wk = WaitKey$
     Input Flush
     If     Len(wk) = 1 Then
       If wk = $Cr Then Exit Function                  'Enter: continue on to next level
       If wk = $Esc Then Function = 1 :Exit Function   'Esc:   quit program
     ElseIf Len(wk) = 2 Then
       Select Case Asc(wk,2)
       Case %leftkey  :ZtiltInc -5                     'arrow key: change view of the current level
       Case %upkey    :XtiltInc -5
       Case %rightkey :ZtiltInc  5
       Case %downkey  :XtiltInc  5
       Case Else      :Iterate Do
       End Select
       Graphic Clear
       Exit Do
     End If
   Loop

 Loop

End Function
'---------------------------------------------------------

%ds = 750                          'graphic window size
$vert = Chr$(179) :$hori = Chr$(196)
$UL = Chr$(218) :$UR = Chr$(191)
$LL = Chr$(192) :$LR = Chr$(217)
$US = Chr$(194) :$LS = Chr$(193)

Function PBMain
 Local winC As Dword               'window handle
 Local i, j As Long                'index of point in Buffer, in Main
 Local nmax As Long                'number of points of last level
 Local lastlevel As Long           'last iteration level
 Local xx, yy, zz As Long          'a point in Buffer

 Cursor Off

 Graphic Window "Hilbert Curve in 3D",  5,0,  %ds,%ds To winC
 Graphic Attach winC,0, ReDraw
 Graphic Width 3             'draw thick lines
 Console Set Loc %ds,200     'move console out of the way
 Console Set Focus           'give console the focus so it gets user's keystrokes

 Graphic Scale (-%ds,%ds)-(%ds,-%ds)    'origin at center, z increases upwards

 pq      = Atn(1)            'pi / 4 radians or 45 degrees
 radianf = pq/45             'factor to convert from degrees to radians
 ' 3D coordinate system
 '     z up
 '     |   y away
 '     |  /
 '     | /
 '     |/______x to right

 xtiltin = 30  :ztiltin = -30              'initial tilt for display
 xtilt = xtiltin*radianf  :ztilt = ztiltin*radianf
 GetTiltMatrix

 ' There are two cubical regions:  "Main" where the
 ' curve is drawn for viewing, and "Buffer" which
 ' contains a half-sized copy of the curve in Main.

 ' We will use a "bottom up" approach and determine
 ' a level from the previous one.

 ' Assume we have the curve in Main at one level.
 ' Move it to Buffer, quarter-sized, and clear Main.
 ' Bisect Main three ways into eight smaller cubes
 ' ordered as follows (each corner in the diagram
 ' is the center of one of the small cubes):
 '
 '         3           6              positions:
 '        /|          /|              (1)   -1 -1 -1
 '       / |         / |              (2)   -1,-1, 1
 '      /  |        /  |              (3)   -1, 1, 1
 '   2 |   |     7 |   |              (4)   -1, 1,-1
 '     |   |       |   |              (5)    1, 1,-1
 '     |   |_______|___|              (6)    1, 1, 1
 '     |  4        |   5              (7)    1,-1, 1
 '     |           |                  (8)    1,-1,-1
 '   1 o         8 o
 '   Begin        End
 '
 ' Copy the contents of Buffer to each of the eight
 ' cubes of Main, oriented various ways, connecting
 ' the end of one copy to the beginning the next.
 '
 ' The copies, 1 through 8, are oriented by rotating
 ' a copy of Buffer about its center.  Numbering the
 ' cubes in order, the rotations are:
 '  (1) rotate about ( 1, 1, 1) by -120 degrees
 '  (2) rotate about ( 1, 1, 1) by +120 degrees
 '  (3) ditto
 '  (4) rotate about ( 1, 0, 0) by  180 degrees
 '  (5) ditto
 '  (6) rotate about (-1, 1, 1) by -120 degrees
 '  (7) ditto
 '  (8) rotate about (-1, 1, 1) by +120 degrees
 ' All these rotations just permute the coordinates.

 ' Start Main at the zero'th level, an isolated point.
 level = 0
 n = 1                              'number of points
 Dim x(1), y(1), z(1)
 x(1) = 0 :y(1) = 0 :z(1) = 0       'center of Main
 stepsize = %stepsize0

 Locate 2,2  :Print "If this small window has the focus then these keys will work:"
 Locate 4,2  :Print $UL;String$(36,$hori);$US+String$(25,$hori);$US;String$(13,$hori);$UR
 Locate  ,2  :Print $vert;"                                    ";$vert;"                         ";$vert;"             ";$vert
 Locate  ,2  :Print $vert;" Enter - quick finish or next level ";$vert;" Arrow key - change view ";$vert;" Esc - quit  ";$vert
 Locate  ,2  :Print $vert;"                                    ";$vert;"                         ";$vert;"             ";$vert
 Locate  ,2  :Print $LL;String$(36,$hori);$LS+String$(25,$hori);$LS;String$(13,$hori);$LR
 Locate 10,2 :Print "Drawing level:  O   (" + Format$(%maxlevel) + " maximum)";
 Locate 24,28 :Print " Algorithms.ARIwatch.com";

 Do

   If DrawMainFinal Then Exit Do    'draw line segments back to front with backlighting

   Incr level :If level > %maxlevel Then Exit Do
   Locate 10,17 :Print level;

   nmax = n*8                       ' = 8^level, the number of points, one less than number of steps
   ReDim Preserve x(nmax), y(nmax), z(nmax)                     'curve points in Main, global
   ReDim Preserve Xbuffer(nmax), Ybuffer(nmax), Zbuffer(nmax)

   Graphic Clear
   stepsize = stepsize / 2

   ' Copy the Main curve to Buffer, shrunken.  The factor
   ' is 1/4 instead of 1/2 to put space between the cubes.
   ' To better see how the copies are connected together,
   ' divide by a number greater than 4, say 6.
   For i = 1 To n
     Xbuffer(i) = x(i) / 4  :Ybuffer(i) = y(i) / 4  :Zbuffer(i) = z(i) / 4
   Next

   ' Fill Main with eight copies of Buffer, rotated per above,
   ' each beginning where the last ended.
   '
   ' For each point in the Buffer cube (index i),
   ' rotate it, apply the translation,
   ' then set the Main point (index j) to it.

   j = 1                                       'each For: same i's, new j's
   For i = 1 To n              'cube (1)
     x(j) = Ybuffer(i) - stepsize              'rotate -120° about (1, 1, 1)
     y(j) = Zbuffer(i) - stepsize              'and translate (-1,-1,-1) * stepsize
     z(j) = Xbuffer(i) - stepsize
     Incr j
   Next
   For i = 1 To n              'cube (2)
     x(j) = Zbuffer(i) - stepsize              'rotate 120° about (1, 1, 1)
     y(j) = Xbuffer(i) - stepsize              'and translate (-1,-1, 1) * stepsize
     z(j) = Ybuffer(i) + stepsize
     Incr j
   Next
   For i = 1 To n              'cube (3)
     x(j) = Zbuffer(i) - stepsize              'rotate same way as last time
     y(j) = Xbuffer(i) + stepsize              'and translate (-1, 1, 1) * stepsize
     z(j) = Ybuffer(i) + stepsize
     Incr j
   Next
   For i = 1 To n              'cube (4)
     x(j) =  Xbuffer(i) - stepsize             'rotate 180° about (1, 0, 0)
     y(j) = -Ybuffer(i) + stepsize             'and translate (-1, 1,-1) * stepsize
     z(j) = -Zbuffer(i) - stepsize
     Incr j
   Next
   For i = 1 To n              'cube (5)
     x(j) =  Xbuffer(i) + stepsize             'rotate same way as last time
     y(j) = -Ybuffer(i) + stepsize             'and translate ( 1, 1,-1) * stepsize
     z(j) = -Zbuffer(i) - stepsize
     Incr j
   Next
   For i = 1 To n              'cube (6)
     x(j) = -Zbuffer(i) + stepsize             'rotate -120° about (-1, 1, 1)
     y(j) = -Xbuffer(i) + stepsize             'and translate ( 1, 1, 1) * stepsize
     z(j) =  Ybuffer(i) + stepsize
     Incr j
   Next
   For i = 1 To n              'cube (7)
     x(j) = -Zbuffer(i) + stepsize             'rotate same way
     y(j) = -Xbuffer(i) - stepsize             'and translate ( 1, -1, 1) * stepsize
     z(j) =  Ybuffer(i) + stepsize
     Incr j
   Next
   For i = 1 To n              'cube (8)
     x(j) = -Ybuffer(i) + stepsize             'rotate 120° about (-1, 1, 1)
     y(j) =  Zbuffer(i) - stepsize             'and translate ( 1,-1,-1) * stepsize
     z(j) = -Xbuffer(i) - stepsize
     Incr j
   Next

   n = j - 1          'n = n*8 = 8^level

   DrawMain           'draw from point to point in order
   Graphic Clear      'only to clear missed points caused by flaky graphic cards
 Loop                 'DrawMainFinal will be executed at start of the next do cycle
End Function