# The Diamond-Square Algorithm for Terrain Generation

The Diamond-Square Algorithm is the natural first stop for generating artificial landscapes.
The algorithm itself is beautifully simple (more details below, and on its
Wikipedia page).
But a casual implementation ended up not working *at all*, prompting me
to look for an existing implementation to learn from. However, most
implementations I found looked hideously complicated (or just hideous),
not necessarily correct, and/or used out-of-date programming languages
and styles. It therefore seemed like a good idea to create a clean,
simple “reference” implementation of this algorithm, using a contemporary
and widely known programming language and style.

## The Diamond-Square Algorithm

The Diamond-Square Algorithm successively populates an n-by-n array, at increasing level of detail. Initially, only the corners need to be initialized; the remaining cells in the array are then populated by averaging over the four nearest populated neighbors, and adding a small random amount.

- The algorithm requires that the edge length of the array is one greater
than an integer power of two:
`n = 1+2**k`

, where`k`

is a positive integer. - At each step, the amplitude of the random noise is reduced by some factor. The smaller the factor, the smoother the resulting landscape.
- There are two possibilities to deal with pixels laying on the edge of the array: using “fixed” boundary conditions, the average is performed only over the populated cells that lie within the array. Using “periodic” boundary conditions, the averaging routine “wraps around” and pulls in values from the opposite site of the array.

The algorithm feels recursive, but isn’t really. Because the value assigned to each cell depends on the cell’s neighbors, the problem does not partition into a smaller version of itself, and it is more natural to process the entire field at each step, but at increasing level of detail.

## Implementation

Here is a Python implementation of the Diamond-Square algorithm.

The entry point is the `make_terrain()`

function, which manages the
iteration over the decreasing cell sizes, until all pixels have been
filled in. It returns a populated np.array.

The `make_terrain()`

function calls the `single_diamond_square_step()`

routine, which performs first a diamond step, followed by a square
step. The square step is broken into two passes: one fills in the
even-numbered rows, the other the even-numbered columns.

The actual averages over the four neighbors is performed by one of
the two functions `fixed()`

or `periodic()`

. In addition to the actual
averaging, they also contain logic to handle the boundary conditions.

The implementation does not initialize the seed values at the four
corners to user-supplied values; instead, they are all set to zero.
(This could be changed easily.) Notice that for periodic boundary
conditions, all four corners represent the *same* pixel and therefore
should be set to the same value.

```
import random
import numpy as np
def fixed( d, i, j, v, offsets ):
# For fixed bdries, all cells are valid. Define n so as to allow the
# usual lower bound inclusive, upper bound exclusive indexing.
n = d.shape[0]
res, k = 0, 0
for p, q in offsets:
pp, qq = i + p*v, j + q*v
if 0 <= pp < n and 0 <= qq < n:
res += d[pp, qq]
k += 1.0
return res/k
def periodic( d, i, j, v, offsets ):
# For periodic bdries, the last row/col mirrors the first row/col.
# Hence the effective square size is (n-1)x(n-1). Redefine n accordingly!
n = d.shape[0] - 1
res = 0
for p, q in offsets:
res += d[(i + p*v)%n, (j + q*v)%n]
return res/4.0
def single_diamond_square_step( d, w, s, avg ):
# w is the dist from one "new" cell to the next
# v is the dist from a "new" cell to the nbs to average over
n = d.shape[0]
v = w//2
# offsets:
diamond = [ (-1,-1), (-1,1), (1,1), (1,-1) ]
square = [ (-1,0), (0,-1), (1,0), (0,1) ]
# (i,j) are always the coords of the "new" cell
# Diamond Step
for i in range( v, n, w ):
for j in range( v, n, w ):
d[i, j] = avg( d, i, j, v, diamond ) + random.uniform(-s,s)
# Square Step, rows
for i in range( v, n, w ):
for j in range( 0, n, w ):
d[i, j] = avg( d, i, j, v, square ) + random.uniform(-s,s)
# Square Step, cols
for i in range( 0, n, w ):
for j in range( v, n, w ):
d[i, j] = avg( d, i, j, v, square ) + random.uniform(-s,s)
def make_terrain( n, ds, bdry ):
# Returns an n-by-n landscape using the Diamond-Square algorithm, using
# roughness delta ds (0..1). bdry is an averaging fct, including the
# bdry conditions: fixed() or periodic(). n must be 1+2**k, k integer.
d = np.zeros( n*n ).reshape( n, n )
w, s = n-1, 1.0
while w > 1:
single_diamond_square_step( d, w, s, bdry )
w //= 2
s *= ds
return d
```

## Using the Algorithm

Here is an example how this algorithm might be called.

```
n = 1 + 2**9 # Edge size of the resulting image in pixes
ds = 0.6 # Roughness delta, 0 < ds < 1 : smaller ds => smoother results
bdry = periodic # One of the averaging routines: fixed or periodic
terrain = make_terrain( n, ds, bdry )
```

The `terrain`

variable contains a populated n-by-n array, which now
can be rendered in some fashion. Here is an example that uses
matplotlib’s `imshow()`

function to create a heatmap version of
the terrain. The resulting image can be saved to file (with
`savefig()`

) or rendered to the screen (using `show()`

).

```
import matplotlib.colors
import matplotlib.pyplot as plt
# Create a colormap
tmp = []
for row in np.loadtxt( "geo-smooth.gpf" ):
tmp.append( [ row[0], row[1:4] ] )
cm = matplotlib.colors.LinearSegmentedColormap.from_list( "geo-smooth", tmp )
# Create an image using the colormap
plt.figure( figsize=( n/100, n/100 ), dpi=100 ) # create n-by-n pixel fig
plt.tick_params( left=False, bottom=False, labelleft=False, labelbottom=False )
plt.imshow( terrain, cmap=cm )
plt.savefig( "terrain.png" ) # Save to file
plt.show() # Show interactive
```

The realistic appearance of the final image depends pretty critically on the colormap chosen to create the image. In fact, using a generic false-color palette, the resulting image does not evoke the impression of a “landscape” at all!

The file `geo-smooth.gpf`

contains a
palette suitable for topographic maps. All color transitions are smooth,
with the exception of the one from “sea” to “land”, which is sharp, hence
giving rise to sharp “coast lines” in the final image.

The roughness of the generated landscape is controlled by the roughness
parameter `ds`

. I have obtained good results with values of `ds`

between
about 0.55 (smooth) and 0.7 (noisy). Below are two typical examples.