home

author: niplav, created: 2020-11-20, modified: 2024-07-09, language: english, status: in progress, importance: 4, confidence: highly likely

The Diamond-Square algorithm is a terrain-generation algorithm for two dimensions (producing a three-dimensional terrain). I generalize the algorithm to any positive number of dimensions, and analyze the resulting algorithm.

Generalizing the Diamond-Square Algorithm to n Dimensions

Libre de la metáfora y del mito
labra un arduo cristal: el infinito
mapa de Aquel que es todas Sus estrellas.

— Jorge Luis Borges, “Spinoza”, 1964

I learned of the diamond-square algorithm by reading through the archives of Chris Wellon's blog, specifically his post on noise fractals and terrain generation. The algorithm is a fairly simple and old one (dating back to the 80s), but not being interested in graphics or game programming, I shelved it as a curiosity.

However, a while later I needed a way to generate high-dimensional landscapes for a simulation, and remembered the algorithm, I felt like I could contribute something here by generalizing the algorithm to produce landscapes in an arbitrary number of dimensions, and that this would be a fun challenge to sharpen my (then fairly weak) Python and numpy skills.

Description

The original (2-dimensional) diamond-square algorithm, in its simplest form, starts with a $2^n+1 \times 2^n+1$ grid of numbers.

It is easiest explained visually:

1. Either a user or the algorithm itself assigns the four corners some values, which can be random.
2. In the diamond step after that, the value in the middle of the grid is determined as the average of the four values in the corners, plus a random value.
3. Next, the middle value every "face" of the grid is determined by the average of the three values in orthogonal directions plus a random value—the square step.
4. The grid is broken down into four sub-grids, and each sub-grid undergoes the diamond step and the square step. The only difference is in the square step: If a point on the grid lies at the face of two sub-grids, it receives the average of all four orthogonal points.
5. The algorithm terminates if each sub-grid is of size $1 \times 1$.

For $n$ dimensions, do that, just higher-dimensional.

We start by initializing an n-dimensional space with zeros, and the corners with random values:

def create_space(dim, size, minval, maxval, factor):
    space=np.zeros([size]*dim)
    corners=(size-1)*get_cornerspos(dim)
    space[*(corners.T)]=np.random.randint(minval, maxval, size=len(corners))

Here, get_cornerspos is just the one-liner return np.array(list(itertools.product([0, 1], repeat=dim))).

We then intialize the variable offsets, and call the recursive diamond-square algorithm:

    offsets=[np.array([0]*dim)]
    return ndim_diamond_square_rec(space, dim, size, offsets, minval, maxval, factor)

Now there are two possible variants of the generalized diamond-square algorithm: the Long Diamond variation and the Long Square variation.

The Long Diamond ⇔ Long Cube Spectrum

Long Diamond

The diamond step of the algorithm starts out with the base case: If the space is only one element big, we return and do nothing (assuming the value has been filled in):

def ndim_diamond_square_rec(space, dim, size, offsets, minval, maxval, factor):
    if size<=1:
        return

We also have to update the size of any axis in the space (not the size of the space itself), we are halving this every recursive call.

    nsize=size//2

Now we come to offsets. Remember above when after the first square step, we moved into a diamond step on the smaller squares? offsets describes where the "left lower corner" of those smaller squares is. We initialized it with zeros, that way we start in a definite corner.

Diamond

Square

Long Square

Code here. I think this is probably the 2nd-most beautiful code I've ever written, just after this absolute smokeshow.

Analysis

Results

One Dimension

Two Dimensions

Three Dimensions