Fractal Landscapes in R: Part Two

My last post showed how, with some simple R code, it is possible to generate a surprisingly realistic mountainous profile, as if seen on the horizon. This was done by taking a two points, calculating their midpoint, offsetting it by a random amount, and then repeating on the two new line segments just created, continuing on for some number of iterations. At each step, the random deviation is reduced, and the length of the line being divided and kinked halves. In this post, I will describe how to apply this algorithm on a rectangular surface, yielding a three-dimensional landscape:

The approach used here is called the diamond-square algorithm, and it is very similar in concept to the line-kinking algorithm. Start with four points arranged in a square in three-dimensional space. Average their values to find the midpoint of the square, and then move this point up or down a random amount. This is the “square step.”

If you draw lines from each corner of the original square to the center, you will now see diamonds (at this stage, they’re actually just half-diamonds) with the midpoints of the original square’s sides at their centers. Set each of these points equal to the average of the diamond’s four corners (for the first iteration, there will only be three points to average instead of four, because the “diamonds” run off the edge of the original square). This is the diamond step. The picture below shows the first iteration through these two steps. New points are in green, and the points they are derived from are connected by dotted lines.

While this algorithm would seem to lend itself to recursion (breaking down the square into smaller and smaller nested squares), this doesn’t actually work that well. The problem is the switch back and forth from averaging the corners of squares to averaging the corners of diamond—they don’t overlap completely, so there’s no good way to pass them down the recursion tree. Instead, we’ll implement the algorithm procedurally (that is, using for loops). The looping is slightly convoluted, since it has to skip certain numbers of cells, and that number has to decrease with each pass through the array. If that’s the sort of bee that will get in your bonnet, you can study the code below (it’s not actually as long as it looks; I put in a lot of comments). Otherwise, just copy and run it to start making awesome fractal landscapes.

# Sam Urmy 10/27/2010
#
# I adapted this implementation of the diamond/square algorithm

# for terrain generation from C++ code by Daniel Beard.
# 
# http://danielbeard.wordpress.com/2010/08/07/terrain-generation-and-smoothing/
# Install/load the rgl package to be able to render the
# terrain in 3-D.  It's worth it.
library(rgl)

mountains <- function(iter, roughness=0.5, m=0, sdev=1) {
  # Function to create mountainous terrain via the 
  # diamond-square algorithm.
  #
  # Arguments:
  # iter: integer. Number of iterations (determines
  #   size of finished terrain matrix).
  # roughness: float, between 0 and 1.  Factor by
  #   which the random deviations are reduced
  #   between iterations.
  # m: optional initial matrix, must be square
  #   with a dimension equal to a power of
  #   two (i.e., 2 ^ number of iterations) plus 1.
  # sdev: optional float. Inital standard deviation
  #   for the random deviations.
  # Value:
  # A square matrix of elevation values.
  
  size <- 2^iter + 1
  # If the user does not supply a matrix,
  # initalize one with zeros.
  if (! m) {
    m <- matrix(0, size, size)
  }
  # Loop through side lengths, starting with the size of the
  # entire matrix, and moving down by factors of 2.
  for (side.length in 2^(iter:1)) {
    half.side <- side.length / 2
    # Square step
    for (col in seq(1, size - 1, by=side.length)) {
      for (row in seq(1, size - 1, by=side.length)) {
        avg <- mean(c(
          m[row, col],        # upper left
          m[row + side.length, col],  # lower left
          m[row, col + side.length],  # upper right
          m[row + side.length, col + side.length] #lower right
        ))
        avg <- avg + rnorm(1, 0, sdev)
        m[row + half.side, col + half.side] <- avg
      }
    }
    
    # Diamond step
    for (row in seq(1, size, by=half.side)) {
      for (col in seq((col+half.side) %% side.length, size, side.length)) {
        # m[row, col] is the center of the diamond
        avg <- mean(c(
          m[(row - half.side + size) %% size, col],# above
          m[(row + half.side) %% size, col], # below
          m[row, (col + half.side) %% size], # right
          m[row, (col - half.side) %% size]  # left
        ))
        m[row, col] <- avg + rnorm(1, 0, sdev)
        # Handle the edges by wrapping around to the
        # other side of the array
        if (row == 0) { m[size - 1, col] = avg }
        if (col == 0) { m[row, size - 1] = avg }
      }
    }
  # Reduce the standard deviation of the random deviation
  # by the roughness factor.
  sdev <- sdev * roughness
  }
  return(m)
}

color.matrix <- function(m, ncolors, palette=terrain.colors) {
  # Takes a matrix, and maps it onto a matrix of the same size
  # in color space.
  # Useful for draping a color map on top of a persp3d map.
  nrows <- dim(m)[1]
  ncols <- dim(m)[2]
  col.index <- m - min(m)
  col.index <- (ncolors - 1) * col.index / max(col.index) + 1
  return(matrix(palette(ncolors)[col.index], nrows, ncols))
}

iter <- 7  # set number of iterations
set.seed(1) # change or delete seed to get different terrains
m <- mountains(iter, .5)

filled.contour(m, col=terrain.colors(36))
# make a color overlay for the 3-D plot
colors <- color.matrix(m, 24)
persp3d(m, aspect=c(1,1,.5), col=colors)

…And that’s it!

It’s worth noting that this is not actually the way landscapes form. The interaction of tectonic uplift and erosion is a bit more complex. But this simple algorithm reproduces mountainous terrain pretty well—at least well enough to fool an untrained eye. The diamond square method found early use in manufacturing background terrain for video and computer games, and it works pretty well for that purpose.

These landscapes, however, do share characteristics with real ones. Most importantly, they are jagged at all spatial scales. Large mountains have crags, and the crags have crags in turn. In the business, this is known as self-similarity, or scale-independent pattern, and it is characteristic of a wide range of natural phenomena and objects. It is this nested roughness that gives them their realistic look.