How to draw fractals in R: Sierpínski Triangles

Fractals are beautiful objects that appeal to both artists and mathematicians alike. This is the first part of my new series showing how to create them in the R programming language. We start with one of the simplest but most iconic fractals, the Sierpínski Triangle.

What is a fractal?

Fractals are not only beautiful but they have strange mathematical properties. They have an infinite amount of detail. For example the Sierpínski Triangle has a infinite number of triangles, which get ever smaller. We will never be able to draw all this detail, but we draw enough to get an impression of the true fractal.

They are also self-similar which means a small part of it contains the same or similar detail of the complete fractal. Looking at the Sierpínski Triangle, we can cut it into three Sierpínski Triangles which are exact copies of the original but half the size.

This leads to the weirdest property of fractals: They have smaller dimensions than you expect. When you draw a square and halve the length of the sides the square you get is four times smaller. When you draw a circle and halve the radius the circle you get is four times smaller. This is what happens with a two dimensional shape but it doesn’t happen with the Sierpínski Triangle and other fractals. Instead of being four times smaller, the three copies we made are a third of the size of the original. This is why mathematicians say the Sierpínski Triangle has a dimension between 1 and 2. Explaining fractional dimensions is beyond the scope of this post, so let’s focus on having fun drawing our first fractal.

Drawing Fractals in R.

R is a programming language widely used for statistical analysis. It’s also a great tool for graphs, graphics and visualisations. It also has a fantastic community who have created a vast library of packages letting you do almost anything. This guide assumes you have installed R and have a little experience with programming. R is completely free and can be downloaded from https://www.r-project.org/. It’s available for Windows, Linux and Mac.

To draw a fractal from this blog, just open R and copy and paste the code into the console. Alternatively you can create a script and tweak and change the fractal to your heart’s content. I’ve included the complete code at the end of the post.

Before We Begin.

The Sierpínski Triangle has been used for decorative patterns for centuries. Best of all, it’s very simple to draw. Starting with an ordinary triangle, we divide into three triangles that are a quarter the size of the original triangle. These triangles then get split into three. Repeat this process a forever and you would get the complete fractal. Just repeating it a few times gives a nice picture like the ones shown here. These pictures may not have the infinite detail of a real fractal, but they look great and they help you understand the fractal.

Before we write any code we need to know what we are actually trying to do. The finished fractal will be a set of triangles so we will use the list structure to store the triangles. The code for a list is simple:

list(.....)

Each triangle will be a data frame. Data frames are how we store tables and spreadsheets in R. The data frame will have three rows, one for each point of the triangle. It will have two columns, one for x and one for y. The code for this will be:

data.frame(x=c(...), y=c(...))

To make and draw the fractal there are four steps;

1. Create a function to shrink a triangle.
2. Create a function to split a triangle into three triangles.
3. Create the fractal starting from one triangle.
4. Draw the fractal we have created.

Shrinking Triangles.

Let’s create a function to shrink a triangle called shrinkTriangle. When we shrink a triangle, we don’t want to shrink it towards the centre of the triangle, but towards one corner. Hence when we call shrinkTriangle it will look like this;

shrinkTriangle(triange, corner)

where triangle is a data frame, and corner is the number of the row that doesn’t change. Let’s build our function step by step. The variable shrunkenTriangle will be our answer. Since one corner of the input triangle doesn’t change, we use it as the basis for calculating shrunkenTriangle.

shrinkTriangle=function(triangle,corner){
   shrunkenTriangle=triangle
   ...
   return(shrunkenTriangle)
}

Next we go through the corners one by one to see which ones change.

shrinkTriangle=function(triangle,corner){
   shrunkenTriangle=triangle
   if(corner!=1){
      ...
   }
   if(corner!=2){
      ...
   }
   if(corner!=3){
      ...
   }
   return(shrunkenTriangle)
}

The new corners are halfway between the old corner and the corner that stays fixed. For example the point half way between corner 1 and corner 2 has x co-ordinate 0.5*(triangle[1,"x"]+triangle[2,"x"]). The y co-ordinate is 0.5*(triangle[1,"y"]+triangle[2,"y"]).

In fact we don’t need to worry about adding x and y here. We can add the rows of triangle together so 0.5(triangle[1,]+triangle[2,]) gives us the new x and y co-ordinates in one package.

shrinkTriangle=function(triangle,corner){
   shrunkenTriangle=triangle
   if(corner!=1){
      shrunkenTriangle[1,]=0.5*(triangle[1,]+triangle[corner,])
   }
   if(corner!=2){
      shrunkenTriangle[2,]=0.5*(triangle[2,]+triangle[corner,])
   }
   if(corner!=3){
      shrunkenTriangle[3,]=0.5*(triangle[3,]+triangle[corner,])
   }
   return(shrunkenTriangle)
}

This code could be shorter using a for loop running over the corners, but this is fine for now. In the next part of this series we will improve and adapt this function so that we can draw more fractals.

Splitting a Triangle into Three.

Our second function takes one triangle and turns it into a list of three smaller triangles. It simply calls shrinkTriangle three times, once for each corner, and stores the new triangles in a list.

sierpinskiSplit=function(inputTriangle){
   return(list(shrinkTriangle(inputTriangle,1),
               shrinkTriangle(inputTriangle,2),
               shrinkTriangle(inputTriangle,3)))
}

Making the Sierpínski Triangle.

We will use co-ordinates that all lie between 0 and 1. We have to calculate the co-ordinates of first or outer triangle ourselves. The base of the triangle goes from x=0,y=0 to x=1,y=0. Using Pythagoras’s theorem gives the y co-ordinate of the third point. I have listed the corners from top to bottom and left to right. Note that the y co-ordinates in R run from bottom to top which is great for drawing graphs, but different to other programs.

firstTriangle=data.frame(x=c(0.5,0,1),
                         y=c(sqrt(3)/2,0,0))

This triangle is completely fine, but it isn’t centred vertically between 0 and 1. Adding 0.05 to all of the y co-ordinates shifts it up slightly to solve this problem.

firstTriangle=data.frame(x=c(0.5,0,1),
                         y=c(sqrt(3)/2+0.05,0.05,0.05))

If we only wanted to draw three triangles, we would write;

myFractal=list(firstTriangle)
myFractal=sierpinskiSplit(myFractal[[1]])

If wanted to draw nine triangles we would then split each triangle again. We create oldFractal to save a copy of the previous stage. We then make myFractal into a list of length nine.

oldFractal=myFractal
myFractal=vector("list",9)
myFractal[1:3]=sierpinskiSplit(oldFractal[[1]])
myFractal[4:6]=sierpinskiSplit(oldFractal[[2]])
myFractal[7:9]=sierpinskiSplit(oldFractal[[3]])

We can also write this with a for loop.

oldFractal=myFractal
myFractal=vector("list",9)
for(iter in seq_along(oldFractal)){
   myFractal[(3*iter-2):(3*iter)]=sierpinskiSplit(oldFractal[[iter]])
}

We have now run this step twice. To run this step seven times, let’s use another for loop. The big difference is that we have to calculate the size of the next iteration when we create a new version of myFractal. This is easy because the list of triangles is always three times larger each time.

for(repeats in 1:7){
   oldFractal=myFractal
   myFractal=vector("list",length(myFractal)*3)
   for(iter in seq_along(oldFractal)){
      myFractal[(3*iter-2):(3*iter)]=sierpinskiSplit(oldFractal[[iter]])
   }
}

Plotting the Finished Fractal.

Let’s start by creating a blank canvas to work on. By default R plots include lots of space for margins and titles. To remove these and create a blank canvas we can use the following code;

oldpar=par(mai=c(0,0,0,0),xaxt="n",yaxt="n")
plot(0:1,0:1,type="n")

To finish we just draw each triangle one by one and our fractal is complete. The last line of code restores the original settings for plotting.

for(iter in seq_along(myFractal)){
   polygon(myFractal[[iter]], col="black")
}

par(oldpar)

Rainbow Variation.

To make a more striking version we will add a black background and colour each triangle a different colour. For my version I also reduced the number of repeated triangle splittings to six.

polygon(firstTriangle,col="black")
colours=rainbow(length(myFractal))
for(iter in seq_along(myFractal)){
   polygon(myFractal[[iter]],col=colours[iter],bg=colours[iter])
}

Full Code

Here is the complete code. Just copy and paste it into R and your first fractal will appear.

shrinkTriangle=function(triangle,corner){
   shrunkenTriangle=triangle
   if(corner!=1){
      shrunkenTriangle[1,]=0.5*(triangle[1,]+triangle[corner,])
   }
   if(corner!=2){
      shrunkenTriangle[2,]=0.5*(triangle[2,]+triangle[corner,])
   }
   if(corner!=3){
      shrunkenTriangle[3,]=0.5*(triangle[3,]+triangle[corner,])
   }
   return(shrunkenTriangle)
}

sierpinskiSplit=function(inputTriangle){
   return(list(shrinkTriangle(inputTriangle,1),
               shrinkTriangle(inputTriangle,2),
               shrinkTriangle(inputTriangle,3)))
}

firstTriangle=data.frame(x=c(0.5,0,1),
                         y=c(sqrt(3)/2+0.05,0.05,0.05))

myFractal=list(firstTriangle)

for(repeats in 1:7){
   oldFractal=myFractal
   myFractal=vector("list",length(myFractal)*3)
   for(iter in seq_along(oldFractal)){
      myFractal[(3*iter-2):(3*iter)]=sierpinskiSplit(oldFractal[[iter]])
   }
}

oldpar=par(mai=c(0,0,0,0),xaxt="n",yaxt="n")

plot(0:1,0:1,type="n")

for(iter in seq_along(myFractal)){
   polygon(myFractal[[iter]], col="black")
}

par(oldpar)