islaApocalypse/Tools/Scripts/DrainageAnalysis.cs

using System;
using System.Collections.Generic;

/// <summary>
/// Drainage-network promotion — C0b part 1 (terrain-water task 21). PURE ANALYSIS:
/// reads the ERODED render heightmap and produces a river PLAN — it changes zero
/// terrain and adds zero water. Standalone numeric (D-035 family; no Godot types).
///
/// Pipeline, built on the task-03 priority-flood family:
///   1. Priority-flood the eroded surface from the map border (Barnes heap+pit
///      variant, 8-connected, same as RunPriorityFloodDiagnostics) — but with a
///      one-ulp epsilon on pit fills, so every filled cell keeps a STRICTLY
///      descending path to its spill. This resolves the ~15,000 erosion pits
///      (task-20 finding) for ROUTING ONLY; the terrain itself is never modified.
///   2. Depressions that are deep AND large enough (the endorheic dials) are NOT
///      filled through: their cells revert to original heights, so flow entering
///      them terminates at the basin minimum. Real closed drainage survives;
///      micro-pits route through.
///   3. D8 flow directions on that routing surface. D8 was reverted as a CARVING
///      technique (task 10 — grid-aligned scratches in the terrain); using it to
///      COMPUTE where water flows is standard hydrology and leaves no mark.
///   4. Flow accumulation by topological (Kahn) propagation — no sort needed.
///   5. Promotion: outlets to the sea ranked by drainage area, top-N (separated)
///      become trunks; main stems traced upstream by max-accumulation; the
///      mountain-exit point found from the along-stem grade; LEAN tributaries and
///      LEAN endorheic terminals marked.
///
/// The plan's lowland courses are provisional: erosion delivered the UPLAND
/// network only (task 18 §3), so below each mountain-exit the traced course is
/// "where the routing surface drains", not a designed river. Part 2 (task 22)
/// routes the lowland reach properly from the mountain-exit points — which is why
/// those points are this analysis's key output.
/// </summary>
public static class DrainageAnalysis
{
	public const float M_PER_UNIT = 251f;

	// Neighbour order is FIXED (it is the deterministic tiebreak).
	private static readonly int[] DX = { -1, -1, -1, 0, 0, 1, 1, 1 };
	private static readonly int[] DY = { -1, 0, 1, -1, 1, -1, 0, 1 };
	private static readonly float[] DIST = {
		1.41421356f, 1f, 1.41421356f, 1f, 1f, 1.41421356f, 1f, 1.41421356f };

	public class Params
	{
		// Endorheic qualification: a depression this deep AND this large is a real
		// closed basin and terminates flow; anything smaller is a pit, filled through.
		public float EndorheicMinDepthM = 2.0f;
		public int EndorheicMinAreaPx = 10000;
		// Endorheic REPORTING is lean: only terminals with at least this much
		// upstream drainage, at most MaxCount of them.
		public int EndorheicMinInflowPx = 50000;
		public int EndorheicMaxCount = 3;

		public int TrunkCount = 3;            // ~3 sea-reaching trunks (developer)
		public int MinOutletSeparationPx = 400; // don't pick 3 mouths of one delta

		public int StemMinAccPx = 1000;       // stem tracing stops below this
		public int TributaryMinAccPx = 30000; // LEAN: a branch must drain this much
		public int TributaryMaxPerTrunk = 4;  // ...and only the top few are marked

		// Mountain-exit: furthest-downstream stem point where the upstream window
		// still sustains this grade (m per px) over ExitWindowPx.
		public float ExitGradeMin = 0.05f;
		public int ExitWindowPx = 100;

		public float SeaLevel = 0.15f;        // flat sea scalar (raw units)
	}

	public class Stream
	{
		public List<(float x, float y)> Course = new();  // downstream-first
		public long DrainageAreaPx;
		public (float x, float y) Head;                  // upstream end
	}

	public class Trunk : Stream
	{
		public (float x, float y) Outlet;                // last land cell before sea
		public (float x, float y) MountainExit;
		public float MountainExitElevM;
		public bool ExitFound;
		public List<Stream> Tributaries = new();
	}

	public class EndorheicTerminal
	{
		public (float x, float y) Terminal;              // basin minimum
		public long DrainageAreaPx;
		public float BasinDepthM;
		public long BasinAreaPx;
	}

	public class Plan
	{
		public List<Trunk> Trunks = new();
		public List<EndorheicTerminal> Endorheics = new();
		public int TerminalBasinCount;                   // basins that qualified as sinks
		public long PitsFilledCount;                     // depressions filled through
		public long LandCells, SeaReachingCells, EndorheicCells, UnroutedCells;
		public List<(float x, float y, long acc)> AllOutletsTop = new(); // top 12, pre-separation
		public Params P;
	}

	/// <param name="isOcean">Row-major mask of THE OCEAN body (WBID == 1) — the
	/// only water that counts as "the sea" for sea-reaching trunks. Below-sea
	/// cells that are NOT ocean (enclosed lagoons, below-datum lake beds) are
	/// ordinary terrain to the router: as depressions they either qualify as
	/// terminal basins (a river legitimately ENDING in a lagoon/lake — reported as
	/// such) or fill and spill onward to the true sea. Without this mask the first
	/// draft called two of its three "sea-reaching" trunks done at enclosed
	/// lagoons, which is exactly the overclaim the gate must not inherit.</param>
	public static Plan Run(float[,] height, int mapSize, bool[] isOcean, Params p)
	{
		int n = mapSize;
		int total = n * n;
		var plan = new Plan { P = p };

		// 1-D row-major copies (idx = x * n + y), same convention as the task-03 pass.
		float[] original = new float[total];
		for (int x = 0; x < n; x++)
			for (int y = 0; y < n; y++)
				original[x * n + y] = height[x, y];

		// --- 1. Priority-flood with one-ulp epsilon (routing surface only) ---
		float[] filled = (float[])original.Clone();
		{
			bool[] visited = new bool[total];
			var heap = new PriorityQueue<int, (float h, int idx)>();
			var pit = new Queue<int>();
			void Seed(int idx)
			{
				if (visited[idx]) return;
				visited[idx] = true;
				heap.Enqueue(idx, (filled[idx], idx)); // idx tiebreak => deterministic
			}
			for (int x = 0; x < n; x++) { Seed(x * n); Seed(x * n + (n - 1)); }
			for (int y = 0; y < n; y++) { Seed(y); Seed((n - 1) * n + y); }

			while (heap.Count > 0 || pit.Count > 0)
			{
				int c = pit.Count > 0 ? pit.Dequeue() : heap.Dequeue();
				float fc = filled[c];
				int cx = c / n, cy = c % n;
				for (int k = 0; k < 8; k++)
				{
					int nx = cx + DX[k], ny = cy + DY[k];
					if (nx < 0 || nx >= n || ny < 0 || ny >= n) continue;
					int ni = nx * n + ny;
					if (visited[ni]) continue;
					visited[ni] = true;
					if (filled[ni] <= fc)
					{
						// One ulp above the parent: strictly descending back out, so
						// D8 never meets an exact flat inside a filled pit.
						filled[ni] = MathF.BitIncrement(fc);
						pit.Enqueue(ni);
					}
					else heap.Enqueue(ni, (filled[ni], ni));
				}
			}
		}

		// --- 2. Depression components; big+deep ones become terminal sinks ---
		// Components of (filled > original), 8-connected — the pools. Qualifying
		// pools revert to ORIGINAL height so flow terminates at their minimum.
		int[] basinId = new int[total];       // 0 = not in a pool
		var basinDepthM = new List<float> { 0f };
		var basinAreaPx = new List<long> { 0L };
		var basinMinCell = new List<int> { -1 };
		{
			var stack = new Stack<int>();
			int nextId = 1;
			for (int i = 0; i < total; i++)
			{
				if (basinId[i] != 0 || filled[i] <= original[i]) continue;
				int id = nextId++;
				long area = 0; float depth = 0f; int minCell = i; float minH = original[i];
				stack.Push(i); basinId[i] = id;
				while (stack.Count > 0)
				{
					int c = stack.Pop();
					area++;
					float d = (filled[c] - original[c]) * M_PER_UNIT;
					if (d > depth) depth = d;
					if (original[c] < minH) { minH = original[c]; minCell = c; }
					int cx = c / n, cy = c % n;
					for (int k = 0; k < 8; k++)
					{
						int nx = cx + DX[k], ny = cy + DY[k];
						if (nx < 0 || nx >= n || ny < 0 || ny >= n) continue;
						int ni = nx * n + ny;
						if (basinId[ni] == 0 && filled[ni] > original[ni])
						{ basinId[ni] = id; stack.Push(ni); }
					}
				}
				basinDepthM.Add(depth); basinAreaPx.Add(area); basinMinCell.Add(minCell);
			}

			bool[] terminal = new bool[nextId];
			for (int id = 1; id < nextId; id++)
			{
				if (basinDepthM[id] >= p.EndorheicMinDepthM && basinAreaPx[id] >= p.EndorheicMinAreaPx)
				{ terminal[id] = true; plan.TerminalBasinCount++; }
				else plan.PitsFilledCount++;
			}
			// Revert terminal pools to the real surface; re-tag basinId to keep only
			// terminal pools (routing needs to know "am I in a terminal basin").
			for (int i = 0; i < total; i++)
			{
				if (basinId[i] == 0) continue;
				if (terminal[basinId[i]]) filled[i] = original[i];
				else basinId[i] = 0;
			}
		}

		// --- 3. D8 flow directions on the routing surface ---
		// dir[i] = 0..7 neighbour, SEA (into a below-sea cell), or NONE (sink).
		const sbyte D_NONE = -1, D_SEA = -2;
		sbyte[] dir = new sbyte[total];
		bool IsSea(int idx) => isOcean[idx];
		for (int i = 0; i < total; i++)
		{
			if (IsSea(i)) { dir[i] = D_NONE; continue; }
			int cx = i / n, cy = i % n;
			float best = 0f; int bestK = -1; bool bestIsSea = false;
			for (int k = 0; k < 8; k++)
			{
				int nx = cx + DX[k], ny = cy + DY[k];
				if (nx < 0 || nx >= n || ny < 0 || ny >= n) continue;
				int ni = nx * n + ny;
				float drop = (filled[i] - filled[ni]) / DIST[k];
				if (drop > best) { best = drop; bestK = k; bestIsSea = IsSea(ni); }
			}
			dir[i] = bestK < 0 ? D_NONE : (bestIsSea ? D_SEA : (sbyte)bestK);
		}

		// --- 4. Flow accumulation (Kahn topological propagation) ---
		int Target(int i)
		{
			if (dir[i] < 0) return -1;
			int cx = i / n, cy = i % n;
			return (cx + DX[dir[i]]) * n + (cy + DY[dir[i]]);
		}
		int[] acc = new int[total];
		{
			byte[] indeg = new byte[total];
			for (int i = 0; i < total; i++)
				if (dir[i] >= 0) indeg[Target(i)]++;
			var q = new Queue<int>();
			for (int i = 0; i < total; i++)
			{
				if (IsSea(i)) continue;
				acc[i] = 1;
				if (indeg[i] == 0) q.Enqueue(i);
			}
			while (q.Count > 0)
			{
				int c = q.Dequeue();
				if (dir[c] < 0) continue;
				int t = Target(c);
				acc[t] += acc[c];
				if (--indeg[t] == 0 && !IsSea(t)) q.Enqueue(t);
			}
		}

		// Bookkeeping: where does each cell's flow END — the sea, WHICH terminal
		// basin, or stuck? Memoised downstream walk. The per-basin totals matter:
		// crediting a terminal basin only with acc at its deepest cell undercounts
		// badly when the basin floor is flat (a lagoon bed scatters inflow across
		// many sub-minima — measured: a 500k-px lagoon system reported under 50k).
		long[] basinInflow = new long[basinMinCell.Count];
		{
			int[] dest = new int[total]; // 0 unknown, -1 sea, -2 stuck, >0 basin id
			var path = new List<int>(4096);
			for (int i = 0; i < total; i++)
			{
				if (IsSea(i) || dest[i] != 0) continue;
				int c = i; path.Clear();
				int result;
				while (true)
				{
					if (dest[c] != 0) { result = dest[c]; break; }
					path.Add(c);
					if (dir[c] == D_SEA) { result = -1; break; }
					if (dir[c] == D_NONE) { result = basinId[c] != 0 ? basinId[c] : -2; break; }
					c = Target(c);
				}
				foreach (int pc in path) dest[pc] = result;
			}
			for (int i = 0; i < total; i++)
			{
				if (IsSea(i)) continue;
				plan.LandCells++;
				if (dest[i] == -1) plan.SeaReachingCells++;
				else if (dest[i] > 0) { plan.EndorheicCells++; basinInflow[dest[i]]++; }
				else plan.UnroutedCells++;
			}
		}

		// --- 5a. Outlets: land cells whose flow enters the sea, ranked by acc ---
		var outlets = new List<(int cell, long acc)>();
		for (int i = 0; i < total; i++)
			if (dir[i] == D_SEA) outlets.Add((i, acc[i]));
		outlets.Sort((a, b) => b.acc.CompareTo(a.acc));

		foreach (var (cell, a) in outlets.GetRange(0, Math.Min(12, outlets.Count)))
			plan.AllOutletsTop.Add((cell / n, cell % n, a));

		// Greedy top-N with separation, so three mouths of one delta can't take
		// all three trunk slots.
		var picked = new List<int>();
		foreach (var (cell, _) in outlets)
		{
			if (picked.Count >= p.TrunkCount) break;
			int cx = cell / n, cy = cell % n;
			bool far = true;
			foreach (int pcell in picked)
			{
				float ddx = cx - pcell / n, ddy = cy - pcell % n;
				if (ddx * ddx + ddy * ddy < (float)p.MinOutletSeparationPx * p.MinOutletSeparationPx)
				{ far = false; break; }
			}
			if (far) picked.Add(cell);
		}

		// upstream max-acc walk shared by trunks and tributaries
		List<int> TraceStem(int fromCell, int minAcc)
		{
			var stem = new List<int> { fromCell };
			int c = fromCell;
			while (true)
			{
				int cx = c / n, cy = c % n;
				int bestN = -1; long bestA = minAcc - 1;
				for (int k = 0; k < 8; k++)
				{
					int nx = cx + DX[k], ny = cy + DY[k];
					if (nx < 0 || nx >= n || ny < 0 || ny >= n) continue;
					int ni = nx * n + ny;
					if (dir[ni] >= 0 && Target(ni) == c && acc[ni] > bestA)
					{ bestA = acc[ni]; bestN = ni; }
				}
				if (bestN < 0) break;
				stem.Add(bestN);
				c = bestN;
			}
			return stem;
		}

		List<(float x, float y)> Decimate(List<int> cells, int step = 4)
		{
			var pts = new List<(float, float)>();
			for (int i = 0; i < cells.Count; i += step)
				pts.Add((cells[i] / n, cells[i] % n));
			if ((cells.Count - 1) % step != 0)
				pts.Add((cells[^1] / n, cells[^1] % n));
			return pts;
		}

		// --- 5b. Trunks: stems, mountain exits, LEAN tributaries ---
		foreach (int outletCell in picked)
		{
			var t = new Trunk
			{
				Outlet = (outletCell / n, outletCell % n),
				DrainageAreaPx = acc[outletCell]
			};
			var stem = TraceStem(outletCell, p.StemMinAccPx);
			t.Course = Decimate(stem);
			t.Head = (stem[^1] / n, stem[^1] % n);

			// Mountain-exit: walk the stem downstream-first; the exit is the
			// furthest-DOWNSTREAM point whose upstream window still sustains the
			// grade — i.e. where the mountains hand the river to the flats.
			// Elevation truth is the ORIGINAL eroded surface, not the fill.
			int w = p.ExitWindowPx;
			for (int i = 0; i + w < stem.Count; i++)
			{
				float rise = (original[stem[i + w]] - original[stem[i]]) * M_PER_UNIT;
				if (rise / w >= p.ExitGradeMin)
				{
					t.ExitFound = true;
					t.MountainExit = (stem[i] / n, stem[i] % n);
					t.MountainExitElevM = original[stem[i]] * M_PER_UNIT;
					break;
				}
			}

			// LEAN tributaries: junction branches off the stem with enough drainage,
			// top few by accumulation.
			var stemSet = new HashSet<int>(stem);
			var cands = new List<(int cell, long acc)>();
			foreach (int sc in stem)
			{
				int cx = sc / n, cy = sc % n;
				for (int k = 0; k < 8; k++)
				{
					int nx = cx + DX[k], ny = cy + DY[k];
					if (nx < 0 || nx >= n || ny < 0 || ny >= n) continue;
					int ni = nx * n + ny;
					if (stemSet.Contains(ni)) continue;
					if (dir[ni] >= 0 && Target(ni) == sc && acc[ni] >= p.TributaryMinAccPx)
						cands.Add((ni, acc[ni]));
				}
			}
			cands.Sort((a, b) => b.acc.CompareTo(a.acc));
			// Dedup: two inflow neighbours at adjacent stem cells are one confluence,
			// not two tributaries — keep only junctions ≥ 30 px apart.
			var taken = new List<int>();
			foreach (var (cell, a) in cands)
			{
				if (taken.Count >= p.TributaryMaxPerTrunk) break;
				int cx2 = cell / n, cy2 = cell % n;
				bool dup = false;
				foreach (int tc in taken)
				{
					float ddx = cx2 - tc / n, ddy = cy2 - tc % n;
					if (ddx * ddx + ddy * ddy < 30f * 30f) { dup = true; break; }
				}
				if (!dup) taken.Add(cell);
			}
			foreach (int cell in taken)
			{
				long a = acc[cell];
				var trib = new Stream { DrainageAreaPx = a };
				var ts = TraceStem(cell, Math.Max(p.StemMinAccPx, (int)(a / 20)));
				trib.Course = Decimate(ts);
				trib.Head = (ts[^1] / n, ts[^1] % n);
				t.Tributaries.Add(trib);
			}
			plan.Trunks.Add(t);
		}

		// --- 5c. LEAN endorheic terminals: terminal basins ranked by TOTAL inflow ---
		{
			var terms = new List<(int id, long inflow)>();
			for (int id = 1; id < basinMinCell.Count; id++)
			{
				int mc = basinMinCell[id];
				if (mc < 0 || basinId[mc] != id) continue; // not a terminal basin
				if (basinInflow[id] >= p.EndorheicMinInflowPx) terms.Add((id, basinInflow[id]));
			}
			terms.Sort((a, b) => b.inflow.CompareTo(a.inflow));
			foreach (var (id, inflow) in terms.GetRange(0, Math.Min(p.EndorheicMaxCount, terms.Count)))
			{
				int mc = basinMinCell[id];
				plan.Endorheics.Add(new EndorheicTerminal
				{
					Terminal = (mc / n, mc % n),
					DrainageAreaPx = inflow,
					BasinDepthM = basinDepthM[id],
					BasinAreaPx = basinAreaPx[id]
				});
			}
		}

		return plan;
	}
}