28
$\begingroup$

From this (now deleted) question I found this site where the author discusses a simple technique for random terrain generation on a sphere.

The method discussed is as follows:

  • start with a collection of points and its center, $c$
  • take some random perturbation $p\in(0, 1)$ and offset, $\lambda $, (I missed this originally) and random vector $v$
  • find the hemispheres created by the plane through $c+\lambda v$ normal to $v$, call them $h_1$ and $h_2$
  • for every point in $h_1$, move it $p$ percent further away from $c$, for every point in $h_2$ move it $p$ percent closer to $c$

How can I implement this? (bonus points for efficiency)

$\endgroup$
  • 1
    $\begingroup$ related: Generate landscape by cutting a plane in 3d $\endgroup$ – Kuba Oct 29 '17 at 10:38
  • $\begingroup$ @Kuba good point. I had the implementation, then thought to ask the question (because who knows what I'm missing in mine). Next time I'll give it a day or so. $\endgroup$ – b3m2a1 Oct 29 '17 at 10:39
  • $\begingroup$ @HighPerformanceMark ah yeah the poster of that question deleted it. I'll note that. $\endgroup$ – b3m2a1 Oct 30 '17 at 8:01
  • 1
    $\begingroup$ This is a classic: Complex planetary surface creation with perlin noise $\endgroup$ – RobC Oct 30 '17 at 14:51
  • $\begingroup$ @RobC interesting. I'd love to see a Mathematica implementation of that. $\endgroup$ – b3m2a1 Oct 30 '17 at 14:55
19
$\begingroup$

Here my two cents. I observed that the major part of the computation is about multiplication. Hence I transformed to logarithms such that we can use summations which can be executed efficiently with Dot. Moreover, I replace the If clause with the listable Sign. Thus, the working horse function looks like this:

getErodedPoints =  Compile[
  {{pt, _Real, 1}, {center, _Real, 1}, {logp, _Real, 1}, {offsets, _Real, 1}, {v, _Real, 2}},
  center + (pt - center) Exp[logp.Sign[v.(pt - center)-offsets]],
  RuntimeAttributes -> {Listable}, Parallelization -> True
];

By the way, here an implementation for a random uniform distribution of points on the 3-sphere:

RandomUnitVector3D[n_] := With[{
   cf = Compile[{{X, _Real, 1}},
     {Cos[X[[2]]] Power[1 - X[[1]]^2, 1/2], Power[1 - X[[1]]^2, 1/2] Sin[X[[2]]], X[[1]]},
     RuntimeAttributes -> Listable, Parallelization -> True
     ]},
  cf[Transpose[{RandomReal[{-1, 1}, n], RandomReal[{-Pi, Pi}, n]}]]
  ]

After these preparations, we can generate our new class-M-planet as follows:

R = DiscretizeRegion[Sphere[], MaxCellMeasure -> 0.0000001];
pts = MeshCoordinates[R];
steps = 25000;
logp = RandomReal[{-.0001, .0001}, steps];
v = RandomUnitVector3D[steps];
center = ConstantArray[0., 3];
offsets = RandomReal[{0., 1.}, steps];
npts = getErodedPoints[pts, center, logp, offsets,v];
r = Sqrt[Dot[npts^2,ConstantArray[1., {3}]]];
Graphics3D[{EdgeForm[],
  GraphicsComplex[
   npts, 
   Polygon[Developer`ToPackedArray[MeshCells[R, 2][[All, 1]]]],
   VertexColors -> ColorData["AlpineColors"] /@ (Rescale[r]^2)
   ]
  },
 Lighting -> "Neutral",
 Boxed -> False
 ]

enter image description here

The call to getErodedPoints with a sphere of about 200000 points and with 25000 steps takes about 10 seconds on my machine.

$\endgroup$
  • $\begingroup$ You can even just do Rescale[r]. The MinMax[r] and {0, 1} are implicit. I think it was common enough in the plotting routines they just made it default. $\endgroup$ – b3m2a1 Oct 30 '17 at 15:05
  • $\begingroup$ @b3m2a1 You're right. This might be quite useful in the future... Thank's for the hint! $\endgroup$ – Henrik Schumacher Oct 30 '17 at 15:10
  • $\begingroup$ As noted here, there is a relatively simple method for generating uniformly-distributed random unit vectors: v = Normalize /@ RandomVariate[NormalDistribution[], {steps, 3}];. (But, since it goes inside a Sign[] anyway, you can omit normalization.) Also, r = Sqrt[Total[npts^2, {2}]]; is much more compact. $\endgroup$ – J. M. will be back soon Oct 31 '17 at 7:44
  • $\begingroup$ @J.M. In a perfect world, Sqrt[Total[npts^2, {2}]] would be faster than Sqrt[Dot[npts^2, ConstantArray[1., {3}]]]. But since this is not the case at the moment (see mathematica.stackexchange.com/questions/158345), I do not rely any more on the depth specs of Total. Normalize is pretty slow, but using the normal distribution for the (unscaled) directions is a very good idea! $\endgroup$ – Henrik Schumacher Oct 31 '17 at 12:19
22
$\begingroup$

I have a basic implementation of this which I use with the output of DiscretizeGraphics:

sphTerrainGenCore =
  Compile[
   {
    {tcoords, _Real, 2},
    {coords, _Real, 2},
    {center, _Real, 1},
    {perturbation, _Real},
    {offset, _Real}
    },
   MapIndexed[
    center +
      (coords[[#2[[1]]]] - center)*
       (1 + perturbation*If[#[[1]] > center[[1]] + offset, 1, -1]) &,
    tcoords
    ]
   ];
sphTerrainGenStep[{coords_, cells_, center_}, {normal_, perturbation_,
     offset_}] :=
  {
   sphTerrainGenCore[
    RotationTransform[{normal, {1, 0, 0}}, center]@coords,
    coords,
    center,
    perturbation,
    offset
    ],
   cells,
   center
   };
sphTerrainGenStep[
   {coords_, cells_, center_},
   steps_Integer,
   perturbationBounds : {_, _} : {.00001, .001},
   offsetBounds : {_, _} : {0, .1}
   ] :=
  Fold[
   sphTerrainGenStep[#, #2] &,
   {coords, cells, center},
   Transpose@{
     RandomReal[{-1, 1}, {steps, 3}],
     RandomReal[perturbationBounds, steps],
     RandomReal[offsetBounds, steps]
     }
   ];
sphTerrainGenStep[{r_?RegionQ, center_},
  steps_Integer,
  perturbationBounds : {_, _} : {.00001, .001},
  offsetBounds : {_, _} : {0, .1}
  ] :=
 With[{ret =
    sphTerrainGenStep[{
      MeshCoordinates[r],
      MeshCells[r, All],
      center
      },
     steps,
     perturbationBounds
     ]
   },
  {
   MeshRegion[ret[[1]], ret[[2]]],
   ret[[3]]
   }
  ]

I then initialize a MeshRegion to work with:

Options[sphTerrainGenInit] =
  Options@DiscretizeGraphics;
sphTerrainGenInit[
   pointNum : _Integer : 10000,
   center : {_?NumericQ, _?NumericQ, _?NumericQ} : {0, 0, 0},
   rad : _Real : 1,
   ops : OptionsPattern[]
   ] :=
  {DiscretizeGraphics[Ball[center, rad], ops], center};
Options[sphTerrainGen] =
  Options@sphTerrainGenInit;
sphTerrainGen[
  steps : _Integer : 100,
  perturbationBounds : {_, _} : {.00001, .001},
  ops : OptionsPattern[]
  ] :=
 With[{base = sphTerrainGenInit[ops]},
  MeshRegion @@ Take[sphTerrainGenStep[base, steps], 2]
  ]

I've then got a basic caching generator and a plotting function:

Options[planetTerrainDataCached] =
  Options[sphTerrainGenInit];
planetTerrainDataCached[0, stepSize_Integer, ops_] :=

  planetTerrainDataCached[0, stepSize, ops] =
   sphTerrainGenInit[FilterRules[{ops}, Options@sphTerrainGenInit]];
planetTerrainDataCached[step_Integer, stepSize_Integer, ops_] :=

  planetTerrainDataCached[step, stepSize, ops] = 
   sphTerrainGenStep[planetTerrainData[step - 1, stepSize, ops], 
    stepSize];
Options[planetTerrainData] =
  Options[planetTerrainDataCached];
planetTerrainData[step_Integer, stepSize : _Integer : 100, 
   ops : OptionsPattern[]] :=

  With[{o = 
     SortBy[Flatten@
       FilterRules[{ops}, Options@planetTerrainDataCached], First]},
   planetTerrainDataCached[step, stepSize, o]
   ];
Options[planetTerrain] =
  Options[planetTerrainData];
planetTerrain[step_Integer, stepSize : _Integer : 100, 
  ops : OptionsPattern[]] :=

 planetTerrainData[step, stepSize, ops][[1]]

Options[planetTerrainPlot] =
  Join[
   Options[SliceDensityPlot3D],
   Options[planetTerrain]
   ];
planetTerrainPlot[i_Integer, stepSize : _Integer : 100, 
   ops : OptionsPattern[]] :=
  planetTerrainPlot[
   planetTerrain[i, stepSize, 
    FilterRules[{ops}, Options[planetTerrain]]],
   ops
   ];
planetTerrainPlot[reg_?RegionQ, ops : OptionsPattern[]] :=

 With[{rb = RegionBounds[reg]},
  SliceDensityPlot3D[
   Norm[{x, y, z}],
   reg,
   {x, rb[[1, 1]], rb[[1, 2]]},
   {y, rb[[2, 1]], rb[[2, 2]]},
   {z, rb[[3, 1]], rb[[3, 2]]},
   Sequence @@
     FilterRules[
      {
       ops,
       ColorFunction -> "AlpineColors",
       Boxed -> False,
       Axes -> False
       },
      Options[SliceDensityPlot3D]
      ] // Evaluate
   ]
  ]

Put this all together:

planetTerrainPlot[1, MaxCellMeasure -> .0001]

planet

We can also animate the steps of the algorithm:

slides1 =
  Map[planetTerrainPlot[#, 1, ImageSize -> 250] &, Range@25];
slides1 // ListAnimate

anim

We can really see how it's a sort-of addition of hemispheres.

Finally we can go to large numbers of steps and see just how craggy it becomes:

planetTerrainPlot[25, MaxCellMeasure -> .0001]

big steps

That's after ~25,000 steps. We could always decrease the cragginess some by playing with the perturbationBounds parameter I defined in my step function, though.

To get more visible terrain contrast, we can take one of our less-craggy meshes and scale its points based on the histogram of Norm values we have:

Rescale[Norm /@ 
   MeshCoordinates[
    planetTerrain[3, MaxCellMeasure -> .0001]]] // SmoothHistogram

hist

This might lead us to do something like this:

regs = {
     #,
     MeshRegion[
      MapThread[
       Which[
         #2 > .99, #3*#,
         #2 < .98, #4*#,
         True, #
         ] &,
       {
        MeshCoordinates[#],
        Rescale[Norm /@ MeshCoordinates[#]],
        RandomReal[{1.01, 1.03}, Length@MeshCoordinates[#]],
        RandomReal[{.97, .98}, Length@MeshCoordinates[#]]
        }],
      MeshCells[#, All]
      ]
     } &@
   planetTerrain[3, MaxCellMeasure -> .0001];
regs // Map[planetTerrainPlot] // Row

craggled

At a minimum, the method described by the author in the link can be implemented like this

$\endgroup$
6
$\begingroup$

What I have is a vastly simpler version, with the slight tweak of possibly using smaller perturbations of the spherical mesh.

To simplify things, start from a discretized (unit) ball (centered on the origin) and extract the points and polygons:

m0 = BoundaryDiscretizeRegion[Ball[], MaxCellMeasure -> {1 -> 0.02}];
pts = MeshCoordinates[m0]; polys = MeshCells[m0, 2];

(You can do scaling and translation later if you need it.)

From here, the code to perform the slicing + scaling is remarkably simple (I take only a few steps n and use an exaggerated value for h for now so that the slicing action is clearly seen):

BlockRandom[SeedRandom[42];
            With[{n = 10, h = 0.01},
                 Do[v = RandomVariate[NormalDistribution[], 3]; p = RandomReal[h];
                    pts += p Sign[pts.v] pts, {n}]];]

Graphics3D[{EdgeForm[], GraphicsComplex[pts, polys]}, Boxed -> False]

first few steps

With n = 1000 and h = 0.002 we see something much more craggy:

some fake planet

Coloring can be done like so:

Graphics3D[{EdgeForm[], 
            GraphicsComplex[pts, polys, 
                            VertexColors -> (ColorData["AlpineColors"] /@
                                             Rescale[Norm /@ pts])]}]

fake alpine planet

Alternatively, replace pts with MeshCoordinates[m0] if you only want to use the perturbations for coloring:

just the coloring


A compiled version goes like this:

makePlanet = Compile[{{msh, _Real, 2}, {h, _Real}, {n, _Integer}},
                     Module[{pts = msh, p, v},
                            Do[v = RandomReal[NormalDistribution[], 3]; p = RandomReal[h];
                               pts += p Sign[pts.v] pts, {n}];
                            pts]];

which should be faster for a larger number of iterations:

BlockRandom[SeedRandom[42]; pp = makePlanet[pts, 0.0005, 10000];
            Graphics3D[{EdgeForm[], 
                        GraphicsComplex[pp, polys, 
                                        VertexColors -> (ColorData["AlpineColors"] /@
                                                         Rescale[Norm /@ pp])]},
                       Boxed -> False, Lighting -> "Neutral"]]

really craggy planet

$\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.