Cloud images using the diamond-square algorithm

(0 comments)

The diamond-square algorithm is a popular method for generating images which resemble terrain, clouds or plasma. Its description on Wikipedia is as good as any and will not be repeated here.

In the code below, the image size is $N = 2^n + 1$. Since plt.imshow interpolates quite well anyway, $n=6$ gives acceptable results for the image size it produces. For high-resolution images, increase n accordingly.

This code is also available on my github page.

import numpy as np
import matplotlib.pyplot as plt

# The array must be square with edge length 2**n + 1
n = 6
N = 2**n + 1
# f scales the random numbers at each stage of the algorithm
f = 1.0

# Initialise the array with random numbers at its corners
arr = np.zeros((N, N))
arr[0::N-1,0::N-1] = np.random.uniform(-1, 1, (2,2))
side = N-1

nsquares = 1
while side > 1:
    sideo2 = side // 2

    # Diamond step
    for ix in range(nsquares):
        for iy in range(nsquares):
            x0, x1, y0, y1 = ix*side, (ix+1)*side, iy*side, (iy+1)*side
            xc, yc = x0 + sideo2, y0 + sideo2
            # Set this pixel to the mean of its "diamond" neighbours plus
            # a random offset.
            arr[yc,xc] = (arr[y0,x0] + arr[y0,x1] + arr[y1,x0] + arr[y1,x1])/4
            arr[yc,xc] += f * np.random.uniform(-1,1)

    # Square step: NB don't do this step until the pixels from the preceding
    # diamond step have been set.
    for iy in range(2*nsquares+1):
        yc = sideo2 * iy
        for ix in range(nsquares+1):
            xc = side * ix + sideo2 * (1 - iy % 2)
            if not (0 <= xc < N and 0 <= yc < N):
                continue
            tot, ntot = 0., 0
            # Set this pixel to the mean of its "square" neighbours plus
            # a random offset. At the edges, it has only three neighbours
            for (dx, dy) in ((-1,0), (1,0), (0,-1), (0,1)):
                xs, ys = xc + dx*sideo2, yc + dy*sideo2
                if not (0 <= xs < N and 0 <= ys < N):
                    continue
                else:
                    tot += arr[ys, xs]
                    ntot += 1
            arr[yc, xc] += tot / ntot + f * np.random.uniform(-1,1)
    side = sideo2
    nsquares *= 2
    f /= 2

plt.imshow(arr, cmap=plt.cm.Blues)
plt.axis('off')
plt.show()

Cloudy sky generated by the diamond-square algorithm

Current rating: 4.5

Comments

Comments are pre-moderated. Please be patient and your comment will appear soon.

There are currently no comments

New Comment

required

required (not published)

optional

required