Diamond-square algorithm for making maps in Mathematica

My previous terrain generator was not very good at making accurate maps. I have recently wanted to revisit my terrain generator and have found the time to write a diamond-square algorithm to make my terrain.

This post provided me with the understanding to create my code and has a good explanation of the concept behind the algorithm.



arraymapmaker[order_, initialjitter1_, initialjitter2_,
jitterfunction1_, jitterfunction2_, read_, filename_] :=
Module[{size, base, mid, map, new, dimlist, squarecorners,
diamondcorners, dim, file, jitter1, jitter2, initialjitter,
cornerx, cornery},

size = 2^order + 1;
initialjitter := RandomReal[{initialjitter1, initialjitter2}];
map = SparseArray[{{1, 1} -> initialjitter, {size, size – 1} ->
0, {1, size – 1} -> 0, {size, 1} -> initialjitter}];
mid[n_] := (n + 1)/2;
dimlist = Reverse[Table[2^i + 1, {i, 1, order}]];
file = filename;

Do[
dim = dimlist[[i]];
squarecorners =
Flatten[Table[{y (dim – 1) + 1, x (dim – 1) + 1}, {y,
0, (size – 1)/(dim – 1) – 1, 1}, {x,
0, (size – 1)/(dim – 1) – 1, 1}], 1];
diamondcorners =
DeleteDuplicates[
Join[Flatten[
Table[{y (dim – 1) + 1, Mod[x (dim – 1) + 1, size – 1, 1]}, {y,
0, (size – 1)/(dim – 1), 1}, {x, 0, (size – 1)/(dim – 1),
1}], 1],
Flatten[Table[{y (dim – 1) + 1 – 1 + mid[dim],
x (dim – 1) + 1 – 1 + mid[dim]}, {y,
0, (size – 1)/(dim – 1) – 1, 1}, {x,
0, (size – 1)/(dim – 1) – 1, 1}], 1]]];
jitter1 = jitterfunction1;
jitter2 = jitterfunction2;

Do[
cornery = squarecorners[[i, 1]];
cornerx = squarecorners[[i, 2]];
new = (map[[cornery, cornerx]] +
map[[cornery, Mod[cornerx – 1 + dim, size – 1, 1]]] +
map[[cornery – 1 + dim, cornerx]] +
map[[cornery – 1 + dim,
Mod[cornerx – 1 + dim, size – 1, 1]]])/4. +
RandomReal[{jitter1, jitter2}];
(map[[cornery – 1 + mid[dim], cornerx – 1 + mid[dim]]] = #) &@
new;

, {i, 1, Length[squarecorners]}];
If[read == 1, Export[file, map, "HarwellBoeing"]];
Print["Completed Squares for dimension " <> ToString[dimlist[[i]]]]

Do[
cornery = diamondcorners[[i, 1]];
cornerx = diamondcorners[[i, 2]];

If[(cornery – 1 + mid[dim]) > size,(*southpole*)
new =
(map[[cornery, cornerx]]
+ map[[cornery, Mod[cornerx – 1 + dim, size – 1, 1]]]
(*+map[[cornery+1-mid[dim],Mod[cornerx-1+mid[dim],size-1,
1]]]*))
/2.(*+RandomReal[{jitter1,jitter2}]*),

If[cornery + 1 – mid[dim] <= 0,(*northpole*)
new =
(map[[cornery, cornerx]]
+ map[[cornery, Mod[cornerx – 1 + dim, size – 1, 1]]]
(*+map[[cornery-1+mid[dim],Mod[cornerx-1+mid[dim],size-1,
1]]]*))
/2.(*+RandomReal[{jitter1,jitter2}]*),

new =
(map[[cornery, cornerx]]
+ map[[cornery, Mod[cornerx – 1 + dim, size – 1, 1]]]
+
map[[cornery – 1 + mid[dim],
Mod[cornerx – 1 + mid[dim], size – 1, 1]]]
+
map[[cornery + 1 – mid[dim],
Mod[cornerx – 1 + mid[dim], size – 1, 1]]])
/4. + RandomReal[{jitter1, jitter2}]]];
(map[[cornery, Mod[cornerx – 1 + mid[dim], size – 1, 1]]] = #) &@
new;
,
{i, 1, Length[diamondcorners]}];
If[read == 1, Export[file, map, "HarwellBoeing"]];
Print["Completed dimension " <> ToString[dimlist[[i]]] <> " of " <>
ToString[dimlist] <> " at " <> ToString[DateString[]]],
{i, 1, Length[dimlist]}];
]

Here is an example of it being run. The options are

  1. order of size, n, determines pixel dimensions of map at 2^n+1 pixels per side
  2. lower bound on the range of initial jitter
  3. upper bound on the range of initial jitter
  4. lower bound on jitter function
  5. upper bound on jitter function, jitter is a random real between the lower and upper bounds
  6. tell the code if you want it to save the file will each iteration, 1 if you do
  7. path to a file you want to save your sparse array, .cua file recommended, if you choose a .txt, .tab, or .csv it will take a while to write big files.
arraymapmaker[8, -1, 1, -100/(i*6)^2,
100/(i*6)^2, 1, "~/Desktop/map7.cua"]


map = Import["~/Desktop/map7.cua"];
normmap = (map – Min[map])/Max[map – Min[map]];
colorf1 =
If[# < 0.25, ColorData["DarkTerrain"][0.06],
If[# < 0.3, ColorData["DarkTerrain"][0.12],
"DarkTerrain"~Blend~(#)]] &;

ArrayPlot[normmap, Mesh -> None, ImageSize -> Medium,
ColorFunction -> colorf1, Frame -> False]

map7

 

 

 

 

 

 

 

 

These maps turned out well for what I wanted though the biggest challenge faced was managing the size. A detailed map would require an order, n, of 10 or 11. 10 was manageable but at 2049×2048 pixels order 11 was taking ages on my computer.

I did look into parallelizing the Do loops with ParallelDo[] but it seems the sparse arrays make this impossible¬†or I couldn’t find the right solution.

I was able to come up with a solution that worker for me and reduced the time significantly. My main purpose for this is to make global maps that get wrapped onto a sphere. A standard projection just needs an image with an aspect ratio of roughly 2:1. A limitation of the diamond-square algorithm is that it needs to have sides of 2^n+1 pixels but actually it can be 2^n+1 in one direction and (2^n)/2+1 in the other meaning I can make 2:1 maps with a slight code modification. Essentially all I needed to do was change every instance where the code calls is changing in the north-south direction and asks for the size (2^n+1) I just have to tell it the half size ((2^n)/2+1).

map8
A 2:1 diamond-square map made in half the time!

 

My code above has the following features that not all diamond-square algorithms have.

  • It is seamless if wrapped east-west, the Mod[]’s take care of this nicely.
  • The north border and south borders where the poles are created talk more closely to the ones nearer the pole, the makes the poles a little more consistent.
  • I have it save out the sparse array every time it completes the diamond and square passes for each order. This makes it easy to pick up where it left off in case something causes it to crash.

 

The color functions are fun to play with on these maps. Here are a few of mine.



cf3 = If[# < 0.2, ColorData["DarkTerrain"][0.08],
If[# < 0.28, ColorData["DarkTerrain"][0.1],
If[# < 0.3, ColorData["DarkTerrain"][0.12],
If[# > 0.95, ColorData["AlpineColors"][0.95],
If[# > 0.51 && # < 0.54, ColorData["DarkTerrain"][0.1],
If[# > 0.50 && # < 0.55, ColorData["DarkTerrain"][0.12],
If[# > 0.71 && # < 0.72, ColorData["DarkTerrain"][0.1],
If[# > 0.7 && # < 0.73, ColorData["DarkTerrain"][0.12],
If[# > 0.56 && # < 0.69, "DarkTerrain"~Blend~(1 – #),
If[# > 0.85 && # < 0.95, "AlpineColors"~Blend~#,
"DarkTerrain"~Blend~#]]]]]]]]]] &;


 map11

cf5 = If[# < 0.15, ColorData["DarkTerrain"][0.08],
If[# < 0.2, ColorData["DarkTerrain"][0.1],
If[# < 0.23, ColorData["DarkTerrain"][0.12],
"CoffeeTones"~Blend~(1 – #)]]] &;

map12

 

 

 

 

 

 

cf8 =
If[# < 0.6, ColorData["DarkTerrain"][0.06],
If[# < 0.7, ColorData["DarkTerrain"][0.12],
"SandyTerrain"~Blend~(#)]] &;

map10

 

 

 

 

 

 

cf9 =
If[# < 0.5, ColorData["DarkTerrain"][0.06],
If[# < 0.6, ColorData["DarkTerrain"][0.12],
"PigeonTones"~Blend~(#)]] &;

map9

 

 

 

 

You’re also able to make bump maps and ocean masks very easily if you need to. A bump map would simply be the Automatic color function to give height (may need to invert depending on your own colors) and an ocean mask would be a color function that sets black to anything under your ocean elevation and white to everything else.