Fractals with R, Part 3: The Hilbert Curve

Imagine a meandering river. Imagine, furthermore, that every meandering has a meandering within it. This is one way that you might visualize the Hilbert Curve, an example of a space-filling curve. The complexity of this “river” may have implications for ecological study design. How, for example, would you measure the distance between the source of the river (lower left) and the ocean (lower right)? For a mathematical description of the Hibert Curve, see Sagan (1994).

Figure 1. Iterations in the construction of the Hilbert Curve.

The graphics in Figure 1 can be produced by cutting and pasting the R script given below into your R console. For more, see Fractals with R, Part 1 and Part 2.

References

Sagan, Hans. 1994. “Space-filling curves”. Springer-Verlag.

R Script

#Written by Allan Roberts, 2012.
hilbert.curve <- function(n){

Double <- function(A){
#The matrix for “H(n)” is equal to “Double(H(n-1))”.

m <- dim(A)[1];
n <- dim(A)[2];
N <- m*n;
B <- A+N;
C <- B+N;
D <- C+N;
E <- cbind(rbind(B,skew.transpose(A)),rbind(C,t(D)));
return(E);
}

Rotate <- function(A){
#Rotates the matrix A clockwise.

m <- dim(A)[1];
n <- dim(A)[2];
N <- m*n;
B <- matrix(0,m,n);
for (i in 1:m) for (j in 1:n) B[j,n+1-i] <- A[i,j]
return(B);
}

skew.transpose <- function(A){
return(Rotate(Rotate(t(A))));
}

rowofx <- function(A,x){

#Returns the row index of the matrix A for entry equal to x.
m <- dim(A)[1];
n <- dim(A)[2];
for (i in 1:m) for (j in 1:n) if (A[i,j]==x) return(i);
}

colofx <- function(A,x){

#Returns the column index of the matrix A for entry equal to x.
m <- dim(A)[1];
n <- dim(A)[2];
for (i in 1:m) for (j in 1:n) if (A[i,j]==x) return(j);
}

Draw <- function(A){
#Draws a graphical representation of the matrix A.
A <- Rotate(A);
m <- dim(A)[1];
n <- dim(A)[2];
N <- m*n;
plot(  (rowofx(A,1)-1)/n, (colofx(A,1)-1)/n, pch=19,cex=0.5,ylim = c(0,1), xlim =c(0,1),     ylab=character(1),xlab=character(1),axes=FALSE);
d <- 1/n;
for (i in 1:(N-1)) lines(c((rowofx(A,i)-1)/n,((rowofx(A,i+1)-1)/n)), c((colofx(A,i)-1)/n,((colofx(A,i+1)-1)/n)),lwd=1);
points((rowofx(A,N)-1)/n, (colofx(A,N)-1)/n, pch=19,cex=0.5);
}

H <- function(n){
#H(1) is shown in Figure 2.
if (n==0) return(matrix(c(2,1,3,4),2,2));
return(Double(H(n-1)));
}

Draw(H(n));
}
hilbert.curve(n=3);