islaApocalypse/Tools/Scripts/RiverCarvePass.cs

using System;
using System.Collections.Generic;

/// <summary>
/// River bed carving — C0b part 2a (terrain-water task 22). The first river pass
/// that MODIFIES terrain: executes the frozen task-21b plan by carving channel
/// beds for the promoted rivers. NO WATER — part 2b puts water into these beds
/// once the routing-style gate picks SHORT or LOWGROUND.
///
/// Standalone numeric (D-035 family). Runs the task-21 DrainageAnalysis in-
/// pipeline (deterministic: same seed → same eroded surface → same plan), then:
///
///   1. LOWLAND ROUTING (the A/B): each routed giant gets a route from its
///      pooling terminal to the nearest OCEAN cell by deterministic Dijkstra.
///      SHORT     — cost ≈ distance, uphill penalised: heads direct, avoids walls.
///      LOWGROUND — cost ≈ elevation above sea: follows the lowest available
///                  ground and wanders like a real lowland river.
///   2. BED CARVING: every promoted course (trunks, routed giants + their lowland
///      reaches, lake-enders, tributaries) is stamped as a parabolic channel with
///      a smoothstep shoulder — a bed for water to sit in, not a canyon. Width and
///      depth grow downstream with drainage. The bed elevation is made MONOTONE
///      NON-INCREASING toward the outlet (water must flow), and is clamped to
///      sea + margin everywhere — the flood-guard discipline erosion established:
///      below-sea cells are read-only, no carve may create inland below-sea cells,
///      so the rendered coastline cannot move. The bed meets the ocean AT the
///      coast, where the terrain itself descends through sea level.
///   3. Crater: nothing inside the protected core is modified (erosion's rule).
///
/// Output-only: the caller applies this to the RENDER map after classify-side
/// data (biomes, WBID) is already computed — the oracle stays byte-identical.
/// </summary>
public static class RiverCarvePass
{
	public const float M_PER_UNIT = 251f;
	public const byte STYLE_SHORT = 0;
	public const byte STYLE_LOWGROUND = 1;

	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 };

	// Routing cost constants. SHORT pays lightly for climbing (8 per metre of rise,
	// so a 10 m wall costs like an 80 px detour — walls are avoided, direction is
	// kept). LOWGROUND pays for BEING high (1 per metre of elevation per px) plus
	// heavily for climbing, so the cheapest corridor is the lowest ground even when
	// that wanders.
	private const float SHORT_UPHILL_PER_M = 8f;
	private const float LOWGROUND_ELEV_PER_M = 1f;
	private const float LOWGROUND_BASE = 0.05f;
	private const float LOWGROUND_UPHILL_PER_M = 50f;

	// Bed geometry: sizes grow downstream from head to mouth, scaled by
	// sqrt(drainage / 1e6) so a 2.3M px giant carves roughly 2.3× deeper/wider at
	// the mouth than a 0.4M px trunk. Kept channel-scale, per the erosion
	// detailing philosophy.
	private const float DEPTH_HEAD_M = 1.0f;
	private const float DEPTH_MOUTH_M = 6.0f;     // × sizeFactor × DepthScale
	private const float HALFWIDTH_HEAD_PX = 2.0f;
	private const float HALFWIDTH_MOUTH_PX = 12.0f; // × sizeFactor × WidthScale
	private const float MIN_BED_SLOPE = 0.002f;   // m per px of enforced descent

	public class Params
	{
		public byte RoutingStyle = STYLE_LOWGROUND;
		public float WidthScale = 1.0f;
		public float DepthScale = 1.0f;
		public float SeaMarginM = 0.2f;           // bed floor above sea, everywhere
		public DrainageAnalysis.Params PlanParams = new();
	}

	public class RiverStat
	{
		public string Name;
		public string Kind;
		public bool SouthernCandidate;
		public long DrainagePx;
		public int CourseLenPx;
		public int RouteLenPx;                    // lowland reach only (routed giants)
		public float RouteStraightPx;
		public float WanderRatio;                 // routeLen / straight-line
		public bool ReachedOcean;
		public float MaxCutM;
		public double VolumeM3;
	}

	public class Stats
	{
		public List<RiverStat> Rivers = new();
		public long CarvedCells;
		public double CarvedVolumeM3;
		public float MaxCutM;
		public double AnalysisSeconds, RoutingSeconds, CarveSeconds;
	}

	public static Stats Apply(float[,] height, int mapSize, bool[] isOcean,
		bool[] isClassifyWater, float southX, float southY,
		float[,] seaMap, float seaFlat,
		float craterCx, float craterCy, float craterCoreRadius,
		Func<double> secondsNow, Params p)
	{
		int n = mapSize;
		var stats = new Stats();
		float SeaAt(int x, int y) => seaMap != null ? seaMap[x, y] : seaFlat;
		float coreSq = craterCoreRadius * craterCoreRadius;

		// --- the frozen plan, recomputed deterministically in-pipeline ---
		double t0 = secondsNow();
		var plan = DrainageAnalysis.Run(height, mapSize, isOcean, isClassifyWater, southX, southY, p.PlanParams);
		stats.AnalysisSeconds = secondsNow() - t0;

		// --- lowland routing for the routed giants (the A/B) ---
		t0 = secondsNow();
		var giantRoutes = new List<List<(float x, float y)>>();
		foreach (var g in plan.Giants)
		{
			if (g.Kind != "routed") { giantRoutes.Add(null); continue; }
			giantRoutes.Add(RouteToOcean(height, n, isOcean,
				(int)g.Terminal.x, (int)g.Terminal.y, p.RoutingStyle, SeaAt));
		}
		stats.RoutingSeconds = secondsNow() - t0;

		// --- carve ---
		t0 = secondsNow();
		int riverIdx = 0;
		foreach (var t in plan.Trunks)
		{
			riverIdx++;
			var course = new List<(float x, float y)>(t.Course);
			course.Reverse();                                  // head → mouth
			var rs = CarveRiver($"trunk{riverIdx}", "ocean-trunk", t.DrainageAreaPx,
				course, null, height, n, SeaAt, coreSq, craterCx, craterCy, p, stats);
			rs.ReachedOcean = true;                            // outlet is on the coast by construction
			foreach (var trib in t.Tributaries)
				CarveTributary(trib, height, n, SeaAt, coreSq, craterCx, craterCy, p, stats);
		}
		int gi = 0;
		foreach (var g in plan.Giants)
		{
			var route = giantRoutes[gi]; gi++;
			var course = new List<(float x, float y)>(g.Course);
			course.Reverse();                                  // head → terminal
			var rs = CarveRiver($"giant{gi}", g.Kind, g.DrainageAreaPx,
				course, route, height, n, SeaAt, coreSq, craterCx, craterCy, p, stats);
			rs.SouthernCandidate = g.SouthernCandidate;
			rs.ReachedOcean = g.Kind != "routed" || (route != null && route.Count > 0);
			foreach (var trib in g.Tributaries)
				CarveTributary(trib, height, n, SeaAt, coreSq, craterCx, craterCy, p, stats);
		}
		stats.CarveSeconds = secondsNow() - t0;
		return stats;
	}

	/// <summary>
	/// Deterministic Dijkstra from the start cell to the nearest ocean cell under
	/// the selected style's cost model. Returns the path start → ocean (1-px steps),
	/// or an empty list if no path exists (reported upstream, never asserted away).
	/// </summary>
	private static List<(float x, float y)> RouteToOcean(float[,] height, int n,
		bool[] isOcean, int sx, int sy, byte style, Func<int, int, float> seaAt)
	{
		int total = n * n;
		var gcost = new float[total];
		var parent = new int[total];
		var closed = new bool[total];
		Array.Fill(gcost, float.MaxValue);
		Array.Fill(parent, -1);

		float ElevM(int x, int y) => MathF.Max(0f, (height[x, y] - seaAt(x, y)) * M_PER_UNIT);

		var pq = new PriorityQueue<int, (float c, int i)>();
		int start = sx * n + sy;
		gcost[start] = 0f;
		pq.Enqueue(start, (0f, start));
		int goal = -1;

		while (pq.Count > 0)
		{
			int c = pq.Dequeue();
			if (closed[c]) continue;
			closed[c] = true;
			if (isOcean[c]) { goal = c; break; }
			int cx = c / n, cy = c % n;
			float hc = height[cx, cy];
			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 (closed[ni]) continue;
				float dhM = MathF.Max(0f, (height[nx, ny] - hc) * M_PER_UNIT);
				float step = style == STYLE_SHORT
					? DIST[k] + dhM * SHORT_UPHILL_PER_M
					: DIST[k] * (LOWGROUND_BASE + ElevM(nx, ny) * LOWGROUND_ELEV_PER_M)
						+ dhM * LOWGROUND_UPHILL_PER_M;
				float nc = gcost[c] + step;
				if (nc < gcost[ni])
				{
					gcost[ni] = nc;
					parent[ni] = c;
					pq.Enqueue(ni, (nc, ni));
				}
			}
		}

		var path = new List<(float x, float y)>();
		if (goal >= 0)
		{
			for (int c = goal; c >= 0; c = parent[c])
				path.Add((c / n, c % n));
			path.Reverse();
		}
		return path;
	}

	private static void CarveTributary(DrainageAnalysis.Stream trib, float[,] height, int n,
		Func<int, int, float> seaAt, float coreSq, float craterCx, float craterCy,
		Params p, Stats stats)
	{
		var course = new List<(float x, float y)>(trib.Course);
		course.Reverse();                                      // head → confluence
		CarveRiver(null, "tributary", trib.DrainageAreaPx, course, null,
			height, n, seaAt, coreSq, craterCx, craterCy, p, stats);
	}

	/// <summary>
	/// Carves one river: densify the course, build a monotone-descending clamped
	/// bed profile, stamp the channel. Returns the per-river stat (also appended
	/// to stats.Rivers unless name is null — tributaries fold into the totals).
	/// </summary>
	private static RiverStat CarveRiver(string name, string kind, long drainagePx,
		List<(float x, float y)> upland, List<(float x, float y)> lowlandRoute,
		float[,] height, int n, Func<int, int, float> seaAt,
		float coreSq, float craterCx, float craterCy, Params p, Stats stats)
	{
		// Full head→mouth polyline: upland stem, then the lowland reach if any.
		var pts = new List<(float x, float y)>(upland);
		if (lowlandRoute != null && lowlandRoute.Count > 1)
			pts.AddRange(lowlandRoute.GetRange(1, lowlandRoute.Count - 1));

		// Densify to ~1-px samples (plan courses are decimated ×4).
		var dense = new List<(float x, float y)>();
		for (int i = 0; i + 1 < pts.Count; i++)
		{
			var a = pts[i]; var b = pts[i + 1];
			float segLen = MathF.Sqrt((b.x - a.x) * (b.x - a.x) + (b.y - a.y) * (b.y - a.y));
			int steps = Math.Max(1, (int)MathF.Ceiling(segLen));
			for (int s2 = 0; s2 < steps; s2++)
				dense.Add((a.x + (b.x - a.x) * s2 / steps, a.y + (b.y - a.y) * s2 / steps));
		}
		if (pts.Count > 0) dense.Add(pts[^1]);
		if (dense.Count < 2) return new RiverStat();

		float sizeFactor = MathF.Sqrt(drainagePx / 1_000_000f);
		var rs = new RiverStat
		{
			Name = name, Kind = kind, DrainagePx = drainagePx,
			CourseLenPx = dense.Count,
			RouteLenPx = lowlandRoute?.Count ?? 0
		};
		if (lowlandRoute != null && lowlandRoute.Count > 1)
		{
			var a = lowlandRoute[0]; var b = lowlandRoute[^1];
			rs.RouteStraightPx = MathF.Sqrt((b.x - a.x) * (b.x - a.x) + (b.y - a.y) * (b.y - a.y));
			// Wander = POLYLINE length over straight-line — cell count undercounts
			// diagonal steps and can read below 1, which is geometrically impossible.
			float polyLen = 0f;
			for (int i = 1; i < lowlandRoute.Count; i++)
			{
				float sdx = lowlandRoute[i].x - lowlandRoute[i - 1].x;
				float sdy = lowlandRoute[i].y - lowlandRoute[i - 1].y;
				polyLen += MathF.Sqrt(sdx * sdx + sdy * sdy);
			}
			rs.RouteLenPx = (int)polyLen;
			rs.WanderRatio = rs.RouteStraightPx > 1f ? polyLen / rs.RouteStraightPx : 1f;
		}

		// Bed profile: raw = terrain − depth(t); then monotone non-increasing
		// downstream; then clamped to sea + margin. The clamp can flatten the tail
		// near the mouth — allowed: non-increasing is what water needs, and the
		// flood guard is absolute.
		int m = dense.Count;
		var bed = new float[m];
		var depth = new float[m];
		var halfW = new float[m];
		for (int i = 0; i < m; i++)
		{
			float t = m > 1 ? (float)i / (m - 1) : 1f;
			depth[i] = (DEPTH_HEAD_M + (DEPTH_MOUTH_M * sizeFactor - DEPTH_HEAD_M) * t) * p.DepthScale;
			if (depth[i] < 0.5f) depth[i] = 0.5f;
			halfW[i] = (HALFWIDTH_HEAD_PX + (HALFWIDTH_MOUTH_PX * sizeFactor - HALFWIDTH_HEAD_PX) * t) * p.WidthScale;
			if (halfW[i] < 1.5f) halfW[i] = 1.5f;
			int cx = (int)dense[i].x, cy = (int)dense[i].y;
			bed[i] = height[cx, cy] - depth[i] / M_PER_UNIT;
		}
		for (int i = 1; i < m; i++)
		{
			float maxAllowed = bed[i - 1] - MIN_BED_SLOPE / M_PER_UNIT;
			if (bed[i] > maxAllowed) bed[i] = maxAllowed;
		}
		for (int i = 0; i < m; i++)
		{
			int cx = (int)dense[i].x, cy = (int)dense[i].y;
			float floor = seaAt(cx, cy) + p.SeaMarginM / M_PER_UNIT;
			if (bed[i] < floor) bed[i] = floor;
		}

		// Stamp: parabolic channel to the rim, smoothstep shoulder back to terrain.
		for (int i = 0; i < m; i++)
		{
			float hw = halfW[i];
			float outer = hw * 2f;
			int cx0 = (int)MathF.Floor(dense[i].x - outer), cx1 = (int)MathF.Ceiling(dense[i].x + outer);
			int cy0 = (int)MathF.Floor(dense[i].y - outer), cy1 = (int)MathF.Ceiling(dense[i].y + outer);
			float rimH = bed[i] + depth[i] / M_PER_UNIT;
			for (int x = cx0; x <= cx1; x++)
			{
				if (x < 0 || x >= n) continue;
				for (int y = cy0; y <= cy1; y++)
				{
					if (y < 0 || y >= n) continue;
					float rx = x - dense[i].x, ry = y - dense[i].y;
					float r = MathF.Sqrt(rx * rx + ry * ry);
					if (r > outer) continue;
					float ddx = x - craterCx, ddy = y - craterCy;
					if (ddx * ddx + ddy * ddy < coreSq) continue;   // crater core protected
					float sea = seaAt(x, y);
					float old = height[x, y];
					if (old < sea) continue;                        // below-sea cells read-only
					float target;
					if (r <= hw)
					{
						float f = r / hw;
						target = bed[i] + (depth[i] / M_PER_UNIT) * f * f;
					}
					else
					{
						float f = (r - hw) / hw;                    // 0..1 across the shoulder
						f = f * f * (3f - 2f * f);                  // smoothstep
						target = rimH + (old - rimH) * f;
					}
					float floor = sea + p.SeaMarginM / M_PER_UNIT;
					if (target < floor) target = floor;
					if (target < old)
					{
						float cutM = (old - target) * M_PER_UNIT;
						height[x, y] = target;
						stats.CarvedCells++;
						stats.CarvedVolumeM3 += cutM;
						if (cutM > stats.MaxCutM) stats.MaxCutM = cutM;
						if (cutM > rs.MaxCutM) rs.MaxCutM = cutM;
						rs.VolumeM3 += cutM;
					}
				}
			}
		}

		if (name != null) stats.Rivers.Add(rs);
		return rs;
	}
}