Procedural Terrain Generation
Diamond Square Algorithm
Introduction
Realistic terrain using four random values and basic math? You bet! That’s what I used to create the terrain above. Procedural generation is a method of creating data algorithmically instead of manually. It is commonly used in video games to create textures or landscapes. The 2016 title, No Man’s Sky, boasted over 18 quintillion planets, each with unique landscapes and creatures, all created by using procedural generation algorithms. The Diamond Square Algorithm is a straight forward algorithm for procedurally generating terrains.
Diamond Square Algorithm
The Diamond Square Algorithm starts with a heightmap. A heightmap is a matrix filled with values. For those not familiar, think of a matrix like an Excel Spreadsheet, it has rows and columns that hold data and allow that data to be manipulated. The matrix’s row and column indices act as the \(X\) and \(Y\) position coordinates. The values placed inside the matrix are heights or the \(Z\) position coordinate. A heightmap maps a 2D coordinate \((X, Y)\) with a height \((Z)\). In other words, a height map can describe the topography of an area.
The diamond square algorithm can be broken out into five main steps:
- Set the Dimensions
- Seed the Corners
- Perform the Diamond Step
- Perform the Square Step
- Rinse and Repeat*
1. Set the Dimensions
To start matrix needs to be square and have the dimensions \((2^n + 1)\) x \((2^n + 1)\).
Where \(n \geqslant 1 \) and \(n\) is an integer. Dimensions like 3x3, 5x5, 9x9 ect. We will look at the 5x5 case for simplicity.
Why \(2^n + 1\)? And what is \(n\) anyway? It’s not the dimensions of the grid. \(n\) is the number of iterations through the diamond square algorithm needed to fill every element in the matrix. The equation, \(2^n + 1\), maintains that each time through the algorithm each ‘square’ has a distinct center element. For example, a 7x7 matrix does not meet the criteria for \(2^n + 1\), where \(n\) is an integer. During the second call to the algorithm, the next ‘square’ does not have a distinct center element.
5x5 Grid : \(\bf{n = 2}\)
\(\bf{5 = 2^2 + 1}\)
7x7 Grid
\(\bf{7 \ne 2^n + 1}\)
9x9 Grid : \(\bf{n = 3}\)
\(\bf{9 = 2^3 + 1}\)
2. Seed the Corners
The matrix starts off empty, with dimensions of \((2^n + 1)\) x \((2^n + 1)\), then any four random values are placed in the corners of the matrix.
Shown as \(A, B, C, D\).
3. Perform the Diamond Step
The diamond step takes the four corner elements of a square, then finds the center element of that square. The value of the center element is set equal to the average of the four corner element values plus a random value \((valRan)\). The center element is shown as \(E\).
The diamond step gets its name because if you were to draw lines showing the step, the lines would be diagonal. Over a large matrix this would generate a diamond pattern.
4. Perform the Square Step
The square step takes the tip elements of the diamond, created in the diamond step, and then finds the center element of that diamond. The value of the center element is set equal to the average of all the tip elements plus a random value. The center elements are shown as \(F, G, H, I\).
The square step gets its name because if you were to draw lines showing the step, the lines would be perpendicular to one another, generating a square pattern.
5. Rinse and Repeat*
The diamond and square steps are repeated over and over until all elements within the matrix are filled with values. The below images show the full progression of the diamond square algorithm over a 5x5 matrix.
Set Dimensions
Seed the Corners
Diamond Step \(\bf{n = 1}\)
Square Step \(\bf{n = 1}\)
Diamond Step \(\bf{n = 2}\)
Square Step \(\bf{n = 2}\)
* There is a catch to this process. If the element is an edge to the matrix then it is only the average of the surrounding three elements. If the element is not an edge then it is the average of the surrounding four elements. Another way to deal with the edge issue is by wrapping the left and the right edges of the matrix together, which I omitted from my coding of this algorithm.
Edge Case
Non-Edge Case
Roughness Coefficient
How should \((valRan)\) be generated and how does it influence the overall look of the terrain?
\((valRan)\) is generated based on something called the Roughness Coefficient \((Rco)\).
\( 0 \leqslant Rco \leqslant 1\). The roughness coefficient determines how varied or rough the terrain is, and is decreased in magnitude every iteration through the algorithm.
At iteration \(n\) : \( -Rcoe^n \leqslant valRan \leqslant Rcoe^n\).
Roughness = 0.10
Roughness = 0.25
Roughness = 0.50
Roughness = 0.90
Resources and Going Further
The diamond square algorithm is only the tip of the iceberg when it comes to procedural content generation. The following links were extremely helpful in the making of this post :
3D Graphic Java: Render Fractal Landscapes
An older article, but the most in-depth as it covers lighting, shadows and view point calculations.
Algorithms for Procedural Content Generation
A listing of many high level concepts related to writing code for procedural content generation.
Terrian Generation with Diamond Square
A full series of posts covering: midpoint displacement, recursive midpoint displacement and the diamond square algorithm.
I took an iterative approach, as opposed to a recursive one, when I coded my version of the Diamond Square Algorithm. I used MatLab because it was the most accessible 3D graphing tool I had on hand. There are three functions: the main script, the plotting sub-function and the random value generation sub-function. The best way to learn something is to try coding it on your own first, but feel free to use and play around with this code.
PTG: Main Script
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
%% Procedural Terrain Generation
% Diamond Square Algorithm
% Sarah Wilson
% https://sarahawilson.github.io/
%% Defined Inputs
n = 3; % Number of Iterations
GRDIM = ((2^n)+1); % Grid Dimensions
RCOE = .2; % Roughness Coefficient
plotNum = 0; % Plot Number
%% Initialize Grid
grid = zeros(GRDIM);
rValUpper = RCOE^n; % Upper Limit for Roughness
rValLower = (-1*rValUpper); % Lower Limit for Roughness
grid(1, 1) = 1; % Pre-seed the corners with a value
grid(1, GRDIM) = 1.4;
grid(GRDIM, 1) = 1;
grid(GRDIM, GRDIM) = 1.2;
stpDim = GRDIM; % Step Dimension to be used in loops
mp = (ceil(stpDim/2)); % Midpoint of each square
%% Main
while (mp > 1)
% Diamond Step
for row = 1:(stpDim - 1):GRDIM
if (row == GRDIM)
break; % If edge of grid by row break out
end
for col = 1:(stpDim - 1):GRDIM
if (col==GRDIM)
break; % If edge of grid by col break out
end
topLeft = grid(row,col);
topRight = grid(row, (col + (stpDim - 1)));
botLeft = grid((row + (stpDim - 1)), col);
botRight = grid((row + (stpDim - 1)),(col + (stpDim - 1)));
valRan = randValGen(rValLower, rValUpper); % Random Value
valD = ((topLeft + topRight + botLeft + botRight)/4) + valRan;
grid((row + (mp - 1)), (col + (mp - 1))) = valD;
end
end
% Square Step
for row = 1:(stpDim - 1):GRDIM
if (row == GRDIM)
break; % If edge of grid by row break out
end
for col = 1:(stpDim - 1):GRDIM
if (col == GRDIM)
break; % If edge of grid by col break out
end
sTopL = grid(row, col); % Get the values for the conrners from the diamond step
sTopR = grid(row, (col + (stpDim - 1)));
sBotL = grid((row + (stpDim - 1)), col);
sBotR = grid((row + (stpDim - 1)),(col + (stpDim - 1)));
sCen = grid((row + (mp - 1)), (col + (mp - 1)));
rValR = randValGen(rValLower, rValUpper);
rValL = randValGen(rValLower, rValUpper);
rValT = randValGen(rValLower, rValUpper);
rValB = randValGen(rValLower, rValUpper);
curRowCen = row + (mp - 1);
curColCen = col + (mp - 1);
if (row == 1)
Top = ((sTopL + sTopR + sCen)/3) + rValR;
grid(row, (col + (mp - 1))) = Top;
end
if (col == 1)
Left = ((sTopL + sBotL + sCen)/3) + rValL;
grid((row + (mp - 1)), col) = Left;
end
if ((row + (stpDim - 1)) == GRDIM)
Bot = ((sBotL + sBotR + sCen)/3) + rValB;
grid(GRDIM, (col + (mp - 1))) = Bot;
else
bnsCen = grid((curRowCen + (stpDim - 1)), curColCen);
% bns Bottom Next Square
Bot = ((sBotL + sBotR + sCen + bnsCen)/4) + rValB;
grid((row + (stpDim - 1)), curColCen) = Bot;
end
if ((col + (stpDim - 1)) == GRDIM)
Right = ((sBotR + sTopR + sCen)/3) + rValR;
grid((row + (mp-1)), GRDIM) = Right;
else
rnsCen = grid(curRowCen, (curColCen + (stpDim - 1)));
Right = ((sTopR + sBotL + sCen + rnsCen)/4) + rValR;
grid(curRowCen, (col + (stpDim - 1))) = Right;
end
end
end
% Get Ready for the Next Iteration
stpDim = mp;
mp = ceil(stpDim/2);
n = n - 1;
rValUpper = RCOE^n;
rValLower = (-1 * rValUpper);
end
%% Final Plot
ptg_plotter( grid, plotNum );
PTG: Random Value Generator Sub-Function
1
2
3
4
5
6
7
%% Random Value Generator
% Sarah Wilson
% https://sarahawilson.github.io/
function [ outRandVal ] = randVal( rghUpper, rghLower )
outRandVal = rghUpper + (rghUpper - rghLower).* rand(1);
end
PTG: Plotting Sub-Function
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
%% PTG Plotter
% Sarah Wilson
% https://sarahawilson.github.io/
function [] = ptg_plotter( grid, plotNum )
figure()
surf(grid);
view([-47 74])
caxis('manual');
demcmap(([0 1.8]));
zlim([0 2.5]);
step = num2str(plotNum);
if (plotNum<10)
fileName = ['grid_0' step]
else
fileName = ['grid_' step]
end
saveas(gcf, fileName,'svg')
close all
end