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:

  1. Multi-Dimensional Arrays
  2. Iteration
  3. Updating the Heightmaps
  4. Slicing Heightmaps
  5. Updating the Algorithm
  6. 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.


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")
        (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)]

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


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


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

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


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.