islaApocalypse-v2/Core/Scripts/DrainageAnalysis.cs

using System;
using System.Collections.Generic;

namespace IslaApocalypse.Core
{
	/// <summary>
	/// ⭐⭐ DRAINAGE ANALYSIS — THE FAITHFUL PORT (chat2/12). From the reference's
	/// <c>Tools/Scripts/DrainageAnalysis.cs</c> at tag <c>pre-rewrite-reference</c> (<c>ab78883</c>), verbatim
	/// in arithmetic and order (D-050). PURE ANALYSIS: it reads the ERODED render heightmap and produces a
	/// river PLAN — it changes zero terrain and adds zero water.
	///
	/// ═══ WHAT THE PORT CHANGES (and nothing else) ═══
	///
	///   • Namespace + location: <c>IslaApocalypse.Core</c> — engine-free, a C++ candidate.
	///   • The yardstick: metres via <see cref="WorldScale.MetresFromRaw"/> (the same <c>× 251f</c>; no literal).
	///   • <see cref="Plan.FullFilled"/> and <see cref="Plan.Filled"/> are EXPOSED (the reference kept the
	///     routing surfaces local) so the caller can prove the routing-fill invariants on them.
	///
	/// ═══ ⚠⚠ THE D8 LANDMINE ═══
	///
	/// D8 flow direction is CORRECT FOR ANALYSIS and is used here for exactly that. It was REVERTED as a
	/// CARVING technique (the prototype's task 10 — straight, grid-aligned grooves; → `Design - Terrain -
	/// D8 Incision Revert`). Use D8 to COMPUTE, never to CARVE. Nothing in this file writes terrain.
	///
	/// ═══ ⚠ ENDORHEIC BASINS ARE EXPECTED, FIRST-CLASS OUTPUT ═══
	///
	/// This terrain's biggest drainages pool inland: erosion delivers the upland network only and cannot
	/// cross the flats. A screen full of endorheic basins is the CORRECT result, not a bug.
	///
	/// ═══ THE REFERENCE'S CLASS DOC (verbatim) ═══
	///
	/// 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
	{
		// 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 const sbyte D_NONE = -1, D_SEA = -2;

		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 GiantCount = 3;            // 21b: top endorheic giants promoted
			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();
		}

		/// <summary>
		/// A promoted endorheic giant (task 21b): one of the island's biggest drainage
		/// systems, which pools inland because erosion could not cross the flats.
		/// Kind "routed" carries a PROVISIONAL route across the flats to the ocean —
		/// the path part 2 would carve, drawn for the gate, not water. Kind
		/// "lake-ender" keeps its lake/lagoon terminal (real geography, developer's
		/// call). Terminal is where the MAIN STEM actually pools (its sub-minimum),
		/// which on a flat basin floor is more truthful than the basin's deepest cell.
		/// ⚠ chat2/12 computes the provisional route as the reference did (it is analysis) but does NOT
		/// promote or draw it — lowland routing is a later task.
		/// </summary>
		public class Giant : Stream
		{
			public (float x, float y) Terminal;
			public (float x, float y) Spill;                 // where the basin overtops
			public float BasinDepthM;
			public long BasinAreaPx;
			public string Kind = "routed";                   // "routed" | "lake-ender"
			public bool SouthernCandidate;
			public bool TerminalInClassifyWater;
			public List<(float x, float y)> ProvisionalRoute; // null for lake-enders
			public bool RouteReachedOcean;
			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 List<Giant> Giants = new();               // 21b: the promoted giants
			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;

			// ---- exposed by the port (the reference kept these local) ----
			/// <summary>D8 direction per cell (row-major x·n+y): 0..7, <see cref="D_SEA"/>, <see cref="D_NONE"/>. Analysis only.</summary>
			public sbyte[] Dir;
			/// <summary>Flow accumulation per cell (row-major); 0 on ocean.</summary>
			public int[] Acc;
			/// <summary>The routing surface after terminal basins reverted (row-major).</summary>
			public float[] Filled;
			/// <summary>The FULL priority-flood fill, before terminal reversion (row-major) — every cell drains to the border on it.</summary>
			public float[] FullFilled;
			/// <summary>Terminal-basin id per cell (row-major), 0 = none.</summary>
			public int[] BasinId;
			/// <summary>Per terminal basin id: total inflow (cells whose flow ends there).</summary>
			public long[] BasinInflow;
		}

		/// <param name="isOcean">Row-major mask of THE OCEAN body — the only water that counts as "the sea" for
		/// sea-reaching trunks (in v2: <c>RegionLabeling.OceanMask</c>, the classify field's border-connected water).
		/// Below-sea cells that are NOT ocean (enclosed lagoons, below-datum lake beds, island-fringe waters) are
		/// ordinary terrain to the router: as depressions they either qualify as terminal basins or fill and spill
		/// onward to the true sea.</param>
		/// <param name="isClassifyWater">Row-major mask of ANY classify water: a giant whose main stem pools inside
		/// classify water is a natural lake-ender; one pooling on dry ground is a route-to-sea candidate.</param>
		/// <param name="southX">Southernmost-town position (or -1 for none): the giant whose terminal lies closest is
		/// flagged the SOUTHERN CANDIDATE and always routed provisionally, per the 21b design — shown, not forced.</param>
		public static Plan Run(float[,] height, int mapSize, bool[] isOcean,
			bool[] isClassifyWater, float southX, float southY, 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];

			float[] plan_fullFilled = null; // set inside step 2, used by 21b routing

			// --- 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 = WorldScale.MetresFromRaw(filled[c] - original[c]);
						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);
				}

				// 21b: the FULL fill (before terminal reversion) is the provisional-
				// routing surface — on it, every basin overtops at its spill and drains
				// to the border, which is exactly "where the water would continue".
				plan_fullFilled = (float[])filled.Clone();

				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).
			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 = WorldScale.MetresFromRaw(original[stem[i + w]] - original[stem[i]]);
					if (rise / w >= p.ExitGradeMin)
					{
						t.ExitFound = true;
						t.MountainExit = (stem[i] / n, stem[i] % n);
						t.MountainExitElevM = WorldScale.MetresFromRaw(original[stem[i]]);
						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]
					});
				}
			}

			// --- 5d. The promoted GIANTS (21b): mixed set, provisional routes ---
			// Top GiantCount terminal basins by TOTAL inflow. Their upland stems are the
			// island's real big rivers; whether each continues to the sea is the gate's
			// decision, previewed here.
			{
				var giantsRanked = 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;
					if (basinInflow[id] >= p.EndorheicMinInflowPx) giantsRanked.Add((id, basinInflow[id]));
				}
				giantsRanked.Sort((a, b) => b.inflow.CompareTo(a.inflow));

				// Does a terminal basin HOLD classify water? The lake-ender test must look
				// at the whole pool, not the stem's single pooling cell — a stem can pool on
				// dry ground a few hundred px short of its lagoon and still be a lagoon river.
				bool[] basinHasLake = new bool[basinMinCell.Count];
				for (int i = 0; i < total; i++)
					if (basinId[i] != 0 && isClassifyWater[i] && !isOcean[i])
						basinHasLake[basinId[i]] = true;

				// The main stem's ENTRY into the basin: the highest-accumulation cell
				// whose flow terminates in this basin. On a flat basin floor the deepest
				// cell sees only local trickles (the task-21 lesson), so the stem is
				// anchored on the strongest feeder instead.
				var bestEntry = new Dictionary<int, int>();
				for (int i = 0; i < total; i++)
				{
					if (dest[i] <= 0) continue;
					if (!bestEntry.TryGetValue(dest[i], out int cur) || acc[i] > acc[cur])
						bestEntry[dest[i]] = i;
				}

				// The giant whose pooling point sits closest to the southernmost town is
				// the SOUTHERN CANDIDATE — always routed provisionally (shown, not forced).
				int southernPick = -1;
				if (southX >= 0f)
				{
					float bestD = float.MaxValue;
					foreach (var (id, _) in giantsRanked.GetRange(0, Math.Min(p.GiantCount, giantsRanked.Count)))
					{
						int mc = basinMinCell[id];
						float ddx = mc / n - southX, ddy = mc % n - southY;
						float d2 = ddx * ddx + ddy * ddy;
						if (d2 < bestD) { bestD = d2; southernPick = id; }
					}
				}

				foreach (var (id, inflow) in giantsRanked.GetRange(0, Math.Min(p.GiantCount, giantsRanked.Count)))
				{
					var g = new Giant { DrainageAreaPx = inflow, BasinDepthM = basinDepthM[id], BasinAreaPx = basinAreaPx[id] };
					if (!bestEntry.TryGetValue(id, out int entry)) entry = basinMinCell[id];

					// Downstream from the strongest feeder to where it actually pools…
					int t2 = entry;
					var down = new List<int> { t2 };
					while (dir[t2] >= 0) { t2 = Target(t2); down.Add(t2); }
					g.Terminal = (t2 / n, t2 % n);
					// …then the full main stem, traced upstream from that pooling point.
					var stem = TraceStem(t2, p.StemMinAccPx);
					g.Course = Decimate(stem);
					g.Head = (stem[^1] / n, stem[^1] % n);
					g.TerminalInClassifyWater = isClassifyWater[t2];

					for (int i = 0; i + p.ExitWindowPx < stem.Count; i++)
					{
						float rise = WorldScale.MetresFromRaw(original[stem[i + p.ExitWindowPx]] - original[stem[i]]);
						if (rise / p.ExitWindowPx >= p.ExitGradeMin)
						{
							g.ExitFound = true;
							g.MountainExit = (stem[i] / n, stem[i] % n);
							g.MountainExitElevM = WorldScale.MetresFromRaw(original[stem[i]]);
							break;
						}
					}

					// Lean tributaries on the giant's stem, same junction rule as trunks.
					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));
					var takenT = new List<int>();
					foreach (var (cell, _) in cands)
					{
						if (takenT.Count >= p.TributaryMaxPerTrunk) break;
						int cx2 = cell / n, cy2 = cell % n;
						bool dup = false;
						foreach (int tc in takenT)
						{
							float ddx = cx2 - tc / n, ddy = cy2 - tc % n;
							if (ddx * ddx + ddy * ddy < 30f * 30f) { dup = true; break; }
						}
						if (!dup) takenT.Add(cell);
					}
					foreach (int cell in takenT)
					{
						var trib = new Stream { DrainageAreaPx = acc[cell] };
						var ts = TraceStem(cell, Math.Max(p.StemMinAccPx, (int)(acc[cell] / 20)));
						trib.Course = Decimate(ts);
						trib.Head = (ts[^1] / n, ts[^1] % n);
						g.Tributaries.Add(trib);
					}

					// Kind: the terminal BASIN holds a classify lake → natural lake-ender;
					// dry pan → route to sea; the southern candidate is always routed.
					g.SouthernCandidate = id == southernPick;
					g.TerminalInClassifyWater = g.TerminalInClassifyWater || basinHasLake[id];
					g.Kind = (basinHasLake[id] && !g.SouthernCandidate) ? "lake-ender" : "routed";

					// PROVISIONAL route (routed giants): walk steepest descent on the FULL
					// fill from the pooling point — the basin overtops at its spill and
					// the walk continues along the terrain's own drainage to the ocean.
					// DRAWN, not carved; part 2 carves along a route like this one.
					// (chat2/12: computed as the reference did; not drawn, not promoted — routing is later.)
					if (g.Kind == "routed")
					{
						var route = new List<int>();
						int c = t2;
						bool spillRecorded = false;
						for (int guard = 0; guard < 4 * n; guard++)
						{
							route.Add(c);
							if (isOcean[c]) { g.RouteReachedOcean = true; break; }
							if (!spillRecorded && basinId[c] != id)
							{ g.Spill = (c / n, c % n); spillRecorded = true; }
							int cx = c / n, cy = c % n;
							float best = float.MaxValue; int bestN = -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 (plan_fullFilled[ni] < best) { best = plan_fullFilled[ni]; bestN = ni; }
							}
							if (bestN < 0 || plan_fullFilled[bestN] >= plan_fullFilled[c]) break; // stuck (report via flag)
							c = bestN;
						}
						g.ProvisionalRoute = Decimate(route);
					}
					plan.Giants.Add(g);
				}
			}

			plan.Dir = dir; plan.Acc = acc; plan.Filled = filled; plan.FullFilled = plan_fullFilled;
			plan.BasinId = basinId; plan.BasinInflow = basinInflow;
			return plan;
		}
	}
}