Erosion noise implementation in C#

ErosionNoise.cs

/*
Erosion noise implementation in C# (with the Unity.Mathematics package),
ported one-to-one from the Shadertoy shader code by Fewes and clayjohn.

Please note that while the majority of the code is provided under the MIT license,
the quite central "erosion" function is adapted from code shared under the
Creative Commons Attribution-NonCommercial-ShareAlike 3.0 Unported License.

Shadertoy comment by Fewes from https://www.shadertoy.com/view/7ljcRW :

	This shader is a fork of clayjohn's awesome "Eroded Terrain Noise":
	https://www.shadertoy.com/view/MtGcWh

	The original function seemed to have some biasing problems which I've attempted to fix.
	I also wanted to improve the visualization to really show how great the result is.
	I take no credit for the actual noise/math. I simply wanted to help showcase it.

Shadertoy comment by clayjohn from https://www.shadertoy.com/view/MtGcWh :

	This shader is the result of a long time dreaming of a noise function that looked like 
	eroded terrain, complete with branching structure, that could be run in a single pass 
	pixel shader. I wanted to avoid anything simulated because then you cannot easily make
	infinite terrains. 

	A word about the method used. I found guil's excellent "Gavoronoise" shader awhile back, 
	it is has a beautiful wavy look. I noticed that the direction of the waves was a based on
	mouse input. So I combined it with iq's wonderful gradient noise with analytical 
	derivatives. First I generate a heightmap with normals using FBM noise, based on the iq's 
	gradient noise. Then I use the curl of the derivatives to choose the direction for input 
	to guil's gavoronoise. This creates the effect of erosion running down the sides of hills.
	Lastly, I compute the analytic derivatives of the erosion noise and add the curl of it to
	the curl of the hills normals for each iteration, that way each layer of the erosion noise
	changes direction based on the previous layer, creating a branching effect.

	The noise and normals are generated in Buffer A. Image is just used to display the output.

	To see what the heightmap looks like as a terrain comment out line 31. I have not done
	anything to prettify the output, it is just a heightmap with simply phong shading.

	Credit to user guil for "Gavoronoise" (https://www.shadertoy.com/view/llsGWl) and to iq 
	for "Noise - Gradient - 2D - Deriv" (https://www.shadertoy.com/view/XdXBRH)
*/

// Copyright 2020 Clay John

// Permission is hereby granted, free of charge, to any person obtaining a copy of this software 
// and associated documentation files (the "Software"), to deal in the Software without restriction, 
// including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, 
// and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do 
// so, subject to the following conditions:

// The above copyright notice and this permission notice shall be included in all copies or 
// substantial portions of the Software.

// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT 
// NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. 
// IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, 
// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE 
// SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

using Unity.Mathematics;

public struct ErosionNoise {

	// Erosion parameters.
	const float EROSION_TILES = 4;
	const int EROSION_OCTAVES = 5;
	const float EROSION_GAIN = 0.5f;
	const float EROSION_LACUNARITY = 2.0f;

	// Scale the input slope, leading to more erosion.
	const float EROSION_SLOPE_STRENGTH = 3.0f;
	// Continuously modify the noise direction based on the previous fractal sample.
	// This is what gives the slopes an interesting "branching" structure.
	// A higher value will give you more branches.
	const float EROSION_BRANCH_STRENGTH = 3.0f;
	// Maximum amount the erosion will modify the base height map.
	const float EROSION_STRENGTH = 0.04f;

	// Base height noise parameters.
	const float HEIGHT_TILES = 3.0f;
	const int HEIGHT_OCTAVES = 3;
	const float HEIGHT_AMP = 0.25f;
	const float HEIGHT_GAIN = 0.1f;
	const float HEIGHT_LACUNARITY = 2.0f;

	const float WATER_HEIGHT = 0.45f;

	public float2 heightmap(float2 uv) {
		float2 p = uv * HEIGHT_TILES;

		// FBM terrain
		float3 n = float3.zero;
		float nf = 1.0f;
		float na = HEIGHT_AMP;
		for (int i = 0; i < HEIGHT_OCTAVES; i++) {
			n += noised(p * nf) * na * new float3(1.0f, nf, nf);
			na *= HEIGHT_GAIN;
			nf *= HEIGHT_LACUNARITY;
		}

		// [-1, 1] -> [0, 1]
		n.x = n.x * 0.5f + 0.5f;

		// Take the curl of the normal to get the gradient facing down the slope
		float2 dir = n.zy * new float2(1.0f, -1.0f) * EROSION_SLOPE_STRENGTH;

		// Now we compute another fbm type noise
		// erosion is a type of noise with a strong directionality
		// we pass in the direction based on the slope of the terrain
		// erosion also returns the slope. we add that to a running total
		// so that the direction of successive layers are based on the
		// past layers
		float3 h = float3.zero;

		float a = 0.5f;
		float f = 1.0f;

		// Smooth valleys
		//a *= (smoothstep(0.0, 1.0, n.x));

		a *= math.smoothstep(WATER_HEIGHT - 0.1f, WATER_HEIGHT + 0.2f, n.x);

		int octaves = EROSION_OCTAVES;

		for (int i = 0; i < octaves; i++) {
			h += erosion(p * EROSION_TILES * f, dir + h.zy * new float2(1.0f, -1.0f) * EROSION_BRANCH_STRENGTH) * a * new float3(1.0f, f, f);
			a *= EROSION_GAIN;
			f *= EROSION_LACUNARITY;
		}

		return new float2(n.x + (h.x - 0.5f) * EROSION_STRENGTH, h.x);
	}
	
	// code adapted from https://www.shadertoy.com/view/llsGWl
	// name: Gavoronoise
	// author: guil
	// license: Creative Commons Attribution-NonCommercial-ShareAlike 3.0 Unported License
	// Code has been modified to return analytic derivatives and to favour 
	// direction quite a bit.
	float3 erosion(in float2 p, float2 dir) {
		float2 ip = math.floor(p);
		float2 fp = math.frac(p);
		float f = 2.0f * math.PI;
		float3 va = float3.zero;
		float wt = 0.0f;
		for (int i = -2; i <= 1; i++) {
			for (int j = -2; j <= 1; j++) {
				float2 o = new float2(i, j);
				float2 h = hash(ip - o) * 0.5f;
				float2 pp = fp + o - h;
				float d = math.dot(pp, pp);
				float w = math.exp(-d * 2.0f);
				wt += w;
				float mag = math.dot(pp, dir);
				// Multiplying pp by zero is a modification in Fewes' version compared to clayjohn's.
				va += new float3(math.cos(mag * f), -math.sin(mag * f) * (pp * 0.0f + dir)) * w;
			}
		}
		return va / wt;
	}

	float2 hash(in float2 xIn) {
		float2 k = new float2(0.3183099f, 0.3678794f);
		float2 x = xIn * k + k.yx;
		return -1.0f + 2.0f * math.frac(16.0f * k * math.frac(x.x * x.y * (x.x + x.y)));
	}

	// from: https://www.shadertoy.com/view/XdXBRH
	// name: Noise - Gradient - 2D - Deriv
	// author: iq
	// license: MIT
	// return gradient noise (in x) and its derivatives (in yz)
	float3 noised(in float2 p) {
		float2 i = math.floor(p);
		float2 f = math.frac(p);

		float2 u = f * f * f * (f * (f * 6.0f - 15.0f) + 10.0f);
		float2 du = 30.0f * f * f * (f * (f - 2.0f) + 1.0f);

		float2 ga = hash(i + new float2(0.0f, 0.0f));
		float2 gb = hash(i + new float2(1.0f, 0.0f));
		float2 gc = hash(i + new float2(0.0f, 1.0f));
		float2 gd = hash(i + new float2(1.0f, 1.0f));

		float va = math.dot(ga, f - new float2(0.0f, 0.0f));
		float vb = math.dot(gb, f - new float2(1.0f, 0.0f));
		float vc = math.dot(gc, f - new float2(0.0f, 1.0f));
		float vd = math.dot(gd, f - new float2(1.0f, 1.0f));

		return new float3(va + u.x * (vb - va) + u.y * (vc - va) + u.x * u.y * (va - vb - vc + vd),
			ga + u.x * (gb - ga) + u.y * (gc - ga) + u.x * u.y * (ga - gb - gc + gd) +
			du * (u.yx * (va - vb - vc + vd) + new float2(vb, vc) - va));
	}
}