content/blog/2016/03/recursive-midpoint-displacement.markdown @ 40accabe3a95 default tip

(
:title "Recursive Midpoint Displacement"
:snip "A cleaner version."
:date "2016-03-07T13:45:00Z"
: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/terrain2.js"></script>

In the [last post][mpd] 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:

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

[mpd]: /blog/2016/02/midpoint-displacement/

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

## 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][disp-arr].  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.

[ndarray]: https://github.com/scijs/ndarray
[slicing]: https://github.com/scijs/ndarray#slicing
[disp-arr]: http://clhs.lisp.se/Body/26_glo_d.htm#displaced_array

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

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

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

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

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

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

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

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

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

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

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

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

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

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][ymir] with completed versions of the various algorithms you
can play with.  [The code for that][ymir-code] should be a lot more readable
than the hacky code for these posts.

[ymir]: http://ymir.stevelosh.com/
[ymir-code]: http://bitbucket.org/sjl/ymir/