21
$\begingroup$

Browsing this community, I have noticed that it seems that the friend Vicente and I are working on something similar, but I have seen a web page in Spanish that addresses the same algorithm, only used for the generation of land.

I do not know if the solution which C.E. provided solves my problem. (Please see How to cut a plane at random points?)

The main idea of this algorithm is on each iteration choose two random points on the plane forming a cut line. The side of the plane on the left of the line is raised and the side on the right of the plane lowered. The result has a very great detail quality, with fractal characteristic; Although if you let the algorithm run too many iterations, it ends up canceling itself and returning a flat terrain.

The following image shows some iterations of this algorithm

iterations

Any clue how to implement such algorithm in Mathematica will be welcome, since I have not been able to do anything to attack that problem, thanks in advance.

$\endgroup$
  • 1
    $\begingroup$ For people interested in terrain generation in general, Bourke discusses a whole lot of methods for this. $\endgroup$ – J. M. will be back soon Jan 31 '17 at 9:11
  • $\begingroup$ @Mr Wizard Thanks for the corrections, so it is better understood my problem $\endgroup$ – bullitohappy Jan 31 '17 at 17:15
31
$\begingroup$

My solution does not follow the algorithm you described to the word (I don't choose two random points and find the line equation that runs through them but instead choose a random line equation) but this should not matter for the result.

step[prev_] := 
With[{rand := RandomReal[{-1, 1}]}, 
prev + rand * Sign[rand x + rand y + rand]]

Note x and y are not scoped so they are global (I choose to leave out proper scoping for the sake of readability)

You can generate a single slice with applying step and plot the result with

0 // step // Plot3D[#, {x, -1, 1}, {y, -1, 1}, ExclusionsStyle -> Gray] &

one slice

Nesting step does yield a increasingly ragged landscape

NestList[step, 0, 40] // 
Map[Plot3D[#, {x, -1, 1}, {y, -1, 1}, ExclusionsStyle -> Gray] &] //ListAnimate

gif

After about 200 iterations (which is fast to generate but takes a while to plot) you get something like this:

200 iterations

Advanced usage:

I can imagine that a uniform distribution of the step hight is probably not very realistic and other distributions might archive better results in terms of how realistic the result looks. My intuition tells me that large hight differences must be far less common than small deviations; so a normal distribution might be a better model. To try this out modify step as follows

step[dist_][prev_] := 
With[{stepsize := RandomVariate[dist], rand := RandomReal[{-1, 1}]}, 
prev + stepsize * Sign[rand x + rand y + rand]]

With NormalDistribution[0, 0.02] as parameter dist and after 1500 iterations (and coloring according to the comment by @J.M.) I got these

final final2

The following animation shows the process up to 300 iterations (with only every fourth frame sampled). The same animation is available in a higher resolution on GIPHY

smallGif

$\endgroup$
  • 5
    $\begingroup$ Coloring the last one with ColorData["GreenBrownTerrain"] or ColorData["SandyTerrain"] ought to look nice... $\endgroup$ – J. M. will be back soon Jan 31 '17 at 8:25
  • $\begingroup$ Very good explanation, step by step, so I learn more, thanks for the time dedicated to implementing it and explaining it, interesting of the random line and even more of the distributions in the section of advadced usage. Review re code detail to improve my skills with mathematica $\endgroup$ – bullitohappy Jan 31 '17 at 17:22
  • 3
    $\begingroup$ @user1935664 You might want to unaccept my answer/ be not so quick to accept, to encourage others to post an answer of their own. People around here usually go bonkers for questions like yours. $\endgroup$ – Sascha Jan 31 '17 at 17:51
10
$\begingroup$

We can take the generated image from my other answer and adjust it

Mathematica graphics

and subsequently binarize it with different gray level thresholds and visualize it as a 3D image:

Image3D@Table[ColorNegate@Binarize[ex, cutoff], {cutoff, 0, 1, 0.01}]

enter image description here

This gives the visual effect.

Rahul provided a simpler way in a comment below (this is based on a different image):

ListPlot3D[ImageData[ex]]

Mathematica graphics

The resolution can be improved by generating a larger image.

$\endgroup$
  • $\begingroup$ +1, for a simpler way to get a 3D plot one could just do ListPlot3D[ImageData[ex]] I think. $\endgroup$ – Rahul Jan 31 '17 at 7:32
  • $\begingroup$ @Rahul Yes, that works very well, I added it to the answer. Thanks. $\endgroup$ – C. E. Jan 31 '17 at 7:35
  • $\begingroup$ @C.E. I think ListPlot3D helps a lot in the representation of your images, with bigger images you could get more details, I'll try and see what happens. $\endgroup$ – bullitohappy Jan 31 '17 at 17:39
  • $\begingroup$ @user1935664 That's what I mean by "the resolution can be improved". The edges are jagged because individual pixels are visible. Make the pixels small enough and it will look smooth. $\endgroup$ – C. E. Jan 31 '17 at 17:55
5
$\begingroup$

Here's an extension on Sascha's method to a) CompiledFunctions and b) arbitrary MeshRegions:

planarTerrainFunction[steps : {{_, _}, ___}, 
   init : Except[_List] : 0] :=
  Block[{terrainX, terrainY},
   With[{core =
      Compile @@
       Hold[
        {{terrainX, _Real}, {terrainY, _Real}},
        Evaluate@
         Fold[
          # + #2[[1]]*
             Sign[#2[[2, 1]]*terrainX + #2[[2, 2]]*terrainY + #2[[2, 
                3]]] &,
          init,
          steps
          ],
        RuntimeOptions -> {"EvaluateSymbolically" -> False}
        ]
     },
    With[{
      multiDim =
       Compile[{{clist, _Real, 2}},
        Map[core[#[[1]], #[[2]]] &, clist],
        RuntimeOptions -> {"EvaluateSymbolically" -> False}
        ]
      },
     Switch[Depth[#],
       3, multiDim@#,
       2, core@#,
       _, $Failed
       ] &
     ]
    ]
   ];
planarTerrainFunction[
   dist : _?DistributionParameterQ,
   range : {_?NumericQ, _?NumericQ},
   steps : _?IntegerQ,
   init : Except[_List]
   ] :=
  planarTerrainFunction[
   Transpose@{
     RandomVariate[dist, steps],
     RandomReal[range, {steps, 3}]
     },
   init
   ];
  planarTerrain[
   fn_Function,
   {coords_, cells_}
   ] :=
  {
   MapThread[
    If[Length[#] == 3, {0, 0, #2} + #, Append[#, #2]] &, {#, fn@#} &@
     coords], cells
   };
planarTerrain[
   {
    dist : _?DistributionParameterQ : NormalDistribution[0, 0.02],
    range : {_?NumericQ, _?NumericQ} : {-1, 1},
    steps : _?IntegerQ : 100,
    init : Except[_List] : 0
    },
   {coords_, cells_}
   ] :=
  planarTerrain[
   planarTerrainFunction[dist, range, steps, init],
   {coords, cells}
   ];
planarTerrain[steps_Integer, {coords_, cells_}] :=
  planarTerrain[
   {
    NormalDistribution[0, 0.02],
    {-1, 1},
    steps,
    0
    },
   {coords, cells}
   ];
Options[planarTerrain] =
  Options@DiscretizeRegion;
planarTerrain[
   fn : _Function | _Integer | {
      Optional[_?DistributionParameterQ, 
       NormalDistribution[0, 0.02]],
      Optional[{_?NumericQ, _?NumericQ}, {-10, 10}],
      Optional[_?IntegerQ, 100],
      Optional[Except[_List], 0]
      },
   reg : _?RegionQ | Automatic : Automatic,
   ops : OptionsPattern[]
   ] :=
  With[{mr =
     Replace[reg,
      Automatic :>
       DiscretizeRegion[Rectangle[], ops, MaxCellMeasure -> .001]
      ]},
   MeshRegion @@
    planarTerrain[fn,
     {MeshCoordinates[mr], MeshCells[mr, All]}
     ]
   ];
Options[planarTerrainPlot] =
  Join[
   Options@SliceDensityPlot3D,
   Options@planarTerrain
   ];
planarTerrainPlot[mr_?RegionQ, ops : OptionsPattern[]] :=

  With[{rb = RegionBounds[mr]},
   SliceDensityPlot3D[
    z,
    mr,
    {x, rb[[1, 1]], rb[[1, 2]]},
    {y, rb[[2, 1]], rb[[2, 2]]},
    {z, rb[[3, 1]], rb[[3, 2]]},
    Evaluate[
     FilterRules[{
       ops,
       ColorFunction -> "AlpineColors",
       Boxed -> False,
       Axes -> False
       },
      Options@SliceDensityPlot3D
      ]
     ]
    ]
   ];
planarTerrainPlot[
   tspec : {_Function | _Integer | {_, _, _, _}, _},
   ops : OptionsPattern[]
   ] :=
  planarTerrainPlot[
   planarTerrain[Sequence @@ tspec, 
    FilterRules[{ops}, Options@planarTerrain]],
   ops
   ];
planarTerrainPlot[t : _Function | _Integer | {_, _, _, _}, 
  ops : OptionsPattern[]] :=
 planarTerrainPlot[{t, Automatic}, ops]

Here's how it'd be used to build a standard mountainous terrain:

planarTerrainPlot[1000,
 MaxCellMeasure -> .0001,
 PerformanceGoal -> "Speed"
 ]

mtn

One nice thing about using MeshRegion is that we can easily plot over as Disk or any arbitrary discretizable region:

disk = DiscretizeRegion[Disk[], MaxCellMeasure -> .0001];

t = planarTerrain[2500, disk];

planarTerrainPlot[t]

disk

Finally, since we're working with CompiledFunctions, the actual terrain generation is fast. The display tends to be the rate-limiting step. Here's a fun illustration of how well the time scales in the number of steps:

annulus = DiscretizeRegion[Annulus[], MaxCellMeasure -> .001];

frames =
  Table[
   With[{gen = AbsoluteTiming@planarTerrain[n, annulus]},
    planarTerrainPlot[
     gen[[2]],
     PerformanceGoal -> "Speed",
     ViewPoint -> {1, 0, 1.5},
     PlotLabel -> 
      "`` Steps | Generated in `` seconds"~
       TemplateApply~{n, gen[[1]]}
     ]
    ],
   {n, Prepend[1]@Join[Range[50, 500, 50], Range[500, 1500, 250]]}
   ];

annulus

One last thing. We can use the various special regions to provide an initial shape for our final region to guide the terrain we want:

hill = ImplicitRegion[
   0 <= x <= \[Pi] && 0 <= y <= \[Pi] && z == Sin[x] Sin[y], {x, y, 
    z}] //
  DiscretizeRegion[#, MaxCellMeasure -> .0001] &

hill

We'll make sure to restrict the sampling toward the center of the hill, too, so that things never get too out of control:

frames2 =
  Table[
   planarTerrainPlot[
    planarTerrain[{NormalDistribution[0, Min@{.005*(1 + 1000/n), .2}],
       n, 0}, hill],
    PerformanceGoal -> "Speed",
    ViewPoint -> {1, 1, 2}
    ],
   {n, {10, 50, 100, 250, 500, 750, 1000}}
   ];
frames2 // ListAnimate

hill perturbed

$\endgroup$

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.