HydraulicErosion.js

// See the paper "Implementation of a method for hydraulic erosion" by Hans Theobald Beyer (http://www.firespark.de/resources/downloads/implementation%20of%20a%20methode%20for%20hydraulic%20erosion.pdf) which has examples for all the parameter values
// particleCount number of particles to drop on the terrain
// particleInertia 0 to 1 defines how much changes in the gradient change the direction of the particle. Values near 0.1 work best
// particleCapacity 0 to 100 ish, defines how much sediment each particle can hold onto. Larger values erodes mountains more
// particleDeposition 0 to 1 limits the sediment that is dropped if the sediment carried by the drop exceeds the drops carry capacity
// particleErosion 0 to 1 determines how much of the free capacity of a drop is filled with sediment in case of erosion
// particleEvaporation 0 to 1 defines how fast the particle's waterAmount decreases. Values less than 0.1 work best
// particleRadius, 1 to 10 ish Used during erosion to decide how far away from the particle sediment should be collected. Increasing this value dramatically decreases performance
// particleMinSlope 0 to 1 values should be less than 0.05 in practice
// particleMaxPath maximum number of iterations a particle can perform
// blend 0 to 1, where 0 discards the erosion and 1 blends it with the map 100%. If the erosion looks right and you just want to soften it, use this with lower value
// particleGravity changing this does very little
function HydraulicErosion(data, size, scale, particleCount, particleInertia = 0.1, particleCapacity = 8, particleDeposition = 0.1, particleErosion = 0.1, particleEvaporation = 0.15, particleRadius = 2, particleMinSlope = 0.005, particleMaxPath = 100, blend = 1, particleGravity = 10)
{
	let oldData = null;
	if (blend != 1)
	{
		oldData = new Float32Array(size**2);
		oldData.set(data);
	}
	
	const weightsDimension = particleRadius * 2;
	const weights = new Float32Array(weightsDimension**2);
	
	for (let index = 0; index < particleCount; ++index)
	{
		// Each particle has { positionX, positionY, directionX, directionY, speed, sediment, water } as an element
		let positionX = Math.random() * size;
		let positionY = Math.random() * size;
		const angle = Math.random() * 2 * Math.PI;
		let directionX = Math.cos(angle);
		let directionY = Math.sin(angle);
		let speed = 0;
		let waterAmount = 1;
		let sediment = 0;
		
		let x = Math.floor(positionX);
		let y = Math.floor(positionY);
		let fractionalX = positionX - x;
		let fractionalY = positionY - y;
		
		// Sample at x + xOffset, y + yOffset called heightXY where the XY are the offsets
		let height00 = SampleNoise(data, size, x, y) * scale;
		let height10 = SampleNoise(data, size, x + 1, y) * scale;
		let height01 = SampleNoise(data, size, x, y + 1) * scale;
		let height11 = SampleNoise(data, size, x + 1, y + 1) * scale;
		// Bilinear interpolate height
		let height = (height00 * (1 - fractionalX) + height10 * fractionalX) * (1 - fractionalY) + (height01 * (1 - fractionalX) + height11 * fractionalX) * fractionalY;
		
		for (let step = 0; step < particleMaxPath; ++step)
		{
			const gradientX = (height10 - height00) * (1 - fractionalY) + (height11 - height01) * fractionalY;
			const gradientY = (height01 - height00) * (1 - fractionalX) + (height11 - height10) * fractionalX;
			directionX = directionX * particleInertia - gradientX * (1 - particleInertia);
			directionY = directionY * particleInertia - gradientY * (1 - particleInertia);
			const length = Math.sqrt(directionX**2 + directionY**2);
			if (length < 0.001)
			{
				// Set a random direction
				const angle = Math.random() * 2 * Math.PI;
				directionX = Math.cos(angle);
				directionY = Math.sin(angle);
			}
			else
			{
				directionX /= length;
				directionY /= length;
			}
			// Move a distance of one in the direction
			positionX += directionX;
			positionY += directionY;
			const xNew = Math.floor(positionX);
			const yNew = Math.floor(positionY);
			const fractionalXNew = positionX - xNew;
			const fractionalYNew = positionY - yNew;
			const height00New = SampleNoise(data, size, xNew, yNew) * scale;
			const height10New = SampleNoise(data, size, xNew + 1, yNew) * scale;
			const height01New = SampleNoise(data, size, xNew, yNew + 1) * scale;
			const height11New = SampleNoise(data, size, xNew + 1, yNew + 1) * scale;
			const heightNew = (height00New * (1 - fractionalXNew) + height10New * fractionalXNew) * (1 - fractionalYNew) + (height01New * (1 - fractionalXNew) + height11New * fractionalXNew) * fractionalYNew;
			const heightDifference = heightNew - height;
			const dropSediment = (amount, active00 = 1, active10 = 1, active01 = 1, active11 = 1) =>
			{
				amount = Math.max(0, amount); // The max is just a check against floating point errors
				// Distributed by bilinear interpolation to the 4 points that the position is inside of
				const y0Amount = amount * (1 - fractionalY);
				const y1Amount = amount * fractionalY;
				const w00 = active00 * y0Amount * (1 - fractionalX);
				const w10 = active10 * y0Amount * fractionalX;
				const w01 = active01 * y1Amount * (1 - fractionalX);
				const w11 = active11 * y1Amount * fractionalX;
				SetNoise(data, size, x, y, (height00 + w00) / scale);
				SetNoise(data, size, x + 1, y, (height10 + w10) / scale);
				SetNoise(data, size, x, y + 1, (height01 + w01) / scale);
				SetNoise(data, size, x + 1, y + 1, (height11 + w11) / scale);
				//SetNoise(depositionMap, size, x, y, Math.min(1, (SampleNoise(depositionMap, size, x, y) + w00) / scale));
				//SetNoise(depositionMap, size, x + 1, y, Math.min(1, (SampleNoise(depositionMap, size, x + 1, y) + w10) / scale));
				//SetNoise(depositionMap, size, x, y + 1, Math.min(1, (SampleNoise(depositionMap, size, x, y + 1) + w01) / scale));
				//SetNoise(depositionMap, size, x + 1, y + 1, Math.min(1, (SampleNoise(depositionMap, size, x + 1, y + 1) + w11) / scale));
				sediment -= amount;
			};
			if (heightDifference >= 0)
			{
				// Calculate the amount of sediment required to raise the height up to heightNew using bilinear interpolation:
				// Below is the equation solved where x and y are fractionalX and fractionalY
				// solve((height00 + sedimentRequiredToReachHeight * (1 - x) * (1 - y)) * (1 - x) * (1 - y) + (height10 + sedimentRequiredToReachHeight * x * (1 - y)) * x * (1 - y) + (height01 + sedimentRequiredToReachHeight * (1 - x) * y) * (1 - x) * y + (height11 + sedimentRequiredToReachHeight * x * y) * x * y = heightNew, sedimentRequiredToReachHeight)
				// const sedimentRequiredToReachHeight = (heightNew + height00 * (-fractionalX * fractionalY + fractionalX + fractionalY - 1) + height10 * fractionalX * (fractionalY - 1) + height01 * fractionalX * fractionalY - height11 * fractionalX * fractionalY - height01 * fractionalY) / ((2 * fractionalX**2 - 2 * fractionalX + 1) * (2 * fractionalY**2 - 2 * fractionalY + 1));
				// const dropAmount = Math.min(sedimentRequiredToReachHeight, sediment);
				// dropSediment(dropAmount);
				// This solution has problems though. It can lift up corners higher than the heightNew. Intuitively if you have a square and you want to add sediment to it at a point inside you don't want one of the corners to go higher than the place you're dropping sediment at. This can be solved by locking corners that are higher than heightNew and stopping their weight from causing their respective height to go beyond heightNew.
				// This algorithm is implemented below:
				
				// This iterates 1 to 4 times
				while (sediment > 0.000001)
				{
					// Any corner above Hei
					const active00 = height00 < heightNew - 0.0001 ? 1 : 0;
					const active10 = height10 < heightNew - 0.0001 ? 1 : 0;
					const active01 = height01 < heightNew - 0.0001 ? 1 : 0;
					const active11 = height11 < heightNew - 0.0001 ? 1 : 0;
					// solve((h_0 + l_0 * w * (1 - x) * (1 - y)) * (1 - x) * (1 - y) + (h_1 + l_1 * w * x * (1 - y)) * x * (1 - y) + (h_2 + l_2 * w * (1 - x) * y) * (1 - x) * y + (h_3 + l_3 * w * x * y) * x * y = h, w) where w is sedimentRequiredToReachHeight and h is heightNew
					const sedimentRequiredToReachHeight = (heightNew + height00 * (-fractionalX * fractionalY + fractionalX + fractionalY - 1) + height10 * fractionalX * (fractionalY - 1) + height01 * fractionalX * fractionalY - height11 * fractionalX * fractionalY - height01 * fractionalY) / ((fractionalY - 1)**2 * (active00 * (fractionalX - 1)**2 + active10 * fractionalX**2) + fractionalY**2 * (active01 * (fractionalX - 1)**2 + active11 * fractionalX**2));
					const sedimentToDrop = Math.min(sedimentRequiredToReachHeight, sediment);
					const y0Amount = sedimentToDrop * (1 - fractionalY);
					const y1Amount = sedimentToDrop * fractionalY;
					const w00 = y0Amount * (1 - fractionalX);
					const w10 = y0Amount * fractionalX;
					const w01 = y1Amount * (1 - fractionalX);
					const w11 = y1Amount * fractionalX;
					if ((!active00 || height00 + w00 < heightNew) &&
						(!active10 || height10 + w10 < heightNew) &&
						(!active01 || height01 + w01 < heightNew) &&
						(!active11 || height11 + w11 < heightNew))
					{
						dropSediment(sedimentToDrop, active00, active10, active01, active11);
						break;
					}
					
					// Calculate the sediment required for one of the active corners to reach the heightNew
					let sedimentRequired = Math.min
					(
						active00 ? (heightNew - height00) / ((1 - fractionalX) * (1 - fractionalY)) : 10000,
						active10 ? (heightNew - height10) / (fractionalX * (1 - fractionalY)) : 10000,
						active01 ? (heightNew - height01) / ((1 - fractionalX) * fractionalY) : 10000,
						active11 ? (heightNew - height11) / (fractionalX * fractionalY) : 10000
					);
					if (sedimentRequired <= sediment)
					{
						// Bring the corner up to the heightNew
						dropSediment(sedimentRequired, active00, active10, active01, active11);
					}
					else
					{
						// Not enough sediment for even one of the corners to reach the heightNew, so we can stop and just drop the sediment
						dropSediment(sediment);
						break;
					}
				}
			}
			else
			{
				const capacity = Math.max(-heightDifference, particleMinSlope) * speed * waterAmount * particleCapacity;
				if (sediment >= capacity)
				{
					// Deposition
					const dropAmount = (sediment - capacity) * particleDeposition;
					dropSediment(dropAmount);
				}
				else
				{
					const gainAmount = Math.min((capacity - sediment) * particleErosion, -heightDifference);
					// Erosion happens over a radius.
					// A radius of 1 represents a cell with 4 corners with the particle inside and the corners are weighted by their distance to the particle.
					// A radius of 2 would thus have 16 corners, 3 would have 36, etc. The general formula is used above when allocating the weights array. corners = (particleRadius * 2)**2
					let totalWeight = 0;
					for (let sampleY = 0; sampleY < weightsDimension; ++sampleY)
					{
						const tempY = ((sampleY - (particleRadius - 1)) - fractionalY)**2;
						for (let sampleX = 0; sampleX < weightsDimension; ++sampleX)
						{
							const weight = Math.max(0, particleRadius - Math.sqrt(((sampleX - (particleRadius - 1)) - fractionalX)**2 + tempY));
							weights[sampleY * weightsDimension + sampleX] = weight;
							totalWeight += weight;
						}
					}
					for (let sampleY = 0; sampleY < weightsDimension; ++sampleY)
					{
						const absoluteY = y + sampleY - (particleRadius - 1);
						for (let sampleX = 0; sampleX < weightsDimension; ++sampleX)
						{
							let weight = weights[sampleY * weightsDimension + sampleX];
							weight /= totalWeight;
							// Note that erosion happens at the old position
							const absoluteX = x + sampleX - (particleRadius - 1);
							const erosionFactor = 1; // Lookup on the map. Rock = 0, and Sand = 1. This defines how much a given material erodes. You can even store layers of material say sediment over rock and wear it away or convert rock to sediment among other ideas
							const erosionAmount = weight * gainAmount * erosionFactor;
							SetNoise(data, size, absoluteX, absoluteY, (SampleNoise(data, size, absoluteX, absoluteY) * scale - erosionAmount) / scale);
							sediment += erosionAmount;
							//SetNoise(erosionMap, size, absoluteX, absoluteY, Math.min(1, SampleNoise(erosionMap, size, absoluteX, absoluteY) + erosionAmount));
						}
					}
				}
			}
			
			waterAmount = waterAmount * (1 - particleEvaporation);
			if (waterAmount < 0.01)
			{
				// Evaported
				break;
			}
			
			speed = Math.sqrt(Math.max(0, speed**2 + heightDifference * -particleGravity)); // Note the -particleGravity. I believe the paper has a typo and left this out
			
			x = xNew;
			y = yNew;
			fractionalX = fractionalXNew;
			fractionalY = fractionalYNew;
			height00 = height00New;
			height10 = height10New;
			height01 = height01New;
			height11 = height11New;
			height = heightNew;
		}
	}
	
	if (blend != 1)
	{
		for (let i = 0; i < data; ++i)
		{
			data[i] = oldData[i] * (1 - blend) + data[i] * blend;
		}
	}
}

// Nice example settings for mountains
HydraulicErosion(noise, size, 250, 500000, 0.4, 4, 0.5, 0.4, 0.15, 2, 0.005, 100, 1, 10);
// For certain effects changing the evaporation lower to like 0.01 and such produces large cuts into the terrain. Combinging multiple passes might be ideal or using a map to control parameters might be ideal.