Normal Calculation

[mattdesl]

Hey Nick, thanks for your blog posts and code. Been learning a lot!

I haven't managed to get any of your code running on my macOS (some issues with g++-13 when trying to run make) but have had some success in porting the code to JavaScript which I'm more familiar with. It runs very slowly in JS; so I might try AssemblyScript/WASM to speed things up, or see if I can do it GPU side with WebGPU/WebGL.

I'm a bit stuck on normal map calculation at the moment, it appears they have a dramatic effect on visuals and performance. I also see you've gone through a variety of them, is there one that is better suited? I also don't fully understand what number to choose for "scale" as it seems pretty arbitrary—is there anything that could be relative to the physical size of my map, or how the camera is viewing it? Is there a reason #3 appears to result in the smoothest normals and mountain peaks, but you've moved away from that in more recent projects?

1. soillib - layermap

2. SimpleHydrology - cellpool

3. SimpleErosion - world (I find this to be the "smoothest" so far)

Here's some differences just when I change the normal function from one to another:

Cheers!

[weigert]

Hi Matt. Great question. Seriously. This might be a long answer.

Normal Calculation

I have indeed gone through lots of iterations of computing the surface normals, because this is a relatively complex subject which might even warrant its own mini blog-post.

Quick Problem Statement

1. A surface is represented by a 2D array of points containing floating point values. These are considered samples of an actual, true, underlying surface which we assume is smooth.
2. The surface normal function takes a 2D vector and returns a normalized 3D vector based on the grid data.
3. The 2D vector can be between grid-points. Therefore, the surface normal function consists of computing the normal at a grid-point and an interpolator.
4. What is the ideal surface normal function in our case?

Observations

1. Ideally, our surface-normal function should reflect the true, underlying smooth surface.
2. To visualize the surface, we would typically resort to some kind of meshing, which is not smooth.
3. A grid-based surface representation cannot be meshed by quads. Instead, it would typically be meshed by triangles with a certain orientation (SE-NW diagonal vs. SW-NE diagonal). This is a source of arbitrary bias.
4. The surface normal will have a very strong influence on an erosion simulation because it dictates the particles movement as influenced by gravity.

Aside: The convex-hull of 4 points on a square with arbitrary vertical displacement will be simplectical (a distorted tetrahedron). You can pick either the top two or bottom two triangles of this simplex to represent your surface, with each yielding the one diagonal (SE-NW vs. SW-NE). Hence the arbitrariness. The surface-grid is actually kind-of meshable by quads if you can rotate them arbitrarily. If you rotate the points along the boundaries of the square linearly, eventually the diagonals will necessarily intersect, leading to a flat quad. Nice. Not useful here though.

So basically, the surface which we visualize (e.g. a triangle mesh) and the surface which we assume is underlying and smooth will be different. I have been wrestling with this problem, even thinking about dropping meshes entirely and using SDF-lets for a linear interpolated surface or something... but I digress.

The differences in the methods are basically changes in how I think about the true underlying surface, how I interpolate, whether I care about accuracy in the visualization, and of course trade-offs between speed and accuracy (and smoothness).

Method 1

The first version you linked is actually a deprecated piece of code from SoilMachine.

The grid-node normal function assumes a triangle-fan around the node point to its (4) neighbors. It computes the normals of all 4 triangles around the point and weights them equally (handling boundaries with the k value). The interpolator is just a pure linear interpolator.

This method assumes a triangle mesh around all grid-nodes. It drops the diagonal bias problem by orienting all triangles around itself. This introduces a worse problem, namely that I am linearly interpolating between four entirely different surfaces. And it's (relatively) expensive.

Method 2

Actually identical to the previous version, just hard to recognize as such.

Method 3

This was very early-days, but was better founded in mathematics. What I do here is I compute the normal as a weighted sum normals to the surface gradient along the various directions.

The issue with this method is that it is unnecessarily expensive (computing lots of normalizations) and with slightly arbitrary weights (they sum to one but aren't accurate otherwise). Also, it doesn't work on map boundaries.

Method 4

The method I recommend you use is actually here: surface.hpp

This method is similar to method 3, in that it computes the normal-vector to the surface gradient. The difference is that I actually compute the surface-gradient using finite-differences, which was better mathematical properties. Also, it is much more efficient. Overall this is best.

The last two methods ignore the triangle mesh entirely and just focus on the grid-points.

Scale

The scale parameter is indeed "seemingly" arbitrary. The actual effect that it has in my system, is that it gives a scaling to the floating point values of the grid, assuming that the x-y distance between each grid-node is exactly 1.0. It thereby relates the map height to its width (visually, and in the simulation). Of course, the floating point values themselves don't necessarily exist in a nice range like [-1, 1], so it becomes more arbitrary that way. You can consider the scale parameter a max-height value for the map.

I made lots of efforts to "de-dimensionalize" my model. It's not perfect, but relating it to an actual physical scale is currently not possible unless you introduce a bunch of assumptions. Instead, the scale parameter interacts with the remaining physics.

In particular, the two aspects most affected by the scale are:

1. Cascading / Talus Angle (ratio between horizontal and vertical distance scale)
2. Rate of Erosion / Deposition / Entrainment (because of exponential approach proportional control offset error)

Of course, it also affects the normal calculation, because each grid-node distance is 1, and the height is scaled.

This de-dimensionalization is something I am actively working on improving for future versions.

My recommendation for flat-ish terrain is a scale of about 40, for mountains + flat terrain 80 and for extreme mountains and deep valleys about 120.

Speeding it up

I think you would have some difficulties doing it GPU side, because the particle-based algorithm doesn't lend itself well to pure horizontal scaling. Access patterns are extremely random, and the tradeoff cost of the atomic random reads and writes might kill the benefit. But honestly, it would have to be benchmarked. The physics could be re-written to a eulerian frame instead of a Lagrangian one like I use, but that makes handling certain other things such as momentum conservation much more difficult to achieve (e.g. if you want to get those rivers to meander nicely). There are other things which can potentially be done, like a concurrency pattern instead with a spatial mutex.

Are you running on a new M1 Mac? I am sure you have already tried install the correct g++ version with brew... brew install gcc@13. I'd be happy to help you get the c++ version running.

btw. this is the repository where I will be developing most of the actual underlying library which runs the simulations from now on. At least that's the plan.

Also I took a look at your website and like what you're doing. Let me know if you ever want to collaborate on something.

[mattdesl]

Hey @weigert thanks so much for this incredibly detailed answer! 🤯

I've modified my code to use your new normal calculation and it seems pretty nice and much faster.

In terms of performance, I've mostly solved this with WASM + AssemblyScript. Right now the entire erosion is happening in WASM which is much faster than JS, and I've optimized your functions a fair bit. Here is my AssemblyScript function based on your above Method 4, for reference, and I've added a scale parameter back in (assuming I did that right?).

function sub_and_length(ax: u32, ay: u32, bx: u32, by: u32): f32 {
  let dx: u32 = bx - ax;
  let dy: u32 = by - ay;
  return Mathf.sqrt(f32(dx * dx + dy * dy));
}

function vec_compute_normal(x: u32, y: u32, scale: f32): void {
  norm_x = norm_y = norm_z = 0;

  // horizontal negative
  let pxa_x: u32 = x <= 0 ? x : x - 1;
  let pxa_y: u32 = y;
  // horizontal positive
  let pxb_x: u32 = x >= columns - 1 ? x : x + 1;
  let pxb_y: u32 = y;
  //vertical negative
  let pya_x: u32 = x;
  let pya_y: u32 = y <= 0 ? y : y - 1;
  // vertical positive
  let pyb_x: u32 = x;
  let pyb_y: u32 = y >= rows - 1 ? y : y + 1;

  norm_x =
    (-(
      load<f32>((pxb_x + pxb_y * columns) << 2) -
      load<f32>((pxa_x + pxa_y * columns) << 2)
    ) /
      sub_and_length(pxa_x, pxa_y, pxb_x, pxb_y)) *
    scale;
  norm_y = 1;
  norm_z =
    (-(
      load<f32>((pyb_x + pyb_y * columns) << 2) -
      load<f32>((pya_x + pya_y * columns) << 2)
    ) /
      sub_and_length(pya_x, pya_y, pyb_x, pyb_y)) *
    scale;

  let len: f32 = norm_x * norm_x + norm_y * norm_y + norm_z * norm_z;
  if (len > 0) {
    let sq: f32 = Mathf.sqrt(len);
    norm_x /= sq;
    norm_y /= sq;
    norm_z /= sq;
  }
}

For reference: this old "Method 3" normal calculation took 250 seconds to calculate with a terrain size of 256x256 (using the same defaults but an optimized version of your code here), and only 14 seconds with your new normal map function. Using a map scale of 40.0 for both.

However, I notice the new function is much more "spikey." Am I doing something wrong? Or have you got something you do to avoid the spikes?

After, with new normal function:

My next steps after this is to try to add climate and your layer maps data structure to my AssemblyScript engine.

I'm on a M2 mac air, will have to do some more digging and maybe open a specific issue in one of your repos to try and get things running.

And thanks! A collaboration would be cool at some point. 😃

[weigert]

Nice nice. Looks good.

It appears that the one thing your simulation is currently missing is the cascading step. I don't mention this in either article on the erosion or the original SimpleHydrology model, but only in my article on SimpleWindErosion.

The latest implementation can be found here, and a description of the process can be found here.

Since I wrote the first two articles on erosion and hydrology, I have done a lot of work and I now consider cascading a key part of the dynamics for the hydrology. Any kind of sediment particle has settling. This is actually what gives lots of mountains their shapes. It's natures smoothing function (like a blur with an activation energy). The effect is extremely obvious for sand, hence why I introduced it then, but it quickly became clear the hydrology needed it too. I just haven't written an update since then. I might do soon.

Note that since the original blog/post, the implementation has since changed a bit. For stability purposes, I no longer have a sorted pass but instead approach an average. This is more "well behaved".

If you do everything right, you should get something that looks more like this:

The image above is a [960, 960] resolution erosion simulation run using this (a run-only + raw data export version of SimpleHydrology, currently the most up-to-date version), visualized using this. The colors are the normal map computed by "method 4".

Note that there are no sharp spikes, and that the slopes all naturally converge to the "angle of rest", leading to that nice jaggedly fractal geometric look.

Relief shaded map generated with this

[weigert]

Also: Yes you have added the scale parameter back in correctly. The reason you do not see it in my normal function implementation, is because I made the normal calculator more abstract. i.e. whichever "height" function the method receives is the height function it uses. It is up to the implementer to scale the values inside the height function appropriately, the normal function doesn't care either way.

The reason why I wrote the code in such a convoluted fashion is to keep generic implementations agnostic towards the system they operate on, by only requiring certain constraints statically at compile time. In the case of the normal function, it requires a map which provides an .oob(vec2) and an .height(vec2). If that is satisfied at compile time, it can run.

Basically just C++ wankery. Still, I personally find it very legible and reuseable.

[mattdesl]

Thanks again for this help! Makes sense about the need for cascading and your approach to coding it all in this style.

I was able to install and run the latest SimpleHydrology code in this repo, and also the relief map tool, and now I'm able to get similar results in C++ and also better understand its performance. Unfortunately the others I'm getting stuck on with linking/installing TinyEngine, but I think I have enough to go forward now and keep experimenting with this.

I'll see where I can go with it, and will try to make any new issues/PRs to help if I find errors or areas for optimization along the way. Thanks! :pray:

[weigert]

FYI: TinyEngine is very difficult to compile on Mac at the moment (even for myself), since OpenGL support is weak and the new ARM chips are also causing some issues.

How would you say that the performance compares? I am curious because I have never written anything in WASM.

Please do open issues if anything comes up! Cheers

[mattdesl]

I just published a comparison between JS and WASM with AssemblyScript and Zig. I haven’t gotten around to benchmarking Rust but imagine it will be just as fast if not faster than Zig. Would be interesting to try C as well. From what I understand, it should be pretty close to native.

https://github.com/mattdesl/wasm-bench

It might give you an idea of performance. One thing I notice is that the soil/erosion code needs a lot of zero-allocation structs for vectors that AssemblyScript doesn’t support (but Rust and Zig do). Another note is that Rust should support SIMD with WASM, which may boost some operations like normal calculation and adds/multiplications.