Fractal Planets

George Woodrow

ABSTRACT

Fractal Landscapes were first demonstrated in 1980 by Loren Carpenter, and he worked on the Fractal landscape seen in the “Genesis Effect” in Star Trek II: The Wrath of Khan. In this paper, I will show how to use the functionality in Wolfram Language to generate fractal landscapes for an entire planet effectively and efficiently. The principles and execution shown here may have application to other interesting problems.

Background and Goals

Background

The graphic above is a frame from the first fractal landscape, generated in 1980.

This is a frame from the planet fly-by in Start Trek II: The Wrath of Khan. Shortly after the movie came out, magazines like Creative Computing published code to generate fractal landscapes, which I implemented in Forth and BASIC

Goals

Implement a function to generate a fractal landscape in Wolfram Language that would use arbitrary points for the seed values.

Implement functionality that would generate a fractal planet surface in Wolfram Language.

I want to use built-in functionality in Wolfram Language (WL) to do as much of the “heavy lifting” as possible.

Code Design

In a presentation such as this, the real nature of code design is often lost. I have created an expanded version of this talk, FractalPlanetsExpanded.nb, that gives a flavor of the interactive design process, as well as going into detail that I can not in a short presentation.

New design for a fractal Landscape

Start with a set of points

Make a grid connecting them.

For each segment of the grid, make a new point: The x- and y-values are the midpoint of the line, and the z-value is a random number higher or lower than the average z-value.

Repeat until the desired level of detail has been achieved.

The two tricky bits are the creation of the grid and the code to select new points. Creating a grid from arbitrary points is a non-trivial task, which was avoided in the old code by pre-creating the final grid. Generting the new z-values is a bit of an art, one that is still being explored.

Enhancing the design to use points on a sphere as starting points.

The main design question for generating a fractal ‘planet’ is whether to work with a flattened surface and then map it to a sphere, or to work with points on (or near) a sphere directly.

Work with flat maps and map the result

Initially, this approach seemed to be the best, since the built-in functions like GeoGraphics seemed to be able to map data to a sphere quite well. However, these functions all used data from a real sphere (Earth, Mars, etc.), and what we wanted to do was to generate data. Dealing with the edges of a flat map (which are not edges on the sphere) generated so many problems, and so much ugly code that the effort was abandoned.

Work with the sphere (and spherical co-ordinates) directly

This approach eliminated all the issues dealing with the ‘edges’, since there are no edges when using spherical co-ordinates. After some initial exploration, the only problem appeared to be that many of the built-in functions in Wolfram Language deal with Cartesian co-ordinates only. There would be a lot of conversion from and to Spherical co-ordinates involved.

Summary

After some experimentation, it was obvious that the best approach was to deal with the surface of the sphere directly. Most of the code used to generate the landscape (on a plane) could be adapted to work with the surface of a ‘planet’.

Fractal Landscape

Getting initial points

Initially, I simply created a table of random triples, with the x-, and y- values ranging from 0 to 1 (emulating the old code).

testData = Table[{RandomReal[{0, 1}], RandomReal[{0, 1}], RandomReal[{-.1, .1}]}, {10}];

I quickly discovered that I needed to add points at the corners to make a square starting point. There is code in the Implementation section that provides examples of the starting values.

Making the grid and extracting the lines

The next step in the process is to find line segments that join each point to its ‘nearest’ neighbors. This is a non-trivial task, so I searched the WL documentation for something that would work. I found DelauneyMesh[], and it mostly worked:

pts=RandomReal[{-1,1},{50,2}];
ℛ=DelaunayMesh[pts];
HighlightMesh[ℛ,Style[0,Black]]

The problem with this function is that if the points are in 3-space, the resulting mesh is a solid. The solution (shown in detail in the enhanced notebook) is to

Flatten the data (eliminate the z-co-ordinate).

Make the mesh using DelauneyMesh[].

Extract the mesh lines, using MeshPrimitives[].

Restore the z-co-ordinate.

This procedure gives a list of line segments that make up the grid we need.

Make a new point

The old program used the grid ({x,y} point only) co-ordinate to get the distance between the two points, and had a parameter that governed the ‘roughness’ of the surface. After some experimentation, I came up with a function that worked:

newPoint[lineEndpoints_List, roughness_, zScale_] := Module[{segmentLength},
  segmentLength = EuclideanDistance[Drop[lineEndpoints[[1]], -1], Drop[lineEndpoints[[2]], -1]];
  Midpoint[lineEndpoints] + {0, 0, zScale * RandomReal[{-segmentLength^roughness, segmentLength^roughness}]}
  ]

Putting it together

I wrote a function that incorporated the above steps and made it so that I can easily use Nest[] to iterate as many times as I want:

landscapeDataOneStep[theData_List, generatingFunction_, theRoughness_, theZScale_: 1.0] := Module[{ theMesh, meshLines, lines3D, i, newpoints},
  (* make a 2-D mesh from just the {x,y} co-ordinates of the data *)
  theMesh = DelaunayMesh[Map[Drop[#, -1] &, theData]];
  (* extract the lines from the mesh *)
  meshLines = MeshPrimitives[theMesh, 1][[All, 1]];
  (* make a table by selecting from the data the corresponding triple *)
  lines3D = ParallelTable[{Select[theData, Drop[#, -1] == meshLines[[i, 1]] &][[1]], Select[theData, Drop[#, -1] == meshLines[[i, 2]] &][[1]]}, {i, 1, Length[meshLines]}];
  (* generate new points by adjusting the midpoint of each segment *)
  newpoints = ParallelMap[generatingFunction[#, theRoughness, theZScale] &, lines3D];
  (* return a set of the old plus new points *)
  Union[theData, newpoints]
  ]

I used ParallelMap[] and ParallelTable[] to speed things up. The result is pretty decent:

testIslandterrain = Nest[landscapeDataOneStep[#,  newPoint, 1.0, 1.0]&, testIslandData,5];

ListPlot3D[testIslandterrain,  PlotRange->{0, 1},  Mesh -> All,InterpolationOrder->1, ImageSize->Large, ColorFunction->"GreenBrownTerrain", ClippingStyle->seaColor , FillingStyle->Black, Filling->Bottom]

I could make the display ‘prettier’, but this is sufficient to show that the maths is ok.

Fractal Planet

The process of generating a fractal planet exactly paralleled the process for making the landscape, with some substitutions, and a couple of wrinkles.

Robust co-ordinate transformation

Dealing with Spherical Co-ordinates is a non-trivial task. According to MathWorld, even the convention of the order of the co-ordinates is not uniform. I kept a graphic of WL’s convention handy:

Note that the value for θ, the “latitiude” cannot be equal to 0 or π. If you are using the built-in functions interactively, this is not a problem. However, making thousands of computations iteratively, it is highly likely that the functions will fail to return an ordered triple. (This really does happen.) I made robust versions of these functions that would always work.

Initial Values

Unlike the plane version, we need to start with enough points to keep the results ‘sphere-like’. The function SpherePoints[] works well, but I ended up using the vertices of the Bucky Ball for most of my work. By (my) convention, the points are in spherical co-ordinates.

Making the grid and extracting the lines

Exactly the same issues were present when making the grid on a Sphere. We used ConvexHullMesh[] to generate the lines. This function used all the points when the points were on a Sphere, so we had to do the same trick of flattening — sphere-izing — the points, extracting the lines and restoring the altitudes.

Make a new point

Making a new point for the fractal planet was more involved than the case for the plane. Just finding the midpoint was a challenge until I remembered vector addition. After that, it was easy: The sum of two equal-length vectors bisects the angle between them:

Manipulate[Column[{Text[Style[StringTemplate["point 1 : `` "][{ pt1[[1]]* 180/π, 90 - pt1[[2]]* 180/π}], Blue]],
Text[Style[StringTemplate["point 2 : `` "][{ pt2[[1]]* 180/π, 90 - pt2[[2]]* 180/π}], Green]],Text[Style[StringTemplate["Result : `` "][
{#[[3]] * 180/π, 90 - #[[2]]* 180/π}&[ToSphericalCoordinates[Normalize[Total[Map[Normalize,{FromSphericalCoordinates[{1, pt1[[2]], pt1[[1]]}], FromSphericalCoordinates[{1, pt2[[2]], pt2[[1]]}]}]]]]]], Red]],Graphics3D[{Opacity[0.3],Sphere[],Gray,
Polygon[{{0,0,0},FromSphericalCoordinates[{1, pt1[[2]], pt1[[1]]}] , FromSphericalCoordinates[{1, pt2[[2]], pt2[[1]]}] }],Opacity[1], Blue, Arrowheads[Appearance-> "Projected"],Arrow[{{0,0,0}, FromSphericalCoordinates[{1, pt1[[2]], pt1[[1]]}] }],Green, Arrow[{{0,0,0}, FromSphericalCoordinates[{1, pt2[[2]], pt2[[1]]}] }],
Red, Arrow[{{0,0,0}, 
Normalize[Total[Map[Normalize,{FromSphericalCoordinates[{1, pt1[[2]], pt1[[1]]}], FromSphericalCoordinates[{1, pt2[[2]], pt2[[1]]}]}]]]}]

}, BoxRatios->{1, 1, 1}, SphericalRegion->True, Axes->False, PlotRange->1.1, ImageSize -> 400, Background->Black, Boxed->False]}, Center],{{pt1, {1, .9},Style[ "point 1", Blue, 18]},{ -3.14159, 3.14159}, {3.14159, 0.000001}, ImageSize-> {200, 100}}, {{pt2, {.5, 2.5}, Style[ "point 2", Green, 18]},{ -3.14159, 3.14159}, {3.14159, 0.000001}, ImageSize-> {200, 100}},ControlPlacement->Left, SaveDefinitions->True]

Computing the new altitude was similar to the planar case, although I needed to add a scaling factor; more on this later.

newPointPlanet[lineEndpoints_List, roughness_, zScale_, planRad_] := Module[{newRadius, segLength, cartCoord, newVector},
  (* The line endpoints are in spherical co-ordinates.
  The new altitude is centered on the mean of the two altitudes. *)
  newRadius = (lineEndpoints[[1, 1]] + lineEndpoints[[2, 1]])/2. ;
  (* convert to Cartesian co-ordinates so we can use the built-in VectorAngle[]
  function. *)
  cartCoord = Map[sphTocart, lineEndpoints];
  (* Vector Angle is in Radians, so the arc length is the angle times the radius. *)
  segLength = planRad * VectorAngle[cartCoord[[1]], cartCoord[[2]]];
  (* Iteration 1: multiply segLength by 3.0 <-- empirically arrived at. *)
  segLength *= 3.0;
  (* get the unit vector that bisects the two unit vector *)
  (*
  	working from inside out, we
  1. normalize the vectors in cartesian co-ordinates
  2. add them together. Since the vector are unit
  vectors,the sum bisects the angle between them.
  3. Normalize the resultant vector.
  4. convert back to spherical co-ordinates.
  *)
  newVector = cartTosph[Normalize[Total[Map[Normalize, cartCoord]]]];
  (* add the fractal part to the altitude. This bit will require some
  work to get right. Right now it is the exact analog of the plane code *)
  newRadius += zScale *RandomReal[ {-(segLength^roughness), segLength^roughness}];
  (* newVector is a unit vector, so we can just replace r (=1) with the new redius *)
  {newRadius, newVector[[2]], newVector[[3]]}
  ]

Putting it together

Once again, I made a function that works with Nest[]. Since ListPlot3D does not work with Sperical Co-ordinates, I had to roll my own function display the planet.

buildPlanetOneStep[startPoints_, newPointFunction_, thePlanetRadius_, theMountainScale_, theRoughness_: 1.0] := Module[ {spheredData, cartData, theSphMesh, meshLines, sphMeshLines, restMeshLines, newPts, i},
  (* set all the radii to thePlanetRadius *)
  spheredData = ParallelMap[(x |-> {thePlanetRadius, x[[2]], x[[3]]}), startPoints];
  (* Transform to cartesian co-ordinates *)
  cartData = ParallelMap[sphTocart, spheredData];
  (* Get the convex hull mesh. Since these are all points on the sphere, they will be included in the mesh *)
  theSphMesh = ConvexHullMesh[cartData];
  (* extract the lines from the mesh *)
  meshLines = MeshPrimitives[theSphMesh, 1][[All, 1]];
  (* convert back to spherical co-ordinates. NOTE levelspec for Map[] *)
  sphMeshLines = ParallelMap[cartTosph, meshLines, {2}];
  (* restore the proper radii values *)
  restMeshLines = ParallelTable[{Select[startPoints, Drop[#, 1] == Drop[sphMeshLines[[i, 1]], 1] &][[1]], Select[startPoints, Drop[#, 1] == Drop[sphMeshLines[[i, 2]], 1] &][[1]]}, {i, 1, Length[sphMeshLines]}];
  (* compute the new points NOTE: this requires a working newPointFunction *)
  newPts = ParallelMap[newPointFunction[#, theRoughness, theMountainScale, thePlanetRadius] &, restMeshLines];
  (* Union them with the starting points and return *)
  Union[startPoints, newPts]
  
  ]

SeedRandom[73]

startPlanetBucky = generateStartingPlanetDataFromBuckyBall[simpleInitialCoords, planetRadius, 0.01];

This code fails in 12.0 with 3 iterations. Try 1.

testPlanetbucky1 =Nest[buildPlanetOneStep[#, newPointPlanet, planetRadius, 0.01, 1.0]&,startPlanetBucky, 1];

simplePlanetDisplay[testPlanetbucky1, planetRadius]

We can load in an example generated in 12.1:

SetDirectory[NotebookDirectory[]];
testPlanetbucky = << "testPlanet3iterations.txt";
SetDirectory[];

simplePlanetDisplay[testPlanetbucky, planetRadius]

In the expanded notebook, I show how to take an ‘interesting’ part of the planet and make a projection of it to a plane so that I can use the ‘plane’ version of the landscape code. Here is an example:

In principle, we are done.

Further Development

It would be really great if we could take the data and make it useful for Functions like ReliefPlot[] and GeoElevation data. It would be necessary to obtain a much more detailed set of points, which would take a long time to generate. The advantage is that the full power of Wolfram Language could be leveraged.

Investigate using Compiled functions with the argument Parallelization->True. This should resolve the problem where SeedRandom[] doesn’t completely determine the order of random numbers.

Resolve the issue with the scaling of the segment length when computing a new point.

Finally construct a fly-by movie.

Write a function that that does triangulations of arbitrary points, in both Cartesian and Spherical co-ordinates. The first part has probably already been written in Wolfram Language, since ListPlot3D[] does this as part of its rendering. Having such a function would greatly simplify the code, and speed it up, since there will be no need for co-ordinate transformations.

Observations and Conclusions

Preamble

People use Mathematica or Wolfram Language for a variety of different reasons. For my part, most of my programming experience was is service to some other goal. Starting in the early 1970s, I wrote code to analyze the results of research experiments, to control scientific instruments, and to develop a real-time process control system. I used Mathematica only for the last of these examples; computers in the 1970s were far too wimpy even for Mathematica 1.0.

I have given talks at the WTC for two of my forays into 'computing for fun' (Eliza and Star Trek for Mathematica), but in both these cases, they were ways for me to learn coding better — in 1973 and now.

My motivation for this project was to learn some mathematics, or rather, to see what kind of Mathematics I would need to learn to do something interesting. At this point, it is doubtful that I will be able to master the entire array of new maths that has arisen since I graduated from college. I found this book: Discrete and Computational Geometry, by Devadoss and O’Rourke, but I was having trouble getting into it. Now that I have finished this project, I understand more about Convex hulls and triangulations. I also see that dealing with trigonometry and Vector algebra is critical, not only for this project, but in general.

Another critical motivations is to see how much reliance I can place on Wolfram Language as a ‘black box’ when doing mathematics. The line between doing the work myself and letting the ‘machine’ do it has shifted considerably in my lifetime, so this is a non-trivial exercise. When I took numerical analysis, I used Friden calculators, and I had intimate knowledge of what it meant to invert a 4 by 4 matrix or solve linear equations using Gauss Jordan. For this project, I computed Delaunay meshes without any knowledge of the mathematics involved, to say nothing of the coding.  If I were to do any more work with Delaunay meshes, it would be necessary to learn the maths, especially if I wanted to make the implementation more efficient.

The most important thing I can share form this project is how to bend Wolfram Language to the needs of your project. For example, the built-in functions to do transformations to and from Spherical co-ordinates would not work for my project, but a small modification did. The mesh generating functions needed pre-processing of the data for it to produce the meshes I needed.

The conclusion

“We don’t do a computation to get a number.
 We do it for insight”.

This idea is found all over. Wolfram  Language is very good at giving us a number (graphic, neural network, etc.) When determining whether or not I make too much reliance on technology, I have decided to draw the line at whether I get insight or not.

Implementation

In this section, I put the ‘final’ code for each section of the project, with comments.

Functions for the plane — Cartesian Co-ordinates

Test data

We have several options: a random set of reals, a random set with a fixed boundary, and just the fixed boundary, which is the traditional starting point. The testIsland set has a negative z for the boundary, so that the outer part of the terrain should be mostly under water. Note that the initial seeds never change their z values.

SeedRandom[23];
testData = Table[{RandomReal[{0,1}], RandomReal[{0,1}], RandomReal[{-.1,.1}]}, {10}];
testData2 = Union[testData, {{0.,0.,0.}, {0.,1.,0.}, {1.,1.,0.}, {1.,0.,0.}}];
testIslandData = Join[testData, {{0.,0.,-0.1}, {0.,1.,-0.1}, {1.,1.,-0.1}, {1.,0.,-0.1}}];
testData3 =  {{0.,0.,0.}, {0.,1.,0.}, {1.,1.,0.}, {1.,0.,0.}};
testIslandData2 =   {{0.,0.,-.1},{0., 0.5, -.2}, {0.,1.,-.1},{.3, .3, .2},{.6, .7, .1},{0.44, 0.37, 0.11},  {1.,1.,-.1}, {1.,0.,-.1}, {1.0, 0.5, -0.2},{ 0.5,0.0, -.2}, {0.5, 1.0, -0.2}};

the code

This is the next iteration of the code. We add a parameter for scale, so that the z values will be more realistic when we start to model larger parts of the terrain. All the functions generating a new point need to have the same argument list!

newPoint[lineEndpoints_List, roughness_, zScale_] := Module[{segmentLength},
segmentLength = EuclideanDistance[Drop[lineEndpoints[[1]], -1], Drop[lineEndpoints[[2]], -1]];
Midpoint[lineEndpoints] + {0,0, zScale * RandomReal[{-segmentLength^roughness, segmentLength^roughness}]}
]

We want to save the ‘bad’ function to show what can go wrong.

newPointBlowsUp[lineEndpoints_List, roughness_, zScale_] := Module[{segmentLength},
segmentLength = EuclideanDistance[lineEndpoints[[1]], lineEndpoints[[2]]];
Midpoint[lineEndpoints] + {0,0, zScale * RandomReal[{-segmentLength^roughness, segmentLength^roughness}]}
]

This function  is the next iteration in the development process. It just has the zScale parameter added

landscapeDataOneStep[theData_List, generatingFunction_, theRoughness_, theZScale_:1.0] := Module[{ theMesh, meshLines,lines3D, i, newpoints},
(* make a 2-D mesh from just the {x,y} co-ordinates of the data *)
theMesh = DelaunayMesh[Map[Drop[#,-1]&, theData]];
(* extract the lines from the mesh *)
meshLines = MeshPrimitives[theMesh,1][[All,1]];
(* make a table by selecting from the data the corresponding triple *)
lines3D =ParallelTable[{Select[theData, Drop[#, -1]== meshLines[[i,1]]&][[1]], Select[theData, Drop[#, -1]== meshLines[[i,2]]&][[1]]}, {i, 1,Length[meshLines]}];
(* generate new points by adjusting the midpoint of each segment *)
newpoints =ParallelMap[generatingFunction[#,theRoughness, theZScale]&, lines3D];
(* return a set of the old plus new points *)
Union[theData, newpoints]
]

Planet functions — Spherical Co-ordinates

Nominal values for the planet.

The following are the constants used for the project. The planetRadius is roughly the radius of the earth in km, and the mountainScale is roughly the ratio of the height of Mt. Everest to the radius of the earth.

planetRadius = 6400.0;
mountainScale = 0.0014;
roughness = 0.9;
seaColor = ;

Test data

We use a function to generate the test data. The function that sets up the initial altitudes is passed to the main function so we can try variants.

Instead of generating a random number between ±scFac*rad, we make a random choice of {-0.5*scFac*rad, 0.5*scFac*rad}. We just want the starting points to be just something other than zero altitude, and (from experience) not at the extremes.

simpleInitialCoords[thePointSphCoords_, rad_, scFac_] := Module[{newAlt},
newAlt = thePointSphCoords[[1]] + RandomChoice[{-0.5*scFac*rad, 0.5*scFac*rad}];
{newAlt, thePointSphCoords[[2]], thePointSphCoords[[3]]}
]

More may be added here.

The following function is set up for ease of understanding and modification. Since it is only used once per computation, we can be less efficient. Specifically, we do not use ParallelMap[] so that the same initial configuration will be produced for the same random seed.

generateStartingPlanetData[nPoints_,setUpAltsFunc_, radius_:1.0, scaleFac_:0.0014] := Module[{startPoints, sphPoints, plPoints, endPts},
(* get the initial points. These are in cartesian co-ordinates aroung {0,0,0} *)
startPoints = SpherePoints[nPoints];
(* transform to spherical co-ordinates *)
sphPoints = Map[ToSphericalCoordinates, startPoints];
(* Set the radius to the passed value *)
plPoints = Map[Times[{radius, 1, 1}, #]&,  sphPoints];
(* adjust the altitude for the initial points *)
endPts = Map[setUpAltsFunc[#, radius, scaleFac]&, plPoints];
endPts
]

NOTE: the above function returns a list in spherical co-ordinates.

The following is the starting point for most of the work after the initial investigation. It does not use ParallelMap[] so that the initial configuration will be the same for the same random seed.

generateStartingPlanetDataFromBuckyBall[setUpAltsFunc_, radius_:1.0, scaleFac_:0.0014] := Module[{startPoints, sphPoints,sphPountsNorm, plPoints, endPts},
(* get the initial points. These are in cartesian co-ordinates aroung {0,0,0} *)
startPoints = Map[Normalize,N[PolyhedronCoordinates[UniformPolyhedron["TruncatedIcosahedron"]]]];
(* transform to spherical co-ordinates *)
sphPoints = Map[ToSphericalCoordinates, startPoints];
(* Set the radius to the passed value *)
plPoints = Map[Times[{radius, 1, 1}, #]&,  sphPoints];
(* adjust the altitude for the initial points *)
endPts = Map[setUpAltsFunc[#, radius, scaleFac]&, plPoints];
endPts
]

Planet building

This function takes a list of starting points (in spherical co-ordinates) and returns a new list that has the original points plus the ones constructed from midpoints of the mesh. It makes extensive use of ParallelMap[] and ParallelTable[]. For single core computers, these computations can take a long time.

In version 12.0, ConvexHullMesh[] will return EmptyRegion if too many points are passed to it. This will happen at the second iteration, typically. To get a useful result, version 12.1 is required.

buildPlanetOneStep[startPoints_, newPointFunction_, thePlanetRadius_, theMountainScale_, theRoughness_:1.0] := Module[ {spheredData, cartData, theSphMesh, meshLines, sphMeshLines,restMeshLines, newPts, i},
(* set all the radii to thePlanetRadius *)
spheredData = ParallelMap[(x |->{thePlanetRadius, x[[2]], x[[3]]}), startPoints];
(* Transform to cartesian co-ordinates *)
cartData = ParallelMap[sphTocart, spheredData];
(* Get the convex hull mesh. Since these are all points on the sphere, they will be included in the mesh *)
theSphMesh = ConvexHullMesh[cartData];
(* extract the lines from the mesh *)
meshLines = MeshPrimitives[theSphMesh,1][[All,1]];
(* convert back to spherical co-ordinates. NOTE levelspec for Map[] *)
sphMeshLines =ParallelMap[cartTosph,meshLines, {2}];
(* restore the proper radii values *)
restMeshLines =ParallelTable[{Select[startPoints, Drop[#, 1]==Drop[sphMeshLines[[i,1]], 1]&][[1]], Select[startPoints, Drop[#, 1]== Drop[sphMeshLines[[i,2]], 1]&][[1]]}, {i, 1,Length[sphMeshLines]}];
(* compute the new points NOTE: this requires a working newPointFunction *)
newPts = ParallelMap[newPointFunction[#, theRoughness, theMountainScale, thePlanetRadius]&, restMeshLines];
(* Union them with the starting points and return *)
Union[startPoints, newPts]

]

Here is the function to get a new point from the endpoints of a segment. Because this function is called by a function using ParallelMap[], it cannot use parallelization, since parallelization cannot be nested. Note that robust co-ordinate transformations must be used to prevent errors.

The scaling factor for segment length, 3.0, was determined empirically so that with a roughness of 1.0, there is sufficient roughness in the terrain. As noted in the main text, the function is sensitive to the value chosen for segment length.

newPointPlanet[lineEndpoints_List, roughness_, zScale_, planRad_] := Module[{newRadius, segLength, cartCoord, newVector},
(* The line endpoints are in spherical co-ordinates.
The new altitude is centered on the mean of the two altitudes. *)
newRadius = (lineEndpoints[[1,1]] + lineEndpoints[[2,1]])/2. ;
(* convert to Cartesian co-ordinates so we can use the built-in VectorAngle[]
function. *)
cartCoord = Map[sphTocart,lineEndpoints];
(* Vector Angle is in Radians, so the arc length is the angle times the radius. *)
segLength = planRad * VectorAngle[cartCoord[[1]], cartCoord[[2]]];
(* Iteration 1: multiply segLength by 3.0 <-- empirically arrived at. *)
segLength *= 3.0;
(* get the unit vector that bisects the two unit vector *)
(*
	working from inside out, we
1. normalize the vectors in cartesian co-ordinates
2. add them together. Since the vector are unit
vectors,the sum bisects the angle between them.
3. Normalize the resultant vector.
4. convert back to spherical co-ordinates.
*)
newVector =cartTosph[Normalize[Total[Map[Normalize,cartCoord]]]];
(* add the fractal part to the altitude. This bit will require some
work to get right. Right now it is the exact analog of the plane code *)
newRadius += zScale *RandomReal[ {-(segLength^roughness), segLength^roughness}];
(* newVector is a unit vector, so we can just replace r (=1) with the new redius *)
{newRadius, newVector[[2]], newVector[[3]]}
]

Robust co-ordinate transformation

The following functions are substitutes for ToSphericalCoordinates[] and FromSphericalCoordinates[]. These functions are robust in that they work with the {0,0,1} vector in Cartesian and {1, 0, ϕ} or {1, π, ϕ} in Spherical coordinates. They are not an exact substitute, since they don’t operate over lists.

cartTosph[{x_, y_, z_}] := Quiet[Check[ToSphericalCoordinates[{x,y,z}], {z,0,0}]]

sphTocart[{r_, θ_, ϕ_}] := CoordinateTransformData["Spherical"->"Cartesian","Mapping"][{r, θ, ϕ}]

Rendering the planet

ListPlot3D does not work with Spherical co-ordinates, so we need to deal with graphics primitives. This function is similar to the one that generates note points in that we use ConvexHullMesh[], which means that we need to ‘sphere-ize’ our co-ordinates and then restore the altitudes later. Unlike buildPlanetOneStep[],  the function gets the mesh primitives at level 2, which returns the co-ordinates for the polygons that make up the mesh. Although I have never seen a mesh polygon that was not a triangle, the code can handle for any n-gon.

The function returns a list of polygons in Cartesian co-ordinates that can be rendered by applying color, texture, etc. as desired. This function is used for rendering only, so the output will be used by Graphics3D[], which operates only in Cartesian co-ordinates.

In version 12.0, ConvexHullMesh[] will return EmptyRegion if too many points are passed to it.  To get a useful result, version 12.1 is required.

generateSurfacePolygons[planetPoints_, thePlanetRadius_] := Module[ {restAlt,spheredData, cartData, theSphMesh, meshPolys, sphMeshPolys,restMeshPolys, PolyPts, i},
(*  make function to restore altitudes in an n-gon*)
restAlt[plPoints_, thePoly_] := Module[ {nVertices, j},
nVertices = Length[thePoly];
Table[Select[plPoints, Drop[#, 1] == Drop[thePoly[[j]],1]&][[1]], {j,  nVertices}]
];
(* set all the radii to thePlanetRadius *)
spheredData = ParallelMap[(x |->{thePlanetRadius, x[[2]], x[[3]]}), planetPoints];
(* Transform to cartesian co-ordinates *)
cartData = ParallelMap[sphTocart, spheredData];
(* Get the convex hull mesh. Since these are all points on the sphere, they will be included in the mesh *)
theSphMesh = ConvexHullMesh[cartData];
(* extract the POLYGONS from the mesh *)
meshPolys = MeshPrimitives[theSphMesh,2][[All,1]];
(* convert back to spherical co-ordinates. NOTE levelspec for Map[] *)
sphMeshPolys =ParallelMap[cartTosph,meshPolys, {2}];
(* restore the proper radii values. The polygons are not necessarily triangles *)
restMeshPolys =ParallelMap[restAlt[planetPoints, #]&, sphMeshPolys];
(* convert to Cartesian co-ordinates, since Graphics3D only works with Cartesian *)
PolyPts = ParallelMap[sphTocart,restMeshPolys, {2}];
(* Convert to Polygon objects *)
ParallelMap[Polygon, PolyPts]
]

This is a helper function for generateSurfacePolygons[].

restoreAltitudes[planetPoints_, thePoly_] := Module[ {nVertices, i},
nVertices = Length[thePoly];
Table[Select[planetPoints, Drop[#, 1] == Drop[thePoly[[i]],1]&][[1]], {i,  nVertices}]
]

This is a function that will render a planet with the polygons rendered in the color specified by setColorForPoly[]. This is a rough approximation of what ListPlot3D[] would do, and can be extensively modifies to make a more realistic surface.

simplePlanetDisplay[planetPoints_, thePlanetRadius_] := Module[{polys, elev, stuff},
polys = generateSurfacePolygons[planetPoints, thePlanetRadius];
elev = Max[ParallelMap[elevForPoly[#, thePlanetRadius] &, polys]];
stuff =  Riffle[ParallelMap[setColorForPoly[#, thePlanetRadius,elev ] &, polys], polys];
Graphics3D[{seaColor, Sphere[{0,0,0},thePlanetRadius], White, EdgeForm [None],stuff}]
]

This function returns only the graphics object.

simplePlanetConstruction[planetPoints_, thePlanetRadius_] := Module[{polys, elev, stuff},
polys = generateSurfacePolygons[planetPoints, thePlanetRadius];
elev = Max[ParallelMap[elevForPoly[#, thePlanetRadius] &, polys]];
stuff =  Riffle[ParallelMap[setColorForPoly[#, thePlanetRadius,elev ] &, polys], polys];
{seaColor, Sphere[{0,0,0},thePlanetRadius], White, EdgeForm [None],stuff}
]

“Fast” only in comparison to the above function.

fastPlanetDisplay[planetPoints_, thePlanetRadius_] := Module[{polys, elev, stuff},
polys = generateSurfacePolygons[planetPoints, thePlanetRadius];
Graphics3D[{seaColor, Sphere[{0,0,0},thePlanetRadius], Brown, polys}]
]

This returns the graphics for later display

fastPlanetConstruction[planetPoints_, thePlanetRadius_] := Module[{polys, elev, stuff},
polys = generateSurfacePolygons[planetPoints, thePlanetRadius];
{seaColor, Sphere[{0,0,0},thePlanetRadius], Brown, polys}
]

Helper functions:

elevForPoly[thePoly_, thePlanetRadius_] := Module[{elev},
EuclideanDistance[{0,0,0}, RegionCentroid[thePoly] ]- thePlanetRadius
]

setColorForPoly[thePoly_, thePlanetRadius_, maxElev_] := Module[{elev},
elev = EuclideanDistance[{0,0,0}, RegionCentroid[thePoly] ]- thePlanetRadius;
ColorData[{"GreenBrownTerrain", {0, maxElev}},elev ]
]

Projection

This code takes a point in spherical co-ordinates and maps it to a point on the plane tangent to the ‘north pole’. Note that the resulting z-value is the altitude above or below ‘sea level’.

getProjectedPoint[sphPt_, planetRad_] := Module[{x, y, z, r, angle},
z = sphPt[[1]] - planetRad;
angle = sphPt[[3]];
r = 2 * planetRad *Sin[sphPt[[2]]]/ (1 + Cos[sphPt[[2]]]);
{x,y} = FromPolarCoordinates[{r, angle}];
{x, y, z}
]