Terrain Generation with Diamond Square
Posted on June 27th, 2016.
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 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.
Overview
The high-level description of Diamond Square goes something like this:
- Initialize the corners to random values.
- Set the center of the heightmap to the average of the corners (plus jitter).
- Set the midpoints of the edges of the heightmap to the average of the four points on the "diamond" around them (plus jitter).
- Repeat steps 2-4 on successively smaller chunks of the heightmap until you bottom out at 3x3 chunks.
Initialization
Initialization is pretty simple:
(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.
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).
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│
└─┴─┴─┴─┴─┘ └─┴─┴─┴─┴─┘ └─┴─┴─┴─┴─┘
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:
(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)))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:
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│
└─┴─┴─┴─┴─┴─┴─┴─┴─┘ └─┴─┴─┴─┴─┴─┴─┴─┴─┘
Once again we'll use a center point and radius to specify the region inside the heightmap:
(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)))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:
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│
└─┴─┴─┴─┴─┘ └─┴─┴─┴─┴─┘
So we just ignore the missing point:
(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).
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:
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│
└─┴─┴─┴─┴─┘ └─┴─┴─┴─┴─┘
That's fine. Then we recurse into the first of the four diamonds:
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│
└─┴─┴─┴─┴─┘ └─┴─┴─┴─┴─┘
Everything's still good. Then we recurse into the square:
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│
└─┴─┴─┴─┴─┘ └─┴─┴─┴─┴─┘
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:
(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 log0 1 2 3 4.(do-stride [x] 0 5 2 (console.log x))will log0 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:
(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:
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│
└─┴─┴─┴─┴─┴─┴─┴─┴─┘ └─┴─┴─┴─┴─┴─┴─┴─┴─┘
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:
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
└─┴─┴─┴─┴─┴─┴─┴─┴─┘
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:
(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:
(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:
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 with completed versions of the various algorithms you can play with, and that code should be a lot more readable than the hacky code for these posts.