# Recursive Midpoint Displacement

Posted on March 7th, 2016.

In the last post we looked at implementing the Midpoint Displacement algorithm. I ended up doing the last step iteratively, which works, but isn't the cleanest way to do it. Before moving on to other algorithms I wanted to clean things up by using a handy library.

The full series of posts so far:

- Multi-Dimensional Arrays
- Iteration
- Updating the Heightmaps
- Slicing Heightmaps
- Updating the Algorithm
- Result

## Multi-Dimensional Arrays

Last week when looking at something unrelated I came across the ndarray
library ("nd" stands for "n-dimensional"). This is a little wrapper around
standard Javascript arrays to add easy multi-dimensional indexing. We're going
to to use `Float64Array`

objects as the underlying storage because they're much
more efficient than the vanilla Javascript arrays and they're fairly well
supported.

It also adds slicing, which is a lot like Common Lisp's displaced arrays. It lets you create a new "array" object with a different "shape" that doesn't have any actual storage of its own, but instead refers back to the original array's data. This is perfect for implementing Midpoint Displacement with recursion.

## Iteration

We're still going to need to iterate over ndarrays at certain points (e.g. when normalizing them) so let's make some helpful macros to do the annoying busywork for us:

`(`*defmacro* do-ndarray [vars array-form & body]
(*let* [array-var (gensym "array")
build
(fn build [vars n]
(*if* (empty? vars)
`(do ~@body)
`(do-times ~(first vars) (aget (.-shape ~array-var) ~n)
~(build (rest vars) (inc n)))))]
`(*let* [~array-var ~array-form]
~(build vars 0))))
(*defmacro* do-ndarray-el [element array-form & body]
(*let* [index (gensym "index")
array (gensym "array")]
`(*let* [~array ~array-form]
(do-times ~index (.-length (.-data ~array))
(*let* [~element (aget (.-data ~array) ~index)]
~@body)))))

Now we can easily iterate over the indices:

```
(do-ndarray [x y] my-ndarray
(console.log "Array[" x "][" y "] is: "
(.get my-ndarray x y)))
```

Or just over the items if we don't need their indices:

```
(do-ndarray-el item my-ndarray
(console.log item))
```

These macros should work for ndarrays of any number of dimensions, and will
compile into ugly but fast Javascript `for`

loops.

## Updating the Heightmaps

To start we'll need to update the heightmap functions to work with ndarrays instead of the normal JS arrays.

Start with a few functions to calculate sizes/incides:

`(`*defn* heightmap-resolution [heightmap]
(aget heightmap.shape 0))
(*defn* heightmap-last-index [heightmap]
(dec (heightmap-resolution heightmap)))
(*defn* heightmap-center-index [heightmap]
(midpoint 0 (heightmap-last-index heightmap)))

Support for reading/writing:

`(`*defn* heightmap-get [heightmap x y]
(.get heightmap x y))
(*defn* heightmap-get-safe [heightmap x y]
(*let* [last (heightmap-last-index heightmap)]
(when (and (<= 0 x last)
(<= 0 y last))
(heightmap-get heightmap x y))))
(*defn* heightmap-set! [heightmap x y val]
(.set heightmap x y val))
(*defn* heightmap-set-if-unset! [heightmap x y val]
(when (== 0 (heightmap-get heightmap x y))
(heightmap-set! heightmap x y val)))

Normalization:

`(`*defn* normalize [heightmap]
(*let* [max (- Infinity)
min Infinity]
(do-ndarray-el el heightmap
(when (< max el) (set! max el))
(when (> min el) (set! min el)))
(*let* [span (- max min)]
(do-ndarray [x y] heightmap
(heightmap-set! heightmap x y
(/ (- (heightmap-get heightmap x y) min)
span))))))

Creation:

`(`*defn* make-heightmap [exponent]
(*let* [resolution (+ (Math.pow 2 exponent) 1)]
(*let* [heightmap (ndarray (new Float64Array (* resolution resolution))
[resolution resolution])]
(set! heightmap.exponent exponent)
(set! heightmap.resolution resolution)
(set! heightmap.last (dec resolution))
heightmap)))

I'm not going to go through all the code here line-by-line because it's mostly just a simple update of the last post.

## Slicing Heightmaps

Part of the Midpoint Displacement is "repeat the process on the four corner squares of this one", and with ndarray we can make getting those corners much simpler:

`(`*defn* top-left-corner [heightmap]
(*let* [center (heightmap-center-index heightmap)]
(-> heightmap
(.lo 0 0)
(.hi (inc center) (inc center)))))
(*defn* top-right-corner [heightmap]
(*let* [center (heightmap-center-index heightmap)]
(-> heightmap
(.lo center 0)
(.hi (inc center) (inc center)))))
(*defn* bottom-left-corner [heightmap]
(*let* [center (heightmap-center-index heightmap)]
(-> heightmap
(.lo 0 center)
(.hi (inc center) (inc center)))))
(*defn* bottom-right-corner [heightmap]
(*let* [center (heightmap-center-index heightmap)]
(-> heightmap
(.lo center center)
(.hi (inc center) (inc center)))))

Each of these will return a "slice" of the underlying ndarray that looks and
acts like a fresh array (e.g. its indices start at `0`

, `0`

), but that uses the
appropriate part of the original array as the data storage.

## Updating the Algorithm

Now we can turn the algorithm into a recursive version. With the slicing functions it's pretty simple. Initializing the corners is still trivial:

`(`*defn* mpd-init-corners [heightmap]
(*let* [last (heightmap-last-index heightmap)]
(heightmap-set! heightmap 0 0 (rand))
(heightmap-set! heightmap 0 last (rand))
(heightmap-set! heightmap last 0 (rand))
(heightmap-set! heightmap last last (rand))))

The meat of the algorithm looks long, but is mostly just calculating all the appropriate numbers with readable names:

`(`*defn* mpd-displace [heightmap spread spread-reduction]
(*let* [last (heightmap-last-index heightmap)
c (midpoint 0 last)
; Get the values of the corners
bottom-left (heightmap-get heightmap 0 0)
bottom-right (heightmap-get heightmap last 0)
top-left (heightmap-get heightmap 0 last)
top-right (heightmap-get heightmap last last)
; Calculate the averages for the points we're going to fill
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)
next-spread (* spread spread-reduction)]
; Set the four edge midpoint values
(heightmap-set-if-unset! heightmap c 0 (jitter bottom spread))
(heightmap-set-if-unset! heightmap c last (jitter top spread))
(heightmap-set-if-unset! heightmap 0 c (jitter left spread))
(heightmap-set-if-unset! heightmap last c (jitter right spread))
; Set the center value
(heightmap-set-if-unset! heightmap c c (jitter center spread))
; Recurse on the four corners if necessary (3x3 is the base case)
(when-not (== 3 (heightmap-resolution heightmap))
(mpd-displace (top-left-corner heightmap) next-spread spread-reduction)
(mpd-displace (top-right-corner heightmap) next-spread spread-reduction)
(mpd-displace (bottom-left-corner heightmap) next-spread spread-reduction)
(mpd-displace (bottom-right-corner heightmap) next-spread spread-reduction))))

The main wrapper function is simple:

`(`*defn* midpoint-displacement [heightmap]
(*let* [initial-spread 0.3
spread-reduction 0.55]
(mpd-init-corners heightmap)
(mpd-displace heightmap initial-spread spread-reduction)
(normalize heightmap)))

## Result

The result looks the same as before, but will generate the heightmaps a lot
faster because it's operating on a `Float64Array`

instead of a vanilla JS array.

The code for these blog posts is a bit of a mess because I've been copy/pasting to show the partially-completed demos. To fix that I've created a little single-page demo with completed versions of the various algorithms you can play with. The code for that should be a lot more readable than the hacky code for these posts.