<< Computer Algorithms | Use Ctrl + or – to enlarge or reduce text size. |

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 0 | Level 1 | Level 2 | Level 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 k

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. | |

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 k

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.

D_{0} means the ray 0 0 0 → 1 1 1, that is, from the center to right rear up.

D_{4} means the ray 0 0 0 → –1 1 1, that is, from the center to left rear up.

D

Here, then, is how to proceed from level k to level k + 1 after making eight half-sized copies of level k:

_{0}Start at octant 1. | Rotate the copy –120° about D_{0}. |

_{0}Go up, to octant 2. | Rotate the copy 120° about D_{0}. |

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. |

_{0}Go up, to octant 6. | Rotate –120° about D_{4}. |

Go forward, to octant 7. | Rotate ditto. |

_{0}Go down, to octant 8. | Rotate 120° about D_{4}. |

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 0 | Level 1 | Level 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.

'--------------------------------------------------------- ' 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