Simulating erosion

I decided to create a procedurally-generated 3D landscape in a CAD program, and wrap it around a globe, which led me to an investigation of erosion algorithms. I include a JavaScript erosion demo at the end of this article for you to play with.

As an aside, let me say that I like OpenSCAD in spite of its idiosyncrasies. Other CAD programs can do many things that OpenSCAD can't, but OpenSCAD is the only 3D modelling software I know of that makes it easy to create procedural or algorithmic parametric designs. It runs on all my computers: Windows, macOS, and Linux — and it's freeware. Its main idiosyncrasy is that it uses a declarative language that requires familiarity with functional programming, a different programming mindset in which all values are evaluated at compile-time instead of run-time, meaning that all "variables" are effectively constants. Even so, one still has conditional branching, looping, and recursive functions, so one can get stuff done. For example, calculating a sum of values over an array isn't done in a loop, it's done in a function calling itself recursively until it reaches a termination condition.

Anyway, to make terrain on a globe, I first need to create a fractal terrain algorithm. It's simple conceptually but it seemed daunting to do it in a declarative language, because one can only initialize array elements (in element order), and not change them afterward. But I managed.

Generating the terrain

The fractal terrain algorithm is simple: Start with a 2×2 grid of random elevations, then subdivide it into a new grid (3×3), filling in the missing grid points as the average of its neighbors plus a random elevation increment that is scaled smaller by $2^{-n}$ for each subdivision level n. Repeat to get a 5×5 grid, then 9×9, 17×17, 33×33, etc. up to $(2^n+1)\times (2^n+1)$.

Here's an example I created in OpenSCAD for n=8.

(Click on any image in this article to see it full size, and click your "Back" button to return to the article.)

A surface generated this way has $2(2^n+1)^2$ triangular facets. At n=8 subdivisions, that's 131,098 facets. Nine subdivisions would create over half a million facets, and a surface with 10 subdivisions would have over 4 million facets. Eight subdivisions is the largest practical value for making an object for 3D printing.

By the way, that image above took several days for me to figure out how to do in OpenSCAD. I was quite pleased with myself when I finally got it to work.

While fractal terrains generated this way look reasonably natural, with hills, valleys, craggy mountain tops, shorelines, islands, and so on, they don't have a weathered look. There is no erosion due to wind or rain. I want to address this.

In this article I cover two techniques for creating eroded landscape. The first one I discarded but it's interesting because it simulates a physical process that has some unpredictability in it. The second one is a more efficient and predictable image processing technique.

Raindrops falling down

About a decade ago, I remember reading a blog series by Shamus Young about procedural terrains, in which he described a simple algorithm for erosion as follows:

For each [x,y,z] value on the map, one at a time, drop an imaginary water droplet on that location.
Find which of the neighboring coordinates is lowest, compared to the droplet's current location. Move the droplet to that new lower location.
Every time the droplet moves, remove a bit of elevation away from where it was.
When the droplet cannot move any lower, abandon it and start at step 1 with a new droplet at a new location. Repeat for all locations over the entire map.

Clearly, the paths taken by successive droplets aren't independent. The path of each new droplet is affected by the erosion from previous droplets. That means the results you get depend on the order in which the droplets fall. Because this is a simulation, it wouldn't be realistic to drop the droplets in a nice orderly fashion by rank and file across the map. A better approach would be to scramble the order of coordinates, so that each point gets one droplet during a pass over the map, but in a random sequence.

Unfortunately, this algorithm is impossible (or at least highly impractical) to implement in OpenSCAD, because multiple droplets can pass over the same point, and you can't change array values once set. To change a value, you must make a copy of the whole array with the changed value in it. Or you'd need to manage unwieldy lists of coordinate values somehow.

But let's explore this anyway. (I'll do this this in JavaScript, using the free Plotly library for 3D surface plots, as well as David Bau's seedrandom library because JavaScript doesn't provide a seedable random number generator for repeatable random sequences. I created an interactive demo at the end of this article to play with.)

From the procedure above, it isn't immediately apparent that removing "a bit of elevation" doesn't mean "a constant bit of elevation." When I tried this, the algorithm dug small deep holes all over my terrain. You see, if a droplet is one step away from its final resting place where it can't go lower, and the elevation step to get there is less than the "bit of elevation" normally removed, then when the droplet moves, it digs out the previous coordinate deeper than the place where it was supposed to stop. Upon finding that the location it just left is now lower, the algorithm then moves the drop back there, digging out the previous location, which becomes lower still... and the droplet ends up digging a hole until a termination condition (minimum z or maximum number of allowable steps). This hole is then widened by subsequent droplets that wander near to it. Here is an example:

So step 3 above should be amended to say that the "bit of elevation" means either a standard tiny elevation increment, or the next elevation step, whichever is less.

Aside from fixing that, this algorithm digs narrow channels. One doesn't get tributaries merging into widening riverbeds like real-world water erosion. Here's a before-and-after illustration:

A nearly trivial improvement is to account for sediment. As a drop flows down, it collects a unit of "sediment" equal to unit of heigh reduction that the drop inflicted on a coordinate. When a drop cannot move anymore, it deposits two units of sediment where it is, and then continues trying to move. It deposits two units because if it can move, it takes one unit away when it moves. Instead of abandoning the drop when it cannot move anymore, I abandon it when it runs out of sediment. The result is interesting:

I like the way that crater at the top left of the image fills up to create a flat bottom. All "material" is preserved this way; any height taken away from one place gets added back in another place lower down.

But let's move on from this. While it's interesting, I can't use this algorithm for my application, as I explained above.

Weathering via image processing

Let's delve instead into image processing. A terrain map is basically a grayscale image, with elevations being represented by shades of gray. An image processing kernel is a 3×3 convolution matrix, performing operations that slide sequentially over each pixel to accomplish blurring, edge, detection, sharpening, and so on. I need something similar: a processing kernel that performs erosion.

The technique known as grayscale erosion looks promising. It works like this:

In the input image, for any grid point [x,y,z], find the minimum z value in a 3×3 grid centered on [x,y].
Set the output image [x,y] coordinate to that minimum z value.
Move on to the next grid point in the input image and repeat until the output image is complete.

That looks simple, fast, and just what I need. Unlike the raindrop simulation, each kernel operation is independent. Although it can't be performed "in place" on one array like the raindrops, that's fine because one can prepare the entire output array by operating on each pixel sequentially without needing to follow slopes — which is compatible with how OpenSCAD requires me to work with arrays.

Here is the reference terrain, the starting configuration used as input for processing:

When I applied this algorithm to the terrain, I didn't notice any significant changes, just a slight softening of some features. So I looped the output back into the input, making 10 passes in total, and I got this:

That sure is ugly! Interestingly, details including some rough spots were preserved. It certainly did get eroded, but in a weird way.

The first thing I need to do is get rid of the square kernel, and make it round, but without making it bigger. This is simple enough to do by interpolation. Imagine inscribing a circle in the 3×3 kernel, moving along a diagonal line from the center to the corner, stopping at the circle, one grid unit along that diagonal. Per the Pythagorean theorem, this is $\frac{1}{\sqrt 2}$ of the diagonal distance, and we use the same fraction of elevation change between the center and the corner. That is, if z₀ is the elevation of the kernel's center, z_c is the elevation of a corner, then the new interpolated corner elevation z_p is:

$$z_p = z_0 + \frac{z_c-z_0}{\sqrt 2} \quad \approx z_0 + 0.707(z_c-z_0)$$

Again using 10 passes over the landscape, this modified round kernel looks better than the square one:

However, it now looks as if someone took a melon baller to the landscape and scooped out divots from it.

Those divots are likely due to the aggressive height reduction imposed by the kernel. Instead of setting the center coordinate to the minimum value of the kernel, how about setting it to partway between the center and minimum values? I can adjust it for more or less weathering simply by changing the reduction proportion, where 0% is no erosion and 100% gives the image above. Here's what I get when I reduce the elevation of the kernel center by 20% of the distance between the kernel's center and lowest elevation, still doing 10 passes over the landscape:

Now that looks pretty good! Compared to the original landscape, it certainly looks well-weathered, not by water erosion, but by wind. Sharp features like mountain peaks and ridges are kept while rough regions are smoothed, just like one would expect.

This has possibilities. I couldn't really improve it further, although I found that if I want a windswept sand-dune look, I make the adjustment percentage variable instead of fixed. That is, use 0% reduction if the kernel center has the lowest height, use 100% reduction if the center has the highest elevation, and proportional values in between. Here is the "windswept sand dune look":

In fact, just two passes of this "windswept sand dune" kernel produces a landscape that is visually indistinguishable from 10 passes of the fixed-percentage reduction method. So I believe I've hit on an efficient algorithm for weathering my fractal landscape.

Back to OpenSCAD

Remember the first image in this article, showing a procedural landscape generated in OpenSCAD? Well, it was fairly easy to apply the processing kernel algorithm to it, although it involves making a new copy of the elevation array for each new pass. Here are the before-and-after pictures (the first being the same as the one at the top of this article). The eroded version uses three passes of the algorithm.

Original landscape

Weathered landscape

So, one big step is done, on the way to creating a globe with a procedural landscape. I plan to make the globe by generating a fractal landscape on each face of a cube, making sure each adjacent face shares a common edge, and project each point out onto a sphere, with the elevation of each point representing a change to the sphere's radius. I expect that doing erosion smoothly over the boundaries between the projected cube faces will present a big challenge. I'll stop here for now, as I have accomplished the first challenge.

Here is a 3D-printed model that I created using my OpenSCAD script, and published on Thingiverse:

Play with it yourself!

As I promised in the first paragraph, here is a JavaScript applet for you to play with. It lets you set a random number seed to generate a fractal terrain and optionally apply erosion passes to it. Have fun.

Random seed (any text):
Number of subdivisions:
Erosion passes (0-20):

Left mouse to rotate, right mouse to pan, scroll wheel to zoom

Search This Blog

NABLU