*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.

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.

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:

- Either a user or the algorithm itself assigns the four corners some values, which can be random.
- 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. - 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**. - 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. - 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)
```

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.

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