In the images below you can see two fractal structures. The first is the famous “Koch snowflake” named after the Swedish mathematician Helge von Koch. The second is based on a “recursive tree”. Both of these structures are “self similar”. These means that they are based on a pattern that keeps on repeating itself in ever greater detail. If you were to zoom in on these structures, you would see the same kind of patterns over and over again. In mathematics such geometric figures are called “fractals”. In this tutorial you will learn how can create and draw such structures yourself using nothing but base R.

.

Lines in space

Essentially, the images above are nothing but a set of lines in two-dimensional space. Therefore, we will start with a function that creates an empty two-dimensional space or simply a canvas to draw an.

# function to create empty canvas
emptyCanvas <- function(xlim, ylim, bg="gray20") {
  par(mar=rep(1,4), bg=bg)
  plot(1, 
       type="n", 
       bty="n",
       xlab="", ylab="", 
       xaxt="n", yaxt="n",
       xlim=xlim, ylim=ylim)
}

# example
emptyCanvas(xlim=c(0,1), ylim=c(0,1))

Admittedly this wasn’t all that fascinating yet. So let’s write a function called drawLine(line) that we can use to draw some lines. Every line will be a simple vector consisting of four coordinates. An x and y coordinate for the starting point and an x and y coordinate for the end point of the line. To draw multiple lines at once, we will create a second function called drawObject(object). This function is really just a wrapper around the first. Its main argument will not be a vector but a matrix. In this matrix each row represents a different line to be plotted.

# function to draw a single line
drawLine <- function(line, col="white", lwd=1) {
  segments(x0=line[1], 
           y0=line[2], 
           x1=line[3], 
           y1=line[4], 
           col=col,
           lwd=lwd)
}

# wrapper around "drawLine" to draw entire objects
drawObject <- function(object, col="white", lwd=1) {
  invisible(apply(object, 1, drawLine, col=col, lwd=lwd))
}

# example
line1 = c(0,0,1,1)
line2 = c(-3,4,-2,-4)
line3 = c(1,-3,4,3)
mat = matrix(c(line1,line2,line3), byrow=T, nrow=3)
mat
##      [,1] [,2] [,3] [,4]
## [1,]    0    0    1    1
## [2,]   -3    4   -2   -4
## [3,]    1   -3    4    3
# draw separately
emptyCanvas(xlim=c(-5,5), ylim=c(-5,5))
drawLine(line1)
drawLine(line2, col="orange", lwd=4)
drawLine(line3, col="cyan", lwd=6)

# draw together
emptyCanvas(xlim=c(-5,5), ylim=c(-5,5))
drawObject(mat, col="yellowgreen", lwd=3)

Creating structures from lines

Next, we define a function called newLine() that we can use to generate new lines based on existing ones. Given an earlier line \(a\), newLine(line, angle, reduce) will create a line \(b\) that originates at the end point of line \(a\). The latter two arguments of newLine() control the angle between the two lines and how much shorter (or longer) \(b\) is in relation to \(a\).

# function to add a new line to an existing one
newLine <- function(line, angle, reduce=1) {
  
  x0 <- line[1]
  y0 <- line[2]
  x1 <- line[3]
  y1 <- line[4]
  
  dx <- unname(x1-x0)                      # change in x direction
  dy <- unname(y1-y0)                      # change in y direction
  l <- sqrt(dx^2 + dy^2)                   # length of the line
  
  theta <- atan(dy/dx) * 180 / pi          # angle between line and origin
  rad <- (angle+theta) * pi / 180          # (theta + new angle) in radians
  
  coeff <- sign(theta)*sign(dy)            # coefficient of direction
  if(coeff == 0) coeff <- -1
  
  x2 <- x0 + coeff*l*cos(rad)*reduce + dx  # new x location
  y2 <- y0 + coeff*l*sin(rad)*reduce + dy  # new y location
  return(c(x1,y1,x2,y2))
  
}
# example
a = c(-1.2,1,1.2,1)
b = newLine(a, angle=-90, reduce=1)
c = newLine(b, angle=45, reduce=.72)
d = newLine(c, angle=90, reduce=1)
e = newLine(d, angle=45, reduce=1.4)

# draw lines
emptyCanvas(xlim=c(-5,5), ylim=c(0,5))
drawLine(a, lwd=3, col="white")
drawLine(b, lwd=3, col="orange")
drawLine(c, lwd=3, col="cyan")
drawLine(d, lwd=3, col="firebrick")
drawLine(e, lwd=3, col="limegreen")

Now that we have our “graphical engine” ready, we will define one final wrapper function called iterate(object, ifun, ...) before we start programming our first fractals. Technically we don’t need this wrapper function to draw our first fractals, but it will come in handy when we start playing around with different kinds of fractals. With iterate() we will be able to quickly apply any sort of fractal function ifun(), like the Koch curve or recursive trees, to any set of lines stored in object.

# function to run next iteration based on "ifun()"
iterate <- function(object, ifun, ...) {
  linesList <- vector("list",0)
  for(i in 1:nrow(object)) {
    old_line <- matrix(object[i,], nrow=1)
    new_line <- ifun(old_line, ...)
    linesList[[length(linesList)+1]] <- new_line
  }
  new_object <- do.call(rbind, linesList)
  return(new_object)
}

The Koch curve

The Koch curve starts from a straight line which is broken into three equally long pieces. Next, two new lines, which are just as long as one the three pieces, are added on top of the initial line. Together with the center piece they form an equilateral triangle. Lastly, the center piece is removed from the initial line. These steps are then recursively applied to every line in the curve. Here is what that looks like after zero, one, two and three iterations:

In terms of programming the Koch curve, we start from the end point of the initial line and make a 180 degree turn. Then we use the fact that all angles of an equilateral triangle have 60 degrees. This means we first turn by negative 60 degrees, then by 120 degrees and finally another negative 60 degrees to reach the initial start point. For the first of the four line segments, we set reduce=1/3, thereafter we set reduce=1, as we do not need to shorten the segments any further. Since we started at the initial end point, we finally need to reorder the coordinates and lines to make sure the resulting object is in the correct order. Here is our final function:

# iterator function: koch curve
koch <- function(line0) {
  
  # new triangle (starting at right)
  line1 <- newLine(line0, angle=180, reduce=1/3)
  line2 <- newLine(line1, angle=-60, reduce=1)
  line3 <- newLine(line2, angle=120, reduce=1)
  line4 <- newLine(line3, angle=-60, reduce=1)
  
  # reorder lines (to start at left)
  line1 <- line1[c(3,4,1,2)]
  line2 <- line2[c(3,4,1,2)]
  line3 <- line3[c(3,4,1,2)]
  line4 <- line4[c(3,4,1,2)]
  
  # store in matrix and return
  mat <- matrix(c(line4,line3,line2,line1), byrow=T, ncol=4)
  return(mat)
  
}

# example: Koch curve (after six iterations)
fractal <- matrix(c(10,0,20,1e-9), nrow=1)
for(i in 1:6) fractal <- iterate(fractal, ifun=koch)
emptyCanvas(xlim=c(10,20), ylim=c(0,3))
drawObject(fractal)

If we start the recursive process with a triangle instead of a simple line, we obtain the Koch snowflake.

# Koch snowflake (after six iterations)
A <- c(0,1e-9)
B <- c(3,5)
C <- c(6,0)
fractal <- matrix(c(A,B,B,C,C,A), nrow=3, byrow=T)
for(i in 1:6) fractal <- iterate(fractal, ifun=koch)
emptyCanvas(xlim=c(-2,8), ylim=c(-2,5))
drawObject(fractal)

Notice that I set some coordinates to 1e9 instead of exactly \(0\). This is to avoid some corner cases but does not make any difference at all when drawing the fractal images.

Recursive trees

Recursive trees are based on the concept of large branches splitting into smaller branches that in turn split into even smaller branches. The whole concept is so intuitive, just have a look at the first couple of iterations and you will immediately understand what is going on.

When programming our new ifun() function called tree(), we need it to draw just two new lines. Both start at the end point of the initial line. The first branches out to the left, the second branches out to the right.

# iterator function: recursive tree
tree <- function(line0, angle=30, reduce=.7, randomness=0) {
  
  # angles and randomness
  angle1 <- angle+rnorm(1,0,randomness)  # left branch
  angle2 <- -angle+rnorm(1,0,randomness) # right branch
  
  # new branches
  line1 <- newLine(line0, angle=angle1, reduce=reduce)   
  line2 <- newLine(line0, angle=angle2, reduce=reduce)
  
  # store in matrix and return
  mat <- matrix(c(line1,line2), byrow=T, ncol=4)
  return(mat)
  
}

# example: recursive tree (after ten iterations)
fractal <- matrix(c(0,0,0,10), nrow=1)
emptyCanvas(xlim=c(-30,30), ylim=c(0,35))
drawObject(fractal)
for(i in 1:10) {
  fractal <- iterate(fractal, ifun=tree, angle=23)
  drawObject(fractal)
}

You will have noticed that the tree() function has an argument called randomness. We can use this argument to introduce slight deviations from the perfect symmetry in the example above. If we also reduce the width of the branches as we move up the tree, we can create more realistically looking trees.

# recursive tree (organic)
set.seed(1234)
fractal <- matrix(c(0,0,0,10), nrow=1)
emptyCanvas(xlim=c(-30,30), ylim=c(0,35))
lwd <- 7
drawObject(fractal, lwd=lwd)
for(i in 1:12) {
  lwd <- lwd*0.75
  fractal <- iterate(fractal, ifun=tree, angle=29, randomness=9)
  drawObject(fractal, lwd=lwd)
}

Lastly, to make things even more interesting, we can also grow multiple trees at once. Here is an example of four trees originating in the same point with different colors for every level of the tree.

# example: recursive tree
Z <- c(0,0)
A <- c(1e-9,5)
B <- c(5,-1e-9)
fractal <- matrix(c(Z,A,Z,B,Z,-A,Z,-B), nrow=4, byrow=T)
emptyCanvas(xlim=c(-20,20), ylim=c(-20,20))
drawObject(fractal)
for(i in 1:11) {
  fractal <- iterate(fractal, ifun=tree, angle=29, reduce=.75)
  drawObject(fractal, col=i+1)
}

Now, you are ready to develop your own iterative functions. Once you defined your rule, simply pass it as ifun to the iterate() function and watch the magic unfold. If you liked this post or fractals in general or if you are just into drawing beautiful pictures with R, you might also enjoy my post on using R to find “Strange Attractors”. You can find this post here.