content/blog/2016/03/recursive-midpoint-displacement.html @ bb75d39cfff3

    {% extends "_post.html" %}

    {% load mathjax %}

    {% hyde
        title: "Recursive Midpoint Displacement"
        snip: "A cleaner version."
        created: 2016-03-07 13:45:00
    %}

    {% block extra_js %}
        <script data-cfasync="false" src="/media/js/three.min.js"></script>
        <script data-cfasync="false" src="/media/js/TrackballControls.js"></script>

        <script data-cfasync="false" src="/media/js/terrain2.js"></script>
    {% endblock extra_js %}

    {% block article %}

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/

[TOC]

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

    {% endblock article %}