Terrain Generation with Midpoint Displacement

Posted on February 19th, 2016.

I'm taking the Game Engine Architecture class at Reykjavík University. My group just finished our midterm project where we played around with procedural terrain generation in Unity. It was a lot of fun and I really enjoyed creating terrain out of numbers, so I thought I'd write up a little introduction to a few of the algorithms.

In this post we'll be looking at Midpoint Displacement.

The full series of posts so far:

  1. Terrain Generation
  2. Resources, Code, and Examples
    1. Wisp
    2. Three.js
  3. Heightmaps
    1. Overview
    2. Code
  4. Random Noise
  5. Midpoint Displacement
    1. Initialize the Corners
    2. Set the Edges
    3. Set the Center
    4. Iterate
  6. Result

Terrain Generation

Procedural generation is a pretty big field. The Wikipedia article gives a pretty brief overview. The Procedural Content Generation Wiki is an entire wiki devoted to procedural generation. There's also a small but active subreddit for interested folks.

Today we're just going to look at a single algorithm for generating realistic-looking terrain called "Midpoint Displacement". It's relatively simple (compared to other algorithms) but it produces terrain that actually looks kind of cool.

Resources, Code, and Examples

There are a lot of resources for learning how terrain generation algorithms work. For midpoint displacement I used the Wikipedia article, the PCG Wiki article, and academic papers like this one.

Unfortunately, while there are a lot of places that describe algorithms like this, there are relatively few that explain it thoroughly. Academic papers in particular tend to have a really bad case of draw-the-rest-of-the-fucking-owl-syndrome. They tend to slap an equation or two on the page and move on without showing any code.

As anyone who's done much programming knows, there's a big difference between saying something like "set the midpoints of the diamonds to the average of the corners" and making a computer actually do that. How do we iterate over the array (row-major or column-major)? Does it matter (yes: if you care about cache coherency)? What about the edges where there are fewer corners? How random (and how fast) is our random number generator (and how much do we care)? How do we write all this code so we can actually read it six months later?

So I'm going to try to do my part to improve the situation by not just describing the algorithm, but showing actual code that runs and generates a terrain.

The full code is here, but we'll be going through the important parts below.

Wisp

I'm using a language called Wisp. It's a small, Clojure-inspired dialect of Lisp that compiles down to vanilla Javascript.

I went with Javascript because it means I can actually put demos inline in this post. Sometimes it's so much easier to understand something when you can actually see it running.

I've tried to write the code so that it's readable even if you don't know Lisp. The point of me showing you this code is not to give you something to copy and paste. The point is that until you actually write and run the code, you don't know all the edge cases and problems that can happen. So by making examples that actually run I feel a bit more confident that I'm explaining things enough.

I haven't focused on performance at all. If you want a fast implementation of this, you'll probably want to use a language with real arrays (i.e. contiguous chunks of memory holding a bunch of unboxed floats or integers).

Three.js

Three.js is a Javascript library for 3D rendering. I'm using it here to get something on the screen you can play with. I'm not going to cover the code to put the terrains on the screen because it's really an entirely separate problem to generating the terrain, and it's mostly uninteresting boilerplate for these little demos.

I haven't been able to get these demos working on my iPhone. Sorry about that, but I'm not much of a frontend developer. You'll have to use a computer to see them. Pull requests are welcome.

Heightmaps

Let's get started. The main data structure we're going to be using is a heightmap.

Overview

Heightmaps are essentially just big two-dimensional arrays of numbers in a given range that represent the height of the ground at a given point. For example, a heightmap with a small "mountain" in the middle might look like:

[[0, 0, 2, 0, 0],
 [0, 2, 5, 1, 0],
 [2, 5, 9, 4, 1],
 [0, 1, 4, 1, 0],
 [0, 0, 1, 0, 0]]

In practice, heightmaps are often represented as a single-dimensional array instead to ensure that all the data is together in one big hunk of memory:

[0, 0, 2, 0, 0,
 0, 2, 5, 1, 0,
 2, 5, 9, 4, 1,
 0, 1, 4, 1, 0,
 0, 0, 1, 0, 0]

The type and range of the numbers in the map varies by program and programmer. Some programs use nonzero integers, some use all the integers, some use floating point, etc.

The heightmaps in this post will contain floats from 0.0 to 1.0. I like using floats between zero and one because it's easy to roughly visualize them in your head (e.g. 0.1 is ten percent of the maximum height).

One common way to visualize heightmaps is by turning them into greyscale images with one pixel per number, with black pixels for lowest value in the map's range and white for the highest. The Wikipedia article has an example of this. This is handy because you can also write some code to read images, and then you can use any image editor to draw your own heightmaps by hand.

We're going to visualize the heightmaps in 3D. It's a lot more fun, and they start to actually look like real terrains.

Code

We're going to start by writing a few functions that will hide away the internal details of how our heightmaps are implemented, so we can talk about them in nice high-level terms for the rest of the program. We'll also need a couple of helper functions along the way.

We'll start with something to create a heightmap:

(defn make-heightmap [exponent]
  (let [resolution (+ 1 (Math.pow 2 exponent))]
    (l (+ "Creating " resolution " by " resolution " heightmap..."))
    (def heightmap
      (new Array (* resolution resolution)))
    (set! heightmap.resolution resolution)
    (set! heightmap.exponent exponent)
    (set! heightmap.last (- resolution 1))
    (zero-heightmap heightmap)))

make-heightmap takes an exponent and creates a heightmap. For reasons that will become clear later, all of our heightmaps must be square, and they'll need to have \(2^n + 1\) rows and columns. This means we can make heightmaps of sizes like 3x3, 5x5, 9x9, 17x17, 33x33, etc. So our function just takes the \(n\) to use and creates a correctly-sized array.

(If you've ever poked around with Unity's Terrain objects you might recognize these "powers of two plus one".)

We create a new array (1-dimensional) and store a few pieces of extra data on it for later use. last is the index of the last row/column, so for a 3x3 array it would be 2.

l is just a simple little logging function:

(defn l [v]
  (console.log v))

We'll need to zero out the heightmap as well:

(defn zero-heightmap [heightmap]
  (do-times i heightmap.length
    (set! (aget heightmap i) 0.0))
  heightmap)

zero-heightmap just iterates from 0 to heightmap.length and sets each element to 0.0. Wisp doesn't have syntax to loop from 0 to some number, but it's a Lisp, so we can make some:

(defmacro when [condition & body]
  `(if ~condition
     (do ~@body)))

(defmacro do-times [varname limit & body]
  (let [end (gensym)]
    `(let [~end ~limit]
       (loop [~varname 0]
         (when (< ~varname ~end)
           ~@body
           (recur (+ 1 ~varname)))))))

If you're not a Lisp person, don't sweat this too much. Just trust me that (do-times i 10 (console.log i)) will print the numbers 0 to 9.

Now that we can create a heightmap, we'll probably want to be able to read its elements:

(defmacro heightmap-get [hm x y]
  `(aget ~hm (+ (* ~y (.-resolution ~hm)) ~x)))

(defn heightmap-get-safe [hm x y]
  (when (and (<= 0 x hm.last)
             (<= 0 y hm.last))
    (heightmap-get hm x y)))

heightmap-get handles the nasty business of translating our human-thinkable x and y coordinates into an index into the single-dimensional array. It's not particularly fast (especially if the compiler won't inline it), but we're not worried about speed for this post.

heightmap-get-safe is a version of heightmap-get that will do bounds checking for us and return nil if we ask for something out of range. Javascript will return undefined if you ask for a non-existent index, but because we're storing the array as one big line of data, using a y value that's too big will actually "wrap around" to the next row, so it's safer to just check it explicitly.

Of course we'll also want to be able to set values in our heightmap:

(defmacro heightmap-set! [hm x y val]
  `(set! (heightmap-get ~hm ~x ~y) ~val))

And finally, we'll want to be able to "normalize" our array. During the course of running our algorithm, it's possible for the values to become less than 0.0 or greater than 1.0. We want our heightmaps to contain only values between zero and one, so rather than change the algorithm we can just loop through the array at the end and "squeeze" the range down to 0 and 1 (or "stretch" it if it ends up smaller):

(defn normalize [hm]
  (let [max (- Infinity)
        min Infinity]
    (do-times i hm.length
      (let [el (aget hm i)]
        (when (< max el) (set! max el))
        (when (> min el) (set! min el))))
    (let [span (- max min)]
      (do-times i hm.length
        (set! (aget hm i)
          (/ (- (aget hm i) min)
             span))))))

Random Noise

Now that we've got a nice little interface to heightmaps we can move on to actually making some terrain! Before we dive into midpoint displacement, let's get a black triangle on the screen. We'll start by just setting every point in the heightmap to a completely random value.

(defn rand []
  (Math.random))

(defn rand-around-zero [spread]
  (- (* spread (rand) 2) spread))

(defn random-noise [heightmap]
  (do-times i heightmap.length
    (set! (aget heightmap i) (rand))))

Let's see it! Hit the "refresh" button to create and render the terrain. You can keep hitting it to get fresh terrains.

Not really a very fun landscape to explore, but at least we've got something running. Onward to the real algorithm.

Midpoint Displacement

Midpoint displacement is a simple, fast way to generate decent-looking terrain. It's not perfect (you'll notice patterns once you start viewing the mesh from different angles) but it's a good place to start.

The general process goes like this:

  1. Initialize the four corners of the heightmap to random values.
  2. Set the midpoints of each edge to the average of the two corners it's between, plus or minus a random amount.
  3. Set the center of the square to the average of those edge midpoints you just set, again with a random jitter.
  4. Recurse on the four squares inside this one, reducing the jitter.

That's the standard description. Now let's draw the fucking owl.

Initialize the Corners

The first step is pretty easy: we just need to take our heightmap:

               Column
          0   1   2   3   4
        ┌───┬───┬───┬───┬───┐
      0 │   │   │   │   │   │
        ├───┼───┼───┼───┼───┤
      1 │   │   │   │   │   │
    R   ├───┼───┼───┼───┼───┤
    o 2 │   │   │   │   │   │
    w   ├───┼───┼───┼───┼───┤
      3 │   │   │   │   │   │
        ├───┼───┼───┼───┼───┤
      4 │   │   │   │   │   │
        └───┴───┴───┴───┴───┘

Find the corners:

               Column
          0   1   2   3   4
        ┌───┬───┬───┬───┬───┐
      0 │ ◉ │   │   │   │ ◉ │
        ├───┼───┼───┼───┼───┤
      1 │   │   │   │   │   │
    R   ├───┼───┼───┼───┼───┤
    o 2 │   │   │   │   │   │
    w   ├───┼───┼───┼───┼───┤
      3 │   │   │   │   │   │
        ├───┼───┼───┼───┼───┤
      4 │ ◉ │   │   │   │ ◉ │
        └───┴───┴───┴───┴───┘

And shove random values in them:

               Column
          0   1   2   3   4
        ┌───┬───┬───┬───┬───┐
      0 │ 2 │   │   │   │ 8 │
        ├───┼───┼───┼───┼───┤
      1 │   │   │   │   │   │
    R   ├───┼───┼───┼───┼───┤
    o 2 │   │   │   │   │   │
    w   ├───┼───┼───┼───┼───┤
      3 │   │   │   │   │   │
        ├───┼───┼───┼───┼───┤
      4 │ 0 │   │   │   │ 3 │
        └───┴───┴───┴───┴───┘

(I'm using 0 to 10 instead of 0 to 1 because it's easier to read in the ASCII-art diagrams.)

The code is about what you'd expect:

(defn mpd-init-corners [heightmap]
  (heightmap-set! heightmap 0 0 (rand))
  (heightmap-set! heightmap 0 heightmap.last (rand))
  (heightmap-set! heightmap heightmap.last 0 (rand))
  (heightmap-set! heightmap heightmap.last heightmap.last (rand)))

(defn midpoint-displacement-d1 [heightmap]
  (mpd-init-corners heightmap))

That was easy.

Set the Edges

For the next step, we need to take our square heightmap with filled-in corners and use the corners to fill in the middle element of each edge:

               Column
          0   1   2   3   4
        ┌───┬───┬───┬───┬───┐
      0 │ 2 │   │ ◉ │   │ 8 │
        ├───┼───┼───┼───┼───┤
      1 │   │   │   │   │   │
    R   ├───┼───┼───┼───┼───┤
    o 2 │ ◉ │   │   │   │ ◉ │
    w   ├───┼───┼───┼───┼───┤
      3 │   │   │   │   │   │
        ├───┼───┼───┼───┼───┤
      4 │ 0 │   │ ◉ │   │ 3 │
        └───┴───┴───┴───┴───┘



               Column
          0   1   2   3   4
        ┌───┬───┬───┬───┬───┐
      0 │ 2 ├───▶ ◉ ◀───┤ 8 │
        ├─┬─┼───┼───┼───┼─┬─┤
      1 │ │ │   │   │   │ │ │
    R   ├─▼─┼───┼───┼───┼─▼─┤
    o 2 │ ◉ │   │   │   │ ◉ │
    w   ├─▲─┼───┼───┼───┼─▲─┤
      3 │ │ │   │   │   │ │ │
        ├─┴─┼───┼───┼───┼─┴─┤
      4 │ 0 ├───▶ ◉ ◀───┤ 3 │
        └───┴───┴───┴───┴───┘


               Column
          0   1   2   3   4
        ┌───┬───┬───┬───┬───┐
      0 │ 2 ├───▶ 5 ◀───┤ 8 │
        ├─┬─┼───┼───┼───┼─┬─┤
      1 │ │ │   │   │   │ │ │
    R   ├─▼─┼───┼───┼───┼─▼─┤
    o 2 │ 1 │   │   │   │5.5│
    w   ├─▲─┼───┼───┼───┼─▲─┤
      3 │ │ │   │   │   │ │ │
        ├─┴─┼───┼───┼───┼─┴─┤
      4 │ 0 ├───▶1.5◀───┤ 3 │
        └───┴───┴───┴───┴───┘

This is why the size of the heightmap has to be \(2^n + 1\) — the \(+ 1\) ensures that the sides are odd lengths, which makes sure they have a middle element for us to set.

Let's make a couple of helper functions:

(defn jitter [value spread]
  (+ value (rand-around-zero spread)))

(defn midpoint [a b]
  (/ (+ a b) 2))

(defn average2 [a b]
  (/ (+ a b) 2))

(defn average4 [a b c d]
  (/ (+ a b c d) 4))

jitter is just a nice way of saying "take this point and move it a random amount up or down".

The other three functions are pretty self-explanatory. Also, two of them are exactly the same function! There's a reason for this.

Lately I've started digging into SICP again and watching the lectures. Abelson and Sussman practice a style of programming that's more about "growing" a language to talk about a problem you're trying to solve. So you start with the base programming language, and then define new terms/syntax that let you talk about your problem at the appropriate level of detail.

So in this case, we've got two functions for conceptually different tasks:

They happen to be implemented in the same way, but I think giving them different names makes them easier to talk about and easier to keep straight in your head.

(defn mpd-displace [heightmap lx rx by ty spread]
  (let [cx (midpoint lx rx)
        cy (midpoint by ty)

        bottom-left (heightmap-get heightmap lx by)
        bottom-right (heightmap-get heightmap rx by)
        top-left (heightmap-get heightmap lx ty)
        top-right (heightmap-get heightmap rx ty)

        top (average2 top-left top-right)
        left (average2 bottom-left top-left)
        bottom (average2 bottom-left bottom-right)
        right (average2 bottom-right top-right)]
    (heightmap-set! heightmap cx by (jitter bottom spread))
    (heightmap-set! heightmap cx ty (jitter top spread))
    (heightmap-set! heightmap lx cy (jitter left spread))
    (heightmap-set! heightmap rx cy (jitter right spread))))

(defn midpoint-displacement-d2 [heightmap]
  (mpd-init-corners heightmap)
  (mpd-displace heightmap
                0 heightmap.last
                0 heightmap.last
                0.1))

mpd-displace is the meat of the algorithm. It takes a heightmap, the spread (how much random "jitter" to apply after averaging the values), and four values to define the square we're working on in the heightmap:

For the first iteration we want to work on the entire heightmap, so we pass in 0 and heightmap.last for the left/right and top/bottom indices.

Run the demo a few times and watch how the midpoints fall between the corner points because they're averaged (with a little bit of jitter thrown in).

Set the Center

The next step is to set the center element of the square to the average of those midpoints we just set, plus a bit of jitter:

               Column
          0   1   2   3   4
        ┌───┬───┬───┬───┬───┐
      0 │ 2 │   │ 5 │   │ 8 │
        ├───┼───┼───┼───┼───┤
      1 │   │   │   │   │   │
    R   ├───┼───┼───┼───┼───┤
    o 2 │ 1 │   │ ◉ │   │5.5│
    w   ├───┼───┼───┼───┼───┤
      3 │   │   │   │   │   │
        ├───┼───┼───┼───┼───┤
      4 │ 0 │   │1.5│   │ 3 │
        └───┴───┴───┴───┴───┘


               Column
          0   1   2   3   4
        ┌───┬───┬───┬───┬───┐
      0 │ 2 │   │ 5 │   │ 8 │
        ├───┼───┼─┬─┼───┼───┤
      1 │   │   │ │ │   │   │
    R   ├───┼───┼─▼─┼───┼───┤
    o 2 │ 1 ├───▶ ◉ ◀───┤5.5│
    w   ├───┼───┼─▲─┼───┼───┤
      3 │   │   │ │ │   │   │
        ├───┼───┼─┴─┼───┼───┤
      4 │ 0 │   │1.5│   │ 3 │
        └───┴───┴───┴───┴───┘


               Column
          0   1   2   3   4
        ┌───┬───┬───┬───┬───┐
      0 │ 2 │   │ 5 │   │ 8 │
        ├───┼───┼─┬─┼───┼───┤
      1 │   │   │ │ │   │   │
    R   ├───┼───┼─▼─┼───┼───┤
    o 2 │ 1 ├───▶3.3◀───┤5.5│
    w   ├───┼───┼─▲─┼───┼───┤
      3 │   │   │ │ │   │   │
        ├───┼───┼─┴─┼───┼───┤
      4 │ 0 │   │1.5│   │ 3 │
        └───┴───┴───┴───┴───┘

We can add this to the displacement function pretty easily:

(defn mpd-displace [heightmap lx rx by ty spread]
  (let [cx (midpoint lx rx)
        cy (midpoint by ty)

        bottom-left (heightmap-get heightmap lx by)
        bottom-right (heightmap-get heightmap rx by)
        top-left (heightmap-get heightmap lx ty)
        top-right (heightmap-get heightmap rx ty)

        top (average2 top-left top-right)
        left (average2 bottom-left top-left)
        bottom (average2 bottom-left bottom-right)
        right (average2 bottom-right top-right)
        center (average4 top left bottom right)]
    (heightmap-set! heightmap cx by (jitter bottom spread))
    (heightmap-set! heightmap cx ty (jitter top spread))
    (heightmap-set! heightmap lx cy (jitter left spread))
    (heightmap-set! heightmap rx cy (jitter right spread))
    (heightmap-set! heightmap cx cy (jitter center spread))))

And now our center is ready to go:

Iterate

To finish things off, we need to iterate so that we can fill in the rest of the values in progressively smaller squares.

A 3x3 square is the base case, because its corners will be filled in to start, and the edges and midpoint get filled in. We just call our displacement function with the top/bottom/left/right values defining the square and we're done.

A 5x5 starts with one iteration that displaces the whole map. Then it does another "pass" that displaces each of the 3x3 "sub-squares":

    ╔═══════════╗───┬───┐      ┌───┬───╔═══════════╗
    ║ 2 │   │ 5 ║   │ 8 │      │ 2 │   ║ 5 │   │ 8 ║
    ║───┼───┼───║───┼───┤      ├───┼───║───┼───┼───║
    ║   │   │   ║   │   │      │   │   ║   │   │   ║
    ║───┼───┼───║───┼───┤      ├───┼───║───┼───┼───║
    ║ 1 │   │3.3║   │5.5│      │ 1 │   ║3.3│   │5.5║
    ╚═══════════╝───┼───┤      ├───┼───╚═══════════╝
    │   │   │   │   │   │      │   │   │   │   │   │
    ├───┼───┼───┼───┼───┤      ├───┼───┼───┼───┼───┤
    │ 0 │   │1.5│   │ 3 │      │ 0 │   │1.5│   │ 3 │
    └───┴───┴───┴───┴───┘      └───┴───┴───┴───┴───┘


    ┌───┬───┬───┬───┬───┐      ┌───┬───┬───┬───┬───┐
    │ 2 │   │ 5 │   │ 8 │      │ 2 │   │ 5 │   │ 8 │
    ├───┼───┼───┼───┼───┤      ├───┼───┼───┼───┼───┤
    │   │   │   │   │   │      │   │   │   │   │   │
    ╔═══════════╗───┼───┤      ├───┼───╔═══════════╗
    ║ 1 │   │3.3║   │5.5│      │ 1 │   ║3.3│   │5.5║
    ║───┼───┼───║───┼───┤      ├───┼───║───┼───┼───║
    ║   │   │   ║   │   │      │   │   ║   │   │   ║
    ║───┼───┼───║───┼───┤      ├───┼───║───┼───┼───║
    ║ 0 │   │1.5║   │ 3 │      │ 0 │   ║1.5│   │ 3 ║
    ╚═══════════╝───┴───┘      └───┴───╚═══════════╝

If we have a bigger heightmap (e.g. 9x9) we'll need three passes:

     Pass 1
    ▩▩▩▩▩▩▩▩▩
    ▩▩▩▩▩▩▩▩▩
    ▩▩▩▩▩▩▩▩▩
    ▩▩▩▩▩▩▩▩▩
    ▩▩▩▩▩▩▩▩▩
    ▩▩▩▩▩▩▩▩▩
    ▩▩▩▩▩▩▩▩▩
    ▩▩▩▩▩▩▩▩▩
    ▩▩▩▩▩▩▩▩▩


     Pass 2
    ▩▩▩▩▩▢▢▢▢  ▢▢▢▢▩▩▩▩▩  ▢▢▢▢▢▢▢▢▢  ▢▢▢▢▢▢▢▢▢
    ▩▩▩▩▩▢▢▢▢  ▢▢▢▢▩▩▩▩▩  ▢▢▢▢▢▢▢▢▢  ▢▢▢▢▢▢▢▢▢
    ▩▩▩▩▩▢▢▢▢  ▢▢▢▢▩▩▩▩▩  ▢▢▢▢▢▢▢▢▢  ▢▢▢▢▢▢▢▢▢
    ▩▩▩▩▩▢▢▢▢  ▢▢▢▢▩▩▩▩▩  ▢▢▢▢▢▢▢▢▢  ▢▢▢▢▢▢▢▢▢
    ▩▩▩▩▩▢▢▢▢  ▢▢▢▢▩▩▩▩▩  ▩▩▩▩▩▢▢▢▢  ▢▢▢▢▩▩▩▩▩
    ▢▢▢▢▢▢▢▢▢  ▢▢▢▢▢▢▢▢▢  ▩▩▩▩▩▢▢▢▢  ▢▢▢▢▩▩▩▩▩
    ▢▢▢▢▢▢▢▢▢  ▢▢▢▢▢▢▢▢▢  ▩▩▩▩▩▢▢▢▢  ▢▢▢▢▩▩▩▩▩
    ▢▢▢▢▢▢▢▢▢  ▢▢▢▢▢▢▢▢▢  ▩▩▩▩▩▢▢▢▢  ▢▢▢▢▩▩▩▩▩
    ▢▢▢▢▢▢▢▢▢  ▢▢▢▢▢▢▢▢▢  ▩▩▩▩▩▢▢▢▢  ▢▢▢▢▩▩▩▩▩


     Pass 3
    ▩▩▩▢▢▢▢▢▢  ▢▢▩▩▩▢▢▢▢  ▢▢▢▢▩▩▩▢▢  ▢▢▢▢▢▢▩▩▩
    ▩▩▩▢▢▢▢▢▢  ▢▢▩▩▩▢▢▢▢  ▢▢▢▢▩▩▩▢▢  ▢▢▢▢▢▢▩▩▩
    ▩▩▩▢▢▢▢▢▢  ▢▢▩▩▩▢▢▢▢  ▢▢▢▢▩▩▩▢▢  ▢▢▢▢▢▢▩▩▩
    ▢▢▢▢▢▢▢▢▢  ▢▢▢▢▢▢▢▢▢  ▢▢▢▢▢▢▢▢▢  ▢▢▢▢▢▢▢▢▢
    ▢▢▢▢▢▢▢▢▢  ▢▢▢▢▢▢▢▢▢  ▢▢▢▢▢▢▢▢▢  ▢▢▢▢▢▢▢▢▢
    ▢▢▢▢▢▢▢▢▢  ▢▢▢▢▢▢▢▢▢  ▢▢▢▢▢▢▢▢▢  ▢▢▢▢▢▢▢▢▢
    ▢▢▢▢▢▢▢▢▢  ▢▢▢▢▢▢▢▢▢  ▢▢▢▢▢▢▢▢▢  ▢▢▢▢▢▢▢▢▢
    ▢▢▢▢▢▢▢▢▢  ▢▢▢▢▢▢▢▢▢  ▢▢▢▢▢▢▢▢▢  ▢▢▢▢▢▢▢▢▢
    ▢▢▢▢▢▢▢▢▢  ▢▢▢▢▢▢▢▢▢  ▢▢▢▢▢▢▢▢▢  ▢▢▢▢▢▢▢▢▢

    ▢▢▢▢▢▢▢▢▢  ▢▢▢▢▢▢▢▢▢  ▢▢▢▢▢▢▢▢▢  ▢▢▢▢▢▢▢▢▢
    ▢▢▢▢▢▢▢▢▢  ▢▢▢▢▢▢▢▢▢  ▢▢▢▢▢▢▢▢▢  ▢▢▢▢▢▢▢▢▢
    ▩▩▩▢▢▢▢▢▢  ▢▢▩▩▩▢▢▢▢  ▢▢▢▢▩▩▩▢▢  ▢▢▢▢▢▢▩▩▩
    ▩▩▩▢▢▢▢▢▢  ▢▢▩▩▩▢▢▢▢  ▢▢▢▢▩▩▩▢▢  ▢▢▢▢▢▢▩▩▩
    ▩▩▩▢▢▢▢▢▢  ▢▢▩▩▩▢▢▢▢  ▢▢▢▢▩▩▩▢▢  ▢▢▢▢▢▢▩▩▩
    ▢▢▢▢▢▢▢▢▢  ▢▢▢▢▢▢▢▢▢  ▢▢▢▢▢▢▢▢▢  ▢▢▢▢▢▢▢▢▢
    ▢▢▢▢▢▢▢▢▢  ▢▢▢▢▢▢▢▢▢  ▢▢▢▢▢▢▢▢▢  ▢▢▢▢▢▢▢▢▢
    ▢▢▢▢▢▢▢▢▢  ▢▢▢▢▢▢▢▢▢  ▢▢▢▢▢▢▢▢▢  ▢▢▢▢▢▢▢▢▢
    ▢▢▢▢▢▢▢▢▢  ▢▢▢▢▢▢▢▢▢  ▢▢▢▢▢▢▢▢▢  ▢▢▢▢▢▢▢▢▢

    ▢▢▢▢▢▢▢▢▢  ▢▢▢▢▢▢▢▢▢  ▢▢▢▢▢▢▢▢▢  ▢▢▢▢▢▢▢▢▢
    ▢▢▢▢▢▢▢▢▢  ▢▢▢▢▢▢▢▢▢  ▢▢▢▢▢▢▢▢▢  ▢▢▢▢▢▢▢▢▢
    ▢▢▢▢▢▢▢▢▢  ▢▢▢▢▢▢▢▢▢  ▢▢▢▢▢▢▢▢▢  ▢▢▢▢▢▢▢▢▢
    ▢▢▢▢▢▢▢▢▢  ▢▢▢▢▢▢▢▢▢  ▢▢▢▢▢▢▢▢▢  ▢▢▢▢▢▢▢▢▢
    ▩▩▩▢▢▢▢▢▢  ▢▢▩▩▩▢▢▢▢  ▢▢▢▢▩▩▩▢▢  ▢▢▢▢▢▢▩▩▩
    ▩▩▩▢▢▢▢▢▢  ▢▢▩▩▩▢▢▢▢  ▢▢▢▢▩▩▩▢▢  ▢▢▢▢▢▢▩▩▩
    ▩▩▩▢▢▢▢▢▢  ▢▢▩▩▩▢▢▢▢  ▢▢▢▢▩▩▩▢▢  ▢▢▢▢▢▢▩▩▩
    ▢▢▢▢▢▢▢▢▢  ▢▢▢▢▢▢▢▢▢  ▢▢▢▢▢▢▢▢▢  ▢▢▢▢▢▢▢▢▢
    ▢▢▢▢▢▢▢▢▢  ▢▢▢▢▢▢▢▢▢  ▢▢▢▢▢▢▢▢▢  ▢▢▢▢▢▢▢▢▢

    ▢▢▢▢▢▢▢▢▢  ▢▢▢▢▢▢▢▢▢  ▢▢▢▢▢▢▢▢▢  ▢▢▢▢▢▢▢▢▢
    ▢▢▢▢▢▢▢▢▢  ▢▢▢▢▢▢▢▢▢  ▢▢▢▢▢▢▢▢▢  ▢▢▢▢▢▢▢▢▢
    ▢▢▢▢▢▢▢▢▢  ▢▢▢▢▢▢▢▢▢  ▢▢▢▢▢▢▢▢▢  ▢▢▢▢▢▢▢▢▢
    ▢▢▢▢▢▢▢▢▢  ▢▢▢▢▢▢▢▢▢  ▢▢▢▢▢▢▢▢▢  ▢▢▢▢▢▢▢▢▢
    ▢▢▢▢▢▢▢▢▢  ▢▢▢▢▢▢▢▢▢  ▢▢▢▢▢▢▢▢▢  ▢▢▢▢▢▢▢▢▢
    ▢▢▢▢▢▢▢▢▢  ▢▢▢▢▢▢▢▢▢  ▢▢▢▢▢▢▢▢▢  ▢▢▢▢▢▢▢▢▢
    ▩▩▩▢▢▢▢▢▢  ▢▢▩▩▩▢▢▢▢  ▢▢▢▢▩▩▩▢▢  ▢▢▢▢▢▢▩▩▩
    ▩▩▩▢▢▢▢▢▢  ▢▢▩▩▩▢▢▢▢  ▢▢▢▢▩▩▩▢▢  ▢▢▢▢▢▢▩▩▩
    ▩▩▩▢▢▢▢▢▢  ▢▢▩▩▩▢▢▢▢  ▢▢▢▢▩▩▩▢▢  ▢▢▢▢▢▢▩▩▩

Let's think about what's happening here.

This 9x9 heightmap was made with an exponent of 3 (\(2^3 + 1 = 9\)) and we need 3 "passes" to finish it. This isn't a coincidence — it's the reason we chose our heightmap resolutions to be \(2^n + 1\).

So at iteration \(i\), we need to displace a \(2^i\)x\(2^i\) collection of squares. Let's go:

(defn midpoint-displacement [heightmap]
  (mpd-init-corners heightmap)
  (loop [iter 0
         spread 0.3]
    (when (< iter heightmap.exponent)
      (let [chunks (Math.pow 2 iter)
            chunk-width (/ (- heightmap.resolution 1) chunks)]
        (do-nested xchunk ychunk chunks
          (let [left-x (* chunk-width xchunk)
                right-x (+ left-x chunk-width)
                bottom-y (* chunk-width ychunk)
                top-y (+ bottom-y chunk-width)]
                (mpd-displace heightmap
                              left-x right-x bottom-y top-y
                              spread))))
      (recur (+ 1 iter) (* spread 0.5))))
  (normalize heightmap))

That's a lot to take in. Let's break it down:

(defn midpoint-displacement [heightmap]
  (mpd-init-corners heightmap)
  (loop [iter 0
         spread 0.3]
    (when (< iter heightmap.exponent)
      ; Process the chunks at this iteration
      (recur (+ 1 iter) (* spread 0.5))))
  (normalize heightmap))

We initialize the corners at the beginning and normalize the heightmap at the end. Then we loop a number of times equal to the exponent we used to make the heightmap, as described above. Each time through the loop we increment iter and cut the spread for those points in half.

For each pass, we figure out what we're going to need to do:

(let [chunks (Math.pow 2 iter)
      chunk-width (/ (- heightmap.resolution 1) chunks)]
  ; Actually do it...
  )

Let's make a quick "nested for loop" syntax. It'll make things easier to read in the main function and we'll need it for later algorithms too:

(defmacro do-nested [xname yname width & body]
  (let [iterations (gensym)]
    `(let [~iterations ~width]
       (do-times ~xname ~iterations
         (do-times ~yname ~iterations
           ~@body)))))

This lets us say something like (do-nested x y 5 (console.log [x y])) to mean:

for (x = 0; x < 5; x++) {
    for (y = 0; y < 5; y++) {
        console.log([x, y]);
    }
}

And now we can see how we process the chunks:

(let [chunks (Math.pow 2 iter)
      chunk-width (/ (- heightmap.resolution 1) chunks)]
  (do-nested xchunk ychunk chunks
    (let [left-x (* chunk-width xchunk)
          right-x (+ left-x chunk-width)
          bottom-y (* chunk-width ychunk)
          top-y (+ bottom-y chunk-width)]
          (mpd-displace heightmap
                        left-x right-x bottom-y top-y
                        spread))))

We loop over the indices of the chunks, calculate the bounds for each one, and displace it. Simple enough, but a bit fiddly to get right. And now we've finished displacing all the points!

We could have done this recursively, but I was in an iterative mood when I wrote this. Feel free to play around with it.

Result

Now that we've got a working algorithm we can play around with the values and see what effect they have on the resulting terrain.

Warning: don't set the exponent higher than 8 unless you want to crash this browser tab. An exponent of 8 results in a heightmap with 66049 elements. Going to 9 is 263169, which most browsers seem to... dislike.



And that's it, we've grown some mountains!

Now that we've got some of the boilerplate down I'm planning on doing a couple more posts like this looking at other noise algorithms, and maybe some erosion algorithms to simulate weathering. Let me know if you've got any specific requests!