content/blog/2016/06/diamond-square.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 Diamond Square"
:snip "Improving on Midpoint Displacement."
:date "2016-06-27T13:35: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/terrain3.js"></script>

In the last two posts we looked at implementing the Midpoint Displacement
algorithm for procedurally generating terrain.  Today we're going to look at
a similar algorithm called Diamond Square that fixes some problems with 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/)

Midpoint Displacement is simple and fast, but it has some flaws.  If you look at
the end result you'll probably start to notice "seams" in the terrain that
appear on perfectly square chunks.  This happens because when you're calculating
the midpoints of the edges of each square you're only averaging two sources of
information.

[Diamond Square][] is an algorithm that works a lot like Midpoint Displacement,
but it adds an extra step to ensure that (almost) every point uses *four*
sources of data.  This reduces the visual artifacts a lot without much extra
effort.

Let's [draw the owl][owl].

[Diamond Square]: https://en.wikipedia.org/wiki/Diamond-square_algorithm
[owl]: https://i.imgur.com/RadSf.jpg

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

## Overview

The high-level description of Diamond Square goes something like this:

1. Initialize the corners to random values.
2. Set the center of the heightmap to the average of the corners (plus jitter).
3. Set the midpoints of the edges of the heightmap to the average of the four
   points on the "diamond" around them (plus jitter).
4. Repeat steps 2-4 on successively smaller chunks of the heightmap until you
   bottom out at 3x3 chunks.

## Initialization

Initialization is pretty simple:

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

Just set the corners to random values, exactly like we did in Midpoint
Displacement.

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

## Square

In the "square" step of Diamond Square, we take a square region of the heightmap
and set the center of it to the average of the four corners (plus jitter).

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

The Wikipedia article calls this the "diamond" step, but that seems backwards to
me.  We're operating on a *square* region of the heightmap, so this should be
the *square* step.

Since we're no longer working with heightmap "slices" like we were in the
previous post, we'll need a way to specify the square chunk of the heightmap
we're working on.  I chose to use a center point (`x` and `y` coordinates) and
a radius:

```clojure
(defn ds-square [heightmap x y radius spread]
  (let [new-height
        (jitter
          (average4
            (heightmap-get heightmap (- x radius) (- y radius))
            (heightmap-get heightmap (- x radius) (+ y radius))
            (heightmap-get heightmap (+ x radius) (- y radius))
            (heightmap-get heightmap (+ x radius) (+ y radius)))
          spread)]
    (heightmap-set! heightmap x y new-height)))
```

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

## Diamond

In the "diamond" step of Diamond Square, we take a diamond-shaped region and set
the center of it to the average of the four corners:

<pre class="lineart">
               Column
         0 1 2 3 4 5 6 7 8        0 1 2 3 4 5 6 7 8
        ┌─┬─┬─┬─┬─┬─┬─┬─┬─┐      ┌─┬─┬─┬─┬─┬─┬─┬─┬─┐
      0 │1│ │ │ │3│ │ │ │3│    0 │1│ │ │ │3│ │ │ │3│
        ├─┼─┼─┼─┼─┼─┼─┼─┼─┤      ├─┼─┼─┼─╱─╲─┼─┼─┼─┤
      1 │ │ │ │ │ │ │ │ │ │    1 │ │ │ │╱│ │╲│ │ │ │
    R   ├─┼─┼─┼─┼─┼─┼─┼─┼─┤      ├─┼─┼─╱─┼─┼─╲─┼─┼─┤
    o 2 │ │ │3│ │ │ │4│ │ │    2 │ │ │3│ │◉│ │4│ │ │
    w   ├─┼─┼─┼─┼─┼─┼─┼─┼─┤      ├─┼─┼─╲─┼─┼─╱─┼─┼─┤
      3 │ │ │ │ │ │ │ │ │ │    3 │ │ │ │╲│ │╱│ │ │ │
        ├─┼─┼─┼─┼─┼─┼─┼─┼─┤      ├─┼─┼─┼─╲─╱─┼─┼─┼─┤
      4 │5│ │ │ │5│ │ │ │5│    4 │5│ │ │ │5│ │ │ │5│
        ├─┼─┼─┼─┼─┼─┼─┼─┼─┤      ├─┼─┼─┼─┼─┼─┼─┼─┼─┤
      5 │ │ │ │ │ │ │ │ │ │    5 │ │ │ │ │ │ │ │ │ │
        ├─┼─┼─┼─┼─┼─┼─┼─┼─┤      ├─┼─┼─┼─┼─┼─┼─┼─┼─┤
      6 │ │ │7│ │ │ │6│ │ │    6 │ │ │7│ │ │ │6│ │ │
        ├─┼─┼─┼─┼─┼─┼─┼─┼─┤      ├─┼─┼─┼─┼─┼─┼─┼─┼─┤
      7 │ │ │ │ │ │ │ │ │ │    7 │ │ │ │ │ │ │ │ │ │
        ├─┼─┼─┼─┼─┼─┼─┼─┼─┤      ├─┼─┼─┼─┼─┼─┼─┼─┼─┤
      8 │9│ │ │ │7│ │ │ │7│    8 │9│ │ │ │7│ │ │ │7│
        └─┴─┴─┴─┴─┴─┴─┴─┴─┘      └─┴─┴─┴─┴─┴─┴─┴─┴─┘

               Column
         0 1 2 3 4 5 6 7 8        0 1 2 3 4 5 6 7 8
        ┌─┬─┬─┬─┬─┬─┬─┬─┬─┐      ┌─┬─┬─┬─┬─┬─┬─┬─┬─┐
      0 │1│ │ │ │3│ │ │ │3│    0 │1│ │ │ │3│ │ │ │3│
        ├─┼─┼─┼─╱┬╲─┼─┼─┼─┤      ├─┼─┼─┼─╱┬╲─┼─┼─┼─┤
      1 │ │ │ │╱│││╲│ │ │ │    1 │ │ │ │╱│││╲│ │ │ │
    R   ├─┼─┼─╱─┼▼┼─╲─┼─┼─┤      ├─┼─┼─╱─┼▼┼─╲─┼─┼─┤
    o 2 │ │ │3├─▶◉◀─┤4│ │ │    2 │ │ │3├─▶4◀─┤4│ │ │
    w   ├─┼─┼─╲─┼▲┼─╱─┼─┼─┤      ├─┼─┼─╲─┼▲┼─╱─┼─┼─┤
      3 │ │ │ │╲│││╱│ │ │ │    3 │ │ │ │╲│││╱│ │ │ │
        ├─┼─┼─┼─╲┴╱─┼─┼─┼─┤      ├─┼─┼─┼─╲┴╱─┼─┼─┼─┤
      4 │5│ │ │ │5│ │ │ │5│    4 │5│ │ │ │5│ │ │ │5│
        ├─┼─┼─┼─┼─┼─┼─┼─┼─┤      ├─┼─┼─┼─┼─┼─┼─┼─┼─┤
      5 │ │ │ │ │ │ │ │ │ │    5 │ │ │ │ │ │ │ │ │ │
        ├─┼─┼─┼─┼─┼─┼─┼─┼─┤      ├─┼─┼─┼─┼─┼─┼─┼─┼─┤
      6 │ │ │7│ │ │ │6│ │ │    6 │ │ │7│ │ │ │6│ │ │
        ├─┼─┼─┼─┼─┼─┼─┼─┼─┤      ├─┼─┼─┼─┼─┼─┼─┼─┼─┤
      7 │ │ │ │ │ │ │ │ │ │    7 │ │ │ │ │ │ │ │ │ │
        ├─┼─┼─┼─┼─┼─┼─┼─┼─┤      ├─┼─┼─┼─┼─┼─┼─┼─┼─┤
      8 │9│ │ │ │7│ │ │ │7│    8 │9│ │ │ │7│ │ │ │7│
        └─┴─┴─┴─┴─┴─┴─┴─┴─┘      └─┴─┴─┴─┴─┴─┴─┴─┴─┘
</pre>

Once again we'll use a center point and radius to specify the region inside the
heightmap:

```clojure
(defn ds-diamond [heightmap x y radius spread]
  (let [new-height
        (jitter
          (safe-average
            (heightmap-get-safe heightmap (- x radius) y)
            (heightmap-get-safe heightmap (+ x radius) y)
            (heightmap-get-safe heightmap x (- y radius))
            (heightmap-get-safe heightmap x (+ y radius)))
          spread)]
    (heightmap-set! heightmap x y new-height)))
```

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

## Edge Cases

You might have noticed that we used `heightmap-get-safe` and `safe-average` in
the `ds-diamond` code.  This is to handle the diamonds on the edge of the
heightmap.  Those diamonds only have three points to average:

<pre class="lineart">
           Column
         0 1 2 3 4        0 1 2 3 4
        ┌─┬─┬─┬─┬─┐      ┌─┬─┬─┬─┬─┐
      0 │1│ │ │ │8│    0 │1│ │ │ │8│
        ├─┼─┼─┼─┼─┤      ├─┼─┼─┼─╱┬╲
      1 │ │ │ │ │ │    1 │ │ │ │╱│││╲
    R   ├─┼─┼─┼─┼─┤      ├─┼─┼─╱─┼▼┤ ╲
    o 2 │ │ │3│ │ │    2 │ │ │3├─▶◉◀──?
    w   ├─┼─┼─┼─┼─┤      ├─┼─┼─╲─┼▲┤ ╱
      3 │ │ │ │ │ │    3 │ │ │ │╲│││╱
        ├─┼─┼─┼─┼─┤      ├─┼─┼─┼─╲┴╱
      4 │0│ │ │ │3│    4 │0│ │ │ │3│
        └─┴─┴─┴─┴─┘      └─┴─┴─┴─┴─┘
</pre>

So we just ignore the missing point:

```clojure
(defn safe-average [a b c d]
  (let [total 0 count 0]
    (when a (add! total a) (inc! count))
    (when b (add! total b) (inc! count))
    (when c (add! total c) (inc! count))
    (when d (add! total d) (inc! count))
    (/ total count)))

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

Technically this means they have less information to work with than the rest of
the points, and you might see some slight artifacts.  But in practice it's not
too bad, and it's a lot nicer than midpoint displacement where *every* line step
uses only *two* points.

Another way to handle this is to also wrap those edge coordinates around and
pull the point from the other side of the heightmap.  This is a bit more work
but gives you an extra benefit: the heightmap will be tileable (**update**: like
everything where someone doesn't show you running code, it's [not actually that
simple](/blog/2016/08/lisp-jam-postmortem/#tiling-diamond-square)).

## Iteration

Unfortunately we can't use recursion to iterate for Diamond Square like we did
for Midpoint Displacement.  We have to interleave the diamond and square steps,
performing each round of squares on the *entire* map before moving on to a round
of diamonds.  If we don't, bad things will happen.

Let's look at an example to see why.  First we perform the initial square step:

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

That's fine.  Then we recurse into the first of the four diamonds:

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

Everything's still good.  Then we recurse into the square:

<pre class="lineart">
            Column
         0 1 2 3 4        0 1 2 3 4
        ┌─┬─┬─┬─┬─┐      ╔═════╗─┬─┐
      0 │1│ │4│ │8│    0 ║1│ │4║ │8│
        ├─┼─┼─┼─┼─┤      ║─◢─◣─║─┼─┤
      1 │ │ │ │ │ │    1 ║ │ │ ║ │ │
    R   ├─┼─┼─┼─┼─┤      ║─◥─◤─║─┼─┤
    o 2 │ │ │3│ │ │    2 ║?│ │3║ │ │
    w   ├─┼─┼─┼─┼─┤      ╚═════╝─┼─┤
      3 │ │ │ │ │ │    3 │ │ │ │ │ │
        ├─┼─┼─┼─┼─┤      ├─┼─┼─┼─┼─┤
      4 │0│ │ │ │3│    4 │0│ │ │ │3│
        └─┴─┴─┴─┴─┘      └─┴─┴─┴─┴─┘
</pre>

But wait, there's a corner missing!  We needed to complete that left-hand
diamond step to get it, but that's still waiting on the stack for this line of
recursion to bottom out.

### A "Strided" Iteration Construct

So we need to use normal iteration.  We can start by writing the code to do
a single round of diamonds or squares on the heightmap.  First we'll make a nice
looping construct called `do-stride`:

```clojure
(defmacro do-stride
    [varnames start-form end-form stride-form & body]
  (let [stride (gensym "stride")
        start (gensym "start")
        end (gensym "end")
        build (fn build [vars]
                (if (empty? vars)
                  `(do ~@body)
                  (let [varname (first vars)]
                    `(loop [~varname ~start]
                       (when (< ~varname ~end)
                         ~(build (rest vars))
                         (recur (+ ~varname ~stride)))))))]
    ; Fix the numbers once outside the nested loops,
    ; and then build the guts.
    `(let [~start ~start-form
           ~end ~end-form
           ~stride ~stride-form]
       ~(build varnames))))
```

This scary-looking macro builds the nested looping structure we'll need to
perform the iteration for our diamonds and squares.  It abstracts away the
tedium of writing out nested `for` loops so we can focus on the rest of the
algorithm.

`do-stride` takes a list of variables to bind, a number to start at, a number to
end before, and a body.  For each variable it starts at `start`, runs `body`,
increments by `stride`, and loops until the value is at or above `end`.  For
example:

* `(do-stride [x] 0 5 1 (console.log x))` will log `0 1 2 3 4`.
* `(do-stride [x] 0 5 2 (console.log x))` will log `0 2 4`.
* `(do-stride [x y] 0 3 2 (console.log [x y]))` will log `[0, 0] [0, 2] [2, 0]
  [2, 2]`.

### Square Iteration

Now that we've got this, we can iterate our square step over the heightmap for
a specific radius:

```clojure
(defn ds-squares [heightmap radius spread]
  (do-stride [x y] radius (heightmap-resolution heightmap) (* 2 radius)
    (ds-square heightmap x y radius spread)))
```

We use `do-stride` to start at the radius and loop until we fall off the
heightmap, stepping by *twice* the radius each time.

Because our heightmaps are always square we can use the same values for both
dimensions.  `do-stride` can take zero or more variables and will just do the
right thing, so we don't need to worry about it.

The stride is double the radius because as we step between squares we move over
the right half of the first *and* the left half of the next:

<pre class="lineart">
               Column
          0 1 2 3 4 5 6 7 8        0 1 2 3 4 5 6 7 8
         ╔═════════╗─┬─┬─┬─┐      ╔═════════╗─┬─┬─┬─┐
       0 ║1│ │ │ │3║ │ │ │3│    0 ║1│ │ │ │3║ │ │ │3│
         ║─┼─┼─┼─┼─║─┼─┼─┼─┤      ║─┼─┼─┼─┼─║─┼─┼─┼─┤
       1 ║ │ │ │ │ ║ │ │ │ │    1 ║ │ │ │ │ ║ │ │ │ │
    R    ║─┼─┼─┼─┼─║─┼─┼─┼─┤     radius ┼─┼─║─┼─┼─┼─┤
    o  2 ║ │ │◉│ │ ║ │◉│ │ │    2 ║◀━━▶◉◀━━━━━▶◉│ │ │
    w    ║─┼─┼─┼─┼─║─┼─┼─┼─┤      ║─┼─┼─ stride ┼─┼─┤
       3 ║ │ │ │ │ ║ │ │ │ │    3 ║ │ │ │ │ ║ │ │ │ │
         ║─┼─┼─┼─┼─║─┼─┼─┼─┤      ║─┼─┼─┼─┼─║─┼─┼─┼─┤
       4 ║5│ │ │ │5║ │ │ │5│    4 ║5│ │ │ │5║ │ │ │5│
         ╚═════════╝─┼─┼─┼─┤      ╚═════════╝─┼─┼─┼─┤
       5 │ │ │ │ │ │ │ │ │ │    5 │ │ │ │ │ │ │ │ │ │
         ├─┼─┼─┼─┼─┼─┼─┼─┼─┤      ├─┼─┼─┼─┼─┼─┼─┼─┼─┤
       6 │ │ │◉│ │ │ │◉│ │ │    6 │ │ │◉│ │ │ │◉│ │ │
         ├─┼─┼─┼─┼─┼─┼─┼─┼─┤      ├─┼─┼─┼─┼─┼─┼─┼─┼─┤
       7 │ │ │ │ │ │ │ │ │ │    7 │ │ │ │ │ │ │ │ │ │
         ├─┼─┼─┼─┼─┼─┼─┼─┼─┤      ├─┼─┼─┼─┼─┼─┼─┼─┼─┤
       8 │9│ │ │ │7│ │ │ │7│    8 │9│ │ │ │7│ │ │ │7│
         └─┴─┴─┴─┴─┴─┴─┴─┴─┘      └─┴─┴─┴─┴─┴─┴─┴─┴─┘
</pre>

### Diamond Iteration

The diamonds are a little more complicated.

First of all, in the `y` dimension we need to start at `0` instead of `radius`
(these will be the top and bottom edge cases we looked at).  Also, in the `x`
direction all the even iterations need to start at `radius`, while the odd
iterations should start at `0`.

Hopefully a picture will make it a bit more clear:

<pre class="lineart">
               Column
         0 1 2 3 4 5 6 7 8
        ┌─┬─┬─┬─┬─┬─┬─┬─┬─┐
      0 │1│ │◉│ │3│ │◉│ │3│ y iteration: 0, start at x = radius
        ├─┼─┼─┼─┼─┼─┼─┼─┼─┤
      1 │ │ │ │ │ │ │ │ │ │
    R   ├─┼─┼─┼─┼─┼─┼─┼─┼─┤
    o 2 │◉│ │3│ │◉│ │4│ │◉│ y iteration: 1, start at x = 0
    w   ├─┼─┼─┼─┼─┼─┼─┼─┼─┤
      3 │ │ │ │ │ │ │ │ │ │
        ├─┼─┼─┼─┼─┼─┼─┼─┼─┤
      4 │5│ │◉│ │5│ │◉│ │5│ y iteration: 2, start at x = radius
        ├─┼─┼─┼─┼─┼─┼─┼─┼─┤
      5 │ │ │ │ │ │ │ │ │ │
        ├─┼─┼─┼─┼─┼─┼─┼─┼─┤
      6 │◉│ │7│ │◉│ │6│ │◉│ y iteration: 3, start at x = 0
        ├─┼─┼─┼─┼─┼─┼─┼─┼─┤
      7 │ │ │ │ │ │ │ │ │ │
        ├─┼─┼─┼─┼─┼─┼─┼─┼─┤
      8 │9│ │◉│ │7│ │◉│ │7│ y iteration: 4, start at x = radius
        └─┴─┴─┴─┴─┴─┴─┴─┴─┘
</pre>

The code for the diamonds is the hairiest bit of the algorithm.  Make sure you
understand it before moving on.  `shift` will be `radius` for even iterations
and `0` for odd ones:

```clojure
(defn ds-diamonds [heightmap radius spread]
  (let [size (heightmap-resolution heightmap)]
    (do-stride [y] 0 size radius
      (let [shift (if (even? (/ y radius)) radius 0)]
        (do-stride [x] shift size (* 2 radius)
          (ds-diamond heightmap x y radius spread))))))
```

### Top-Level Iteration

Finally we can put it all together.  We initialize the corners and starting
values, then loop over smaller and smaller radii and just call `ds-squares` and
`ds-diamonds` to do the heavy lifting:

```clojure
(defn diamond-square [heightmap]
  (let [initial-spread 0.3
        spread-reduction 0.5
        center (heightmap-center-index heightmap)
        size (aget heightmap.shape 0)]
    (ds-init-corners heightmap)
    (loop [radius center
           spread initial-spread]
      (when (>= radius 1)
        (ds-squares heightmap radius spread)
        (ds-diamonds heightmap radius spread)
        (recur (/ radius 2)
               (* spread spread-reduction))))
    (sanitize heightmap)))
```

## Result

The end result looks like this:

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

It looks a lot like the result from Midpoint Displacement, but without the ugly
seams around the square chunks.

The code for these blog posts is a bit of a mess because I've been copy/pasting
to show the partially-completed demos after each step.  I've created a little
[single-page demo][ymir] with completed versions of the various algorithms you
can play with, and [that code][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/