Fractal Snowflakes: A Recursive Algorithm

When I was out for a walk in the snow the other day, I thought of a simple, recursive way to draw a snowflake.

Start with a number of line segments, all starting at a single point, radiating equally-spaced in each direction like an asterisk.
At some fractional distance of each segment, add two more segments, splitting off from the first at some defined angle, in either direction.
Repeat the previous step n times, with each iteration of new segments splitting off from the segments that were generated during the previous iteration.

Ok! So let’s build a function do do it. We’ll make it fairly flexible, to allow the user to control the number of sides it starts with, the angle each split makes from the previous segment, the fractional distance each split occurs at, and, of course, the number of iterations.

Having made this plan, here’s a skeleton of the function. We’re going to leave the number of iterations unspecified, but will give the number of sides a default value of 6, the angle a default value of 60 degrees, and each split occurring halfway down the previous segment.

snowflake <- function(n, nsides=6, angle=60, wherebranch=0.5) {
  
}

The algorithm above had two main components: an initial step setting things up, and a recursive step to be repeated. Let’s start with the initial step. First, we’ll make an empty plot. The plot() call below makes an empty plot from -1 to 1 in the x- and y-directions, with no axes labels and no axis ticks. Specifying asp=1 forces an aspect ratio of 1, meaning that shapes we draw will not be stretched in the x- or y-directions.

We’ll plan on using the segments() function to draw the line segments. This takes vectors of x- and y-coordinates for the start and end of each line segment, so we know we need four vectors in all. The starting coordinate vectors are easy: since we’ll start at (0,0), these are just vectors of zeroes. The ending coordinates can be derived from basic trigonometry if we know the angles the segments make from a starting point. Dividing 360 degrees by the number of sides gives us the angle between segments, and it’s not hard to generate a vector giving the sequence of angles.

Trig functions in R use radians instead of degrees, so we need to convert the between-segment angle to radians by multiplying by pi/180. While we’re at it, we’ll convert the user-specified split angle to radians.

snowflake <- function(n, nsides=6, angle=60, wherebranch=0.5) {
  plot(NA, xlim=c(-1,1), ylim=c(-1,1), xlab="", ylab="", xaxt='n', yaxt='n', asp=1)
  
  xstart <- rep(0,nsides)
  ystart <- rep(0,nsides)
  
  firstangle <- 360/nsides*pi/180
  angle <- angle*pi/180
  
  xend <- sin(firstangle*(1:nsides))
  yend <- cos(firstangle*(1:nsides))
  
  segments(xstart, ystart, xend, yend)
}

Let’s try it.

snowflake()

Seems to work! Now for the fun part: the recursive step. We already did one iteration, so we’ll create a loop that runs from step 2 to n, with each iteration using the information from the previous one to calculate its values.

The segments(xstart, ystart, xend, yend) call worked fine for us to draw the segments, so let’s just use it again. This means we need x- and y- coordinate vectors of the new segments we want to draw for each step. The x- and y-coordinates of the start points are calculated by going the fractional distance specified by wherebranch between the previous starting and ending points.

In each step, we need to generate more new segments than we did last time. We need branch off in either direction from each previous segment, plus keep the original branch active - therefore, we know we need to generate three times as many new segments as last time. Note: we will also need to make sure that the starting and ending coordinate vectors are in the same order!

New lines are highlighted with ##########.

snowflake <- function(n, nsides=6, angle=60, wherebranch=0.5) {
  plot(NA, xlim=c(-1,1), ylim=c(-1,1), xlab="", ylab="", xaxt='n', yaxt='n', asp=1)
  xstart <- rep(0,nsides)
  ystart <- rep(0,nsides)
  firstangle <- 360/nsides*pi/180
  angle <- angle*pi/180
  xend <- sin(firstangle*(1:nsides))
  yend <- cos(firstangle*(1:nsides))
  segments(xstart, ystart, xend, yend)
  
  for(i in 2:n) {                                             ##########
    xstart <- rep(xstart + wherebranch*(xend-xstart),3)       ##########
    ystart <- rep(ystart + wherebranch*(yend-ystart),3)       ##########
  }                                                           ##########
}

Having just defined NEW values for xstart and ystart, we can define values for the ending coordinates by going the right distance in the right direction from the start points. The distance is the same as the distance from the branching location to the previous endpoint, effectively 1-wherebranch multiplied by itself, the number of times the loop has previously repeated. The directions can be calculated by the angles formed by the previous segments, plus the branching angle. Since we repeated the starting coordinate vectors three times, we’ll create coordinate vectors calculated by the previous angles minus the branching angle, then the previous angles themselves, then the previous angles plus the branching angle.

The branching angle gets specified by the user, but we need to make sure we store the angles this iteration uses so that they can be used in the next iteration. And we also need to add a vector of angles that the first step used!

snowflake <- function(n, nsides=6, angle=60, wherebranch=0.5) {
  plot(NA, xlim=c(-1,1), ylim=c(-1,1), xlab="", ylab="", xaxt='n', yaxt='n', asp=1)
  xstart <- rep(0,nsides)
  ystart <- rep(0,nsides)
  firstangle <- 360/nsides*pi/180
  angle <- angle*pi/180
  xend <- sin(firstangle*(1:nsides))
  yend <- cos(firstangle*(1:nsides))
  
  last_ang <- firstangle*(1:nsides)                                                     ########
  
  segments(xstart, ystart, xend, yend)
  for(i in 2:n) {
    end_l <- (1-wherebranch)^(i-1)                                                      ########
    
    xstart <- rep(xstart + wherebranch*(xend-xstart),3)
    ystart <- rep(ystart + wherebranch*(yend-ystart),3)
    
    xend <- xstart + end_l*c(sin(last_ang-angle), sin(last_ang), sin(last_ang+angle))   ########
    yend <- ystart + end_l*c(cos(last_ang-angle), cos(last_ang), cos(last_ang+angle))   ########
    
    last_ang <- c(last_ang-angle, last_ang, last_ang+angle)                             ########
    
    segments(xstart, ystart, xend, yend)                                                ########
  }
}

Let’s try it out.

snowflake(9, wherebranch=.55)

Here it is again with some additional pieces, an logical color argument that prints lines on a rainbow color ramp, and a ... so you can pass in additional plotting arguments. Feel free to add some festive holiday cheer to your office.

snowflake <- function(n, nsides=6, angle=60, wherebranch=0.5, color=T, ...) {
  plot(NA, xlim=c(-1,1), ylim=c(-1,1), xlab="", ylab="", xaxt='n', yaxt='n', asp=1, ...=...)
  xstart <- rep(0,nsides)
  ystart <- rep(0,nsides)
  firstangle <- 360/nsides*pi/180
  angle <- angle*pi/180
  xend <- sin(firstangle*(1:nsides))
  yend <- cos(firstangle*(1:nsides))
  last_ang <- firstangle*(1:nsides)
  if(color) cols <- rainbow(n)
  else cols <- rep(1,n)
  segments(xstart, ystart, xend, yend, col=cols[1])
  for(i in 2:n) {
    end_l <- (1-wherebranch)^(i-1)
    xstart <- rep(xstart + wherebranch*(xend-xstart),3)
    ystart <- rep(ystart + wherebranch*(yend-ystart),3)
    xend <- xstart + end_l*c(sin(last_ang-angle), sin(last_ang), sin(last_ang+angle))
    yend <- ystart + end_l*c(cos(last_ang-angle), cos(last_ang), cos(last_ang+angle))
    last_ang <- c(last_ang-angle, last_ang, last_ang+angle)
    segments(xstart, ystart, xend, yend, col=cols[i])
  }
}

snowflake(9, angle=70, wherebranch=.45, main="Merry Christmas!")

Check this out in webapp form, at https://mbtyers.shinyapps.io/snowflake/!

Fractal Snowflakes: A Recursive Algorithm

Matt Tyers