content/blog/2016/02/midpoint-displacement.markdown @ aa7aebed85bc

Remove the last few remaining Source Hut links
author Steve Losh <steve@stevelosh.com>
date Thu, 23 Jan 2020 00:13:09 -0500
parents f5556130bda1
children (none)
(
:title "Terrain Generation with Midpoint Displacement"
:snip "A first step toward growing worlds with computers."
:date "2016-02-19T19:45:00Z"
:mathjax t
:draft nil

)

<script defer src="/static/js/terrain/jquery.js" type="text/javascript"></script>
<script defer src="/static/js/terrain/three.min.js"></script>
<script defer src="/static/js/terrain/TrackballControls.js"></script>
<script defer src="/static/js/terrain/terrain1.js"></script>

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:

* [Midpoint Displacement](/blog/2016/02/midpoint-displacement/)
* [Recursive Midpoint Displacement](/blog/2016/03/recursive-midpoint-displacement/)
* [Diamond Square](/blog/2016/06/diamond-square/)

[project]: https://www.youtube.com/watch?v=G5u79w4qiAA
[Unity]: http://unity3d.com/

<div id="toc"></div>

## Terrain Generation

Procedural generation is a pretty big field. The [Wikipedia
article][wikipedia-pg] gives a pretty brief overview.  The [Procedural Content
Generation Wiki][pcg-wiki] is an entire wiki devoted to procedural generation.
There's also a small but active [subreddit][r-pg] 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.

[wikipedia-pg]: https://en.wikipedia.org/wiki/Procedural_generation
[pcg-wiki]: http://pcg.wikidot.com/
[r-pg]: https://www.reddit.com/r/proceduralgeneration/

## 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][wikipedia-ds],
the [PCG Wiki article][pcgw-midpoint], and academic papers like [this
one][paper].

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][owl].  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][code], but we'll be going through the important parts
below.

[wikipedia-ds]: https://en.wikipedia.org/wiki/Diamond-square_algorithm
[pcgw-midpoint]: http://pcg.wikidot.com/pcg-algorithm:midpoint-displacement-algorithm
[paper]: http://micsymposium.org/mics_2011_proceedings/mics2011_submission_30.pdf
[owl]: https://i.imgur.com/RadSf.jpg
[code]: https://bitbucket.org/sjl/stevelosh/src/default/static/js/terrain1.wisp

### Wisp

I'm using a language called [Wisp][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).

[wisp]: https://github.com/Gozala/wisp

### Three.js

[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][slc] are welcome.

[slc]: http://bitbucket.org/sjl/stevelosh/
[three.js]: http://threejs.org/

## 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:

```text
[[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:

```text
[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][heightmap] 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.

[heightmap]: https://en.wikipedia.org/wiki/Heightmap

### 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:

```clojure
(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:

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

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

```clojure
(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:

```clojure
(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:

```clojure
(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:

```clojure
(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):

```clojure
(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.

```clojure
(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.

<div id="demo-random" class="threejs"></div>

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

[black triangle]: http://rampantgames.com/blog/?p=7745

## 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:

<pre class="lineart">
               Column
          0   1   2   3   4
        ┌───┬───┬───┬───┬───┐
      0 │   │   │   │   │   │
        ├───┼───┼───┼───┼───┤
      1 │   │   │   │   │   │
    R   ├───┼───┼───┼───┼───┤
    o 2 │   │   │   │   │   │
    w   ├───┼───┼───┼───┼───┤
      3 │   │   │   │   │   │
        ├───┼───┼───┼───┼───┤
      4 │   │   │   │   │   │
        └───┴───┴───┴───┴───┘
</pre>

Find the corners:

<pre class="lineart">
               Column
          0   1   2   3   4
        ┌───┬───┬───┬───┬───┐
      0 │ ◉ │   │   │   │ ◉ │
        ├───┼───┼───┼───┼───┤
      1 │   │   │   │   │   │
    R   ├───┼───┼───┼───┼───┤
    o 2 │   │   │   │   │   │
    w   ├───┼───┼───┼───┼───┤
      3 │   │   │   │   │   │
        ├───┼───┼───┼───┼───┤
      4 │ ◉ │   │   │   │ ◉ │
        └───┴───┴───┴───┴───┘
</pre>

And shove random values in them:

<pre class="lineart">
               Column
          0   1   2   3   4
        ┌───┬───┬───┬───┬───┐
      0 │ 2 │   │   │   │ 8 │
        ├───┼───┼───┼───┼───┤
      1 │   │   │   │   │   │
    R   ├───┼───┼───┼───┼───┤
    o 2 │   │   │   │   │   │
    w   ├───┼───┼───┼───┼───┤
      3 │   │   │   │   │   │
        ├───┼───┼───┼───┼───┤
      4 │ 0 │   │   │   │ 3 │
        └───┴───┴───┴───┴───┘
</pre>

(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:

```clojure
(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.

<div id="demo-mpd-1" class="threejs"></div>

### 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:

<pre class="lineart">
               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 │
        └───┴───┴───┴───┴───┘
</pre>

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:

```clojure
(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:

* `midpoint` finds the middle point between two points (indexes in an array in
  this case).
* `average2` finds the average of two values (floating point height values
  here).

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.

```clojure
(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:

* `lx`: the x coordinate for the left-hand corners
* `rx`: the x coordinate for the right-hand corners
* `by`: the y coordinate for the bottom corners
* `ty`: the y coordinate for the top corners

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).

<div id="demo-mpd-2" class="threejs"></div>

[SICP]: http://www.amazon.com/dp/0262510871/?tag=stelos-20
[lectures]: https://youtu.be/2Op3QLzMgSY

### 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:

<pre class="lineart">
               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 │
        └───┴───┴───┴───┴───┘
</pre>

We can add this to the displacement function pretty easily:

```clojure
(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:

<div id="demo-mpd-3" class="threejs"></div>

### 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":

<pre class="lineart">
    ╔═══════════╗───┬───┐      ┌───┬───╔═══════════╗
    ║ 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 ║
    ╚═══════════╝───┴───┘      └───┴───╚═══════════╝
</pre>

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

<pre class="lineart">
     Pass 1
    ▩▩▩▩▩▩▩▩▩
    ▩▩▩▩▩▩▩▩▩
    ▩▩▩▩▩▩▩▩▩
    ▩▩▩▩▩▩▩▩▩
    ▩▩▩▩▩▩▩▩▩
    ▩▩▩▩▩▩▩▩▩
    ▩▩▩▩▩▩▩▩▩
    ▩▩▩▩▩▩▩▩▩
    ▩▩▩▩▩▩▩▩▩


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


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

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

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

    ▢▢▢▢▢▢▢▢▢  ▢▢▢▢▢▢▢▢▢  ▢▢▢▢▢▢▢▢▢  ▢▢▢▢▢▢▢▢▢
    ▢▢▢▢▢▢▢▢▢  ▢▢▢▢▢▢▢▢▢  ▢▢▢▢▢▢▢▢▢  ▢▢▢▢▢▢▢▢▢
    ▢▢▢▢▢▢▢▢▢  ▢▢▢▢▢▢▢▢▢  ▢▢▢▢▢▢▢▢▢  ▢▢▢▢▢▢▢▢▢
    ▢▢▢▢▢▢▢▢▢  ▢▢▢▢▢▢▢▢▢  ▢▢▢▢▢▢▢▢▢  ▢▢▢▢▢▢▢▢▢
    ▢▢▢▢▢▢▢▢▢  ▢▢▢▢▢▢▢▢▢  ▢▢▢▢▢▢▢▢▢  ▢▢▢▢▢▢▢▢▢
    ▢▢▢▢▢▢▢▢▢  ▢▢▢▢▢▢▢▢▢  ▢▢▢▢▢▢▢▢▢  ▢▢▢▢▢▢▢▢▢
    ▩▩▩▢▢▢▢▢▢  ▢▢▩▩▩▢▢▢▢  ▢▢▢▢▩▩▩▢▢  ▢▢▢▢▢▢▩▩▩
    ▩▩▩▢▢▢▢▢▢  ▢▢▩▩▩▢▢▢▢  ▢▢▢▢▩▩▩▢▢  ▢▢▢▢▢▢▩▩▩
    ▩▩▩▢▢▢▢▢▢  ▢▢▩▩▩▢▢▢▢  ▢▢▢▢▩▩▩▢▢  ▢▢▢▢▢▢▩▩▩
</pre>

Let's think about what's happening here.

* On pass `0`, we break up the grid into a 1x1 set of chunks and displace them (well, "it").
* On pass `1`, we break up the grid into a 2x2 set of chunks and displace them.
* On pass `2`, we break up the grid into a 4x4 set of chunks and displace them.

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:

```clojure
(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:

```clojure
(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:

```clojure
(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:

```clojure
(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:

```javascript
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:

```clojure
(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.

<div id="demo-mpd-4" class="threejs"></div>

## 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.

<div id="demo-final" class="threejs"></div>
<div id="demo-final-controls">
    <p>
        <label>Exponent
            <input type="number" id="input-exponent" value="5"></input>
        </label>
        <br/>
        <label>Starting Spread
            <input type="number" id="input-starting-spread" value="0.3"></input>
        </label>
        <br/>
        <label>Spread Reduction Constant
            <input type="number" id="input-spread-reduction" value="0.5"></input>
        </label>
    </p>
</div>

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!