﻿ The Diamond Square Algorithm
HOME Mathematics Computer Science
The Diamond Square Algorithm
Prerequisites
Show All
Hide All
Heightmaps
Show/Hide
﻿
Introduction
• Heightmaps are a way of representing a 3D surface as a 2D image.
• Each pixel in the 2D image corresponds to a specific point on the 3D surface.
• The brightness of the pixel determines the height of the point:
• A dark pixel is a low point
• A bright pixel is a high point
Image rendered in SketchUp
Code
from PIL import Image

im = Image.new("L", (256, 256))

for x in range(256):
for y in range(256):
pix[x,y] = 256 - abs(x - 128) - abs(y - 128)

im.save("test.png", "PNG")
im.show()
Min-Max Normalization
Show/Hide
﻿
Introduction
• Min-max normalization is an operation which rescales a set of data.
• This can be useful when:
• Comparing data from two different scales
• Converting data to a new scale
• In most situations, data is normalized to a fit a target range of [0, 1]
• The smallest value in the original set would be mapped to 0
• The largest value in the original set would be mapped to 1
• Every other value would be mapped to a value somewhere between these two bounds
• It is also called:
• Feature Scaling
• Min-Max Scaling
• Rescaling
• Normalization
Normalizing to [0, 1]
• A set of numbers will have:
• A smallest value:
• This is also called the lower bound or least element
• It is denoted: min(x)
• A largest value:
• This is also called the upper bound or greatest element
• It is denoted: max(x)
• A range:
• The difference between the smallest and largest values
• It is denoted: max(x) - min(x)
• Normalization is the process of changing the lower and upper bounds to be 0 and 1 respectively
Algorithm
1. First we modify the data to have a lower bound of 0. To do this we subtract min(x) from each value:
2. Then we modify the data to have an upper bound of 1. We do this by dividing each value by the original range:
3. Finally, if we combine these two steps we get:
Example
Normalize the following data:
1. First we calculate the lower bound, upper bound, and range:
• min(x) = 7
• max(x) = 21
• max(x) - min(x) = 14
2. Next we subtract the lower bound from each value:
3. Finally, we divide by the range:
Code
import numpy as np

def normalize(x):
min = np.min(x)
max = np.max(x)
range = max - min

return [(a - min) / range for a in x]

x = [7, 21, 13, 15]
normalizedX = normalize(x)

print(normalizedX) # prints: [0.0, 1.0, 0.42857142857142855, 0.5714285714285714]
Normalizing from [0, 1]
• If our numbers are in the range [0, 1] then we can scale them to have a different lower and upper bound
• To achieve this, we simply do the reverse of normalization:
1. Find the new range by subtracting the lower bound from the upper bound
2. Multiply each value by the new range
3. and add the new lower bound to each value:
Example
Normalize the following data to have a lower bound of 3 and an upper bound of 24:
1. first we calculate the range:
2. Then we multiply each value by the range:
3. Finally, we add the lower bound to each value:
Code
def normalize(normalizedX, newLowerBound, newUpperBound):
range = newUpperBound - newLowerBound

return [a * range + newLowerBound for a in normalizedX]

normalizedX = [0.0, 1.0, 3/7, 4/7]
x = normalize(normalizedX, 3, 24)

print(x) # prints: [3.0, 24.0, 12.0, 15.0]
Normalizing from one range to another
• Sometimes we need to normalize data in which neither the source range nor the target range is [0, 1]
• In these situations, we first normalize the data to range of [0, 1], and then normalize it again to the true target range.
• These two steps can be combined:
Code
import numpy as np

def normalize(x, newLowerBound, newUpperBound):
min = np.min(x)
max = np.max(x)
range = max - min
newRange = newUpperBound - newLowerBound

return [((a - min) / range) * newRange + newLowerBound for a in x]

x = [7, 21, 13, 15]
y = normalize(x, 3, 24)

print(y) # prints: [3.0, 24.0, 12.0, 15.0]
2n + 1
Show/Hide
﻿
Introduction
• We can generate a sequence of integers using the following equation:
• These numbers have an important property when it comes to calculating their midpoints:
• To calculate the midpoint between two numbers, you can use this equation:
• For example: the midpoint between 1 and 17 is:
• The midpoint between 1 and any number in the sequence 2n + 1 is equal to the previous number in the sequence
• For example:
• The midpoint between 1 and 33 is 17
• The midpoint between 1 and 17 is 9
• The midpoint between 1 and 9 is 5
• Proof:
• the midpoint between 1 and 2n + 1 is:
• A number which is in the sequence 2n + 1 can be recursively split into two segments until each segment has a size of 1:
Code
from collections import deque

# The following code creates an array of numbers starting at 0 and smoothly increasing to 255

n = 3

# creates an array of 0's of size 2^n + 1
x = [0]*(2**n + 1)

# set the value for the last element
x[len(x) - 1] = 255

# we create a queue of the line segments we need to subdivide
q = deque()

# and add the first segment to it
q.append((0, len(x) - 1))

# now we go through and subdivide each segment
while len(q) != 0:
left, right = q.popleft()
midpoint = (left + right) // 2

# set the midpoint to the average of the two ends
x[midpoint] = (x[left] + x[right]) // 2

# if the width of the segment is greater than 2 then it can be subdivided
if right - left > 2:
q.append((left, midpoint))
q.append((midpoint, right))

print(x) # prints: [0, 31, 63, 95, 127, 159, 191, 223, 255]
Noise Functions
Show/Hide
﻿
Random Values
We can generate random values using a pseudorandom number generator:
Example
 X 1 2 3 4 ... rand(X) 28 21 96 96 ...
Visualization
Code
import random

randX = [random.randint(0, 100) for x in range(0, 50)]
print(randX)
Noise
Many areas of procedural generation require random values, but the randomness needs to look more organic:
Visualization
This is the job of a noise function. It produces random data, but the data has an underlying organic nature.
Higher Dimensions
• The above examples were in one dimension
• X is an array of numbers, for each number in X there is a corresponding value.
• However it is often useful to generate noise for 2 or more dimensions:

Heightmap

3D Terrain

Image rendered in SketchUp
• In the above example a random value is assigned to each (x, y) coordinate pair:
• 2D noise functions are often used to generate terrains in which the height values are the result of the noise function applied to each point.
• (x, y, z) = (x, y, noise(x, y))
• Noise functions can be extended into 3 or more dimensions:
• For example, a 3D noise function could be used for modelling gas. The random value at each 3D point could describe its density.
• 3D noise functions can also be used to animate 2D noise functions:
• A 3D array of data is produced and then cut into a series of 2D slices.
• You can then convert each slice into an individual frame of an animation.
• This can then be used for when you need 2D noise that changes over time, such as when animating the formation of clouds.
Midpoint Displacement in one dimension
Show/Hide
﻿
Introduction
• Midpoint Displacement is a simple algorithm for generating realistic looking height ranges.
• It is a recursive algorithm that subdivides a line into segments, adding a decreasing amount of randomness at each round.
Algorithm
1. We start with a set of heights initialized to 0. The length of the set must conform to 2n + 1.
2. Then we set the first and last heights to random value. In this case they are 5 and 2:
3. Next we set the midpoint to the average of the two endpoints:
4. And move (or displace) the midpoint up or down by a random amount. In this case we move it down by 1:
5. And repeat the midpoint displacement recursively for the newly created line segments:
6. Each time we subdivide the line segments, the amount we displace the midpoint should also decrease:
Code
from PIL import Image
import PIL.ImageDraw as ImageDraw
import random
import numpy as np
from math import floor
from collections import deque

def main():
imageWidth = 513
imageHeight = 256
roughness = 200

# we initialize all the heights to 0
heights = [0]*imageWidth

# and then initialize the first and last heights to random values
heights[0] = random.randint(0, 256)
heights[imageWidth - 1] = random.randint(0, 256)

# we create a queue of the line segments we need to subdivide
q = deque()

# and add the first segment to it
q.append((0, imageWidth - 1, roughness))

# now we go through and subdivide each segment
while len(q) != 0:
left, right, randomness = q.popleft()
center = (left + right + 1) // 2

# set the midpoint to the average of the two ends
heights[center] = (heights[left] + heights[right]) // 2

# and then add some randomness
heights[center] = heights[center] + random.randint(-randomness, randomness)

# if the width of the segment is greater than 2 then it can be subdivided
if right - left > 2:
# when we add the new line segments to the queue, we reduce the amount of randomness
q.append((left, center, floor(randomness // 2)))
q.append((center, right, floor(randomness // 2)))

# finally we render these heights
im = renderHeights(heights, imageWidth, imageHeight)
im.save("test.png", "PNG")
im.show()

# draws the given heights on the image
# the heights should be normalized to be between 0 and image.height
def renderHeights(heights, imageWidth, imageHeight):

# we normalize the heights so that its minimum and maximum values
# fit nicely on the image we are about to draw
heights = normalize(heights, imagePadding, imageHeight - imagePadding)

# we need to convert the heights to a series of points to be plotted
points = [(i, heights[i]) for i in range(0, imageWidth)]

# and add the corners so that tha polygon is nice and closes properly
points.insert(0, (0, 0))
points.append((imageWidth - 1, 0))

im = Image.new("L", (imageWidth, imageHeight))
draw = ImageDraw.Draw(im)
draw.polygon(points, fill=200)

return im

# normalizes an array of numbers such that the smallest number will be
# the newLowerBound and the largest number will be the newUpperBound
# every other number will be scaled to be between those points
def normalize(data, newLowerBound, newUpperBound):
min = np.min(data)
max = np.max(data)
range = max - min
newRange = newUpperBound - newLowerBound

return [(a - min) * newRange / range + newLowerBound for a in data]

main()
Midpoint Displacement in two dimensions
Show/Hide
﻿
Introduction
• The Midpoint Displacement algorithm can be extended into two dimensions
• This allows us to generate 2D heightmaps that can be used for terrain generation
Image rendered in SketchUp
Algorithm
1. We start by setting the corners to random values:
2. Then we displace the midpoint of each edge:
3. And set the center to be the average of the 4 midpoints + a random displacement:
4. Finally, we apply the above two steps recursively:
Code
import random
from math import floor
from collections import deque

def main():
size = 129
heightmap = [[0]*size for i in range(size)]

heightmap[0][0] = random.randint(0, 256)
heightmap[size - 1][0] = random.randint(0, 256)
heightmap[0][size - 1] = random.randint(0, 256)
heightmap[size - 1][size - 1] = random.randint(0, 256)

q = deque()
q.append((0, 0, size - 1, size - 1, 200))

while len(q) != 0:
top, left, bottom, right, randomness = q.popleft()

centerX = (left + right) // 2
centerY = (top + bottom) // 2

heightmap[centerX][top] = (heightmap[left][top] + heightmap[right][top]) // 2 + random.randint(-randomness, randomness)
heightmap[centerX][bottom] = (heightmap[left][bottom] + heightmap[right][bottom]) // 2 + random.randint(-randomness, randomness)
heightmap[left][centerY] = (heightmap[left][top] + heightmap[left][bottom]) // 2 + random.randint(-randomness, randomness)
heightmap[right][centerY] = (heightmap[right][top] + heightmap[right][bottom]) // 2 + random.randint(-randomness, randomness)

heightmap[centerX][centerY] = (heightmap[left][top] +
heightmap[right][top] +
heightmap[left][bottom] +
heightmap[right][bottom]) // 4 + \
random.randint(-randomness, randomness)

if right - left > 2:
q.append((top, left, centerY, centerX, randomness // 2))
q.append((top, centerX, centerY, right, randomness // 2))
q.append((centerY, left, bottom, centerX, randomness // 2))
q.append((centerY, centerX, bottom, right, randomness // 2))

main()

The Diamond Square Algorithm
﻿
Introduction
• The diamond square algorithm allows us to generate 2D heightmaps that can be used for terrain generation.
• It is an improvement on the 2D midpoint displacement algorithm, which can produce noticeable creases.

Heightmap

3D Terrain

Image rendered in SketchUp
Algorithm
1. We start by setting all the corners to the same random value:
2. Then we set the center to be the average of the 4 corners + a random displacement. This is called the diamond step.
3. And set the midpoint of each edge to be the average of its corner points + a random displacement. This is called the square step.
• Some points are outside of the bounds of the heightmap. I.e. there is no point to the left of the left edge.
• In this case, we wrap around to the other side of the heightmap: the point left of the left edge is the same as the point to the left of the right edge.
4. Finally, we apply the above steps recursively, reducing the random displacement with each pass:
Code
import random
import numpy as np
from math import floor
from PIL import Image
from collections import deque

heightmapWidth = 65

# initialize the heightmap to 0's
heightmap = [[0]*heightmapWidth for i in range(heightmapWidth)]

# set the corner points to the same random value
rand = random.randint(0, 256)
heightmap[0][0] = rand
heightmap[heightmapWidth - 1][0] = rand
heightmap[0][heightmapWidth - 1] = rand
heightmap[heightmapWidth - 1][heightmapWidth - 1] = rand

# set the randomness bounds, higher values mean rougher landscapes
randomness = 128
tileWidth = heightmapWidth - 1

# we make a pass over the heightmap
# each time we decrease the side length by 2
while tileWidth > 1:
halfSide = tileWidth // 2

# set the diamond values (the centers of each tile)
for x in range(0, heightmapWidth - 1, tileWidth):
for y in range(0, heightmapWidth - 1, tileWidth):
cornerSum = heightmap[x][y] + \
heightmap[x + tileWidth][y] + \
heightmap[x][y + tileWidth] + \
heightmap[x + tileWidth][y + tileWidth]

avg = cornerSum / 4
avg += random.randint(-randomness, randomness)

heightmap[x + halfSide][y + halfSide] = avg

# set the square values (the midpoints of the sides)
for x in range(0, heightmapWidth - 1, halfSide):
for y in range((x + halfSide) % tileWidth, heightmapWidth - 1, tileWidth):
avg = heightmap[(x - halfSide + heightmapWidth - 1) % (heightmapWidth - 1)][y] + \
heightmap[(x + halfSide) % (heightmapWidth - 1)][y] + \
heightmap[x][(y + halfSide) % (heightmapWidth - 1)] + \
heightmap[x][(y - halfSide + heightmapWidth - 1) % (heightmapWidth - 1)]

avg /= 4.0
avg += random.randint(-randomness, randomness)

heightmap[x][y] = avg

# because the values wrap round, the left and right edges are equal, same with top and bottom
if x == 0:
heightmap[heightmapWidth - 1][y] = avg
if y == 0:
heightmap[x][heightmapWidth - 1] = avg

# reduce the randomness in each pass, making sure it never gets to 0
randomness = max(randomness // 2, 1)
tileWidth //= 2