Tag Archives: geometry

How to draw the Hilbert Curve using a computer program

In this note we will learn how to draw the Hilbert Curve using a computer program written in Scheme.

There are a number of descriptions in books, papers, and online that talk about how to draw the curve. (For a list of the materials that I reviewed while trying to learn this, see the References section below.) Unfortunately, I found most of the descriptions very difficult to understand!

Most descriptions of how to draw the Hilbert curve seem to fall into one of the following categories:

  1. Using magic tables of numbers to generate coordinates, which is fine if you already understand how to construct the curve and then work backwards from there to figure out how to generate those magic numbers. But it doesn’t really help you understand how to draw it in the first place.
  2. A “grammatical” approach using L-systems (about which see the wiki article linked above), which is fine if you want to take a detour into writing an interpreter for L-systems. Unfortunately, I did not — or at least, I did not think it would teach me much about how to actually draw the Hilbert curve. But it’s definitely a cool approach!
  3. Iterating over C arrays and doing magical bitshifts to represent both the points on the curve and the various rotations/reflections that are applied to those points.

In my opinion, none of these approaches are really useful for learning. They are fine if you are looking for implementation techniques for something you already understand.

Thankfully all was not lost. I spent a lot of time reading, and scratching my head, and drawing curves in my notebook. Then, I found the excellent article Crinkly Curves by Brian Hayes. In addition to the great information in the article, Hayes linked to a really helpful visualization on his website called “Constructing the Hilbert Curve”. I urge you to try it, since seeing the visualization is important to what follows.

Seeing this is what helped me finally understand. Thank you, Brian Hayes!

The algorithm displayed in Hayes’ visualization is approximately what I ended up implementing, and it looked like this in Scheme:

(define (hilbert order width)
  (let loop ((count 0)
             (shape (h0 width)))
    (if (= count order)
        shape
        (let* ((next-shape (shrink-shape shape 0.5))
               (upper-left (reflect-shape (rotate-shape next-shape (* 3.14159 1.5))))
               (upper-right (translate-shape next-shape 'east width))
               (lower-left (translate-shape next-shape 'south width))
               (lower-left* (reflect-shape (rotate-shape lower-left (* 3.14159 .5))))
               (lower-right (translate-shape next-shape 'southeast width)))
          (loop (+ count 1)
                (merge-shapes upper-left upper-right lower-right lower-left*))))))

It works like this:

  1. Given a “shape” (starting with the zeroth order Hilbert curve, which is just point in the center of a square), make 4 copies of that shape, each of which fits in a square 1/2 the width of the original square. Then, “move” the shape to one of the 4 corners of the original square. (The “move” is accomplished by a call to translate-shape as shown above.)
  2. Depending on which corner of the full-sized square the shape is now in, it may need to be rotated (spun) around its axis so it points in the right direction. In some cases it may also need to be reflected so that the locality of the points in the curve is correctly preserved. Note that the calls to rotate-shape take a radians argument, a unit for measuring angles. 360 degrees is equal to (* 3.14159 2) in radians. For more information, see this article.
  3. Keep repeating the previous two steps, bumping a counter and accumulating a list of points until we have all of the points for an Nth order curve. We accumulate the list of points using merge-shapes, which is really just append.
  4. Base case: The counter is now equal to the order argument, telling us we have an Nth order Hilbert curve. So we return the accumulated list of points.

These “points” are really just 2-element (X Y) lists, as you can see below:

> (hilbert 2 600))

((74.99970147174234 75.00029852944591) (224.9997014705541 74.99970147174234)
 (225.00029852825764 224.9997014705541)
 (75.00029852944591 225.00029852825764)
 (75 375)
 (75 525)
 (225 525)
 (225 375)
 (375 375)
 (375 525)
 (525 525)
 (525 375)
 (525.0000995095512 224.99990049031675)
 (375.00009950968325 225.00009950955126)
 (374.99990049044874 75.00009950968327)
 (524.9999004903167 74.99990049044875))

Using some not-so-interesting utility code that builds an SVG (XML) file and does some geometry things like shape rotation, etc., we can draw this list of points by outputting it into an SVG polyline:

(let ((content (tag->string (svg (shape->polyline (rotate-shape (hilbert 5 600) (* 3.14159 0.5)))))))
  (with-output-to-file "hilbert.svg"
    (lambda ()
      (display content))))

The file hilbert.svg output by the above code looks like this:

We can make an even prettier picture if we stack the curves atop each other in the same image as follows:

(let ((content (tag->string
                (svg (reverse (list
                               (shape->polyline (rotate-shape (hilbert 1 600) (* 3.14159 0.5)) 'black)
                               (shape->polyline (rotate-shape (hilbert 2 600) (* 3.14159 0.5)) 'red)
                               (shape->polyline (rotate-shape (hilbert 3 600) (* 3.14159 0.5)) 'blue)
                               (shape->polyline (rotate-shape (hilbert 4 600) (* 3.14159 0.5)) 'orange)
                               (shape->polyline (rotate-shape (hilbert 5 600) (* 3.14159 0.5)) 'yellow)))))))
  (with-output-to-file "multi-hilbert.svg" (lambda () (display content))))

Here is the image in multi-hilbert.svg:

References

Just for Fun: Estimating pi with Scheme

../img/double-star-from-pentagon-backlit.jpg

A while back I shared some Perl code for calculating the circumference of a circle without knowing 𝛑. Just for fun, and due to my longtime infatuation with all things Schemish, I’ve written a little pi approximator in Scheme. It uses the idea that we can approximate a circle using smaller and smaller triangles stacked on top of each other. (See previously for a better explanation with a picture.)

And now, the code!

;;;; pi.scm -- Estimate the value of 𝛑 using smaller and smaller
;;;; triangles.

;;; Call it like so: (pi-estimate n), where n is the number of
;;; iterations you'd like to go through. It doesn't take many to get
;;; pretty accurate.

(define reference-pi 3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679)

(define (circumference radius)
  (* 3.14159 2 radius))

(define (square x)
  (* x x))

(define (hypotenuse a b)
  (sqrt (+ (expt a 2) (expt b 2))))

(define (pi-iter radius a b count maxcount)
  (let* ((hyp (hypotenuse a b))
         (newbase (- radius (sqrt (- (square radius)
                                     (square (/ hyp 2)))))))
    (if (= count maxcount)
        (/
         (* 2 hyp (expt 2 (+ 1 count)))
         (* 2 radius))
        (pi-iter radius newbase (/ hyp 2) (+ count 1) maxcount))))

(define (pi-estimate iterations)
  (pi-iter 128 128 128 0 iterations))

(Origami image courtesy Melisande under Creative Commons license.)

How to calculate the circumference of a circle (without knowing pi)

Some time ago I posted the following little challenge to www.reddit.com/r/programming:

Dear Proggit: Imagine that you must devise a method to calculate the circumference of a circle of any radius. The catch? You are unaware of the existence of Pi.

The other catch? You must begin right now. You may use only what tools you have in front of you (text editor, calculator, pencil and paper with compass(!)), without looking anything up. That means no Google searches, no Wikipedia, no formulas in a book, nothing.

Once you have devised a method, please write a small program in the language of your choice that performs the calculation (the more accurate, the better). Bonus points will be awarded for brevity, cleverness, etc.

I will post my little solution sometime tomorrow. Enjoy!

One final note: If you are a total geometry/math wizard who thinks this is lame, please refrain from posting; I want this to be a fun problem for those of us who will get something out of it.

https://logicgrimoire.files.wordpress.com/2012/09/wpid-circle-diagram.png

Below is my method for calculating (or rather, estimating) the circumference of a circle without using pi. This method involves summing the lengths of the hypotenuses of smaller and smaller triangles, as shown in the diagram above.

Begin by cutting the given circle into quarters. Using the radius r, find the hypotenuse of triangle A using Pythagoras’ method. Cut the hypotenuse in half; using that, and the radius, you can calculate the length of line segment Q. Subtract Q from r, and you have the length of line segment P, which is one side of the right triangle B. That, combined with half the hypotenuse of A, will allow you to calculate the hypotenuse of B. Continue, calculating the hypotenuses of smaller and smaller triangles. Sum as many lengths as are required to cover the quarter-circle arc, and multiply by 4 (or multiply the current length by 2^(n*+1), where *n is the number of iterations of the method you’ve performed).

Below I’ve provided an implementation in Perl. Note that the Perl interpreter doesn’t optimize away recursive function calls (that I know of), so you could in theory blow up the stack. In practice, it only takes a few iterations to arrive at a reasonable estimate, so it’s not a problem.

#!/usr/bin/env perl

use strict;
use warnings;
use autodie;
use Math::Trig ':pi';

my ($radius, $maxiter) = @ARGV;
my $iterations = 0;

sub main {
  my ($radius, $a, $b) = @_;
  my $hyp = hypotenuse($a, $b);
  my $newbase = $radius - sqrt($radius**2 - ($hyp / 2)**2);
  $iterations++;
  die "Current hypotenuse is $hyp after $iterations iterations\n"
    . "Estimated circumference: " . $hyp * (2**($iterations+1)) . "\n"
      . "Actual circumference: " . ($radius*2) * pi . "\n"
        if $iterations == $maxiter;
  main($radius, $newbase, ($hyp / 2));
}

sub hypotenuse {
  my ($a, $b) = @_;
  return sqrt($a**2 + $b**2);
}

main($radius, $radius, $radius);

As you can see, we pretty much mirror the textual description given above, but in code. We then compare our work to the “real” value of pi using the Perl core’s Math::Trig constant. Save the above into a file named circ.pl, and you can run it like so:

mel@foo:~$ ./circ.pl 4 12 # Given a radius of 4, run through 12 iterations
Current hypotenuse is 0.00306796150057116 after 12 iterations
Estimated circumference: 25.132740612679
Actual circumference: 25.1327412287183

I might not want to try to calculate planetary orbits using this kind of rough estimation, but for 99% of real-world cases, it’s good enough.