Generating Hilbert curves

Posted by Steve Eddins, January 25, 2012

64 views (last 30 days) | 0 Likes | 5 comments

This week I came across some files I wrote about 16 years ago to compute Hilbert curves. A Hilbert curve is a type of fractal curve; here is a sample:

I can't remember why I was working on this. Possibly I was anticipating that 16 years in the future, during an unusually mild New England winter, I would be looking for a blog topic.

Anyway, there are several interesting ways to code up a Hilbert curve generator. My old code for generating the Hilbert curve followed the J. G. Griffiths, "Table-driven Algorithms for Generating Space-Filling Curves," Computer-Aided Design, v.17, pp. 37-41, 1985. Here's the basic procedure:

First, initialize four curves, \( A_0 \), \( B_0 \), \( C_0 \), and \( D_0 \), to be empty (no points).

Then, build up the Hilbert curve iteratively as follows:

\[ A_{n+1} = [B_n, N, A_n, E, A_n, S, C_n] \]

\[ B_{n+1} = [A_n, E, B_n, N, B_n, W, D_n] \]

\[ C_{n+1} = [D_n, W, C_n, S, C_n, E, A_n] \]

\[ D_{n+1} = [C_n, S, D_n, W, D_n, N, B_n] \]

where \( N \) represents a unit step up, \( E \) is a unit step to the right, \( S \), is a unit step down, and \( W \) is a unit step left.

One way to code this procedure is to incrementally build up a set of vectors that define the step from one point on the path to the next, and then to use a cumulative summation at the end to turn the steps into x-y coordinates. Here's how you might do it.

A = zeros(0,2);
B = zeros(0,2);
C = zeros(0,2);
D = zeros(0,2);

north = [ 0  1];
east  = [ 1  0];
south = [ 0 -1];
west  = [-1  0];

order = 3;
for n = 1:order
  AA = [B ; north ; A ; east  ; A ; south ; C];
  BB = [A ; east  ; B ; north ; B ; west  ; D];
  CC = [D ; west  ; C ; south ; C ; east  ; A];
  DD = [C ; south ; D ; west  ; D ; north ; B];

  A = AA;
  B = BB;
  C = CC;
  D = DD;
end

A = [0 0; cumsum(A)];

plot(A(:,1), A(:,2), 'clipping', 'off')
axis equal, axis off

I was curious to see what might be on the MATLAB Central File Exchange, so I searched for "hilbert curve" and found several interesting contributions. There are a couple of 3-D Hilbert curve generators, and several different ways of coding up a 2-D Hilbert curve generator. I was particularly interested in the Fractal Curves contribution by Jonas Lundgren.

Jonas shows a much more compact implementation of the ideas above using complex arithmetic. It looks like this:

a = 1 + 1i;
b = 1 - 1i;

% Generate point sequence
z = 0;
order = 5;
for k = 1:order
    w = 1i*conj(z);
    z = [w-a; z-b; z+a; b-w]/2;
end

plot(z, 'clipping', 'off')
axis equal, axis off