Generate landscape by cutting a plane in 3d

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

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.

• For people interested in terrain generation in general, Bourke discusses a whole lot of methods for this. – J. M.'s technical difficulties Jan 31 '17 at 9:11
• @Mr Wizard Thanks for the corrections, so it is better understood my problem – bullitohappy Jan 31 '17 at 17:15

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


Nesting step does yield a increasingly ragged landscape

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


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

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

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

• Coloring the last one with ColorData["GreenBrownTerrain"] or ColorData["SandyTerrain"] ought to look nice... – J. M.'s technical difficulties Jan 31 '17 at 8:25
• 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 – bullitohappy Jan 31 '17 at 17:22
• @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. – Sascha Jan 31 '17 at 17:51

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

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


This gives the visual effect.

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

ListPlot3D[ImageData[ex]]


The resolution can be improved by generating a larger image.

• +1, for a simpler way to get a 3D plot one could just do ListPlot3D[ImageData[ex]] I think. – user484 Jan 31 '17 at 7:32
• @Rahul Yes, that works very well, I added it to the answer. Thanks. – C. E. Jan 31 '17 at 7:35
• @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. – bullitohappy Jan 31 '17 at 17:39
• @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. – C. E. Jan 31 '17 at 17:55

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_}
] :=
{
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"
]


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]


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


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


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