Improve my implementation of the Diamond square algorithm

Some time ago, I began to explore the magnificent world of the terrain generation using the computer. Investigating on the Internet know about the existence of the diamond square algorithm. Try to program that algorithm in Wolfram Mathematica with my little knowledge that I have, getting the following code that I share with you. I'm getting images with some very visible discontinuities but I do not know why, I like to ask for your help to improve the images taken. Any help is welcome.

Here is my code:

BUILT FUNCTIONS

offset[gridsize_] := Module[{val},
];
val]

centro[matrix_List] :=
x = Ceiling[Length[matrix]/2];
media = Mean[{matrix[[1, 1]], matrix[[Length[matrix], 1]],
matrix[[1, Length[matrix]]],
matrix[[Length[matrix], Length[matrix]]]}];

diamante[matrix_List] :=
Module[{a, b, c, d, e, f, g, h, i, k,
modificada = matrix}, {a, b} = {1, 1}; {c, e} = {Length[matrix],
1}; {f, g} = {1, Length[matrix]}; {h, i} =
Dimensions[matrix]; {d, k} = Ceiling[Dimensions[matrix]/2];
Ceiling[Mean[{matrix[[a, b]], matrix[[c, e]], matrix[[d, k]]}] +
offset[c]];
Ceiling[Mean[{matrix[[a, b]], matrix[[f, g]], matrix[[d, k]]}] +
offset[c]];
Ceiling[Mean[{matrix[[h, i]], matrix[[c, e]], matrix[[d, k]]}] +
offset[c]];;
Ceiling[Mean[{matrix[[h, i]], matrix[[f, g]], matrix[[d, k]]}] +
offset[c]];

subgrid[prb_] :=
Module[{upperleft, upperight, half, lowerleft, loweright},
half = Ceiling[Length[prb]/2.];
upperleft = prb[[1 ;; half, ;;]];
lowerleft = prb[[1 ;; half, ;;]];
upperleft = #[[;; half]] & /@ upperleft;
lowerleft = #[[half ;;]] & /@ lowerleft;
upperight = prb[[half ;;, ;;]];
loweright = prb[[half ;;, ;;]];
upperight = #[[;; half]] & /@ upperight;
loweright = #[[half ;;]] & /@ loweright;
{upperleft, upperight, lowerleft, loweright}
]

insertar[chica_List, mayor_List, {xo_, yo_}, {xf_, yf_}] :=
Module[{grande}, grande = mayor;
grande[[xo ;; yo, xf ;; yf]] = chica; grande]

fractdc[red_List] :=
Module[{matr, inter, miniaturas}, matr = centro[red];
matr = diamante[matr];
miniaturas = subgrid[matr];
If[Length[First[miniaturas]] > 2, inter = fractdc[miniaturas[[1]]];
matr = insertar[inter,
matr, {1, Length[inter]}, {1, Length[inter]}]];
If[Length[miniaturas[[2]]] > 2, inter = fractdc[miniaturas[[2]]];
matr = insertar[inter,
matr, {Length[inter], Length[matr]}, {1, Length[inter]}]];
If[Length[miniaturas[[3]]] > 2, inter = fractdc[miniaturas[[3]]];
matr = insertar[inter,
matr, {1, Length[inter]}, {Length[inter], Length[matr]}]];
If[Length[miniaturas[[4]]] > 2, inter = fractdc[miniaturas[[4]]];
matr = insertar[inter,
matr, {Length[inter], Length[matr]}, {Length[inter],
Length[matr]}]];
matr];


INITIAL PARAMETERS

n = 2^7+ 1; mat = ConstantArray[0, {n, n}];

(*Initialize mat*)
altitud = 100; rugosidad = 1; ladera = 1; mat[[1, 1]] = altitud;
mat[[1, Length[mat]]] = altitud;
mat[[Length[mat], 1]] = altitud;
mat[[Length[mat], Length[mat]]] = altitud;

(*Start iterations*)
mat = fractdc[mat];

ListPlot3D[mat, Axes -> False,
Boxed -> False, Mesh -> None,
ColorFunction -> ColorData["SouthwestColors"] , Background -> Black ]


EXAMPLES OF IMAGES OBTAINED WITH THE CODE THAT HAS BEEN SUPPLIED

• So it is not really mathematica problem but a general algorithm implementation review?
– Kuba
Commented Mar 3, 2018 at 6:21
• @Kuba you are right, the doubt arises because of the images obtained, as you have already appreciated Commented Mar 3, 2018 at 6:26
• But isn't it off topic then?
– Kuba
Commented Mar 3, 2018 at 6:31
• @Kuba I think the problem is interesting and we had also similar questions. Commented Mar 3, 2018 at 6:51
• @bullitohappy If I understand your code correctly, you recursively chop the square in finer and finer subblocks and perform fractdc on them. Despite being quite elegant, this is not the diamond-square algorithm. Look here, the last image in the Visualization section: There is also averaging between the subblocks. Commented Mar 3, 2018 at 7:10

I found the algorithm interesting, so I played around with it. This is my implementation of the algorithm. I try to leverage in-place modification of the global matrix and as many vectorized operations as possible. That's why further compilation into a CompiledFunction won't lead to any improvements. (Since CompiledFunctions are not allowed to perform in-place modification, Compile would lead to a loss in performance.)

First, a fused variant of a diamond-square iteration; note the attribute HoldAll that allows us to call diamondsquare by reference. Moreover, we submit an iteration counter and a reservoir of random numbers to this function. The latter allows the user to modify the noise model easily. I cheat a little bit at the boundaries by reflecting the matrix along its boundary edges; that's why I construct a helper list.

ClearAll[diamondsquare]
SetAttributes[diamondsquare, HoldAll];
diamondsquare[A_, iter_, rand_] :=
Module[{d, d2, n, r, shift, ilist, jlist, list, m},
n = Length[A];
m = Ceiling[Log2[n - 1]];
r = 2^(iter - 1);
d = 2^(m - iter);
shift = 1 + d;
d2 = 2 d;
(*diamond step*)
ilist = Range[shift, n, d2];
jlist = Range[shift, n, d2];
A[[ilist, jlist]] = 0.25 Plus[
A[[ilist + d, jlist + d]],
A[[ilist - d, jlist + d]],
A[[ilist - d, jlist - d]],
A[[ilist + d, jlist - d]]
] + rand[[1 ;; r, 1 ;; r]];
(*square step*)
ilist = Range[shift, n, d2];
jlist = Range[1, n, d2];
A[[ilist, jlist]] = 0.25 Plus[
A[[ilist - d, jlist]],
A[[ilist + d, jlist]],
list = jlist - d; list[[1]] = list[[2]]; A[[ilist, list]],
list = jlist + d; list[[-1]] = list[[-2]]; A[[ilist, list]]
] + rand[[r + 1 ;; 2 r, 1 ;; r + 1]];
ilist = Range[1, n, d2];
jlist = Range[shift, n, d2];
A[[ilist, jlist]] = 0.25 Plus[
A[[ilist, jlist - d]],
A[[ilist, jlist + d]],
list = ilist - d; list[[1]] = list[[2]]; A[[list, jlist]],
list = ilist + d; list[[-1]] = list[[-2]]; A[[list, jlist]]
] + rand[[2 r + 1 ;; 3 r + 1, 1 ;; r]];
];


Here is a usage example. First, we set up a noise level σ, an altitute and the starting matrix A:

m = 8;
n = 2^m + 1;
σ = 1.;
altitude = 1.;
SeedRandom[123];
AbsoluteTiming[
A = ConstantArray[0., {n, n}];
A[[{1, n}, {1, n}]] = RandomReal[{-altitude, altitude}, {2, 2}];
][[1]]


0.000119

Do[
rand = RandomReal[
{-σ 2^-iter, σ 2^-iter},
{3 2^(iter - 1) + 1, 2^(iter - 1) + 1}
];
diamondsquare[A, iter, rand];
, {iter, 1, m}]; // AbsoluteTiming // First


0.004468

And here is the result:

ArrayPlot[Rescale[A], DataReversed -> {True, False}, ColorFunction -> "AlpineColors"]
ListPlot3D[A, Axes -> False, Boxed -> False, Mesh -> None,
ColorFunction -> ColorData["AlpineColors"], Background -> Black,
PlotRange -> All]


• I was surprised by the code that you share besides being much shorter in faster. I also mean that the use of HoldAll, this leads me to review your code in more detail to understand it better, thank you for your explanation and for the patience to help me with this algorithm Commented Mar 4, 2018 at 1:16
• It was a pleasure! Feel free to contact me for further questions. Commented Mar 4, 2018 at 1:21

Another implementation (of course less efficient than that of the king of optimisation, Henrik Schumacher):

init := {mat[[1, 1]], mat[[-1, 1]], mat[[-1, -1]], mat[[1, -1]]} =
RandomReal[1, 4]
diamond[p_?EvenQ, e_] := Block[{},
coordsCenters =
Flatten[Table[{p/2 + 1 + i*p, p/2 + 1 + j*p}, {i, 0, nn/p - 1}, {j,
0, nn/p - 1}], 1];
coordsCorners =
Table[coordsCenters[[i]] + p/2*# & /@ {{1, 1}, {-1,
1}, {-1, -1}, {1, -1}}, {i, Length@coordsCenters}];
Table[mat[[Sequence @@ coordsCenters[[i]]]] =
Mean[mat[[Sequence @@ #]] & /@ coordsCorners[[i]]] +
RandomReal[{-e, e}], {i, Length@coordsCenters}];
]
square[p_?EvenQ, e_] := Block[{},
coordsCenters =
Select[Flatten[
Table[{1 + p/2*i, 1 + (1 + (-1)^(i))*p/4 + p*j}, {i, 0,
2 nn/p}, {j, 0, nn/p}], 1], Max[#] <= nn &];
coordsVertices =
Table[coordsCenters[[i]] + p/2*# & /@ {{-1, 0}, {0, 1}, {1,
0}, {0, -1}}, {i, Length@coordsCenters}];
Table[mat[[Sequence @@ coordsCenters[[i]]]] =
Mean[mat[[Sequence @@ #]] & /@
Select[coordsVertices[[i]], 1 <= Min@# && Max@# <= nn &]] +
RandomReal[{-e, e}], {i, Length@coordsCenters}];
]

show[i_] :=
Block[{}, diamond[2^i, e*2.^(i - n)]; square[2^i, 0.5*2^(i - n)];
ListPlot3D[mat, Axes -> False, Boxed -> False, Mesh -> None,
ColorFunction -> ColorData["AlpineColors"], Background -> Black,
PlotRange -> All]]


Use as

n = 8;
e = 1;
nn = 2^n + 1;
mat = ConstantArray[0, {nn, nn}];
init
tab = Table[show[i], {i, n, 1, -1}]


Note: I used's Henrik's code for the visualisation.

Note 2: The only non-obvious thing is that the noise amplitude e should decrease with the step. I've taken this idea from Henrik's answer (again) but I don't know if this is part of the regular algorithm.

• @anderstoood I liked your way of using the Table command in your code, I think I do not take it into account when making my code, so my surprise. Up to where e must decrease in each iteration since the dimensions of each sub-block are reduced. But I must thoroughly review the wikipedia page that Henrik suggested to tell you what should happen to e, many thanks for your code Commented Mar 4, 2018 at 1:28
• Just a thing: the way I generate the points in coordsCorner, coordsCenter and coordsVertices can certainly be simplified a lot and more efficiently. I just did not take the time to improve that. Commented Mar 4, 2018 at 1:30
• Very nice animation! I looked only at MatrixPlots during development so I was quite surprised when I saw the animation for the first time. Btw.: This shrinking of the noise level is quite typical for fractal terrain algorithms (e.g., Perlin noise). Actually, one can obtain quite many different types of terrain by playing around with the noise levels; no need to scale it in this particular way... Commented Mar 4, 2018 at 2:24