Chunlei Chen: Hilbert Curve Concepts & Implementation

Hilbert Curve
Concepts & Implementation

Introduction

The web notes by Andrew Cumming of Napier University Edinburgh, provide a very good introduction to the Hilbert curve. His pseudo code, slightly modified, is shown in listing 1.

Listing 1 (from Andrew Cumming)

procedure hilbert(x, y, xi, xj, yi, yj, n)
/* x and y are the coordinates of the bottom left corner */
/* xi & xj are the i & j components of the unit x vector of the frame */
/* similarly yi and yj */
if (n <= 0) then
   LineTo(x + (xi + yi)/2, y + (xj + yj)/2);
else
   {
   hilbert(x,           y,           yi/2, yj/2,  xi/2,  xj/2, n-1);
   hilbert(x+xi/2,      y+xj/2 ,     xi/2, xj/2,  yi/2,  yj/2, n-1);
   hilbert(x+xi/2+yi/2, y+xj/2+yj/2, xi/2, xj/2,  yi/2,  yj/2, n-1);
   hilbert(x+xi/2+yi,   y+xj/2+yj,  -yi/2,-yj/2, -xi/2, -xj/2, n-1);
   }
end procedure;

To fully understand how Andrews recursive function operates is not easy and in order to do so it is necessary to break-down his Hilbert procedure to that the role that each "part" plays in the construction of the curve. We will begin by looking at the inputs of the Hilbert function.

Inputs

In addition to the counter (n) that tracks the level from which the function begins to recursively call itself, there are six other inputs to Andrew's procedure. The first two define the x and y coordinates of an input point that will be used in the calculation of a point on the Hilbert curve. The next 4 values define two vectors. Before we see how they are used we must understand what is meant by the inputs labelled xi, xj, yi and yj. Hopefully, figure 1 will make them clearer.

The large black dot, labelled P, in figure 1 is defined by its x and y coordinates - shown as small dots on the two axes. We can see, from the two arrows, that each dot marks the head of a vector. We might call these the x-vector and the y-vector.

Figure 1

Although these vectors represent the x and y coordinates of P they, like any vector, have their own coordinates. For example, the coordinates of the x-vector might be 1.0, 0.0, while the coordinates of the y-vector might be 0.0, 0.5. To prevent the confusion of labelling the coordinates of these vectors as "x" and "y" we, instead, refer to them by the labels "i" and "j". For example,

x-vector              y-vector
    xi = 1.0 and          yi = 0.0 and
    xj = 0.0              yj = 0.5

We can make use of these vectors to move (ie. translate) a copy of point P. For example, in Andrew's procedure we note this statement.

LineTo(x + (xi + yi)/2, y + (xj + yj)/2);

Rewritten for point P we have,

Px = Px + (xi + yi)/2
       = 1 + (1 + 0)/2
       = 1.5

    Py = Py + (xj + yj)/2
       = 0.5 + (0 + 0.5)/2
       = 0.75

The effect of using the i's and j's on point P can be seen in figure 2.

Figure 2

Point P could also have been translated by directly changing its x and y values. However, using the "i" and "j" coordinates (x-vector and y-vectors) moves the point in a structured and proportional fashion. In the context of their use in Andrew's procedure, the "i" and "j" components provide, what are in effect, additional "handles" to control the placement of points along the Hilbert curve.

Modifying the Inputs

Looking at the hilited code shown in listing 2 we see the Hilbert procedure calls itself four times. On each call it passes slightly modified versions of its original input values to itself.

Listing 2

procedure hilbert(x, y, xi, xj, yi, yj, n)
if (n <= 0) then
   LineTo(x + (xi + yi)/2, y + (xj + yj)/2);
else
   {
   hilbert(x,           y,           yi/2, yj/2,  xi/2,  xj/2, n-1);
   hilbert(x+xi/2,      y+xj/2 ,     xi/2, xj/2,  yi/2,  yj/2, n-1);
   hilbert(x+xi/2+yi/2, y+xj/2+yj/2, xi/2, xj/2,  yi/2,  yj/2, n-1);
   hilbert(x+xi/2+yi,   y+xj/2+yj,  -yi/2,-yj/2, -xi/2, -xj/2, n-1);

   }
end procedure;

//      x  y  xi  xj  yi  yj  n
    hilbert(0, 0, 1,  0,  0,  1,  1);

Tracking the results of the recursions will show what happens to original x,y coordinates - shown in above in blue. The i and the j vectors will be ignored. The four recursive calls will be labelled A, B, C and D ie.

A    hilbert(...);
    B    hilbert(...);
    C    hilbert(...);
    D    hilbert(...);

Performing the arithmetic on the x and y's for the first level of recursion for A, B, C and D generates the following x,y values,

x    y
    A    hilbert(0.0  0.0, ...);
    B    hilbert(0.5  0.0, ...);
    C    hilbert(0.5  0.5, ...);
    D    hilbert(0.5  1.0, ...);

The value of x and y can be plotted using the LineTo - listing 3.

Listing 3

procedure hilbert(x, y, xi, xj, yi, yj, n)
if (n <= 0) then
   LineTo(x  + (xi + yi)/2, y  + (xj + yj)/2");
else
   {
   hilbert(x,           y,           yi/2, yj/2,  xi/2,  xj/2, n-1);
   hilbert(x+xi/2,      y+xj/2 ,     xi/2, xj/2,  yi/2,  yj/2, n-1);
   hilbert(x+xi/2+yi/2, y+xj/2+yj/2, xi/2, xj/2,  yi/2,  yj/2, n-1);
   hilbert(x+xi/2+yi,   y+xj/2+yj,  -yi/2,-yj/2, -xi/2, -xj/2, n-1);
   }
end procedure;

Figures 3 to 6 show the development of a Hilbert curve from 1 to 4 iterations. The colored dots correspond to the four recursions ie.

A    hilbert(...); red dot
    B    hilbert(...); green dot
    C    hilbert(...); blue dot
    D    hilbert(...); purple dot

Figure 3 shows the basic building block of the Hilbert curve is a open square and that increasingly complex patterns, for example figure 4, are constructed from connected repetitions of this shape. In figure 4 the connections that join the four repititions (2 iterations) are shown as dotted lines.

	Figure 3 1 interation		Figure 4 2 interations
	Figure 5 3 interations		Figure 6 4 interations

Basic Code

Listings 4 and 5 provide code for so-called RenderMan procedural primitives, one written in Python and the other in Tcl. For information about procedural primitives, or helper apps as they are often called, refer to the tutorial "RenderMan Procedural Primitives: Implementations in Python, Tcl and 'C'". A sample rib file to run the implementations is given by listing 6.

Listing 4 (hilbert.py)

import sys, math
  
def hilbert(x0, y0, xi, xj, yi, yj, n):
    if n <= 0:
        X = x0 + (xi + yi)/2
        Y = y0 + (xj + yj)/2
        
        # Output the coordinates of the cv
        print '%s %s 0' % (X, Y)
    else:
        hilbert(x0,               y0,               yi/2, yj/2, xi/2, xj/2, n - 1)
        hilbert(x0 + xi/2,        y0 + xj/2,        xi/2, xj/2, yi/2, yj/2, n - 1)
        hilbert(x0 + xi/2 + yi/2, y0 + xj/2 + yj/2, xi/2, xj/2, yi/2, yj/2, n - 1)
        hilbert(x0 + xi/2 + yi,   y0 + xj/2 + yj,  -yi/2,-yj/2,-xi/2,-xj/2, n - 1)
        
def main():
    args = sys.stdin.readline()
    # Remain the loop until the renderer releases the helper...
    while args:
        arg = args.split()
        # Get the inputs
        pixels = float(arg[0])
        ctype = arg[1]
        reps = int(arg[2])
        width = float(arg[3])
        
        # Calculate the number of curve cv's
        cvs = int(math.pow(4, reps))
            
        # Begin the RenderMan curve statement
        print 'Basis \"b-spline\" 1 \"b-spline\" 1'
        print 'Curves \"%s\" [%s] \"nonperiodic\" \"P\" [' % (ctype, cvs)
    
        # Create the curve
        hilbert(0.0, 0.0, 1.0, 0.0, 0.0, 1.0, reps)
    
        # End the curve statement
        print '] \"constantwidth\" [%s]' % width
      
        # Tell the renderer we have finished   
        sys.stdout.write('\377')
        sys.stdout.flush()
        
        # read the next set of inputs
        args = sys.stdin.readline()
if __name__ == "__main__":
    main()

Listing 5 (hilber.tcl)

fconfigure stdout -translation binary
  
# This is where the control vertices of the Hilbert curve are generated    
proc hilbert { x0 y0 xi xj yi yj n } {
    if { $n <= 0 } {
        set X [expr $x0 + ($xi + $yi)/2]
        set Y [expr $y0 + ($xj + $yj)/2]
        
        # Output the coordinates of the cv
        puts "$X $Y 0 "
    } else {
        set XI [expr $xi/2]
        set XJ [expr $xj/2]
        set YI [expr $yi/2]
        set YJ [expr $yj/2]
        
        # Begin recursion
        hilbert $x0                 $y0                $YI  $YJ  $XI  $XJ [expr $n - 1]
        hilbert [expr $x0+$XI]     [expr $y0+$XJ]      $XI  $XJ  $YI  $YJ [expr $n - 1]
        hilbert [expr $x0+$XI+$YI] [expr $y0+$XJ+$YJ]  $XI  $XJ  $YI  $YJ [expr $n - 1]
        hilbert [expr $x0+$XI+$yi] [expr $y0+$XJ+$yj] -$YI -$YJ -$XI -$XJ [expr $n - 1]
        }
    }
  
# Remain the loop until the renderer releases the helper...
while { [gets stdin args] != -1 } {
    # Get the inputs
    set pixels [lindex $args 0]
    set ctype  [lindex $args 1]
    set reps   [lindex $args 2]
    set width  [lindex $args 3]
  
    # Calculate the number of curve cv's
    set cvs [expr int(pow(4, $reps))]
        
    # Begin the RenderMan curve statement
    puts "Basis \"b-spline\" 1 \"b-spline\" 1"
    puts "Curves \"$ctype\" \[$cvs\] \"nonperiodic\" \"P\" \["
    
    # Create the curve
    hilbert 0.0  0.0  1.0  0  0  1.0 $reps
    
    # End the curve statement
    puts "\] \"constantwidth\" \[ $width \]"
  
    # Tell the renderer we have finished   
    puts "\377" 
    flush stdout
    }

The inputs (arguments) to the helper apps are,
curve type - either linear or cubic
number of iterations
width of the curve

Figure 7
cubic - 4 iterations

Listing 6 (hilbert.rib)

Display "iterations" "framebuffer" "rgba"
Format 300 300 1
Projection "perspective" "fov" 12
ShadingRate 1
  
Translate  -0.5 -0.5 5
Rotate 0 1 0 0
Rotate 0   0 1 0
Scale 1 1 -1
WorldBegin
    TransformBegin
        Surface "constant"
        Color 1 1 1
        # For linux and MacOSC the paths to tclsh and
        # python must be fully specified ie,
        # /usr/bin/tclsh
        # /usr/bin/python
        # Args: curve type, iterations and curve width
        Procedural "RunProgram" ["tclsh  FULL_PATH/hilbert_helper.tcl" "cubic 4 0.005"] 
                    [-1 1 -1 1 -1 1]
        #Procedural "RunProgram" ["python  FULL_PATH/hilbert_helper.py" "linear 4 0.005"] 
        #           [-1 1 -1 1 -1 1]
    TransformEnd
WorldEnd

Mistakes can be Interesting

When entering the code for the Python and Tcl scripts I made a couple of mistakes that produced some interesting results. The first error occured in a python script.

hilbert(x0,               y0,               yi/2, yj/2, xi/2, xj/2, n - 1)
hilbert(x0 + xi/2,        y0 + xj/2,        xi/2, xj/2, yi/2, xj/2, n - 1)
hilbert(x0 + xi/2 + yi/2, y0 + xj/2 + yj/2, xi/2, xj/2, yi/2, xj/2, n - 1)
hilbert(x0 + xi/2 + yi,   y0 + xj/2 + yj,  -yi/2,-yj/2,-xi/2,-xj/2, n - 1)

Figure 8

The second error effected a tcl script.

hilbert $x0                    $y0                     $YI  $YJ  $XI  $XJ [expr $n - 1]
hilbert [expr $x0 + $XI]       [expr $y0 + $XJ]        $XI  $XJ  $YI  $YJ [expr $n - 1]
hilbert [expr $x0 + $XI + $YI] [expr $y0 + $XJ + $YJ]  $XI  $XJ  $YI  $YJ [expr $n - 1]
hilbert [expr $x0 + $XI + $YI] [expr $y0 + $XJ + $yj] -$YI -$YJ -$XI -$XJ [expr $n - 1]

Figure 9

Chunlei Chen

Search This Blog

Monday, March 14, 2011

Hilbert Curve Concepts & Implementation