islaApocalypse/Tools/Scripts/HydraulicErosion.cs

using System;

/// <summary>
/// Droplet-based hydraulic erosion (terrain-water task 17, Phase C0) — the organic
/// carve-AND-deposit pass, Lague/Beyer lineage. Pure numeric over the height array
/// (D-035; a named future C++ candidate, kept standalone — no Godot types at all,
/// System.MathF only, own deterministic PCG32 RNG).
///
/// Each droplet spawns on land (spawn probability weighted toward high ground),
/// then walks downhill with inertia, carrying water and sediment. Where the ground
/// is steep and it moves fast it ERODES (up to capacity, spread over a small brush
/// so no single-cell spikes — the anti-artifact that killed the D8 predecessor);
/// where it flattens out it DEPOSITS, building valley floors and fans. Water
/// evaporates each step; the droplet dies at its lifetime, at the map edge, or on
/// reaching the sea (its remaining sediment is lost to the ocean).
///
/// OUTPUT-ONLY: this pass is applied to the RENDER height map only; the classify
/// map never sees it (the caller owns that split — see MapGenerator).
///
/// The three hard governors (the pass provably cannot run away):
///   1. DropletCount   — total droplets (the main detail/cost dial).
///   2. Lifetime       — max steps per droplet; no infinite wandering.
///   3. CarveCapM      — max cumulative erosion depth per cell, in metres,
///                       enforced against a per-cell ledger. The runaway-trench
///                       guard, and what keeps this a DETAILING pass.
///
/// The sea clamp (the "don't over-flood" guard): erosion never lowers any cell
/// below its local sea level + SeaMarginM, and cells already below sea are
/// read-only — never eroded, never deposited on. Land stays land, sea stays sea;
/// the rendered coastline cannot move. Deposition only raises land cells.
///
/// The crater exclusion: no cell within CraterExclRadius of the impact centre is
/// modified (droplets may traverse). The carve remains the final authority on its
/// own terrain.
///
/// Heights in the array are raw blueprint units (1 unit = 251 m). All sediment
/// accounting below is done in METRES and converted only when a delta is applied,
/// so untouched cells keep their exact bit pattern — the invariants above are
/// exact, not statistical.
/// </summary>
public static class HydraulicErosion
{
	public const ushort VERSION = 1;

	// Deterministic RNG stream: seeded from resolvedSeed + this offset, so a seed
	// reproduces exactly and the stream is decorrelated from every noise field
	// (7409/8117/… are taken; see MakeModulationNoise call sites).
	public const int SEED_OFFSET = 9271;

	public const float M_PER_UNIT = 251f;

	// Crater exclusion factor: erosion stays outside 1.2 × CraterRadius — fully
	// clear of both the physical carve (0.80×) and the detail feather (1.05×).
	public const float CRATER_EXCL_FACTOR = 1.2f;

	// Spawn: droplets source in the mountains, never the ocean. A land point is
	// accepted with probability SPAWN_FLOOR + (1-SPAWN_FLOOR) · relative elevation,
	// after at most SPAWN_TRIES rejection-sampling attempts (then the droplet is
	// skipped and counted — on any real island this is vanishingly rare).
	private const int SPAWN_TRIES = 16;
	private const float SPAWN_FLOOR = 0.15f;

	private const float MIN_WATER = 0.005f; // droplet dies when effectively dry
	private const float MIN_DIR = 1e-10f;   // below this, direction is re-drawn at random

	public struct Params
	{
		public int DropletCount;      // governor 1
		public int Lifetime;          // governor 2
		public float CarveCapM;       // governor 3 (metres)
		public float SeaMarginM;      // sea clamp margin (metres)
		public int BrushRadius;       // erosion brush radius, px
		public float Inertia;         // 0 = pure gradient descent, 1 = never turns
		public float CapacityFactor;  // sediment capacity multiplier
		public float MinSlopeM;       // capacity slope floor, metres per px
		public float ErodeRate;       // fraction of remaining capacity eroded per step
		public float DepositRate;     // fraction of surplus sediment dropped per step
		public float Evaporation;     // water lost per step (fraction)
		public float Gravity;         // speed gain per metre of drop
		public int Seed;              // resolvedSeed + SEED_OFFSET
	}

	public class Stats
	{
		public int Spawned;
		public int SkippedNoLand;
		public long Steps;
		public int DiedLifetime, DiedEdge, DiedSea, DiedDry;
		public double ErodedVolumeM3;     // 1 px = 1 m², so metres of depth sum to m³
		public double DepositedVolumeM3;
		public float MaxCellErosionM;     // must end ≤ CarveCapM
		public long ErodedCells;          // cells with any net ledger erosion
	}

	// PCG32 (O'Neill) — tiny, deterministic, trivially portable to C++.
	private struct Pcg32
	{
		private ulong _state;
		public Pcg32(int seed) { _state = 0; NextU(); _state += (ulong)(uint)seed; NextU(); }
		public uint NextU()
		{
			ulong old = _state;
			_state = old * 6364136223846793005UL + 1442695040888963407UL;
			uint xorshifted = (uint)(((old >> 18) ^ old) >> 27);
			int rot = (int)(old >> 59);
			return (xorshifted >> rot) | (xorshifted << (-rot & 31));
		}
		public float NextF() => (NextU() >> 8) * (1f / 16777216f); // [0,1)
	}

	/// <summary>
	/// Runs the pass in place on <paramref name="height"/>. Sea level per cell is
	/// <paramref name="seaMap"/>[x,y] when non-null, else the flat scalar
	/// <paramref name="seaFlat"/>. Throws (refusing the generation) if a governor
	/// bound is violated on exit — the caller treats that as a build failure.
	/// </summary>
	public static Stats Apply(float[,] height, int mapSize, float[,] seaMap, float seaFlat,
		float craterCx, float craterCy, float craterExclRadius, Params p)
	{
		var stats = new Stats();
		var rng = new Pcg32(p.Seed);
		float capUnits = p.CarveCapM / M_PER_UNIT;
		if (p.DropletCount <= 0 || capUnits <= 0f) return stats;

		// Per-cell cumulative-erosion ledger — governor 3's enforcement record.
		float[,] eroded = new float[mapSize, mapSize];

		// Spawn weighting needs the seed's top height.
		float hTop = float.MinValue;
		for (int x = 0; x < mapSize; x++)
			for (int y = 0; y < mapSize; y++)
				if (height[x, y] > hTop) hTop = height[x, y];

		// Erosion brush: all offsets within BrushRadius, cone-weighted (1 - d/r),
		// normalized. Radius 0 degrades to the single cell.
		int r = Math.Max(p.BrushRadius, 0);
		int brushN = 0;
		for (int dx = -r; dx <= r; dx++)
			for (int dy = -r; dy <= r; dy++)
				if (MathF.Sqrt(dx * dx + dy * dy) <= r + 1e-4f) brushN++;
		int[] brushDx = new int[brushN], brushDy = new int[brushN];
		float[] brushW = new float[brushN];
		{
			int i = 0; float wSum = 0f;
			for (int dx = -r; dx <= r; dx++)
				for (int dy = -r; dy <= r; dy++)
				{
					float d = MathF.Sqrt(dx * dx + dy * dy);
					if (d > r + 1e-4f) continue;
					brushDx[i] = dx; brushDy[i] = dy;
					brushW[i] = r > 0 ? 1f - d / (r + 1f) : 1f;
					wSum += brushW[i]; i++;
				}
			for (int j = 0; j < brushN; j++) brushW[j] /= wSum;
		}

		float exclSq = craterExclRadius * craterExclRadius;
		float SeaAt(int cx, int cy) => seaMap != null ? seaMap[cx, cy] : seaFlat;
		bool Excluded(int cx, int cy)
		{
			float ddx = cx - craterCx, ddy = cy - craterCy;
			return ddx * ddx + ddy * ddy < exclSq;
		}

		for (int drop = 0; drop < p.DropletCount; drop++)
		{
			// --- spawn (land only, elevation-weighted) ---
			float px = -1f, py = -1f;
			for (int attempt = 0; attempt < SPAWN_TRIES; attempt++)
			{
				float sx = 1f + rng.NextF() * (mapSize - 3);
				float sy = 1f + rng.NextF() * (mapSize - 3);
				int cx = (int)sx, cy = (int)sy;
				float h = height[cx, cy];
				float sea = SeaAt(cx, cy);
				if (h < sea) { continue; }
				float rel = hTop > sea ? Math.Clamp((h - sea) / (hTop - sea), 0f, 1f) : 0f;
				if (rng.NextF() < SPAWN_FLOOR + (1f - SPAWN_FLOOR) * rel) { px = sx; py = sy; break; }
			}
			if (px < 0f) { stats.SkippedNoLand++; continue; }
			stats.Spawned++;

			float dirX = 0f, dirY = 0f, speed = 1f, water = 1f, sedimentM = 0f;

			for (int step = 0; step < p.Lifetime; step++)
			{
				stats.Steps++;
				int xi = (int)px, yi = (int)py;
				float fx = px - xi, fy = py - yi;

				// Bilinear height + gradient at the current position.
				float h00 = height[xi, yi], h10 = height[xi + 1, yi];
				float h01 = height[xi, yi + 1], h11 = height[xi + 1, yi + 1];
				float gradX = (h10 - h00) * (1f - fy) + (h11 - h01) * fy;
				float gradY = (h01 - h00) * (1f - fx) + (h11 - h10) * fx;
				float hOld = h00 * (1f - fx) * (1f - fy) + h10 * fx * (1f - fy)
						   + h01 * (1f - fx) * fy + h11 * fx * fy;

				// Inertia blend, then one unit step.
				dirX = dirX * p.Inertia - gradX * (1f - p.Inertia);
				dirY = dirY * p.Inertia - gradY * (1f - p.Inertia);
				float len = MathF.Sqrt(dirX * dirX + dirY * dirY);
				if (len < MIN_DIR)
				{
					float ang = rng.NextF() * 2f * MathF.PI;
					dirX = MathF.Cos(ang); dirY = MathF.Sin(ang); len = 1f;
				}
				dirX /= len; dirY /= len;
				px += dirX; py += dirY;

				if (px < 1f || px >= mapSize - 2 || py < 1f || py >= mapSize - 2)
				{ stats.DiedEdge++; break; }

				int nxi = (int)px, nyi = (int)py;
				float nfx = px - nxi, nfy = py - nyi;
				float n00 = height[nxi, nyi], n10 = height[nxi + 1, nyi];
				float n01 = height[nxi, nyi + 1], n11 = height[nxi + 1, nyi + 1];
				float hNew = n00 * (1f - nfx) * (1f - nfy) + n10 * nfx * (1f - nfy)
						   + n01 * (1f - nfx) * nfy + n11 * nfx * nfy;

				// Reached the sea: die; the sediment is the ocean's now.
				if (hNew < SeaAt(nxi, nyi)) { stats.DiedSea++; break; }

				float dhM = (hNew - hOld) * M_PER_UNIT;
				float capacityM = MathF.Max(-dhM, p.MinSlopeM) * speed * water * p.CapacityFactor;

				if (dhM > 0f || sedimentM > capacityM)
				{
					// Moving uphill (fill the pit behind us, at most the rise) or
					// over capacity (drop a fraction of the surplus): DEPOSIT at
					// the OLD position, bilinear over its 4 cells.
					float amountM = dhM > 0f ? MathF.Min(dhM, sedimentM)
											 : (sedimentM - capacityM) * p.DepositRate;
					if (amountM > 0f)
					{
						float w00 = (1f - fx) * (1f - fy), w10 = fx * (1f - fy);
						float w01 = (1f - fx) * fy, w11 = fx * fy;
						sedimentM -= DepositCell(height, eroded, xi, yi, amountM * w00, stats, SeaAt, Excluded)
								   + DepositCell(height, eroded, xi + 1, yi, amountM * w10, stats, SeaAt, Excluded)
								   + DepositCell(height, eroded, xi, yi + 1, amountM * w01, stats, SeaAt, Excluded)
								   + DepositCell(height, eroded, xi + 1, yi + 1, amountM * w11, stats, SeaAt, Excluded);
					}
				}
				else
				{
					// Under capacity on a downhill move: ERODE, spread over the
					// brush, never more than the drop itself (no digging pits).
					float amountM = MathF.Min((capacityM - sedimentM) * p.ErodeRate, -dhM);
					if (amountM > 0f)
					{
						for (int b = 0; b < brushN; b++)
						{
							int cx = xi + brushDx[b], cy = yi + brushDy[b];
							if (cx < 0 || cx >= mapSize || cy < 0 || cy >= mapSize) continue;
							if (Excluded(cx, cy)) continue;
							float sea = SeaAt(cx, cy);
							float hCell = height[cx, cy];
							if (hCell < sea) continue; // below-sea cells are read-only
							float want = amountM * brushW[b];
							float bySea = MathF.Max(0f, (hCell - (sea + p.SeaMarginM / M_PER_UNIT)) * M_PER_UNIT);
							float byCap = MathF.Max(0f, p.CarveCapM - eroded[cx, cy]);
							float take = MathF.Min(want, MathF.Min(bySea, byCap));
							if (take <= 0f) continue;
							height[cx, cy] = hCell - take / M_PER_UNIT;
							if (eroded[cx, cy] == 0f) stats.ErodedCells++;
							eroded[cx, cy] += take;
							if (eroded[cx, cy] > stats.MaxCellErosionM) stats.MaxCellErosionM = eroded[cx, cy];
							sedimentM += take;
							stats.ErodedVolumeM3 += take;
						}
					}
				}

				speed = MathF.Sqrt(MathF.Max(0f, speed * speed - dhM * p.Gravity));
				water *= 1f - p.Evaporation;
				if (water < MIN_WATER) { stats.DiedDry++; break; }
				if (step == p.Lifetime - 1) stats.DiedLifetime++;
			}
		}

		// Governor 3, proven on exit rather than assumed: the ledger's maximum must
		// respect the cap (float addition of clamped takes cannot exceed it by more
		// than rounding; allow one ulp-scale epsilon).
		if (stats.MaxCellErosionM > p.CarveCapM * (1f + 1e-5f))
			throw new InvalidOperationException(
				$"[HydraulicErosion] CARVE-CAP VIOLATION: a cell accumulated {stats.MaxCellErosionM} m against cap {p.CarveCapM} m. Refusing to generate.");

		return stats;
	}

	/// <summary>
	/// Deposits up to <paramref name="amountM"/> metres on one cell; returns what was
	/// actually placed. Below-sea cells and crater-excluded cells take nothing —
	/// deposition only ever raises LAND, so the coastline cannot move and the sea
	/// cannot shallow. A cell's ledgered erosion is paid back first, so erode-then-
	/// deposit at one cell frees cap headroom instead of double-counting.
	/// </summary>
	private static float DepositCell(float[,] height, float[,] eroded, int cx, int cy,
		float amountM, Stats stats, Func<int, int, float> seaAt, Func<int, int, bool> excluded)
	{
		if (amountM <= 0f) return 0f;
		float hCell = height[cx, cy];
		if (hCell < seaAt(cx, cy)) return 0f;
		if (excluded(cx, cy)) return 0f;
		height[cx, cy] = hCell + amountM / M_PER_UNIT;
		eroded[cx, cy] = MathF.Max(0f, eroded[cx, cy] - amountM);
		stats.DepositedVolumeM3 += amountM;
		return amountM;
	}
}