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
kis 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.
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.
make_terrain() function calls the
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
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 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 - 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 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 )
terrain variable contains a populated n-by-n array, which now
can be rendered in some fashion. Here is an example that uses
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
import matplotlib.colors import matplotlib.pyplot as plt # Create a colormap tmp =  for row in np.loadtxt( "geo-smooth.gpf" ): tmp.append( [ row, 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!
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
ds. I have obtained good results with values of
about 0.55 (smooth) and 0.7 (noisy). Below are two typical examples.