6
$\begingroup$

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},
  If[ladera == 0, val = RandomInteger[{0, rugosidad*gridsize}], 
   val = RandomInteger[{-rugosidad*gridsize, rugosidad*gridsize}]
   ];
  val]

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

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]; 
  modificada[[d, 1]] = 
   Ceiling[Mean[{matrix[[a, b]], matrix[[c, e]], matrix[[d, k]]}] + 
     offset[c]];
  modificada[[1, d]] = 
   Ceiling[Mean[{matrix[[a, b]], matrix[[f, g]], matrix[[d, k]]}] + 
     offset[c]];
  modificada[[c, d]] = 
   Ceiling[Mean[{matrix[[h, i]], matrix[[c, e]], matrix[[d, k]]}] + 
     offset[c]];;
  modificada[[d, c]] = 
   Ceiling[Mean[{matrix[[h, i]], matrix[[f, g]], matrix[[d, k]]}] + 
     offset[c]];
  modificada]

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

gruno

grcinco

grtres

$\endgroup$
  • 1
    $\begingroup$ So it is not really mathematica problem but a general algorithm implementation review? $\endgroup$ – Kuba Mar 3 '18 at 6:21
  • $\begingroup$ @Kuba you are right, the doubt arises because of the images obtained, as you have already appreciated $\endgroup$ – bullitohappy Mar 3 '18 at 6:26
  • $\begingroup$ But isn't it off topic then? $\endgroup$ – Kuba Mar 3 '18 at 6:31
  • 1
    $\begingroup$ @Kuba I think the problem is interesting and we had also similar questions. $\endgroup$ – Henrik Schumacher Mar 3 '18 at 6:51
  • 1
    $\begingroup$ @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. $\endgroup$ – Henrik Schumacher Mar 3 '18 at 7:10
8
$\begingroup$

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]

enter image description here

enter image description here

$\endgroup$
  • $\begingroup$ 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 $\endgroup$ – bullitohappy Mar 4 '18 at 1:16
  • $\begingroup$ It was a pleasure! Feel free to contact me for further questions. $\endgroup$ – Henrik Schumacher Mar 4 '18 at 1:21
  • $\begingroup$ @ Henrik Schumacher Hello Henrik, I would like to communicate with you to present my next project that will also be about terrain generation. Please write me to my email that is 22padrinos@gmail.com to explain more details of this project. I am looking forward to your email, an apology for using this community, but I did not find a better way to communicate with you. $\endgroup$ – bullitohappy Mar 10 '18 at 21:24
5
$\begingroup$

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}]

enter image description here

enter image description here

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.

$\endgroup$
  • $\begingroup$ @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 $\endgroup$ – bullitohappy Mar 4 '18 at 1:28
  • $\begingroup$ 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. $\endgroup$ – anderstood Mar 4 '18 at 1:30
  • 1
    $\begingroup$ 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... $\endgroup$ – Henrik Schumacher Mar 4 '18 at 2:24

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.