Hilbert curve in R

hilbert.r

# Hilbert curve in R

# Make some matrices to do rotation and flip
rotmat = function(theta) matrix(c(cos(theta),-sin(theta),sin(theta), cos(theta)),nrow=2)
rotclock = rotmat(pi/2)
rotaclock = rotmat(-pi/2)
flip=matrix(c(-1,0,0,1),nrow=2)
ident=matrix(c(1,0,0,1),nrow=2)

# Define the recursive function to generate the points in the right order
hilbert <- function(n,x,y,size,mat){
  pts = (cbind(c(-1,-1,1,1),c(-1,1,1,-1)) %*% mat)*size/4
  #if(runif(1)<0.1 & n>1) n=n-1 
  if(n==1){
    xs <<- c(xs,x +size/2+ pts[,1] )
    ys <<- c(ys,y +size/2+ pts[,2] )
   }
  else{
    x= x+size/4
    y= y+size/4
    hilbert(n-1,x+pts[1,1],y+pts[1,2],size/2,mat%*%flip%*%rotaclock) # bottom left
    hilbert(n-1,x+pts[2,1],y+pts[2,2],size/2,mat)                    # top left
    hilbert(n-1,x+pts[3,1],y+pts[3,2],size/2,mat)                    # top right
    hilbert(n-1,x+pts[4,1],y+pts[4,2],size/2,mat%*%flip%*%rotclock)  # bottom right
  }
}

# Generate the curve and draw the lines
doHilbert <- function(o,x=0,y=0,size=1,mat=ident,...){
  xs <<- numeric(0)
  ys <<- numeric(0)
  hilbert(o,x,y,size,mat)  
  lines(xs,ys,...)
}

# Make a png with two Hilbert curves overlaid
png("hilbert1.png",width=1000,height=1000,type="cairo")
par(mfrow=c(1,1),mar=3*c(1,1,1,1))
plot(NA,xlim=c(0,1),ylim=c(0,1),axes=F,ann=F,asp=1)
doHilbert(6,col="red",lwd=1.5)
doHilbert(5,col="black",lwd=2)
dev.off()