Fractal Landscapes in R

The past few days I’ve gone on a bit of a fractals kick. It’s partly been motivated by the recent passing of Benoit Mandelbrot, the iconoclastic mathematical genius who coined the term “fractal” and profoundly shook up a number of fields by introducing fractal concepts into them. Many of you have no doubt seen the Mandelbrot Set before—it’s probably the most famous mathematical fractal shape out there. If you haven’t, check out this video (via kottke.org):

Mandelbrot Fractal Set Trip To e214 HD from teamfresh on Vimeo.

The day he died, I received a link to his obituary in an email from a certain humorist I know, with the following in the subject line: “founder of fractal geometry dead, stoners everywhere in mourning.”

But aside from psychedelic entertainment value, fractal geometry turns out to be tremendously useful in describing patterns found in nature. It is especially useful compared with Euclidean geometry, the familiar system of straight lines, squares, triangles, circles, etc. As anyone who has ever looked at nature can tell you, nothing is straight. Crookedness, roughness, and convoluted patterns are the norm. Using euclidean geometry, there is no simple way to describe a coastline, or mountain range, or the fluctuations of animal populations. Fractal geometry, on the other hand, allows us to approximate many of these shapes surprisingly well using very simple equations.

One neat fractal trick is generating realistic-looking mountainscapes using only a few simple operations repeated over and over on themselves (generally speaking, repeating an operation over and over on itself is the way all fractal patterns arise). Consider a two-dimensional case first, for simplicity: we’ll generate the profile of a mountain range, as seen on the horizon. Code to generate this pattern is below. I’ve written my examples here in R, since I’m familiar with it’s graphics system and since my mostly ecological audience is more likely to be familiar with this language than anything else (I presume mostly ecological…who are you people, anyway?).

mountainize <- function(a, roughness, sdev) {
    # Recursive function to create a one-dimensional "mountain range."
    # Arguments:
    #    a: numeric vector
    #    roughness: number between 0 and 1, controlling how 
    #    quickly the random perturbations tail off as the scale decreases.
    #    sdev: standard deviation of the random perturbations.
    # Value:
    #    A numeric vector the same length as a.
    size <- length(a)
    if (size < 2) {
        return(a)
    } else {
        middle <- (size + 1) / 2
        a[middle] <- mean(c(a[1], a[size])) + rnorm(1, 0, sdev)
        a[1:middle] <- mountainize(a[1:middle], 
                                    roughness, sdev * roughness)
        a[middle:size] <- mountainize(a[middle:size], 
                                    roughness, sdev * roughness)
        return(a)
    }
}

iterations <- 8  # set the number of times to repeat the algorithm
mountains <- rep(0, 2 ^ iterations + 1)  # make an array of the correct size
set.seed(1)   # delete or change the seed to get other mountain ranges
mountains <- mountainize(mountains, 0.5, 1)
plot(mountains, ty='l')

The function mountainize takes a vector of numbers (in this case, just a bunch of zeros to start) and makes the middle number equal to the mean of the two ends, plus a random number from the normal distribution. It then divides the vector into two parts—from the beginning to the middle, and the middle to the end—and repeats the procedure on each half, but with a smaller random deviation. This recursion continues until the sub-section of the vector is only one element long, at which point it unwinds and returns the whole (now jagged) array. The picture below shows what this looks like step-by-step. With each additional iteration, the number of points doubles. After the seventh or eighth time, it looks pretty good. Voila! A mountain range!

Some technical notes. First, the array has to be an appropriate length for the divide-and-conquer algorithm to work—namely, 2ⁱ + 1, where i is the desired number of iterations (i.e. the depth of recursion). Second, there is an argument to the mountainize function called roughness. This parameter controls how much the standard deviation of the random numbers decreases from iteration to iteration, and hence how rough or smooth the mountain range ends up. Run the code yourself and play around with it to see how it affects the finished product.

I need to go to sleep, but I will leave off for now with a teaser: With just a little more effort, it is possible to use the same principle to create awesome looking three-dimensional terrain maps.

Don’t you want to know how to make something like this? Tune in tomorrow to find out how…