islaApocalypse/Tools/Scripts/RiverCarvePass.cs
beezm c1f62f6899 feat: tributaries wet full length; water fills the channel; banks flare as shores (terrain-water task 25)
The gate's two river notes.

FILL THE TRIBUTARIES. RiverTribWaterMinFlow defaults to 0, so every promoted
tributary carries water its whole length (800 of 1574 reaches are now on
tributaries, 65,364 of 268,863 wet px). The graceful tip survives: the fade zone
is now always TaperPx long starting wherever water begins, so with a threshold it
spans (WetFrom-Taper -> WetFrom) and with full watering it spans the first
TaperPx from the HEAD. Fill is about LENGTH, taper is about the END.

FILL THE CHANNEL, NOT THE FLOOR. The water surface is now RiverFillFraction (80%)
of the LOCAL bed depth rather than a fixed height, floored by RiverWaterDepthM as
an absolute minimum. A fixed height is wrong at both ends — a trickle in a 13 m
mouth and over the rim at a 1.5 m head — where a fraction is right at every
scale and cannot spill onto the plain.

BANKS AS SHORES. The shoulder flares RiverBankFlare (3.2, was a hard 2x)
half-widths with a smootherstep instead of smoothstep, so the rim leaves nearly
tangent to the water plane. That alone cut too much where a course crosses a
ridge (cells deeper than 20 m went 339 -> 1917, max 26 -> 37 m), so the SHOULDER
may not lower a cell more than RiverBankMaxCutM (3 m) — measured against the
PRE-PASS surface, because a per-write cap let overlapping stamps each take
another 3 m and changed the carved volume by 1 m3 out of 1.29 M. Gentle ground
still flares into a shore; a ridge crossing keeps steep walls, which is what a
gorge looks like. Net result is better than BOTH predecessors: max cut 16.8 m
(task 24: 26.0), zero cells past 20 m (task 24: 339), carve p95 5.26 m.

Guards unchanged: BIOME oracle md5-identical to the task-22 baseline, 0
newly-below-sea, 0 below-sea cells modified, lowest carved cell exactly
sea+margin, island top 457.65 m exact, crater core excluded.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
2026-08-11 21:11:41 -04:00

659 lines
28 KiB
C#
Raw Blame History

This file contains ambiguous Unicode characters

This file contains Unicode characters that might be confused with other characters. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.

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
// Task 24: a tributary reach is WET where its along-course flow exceeds this;
// upstream of that it TAPERS to dry over TribTaperPx rather than ending in a
// wall of water. 0 = water tributaries end to end.
public int TribWaterMinFlowPx = 0; // task 25: 0 = tributaries wet full length
public int TribTaperPx = 120;
// Task 25 (the gate's "trickle at the bottom of a ditch" note): the reach's
// water surface is a FRACTION of the local bed depth rather than a fixed
// height — the channel reads FILLED at every scale and cannot overfill onto
// the plain, which a fixed depth does at the shallow heads. WaterDepthM
// survives as the absolute minimum so tiny channels still hold water.
public float FillFraction = 0.80f;
// Bank shoulders flare this many half-widths beyond the channel (was a hard
// 2×) with a gentler-than-smoothstep curve, so water meets land as a shore.
public float BankFlare = 3.2f;
// ...but the SHOULDER may never lower a cell by more than this. A gentle flare
// across a ridge would otherwise cut a big notch: widening the flare alone took
// cells deeper than 20 m from 339 to 1917 (measured). With the cap, gentle
// ground still flares into a shore — which is where the "water in a groove"
// complaint lives — while a ridge crossing keeps steep walls, which is what a
// gorge actually looks like. The channel bed itself is not affected by this.
public float BankMaxCutM = 3.0f;
// Lake-enders route to the nearest water body of at least this size — the
// nearest wet PIXEL was a puddle (task-23 gate finding).
public int LakeMinTargetPx = 20_000;
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;
}
/// <summary>The carved geometry the water stage consumes (main rivers only).</summary>
public class CarvedRiver
{
public string Name, Kind;
public bool Southern;
public long DrainagePx;
public List<(float x, float y)> Dense; // head → mouth, ~1-px samples
public float[] Bed; // raw units, monotone non-increasing
public float[] HalfW; // px
public float[] DepthM; // local bed depth, metres (task 25 fill)
public bool ReachedWaterTerminal; // lake-enders: extension reached classify water
// Task 24: along-course flow (px of drainage) per sample, and the first index
// that carries water. Between WetFrom-TaperPx and WetFrom the water tapers
// (narrowing and shallowing to the bed) so a stream head fades out.
public float[] Flow;
public int WetFrom;
public int TaperPx;
}
public class Stats
{
public List<RiverStat> Rivers = new();
public List<CarvedRiver> Carved = new(); // for the stepped-water stage (task 23)
internal float[] PrePass; // cumulative-cut baseline
public long CarvedCells;
public double CarvedVolumeM3;
public float MaxCutM; // CUMULATIVE vs pre-pass heights (task-22 nit 1 fixed)
public double AnalysisSeconds, RoutingSeconds, CarveSeconds;
}
/// <param name="isSignificantWater">Row-major mask of classify water belonging to
/// bodies of at least LakeMinTargetPx cells (task 24). Lake-enders route to THIS,
/// not to any wet pixel: the task-23 build routed one into a puddle a few hundred
/// px short of the obvious lagoon, because "nearest classify water" is satisfied
/// by a 3-cell pond.</param>
public static Stats Apply(float[,] height, int mapSize, bool[] isOcean,
bool[] isClassifyWater, bool[] isSignificantWater, 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();
// Pre-pass snapshot: max-cut is measured CUMULATIVELY against the heights
// this pass found, not per-write — overlapping stamps re-cut a cell and the
// per-write number understated the true deepest cut ~4× (task-22 nit 1).
float[] pre = new float[n * n];
for (int x = 0; x < n; x++)
for (int y = 0; y < n; y++)
pre[x * n + y] = height[x, y];
stats.PrePass = pre;
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
// Lake-enders (task 23): the stem pools on dry ground short of its lake
// BECAUSE its pooling point is a local minimum — a blind descent walk
// dead-ends there immediately (measured: 0 steps). Route to the nearest
// classify-water cell with the same lowground Dijkstra the routed giants
// use, so the bed (and then the water) actually joins the lake.
bool reachedLake = false;
if (g.Kind == "lake-ender")
{
// Target SIGNIFICANT water (task 24). Fall back to any classify water
// only if no significant body is reachable, so a seed whose lake-ender
// genuinely has only small ponds still connects rather than dead-ending.
var ext = RouteToOcean(height, n, isSignificantWater,
(int)g.Terminal.x, (int)g.Terminal.y, STYLE_LOWGROUND, SeaAt);
if (ext.Count == 0)
ext = RouteToOcean(height, n, isClassifyWater,
(int)g.Terminal.x, (int)g.Terminal.y, STYLE_LOWGROUND, SeaAt);
if (ext.Count > 0) { route = ext; reachedLake = true; }
}
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) : reachedLake;
if (stats.Carved.Count > 0) stats.Carved[^1].ReachedWaterTerminal = reachedLake;
foreach (var trib in g.Tributaries)
CarveTributary(trib, height, n, SeaAt, coreSq, craterCx, craterCy, p, stats);
}
stats.CarveSeconds = secondsNow() - t0;
return stats;
}
// ---- Route smoothing (task 23): the road pass's AAA pipeline, numerically ----
// RDP(4.0) decimation + 4 Chaikin corner-cutting passes, endpoints pinned —
// the same constants and structure as MapGenerator.SmoothPath, mirrored here
// because this pass is Godot-free. Kills the 8-connected Dijkstra 45° kinks;
// the bed then carves along the smoothed centreline.
private static List<(float x, float y)> SmoothCourse(List<(float x, float y)> raw)
{
if (raw.Count < 3) return raw;
var dec = Rdp(raw, 0, raw.Count - 1, 4.0f);
if (dec.Count < 3) return raw;
var sm = dec;
for (int pass = 0; pass < 4; pass++)
{
var nxt = new List<(float x, float y)>(sm.Count * 2) { sm[0] };
for (int i = 0; i + 1 < sm.Count; i++)
{
var a = sm[i]; var b = sm[i + 1];
nxt.Add((a.x * 0.75f + b.x * 0.25f, a.y * 0.75f + b.y * 0.25f));
nxt.Add((a.x * 0.25f + b.x * 0.75f, a.y * 0.25f + b.y * 0.75f));
}
nxt.Add(sm[^1]);
sm = nxt;
}
return sm;
}
private static List<(float x, float y)> Rdp(List<(float x, float y)> pts, int i0, int i1, float tol)
{
if (i1 - i0 <= 1) return new List<(float x, float y)> { pts[i0], pts[i1] };
var a = pts[i0]; var b = pts[i1];
float abx = b.x - a.x, aby = b.y - a.y;
float abLen = MathF.Sqrt(abx * abx + aby * aby);
float maxD = 0f; int maxI = i0;
for (int i = i0 + 1; i < i1; i++)
{
float d = abLen < 1e-6f
? MathF.Sqrt((pts[i].x - a.x) * (pts[i].x - a.x) + (pts[i].y - a.y) * (pts[i].y - a.y))
: MathF.Abs(abx * (a.y - pts[i].y) - (a.x - pts[i].x) * aby) / abLen;
if (d > maxD) { maxD = d; maxI = i; }
}
if (maxD <= tol) return new List<(float x, float y)> { pts[i0], pts[i1] };
var left = Rdp(pts, i0, maxI, tol);
var right = Rdp(pts, maxI, i1, tol);
left.RemoveAt(left.Count - 1);
left.AddRange(right);
return left;
}
/// <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[] targets, 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 (targets[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
// Task 24: tributaries are registered as carved rivers (name "trib") so the
// water stage sees them — they carved but stayed dry in task 23. They are NOT
// added to stats.Rivers, so the per-river console table stays the 6 mains.
CarveRiver("trib", "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.
// ONLY the lowland reach is smoothed: the Dijkstra 45° kinks live there, on
// near-flat ground where a rounded corner costs nothing. The upland stems
// already thread the carved valley FLOORS — smoothing them off-line cut
// valley walls (measured: max cut 14.6 → 27.3 m before this was split).
var pts = new List<(float x, float y)>(upland);
if (lowlandRoute != null && lowlandRoute.Count > 1)
{
var smoothedRoute = SmoothCourse(lowlandRoute);
pts.AddRange(smoothedRoute.GetRange(1, smoothedRoute.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 * p.BankFlare;
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
{
// Shoulder: 0 at the rim → 1 at natural ground, over a flare of
// (BankFlare-1) half-widths. Smootherstep (6t⁵15t⁴+10t³) leaves
// the rim nearly tangent to the water plane, so the bank reads as
// a shore rather than the wall a plain smoothstep left.
float f = (r - hw) / MathF.Max(1e-3f, hw * (p.BankFlare - 1f));
if (f > 1f) f = 1f;
f = f * f * f * (f * (6f * f - 15f) + 10f);
target = rimH + (old - rimH) * f;
// Clamp against the PRE-PASS surface, not the current height:
// overlapping stamps re-visit a cell, so a per-write cap lets each
// pass take another BankMaxCut (measured: capping against `old`
// changed the carved volume by 1 m³ out of 1.29 M — i.e. nothing).
float shoulderFloor = stats.PrePass[x * n + y] - p.BankMaxCutM / M_PER_UNIT;
if (target < shoulderFloor) target = shoulderFloor;
}
float floor = sea + p.SeaMarginM / M_PER_UNIT;
if (target < floor) target = floor;
if (target < old)
{
float cutM = (old - target) * M_PER_UNIT;
// Cumulative depth vs the PRE-PASS surface (nit 1): the
// honest "how deep did we cut here in total" number.
float cumM = (stats.PrePass[x * n + y] - target) * M_PER_UNIT;
height[x, y] = target;
stats.CarvedCells++;
stats.CarvedVolumeM3 += cutM;
if (cumM > stats.MaxCutM) stats.MaxCutM = cumM;
if (cumM > rs.MaxCutM) rs.MaxCutM = cumM;
rs.VolumeM3 += cutM;
}
}
}
}
if (name != null)
{
if (kind != "tributary") stats.Rivers.Add(rs);
// Along-course flow: drainage grows from a headwater trickle to the full
// figure at the mouth. Quadratic in t so the substantial lower half
// dominates — a first-order stand-in for real accumulation, which the
// plan only records per-river.
var flow = new float[m];
for (int i = 0; i < m; i++)
{
float t = m > 1 ? (float)i / (m - 1) : 1f;
flow[i] = drainagePx * (0.05f + 0.95f * t * t);
}
int wetFrom = 0;
if (kind == "tributary" && p.TribWaterMinFlowPx > 0)
{
wetFrom = m; // dry unless the threshold is met
for (int i = 0; i < m; i++)
if (flow[i] >= p.TribWaterMinFlowPx) { wetFrom = i; break; }
}
stats.Carved.Add(new CarvedRiver
{
Name = name, Kind = kind, DrainagePx = drainagePx,
Dense = dense, Bed = bed, HalfW = halfW, DepthM = depth,
Flow = flow, WetFrom = wetFrom, TaperPx = p.TribTaperPx
});
}
return rs;
}
/// <summary>
/// The stepped-water builder (task 23, part 2b): segments each carved main
/// river into REACHES — flat water bodies stepping down the bed toward the
/// outlet — and stamps their ids into the WBID grid. Reuses the existing
/// levels-not-cells water model exactly: one body per reach, one flat level
/// each; the writer derives WSRF from body levels as it always has. The step
/// drops are the smoothing dial (smaller drop = more, finer steps); tilted
/// water is the deferred model B and is NOT built here.
///
/// Emission is plain data (no engine types): the caller turns reaches into
/// WBTB entries. Wet cells: inside the channel half-width, currently dry in
/// WBID, at/above sea (below-sea cells belong to the ocean/crater-seam rule),
/// bed below the reach level. Existing water bodies are never overwritten —
/// a river MEETS a lake or the sea, it does not repaint them.
/// </summary>
public class Reach
{
public ushort Id;
public string River;
public float Level; // raw units
public int PixelCount;
public double Cx, Cy; // centroid accumulators → mean
}
public static List<Reach> AddSteppedWater(float[,] height, int mapSize,
ushort[,] wbid, ushort firstId, List<CarvedRiver> rivers,
float[,] seaMap, float seaFlat, float craterCx, float craterCy,
float craterCoreRadius, float stepDropM, float waterDepthM, float fillFraction)
{
int n = mapSize;
float SeaAt(int x, int y) => seaMap != null ? seaMap[x, y] : seaFlat;
float coreSq = craterCoreRadius * craterCoreRadius;
var reaches = new List<Reach>();
ushort nextId = firstId;
foreach (var r in rivers)
{
int m = r.Dense.Count;
if (m < 2) continue;
// Task 24: start at the wet-from index (tributary threshold; 0 for mains),
// but back up by the taper length so the transition is a FADE, not a wall.
int i = Math.Max(0, r.WetFrom - r.TaperPx);
if (i >= m) continue; // entirely below threshold: dry
// The fade zone is always TaperPx long starting at i. With a threshold it
// spans (WetFrom-Taper → WetFrom); with tributaries fully watered
// (WetFrom = 0, task 25) it spans the first TaperPx from the HEAD, so a
// full-length tributary still fades in at its tip instead of starting as a
// wall of water.
int taperStart = i, taperEnd = Math.Min(m - 1, i + r.TaperPx);
float lastLevel = float.MaxValue;
while (i < m)
{
// Reach spans from i while the bed stays within stepDropM of the
// reach's starting bed; its flat level sits waterDepthM above that
// start (deepening toward the next step — the pool behind a riffle).
float startBed = r.Bed[i];
// Fill the channel: the surface sits at FillFraction of the LOCAL bed
// depth, never below the absolute minimum. A fixed height made deep
// mouths a trickle and overtopped shallow heads; a fraction is right at
// both ends and cannot spill onto the plain.
float localDepthM = r.DepthM != null ? r.DepthM[i] : waterDepthM;
float fillM = MathF.Max(waterDepthM, localDepthM * fillFraction);
if (fillM > localDepthM) fillM = localDepthM; // never above the rim
float level = startBed + fillM / M_PER_UNIT;
if (level >= lastLevel) // enforce strict descent
level = lastLevel - 0.01f / M_PER_UNIT;
int j = i;
while (j < m && r.Bed[j] > startBed - stepDropM / M_PER_UNIT) j++;
var reach = new Reach { Id = nextId, River = r.Name, Level = level };
for (int k2 = i; k2 < j; k2++)
{
// Taper: 0 at the dry end of the fade zone → 1 at full water. The
// water narrows AND its surface drops to the bed, so a stream head
// thins out and disappears instead of ending in a flat wall.
float taper = 1f;
if (taperEnd > taperStart && k2 < taperEnd)
taper = (float)(k2 - taperStart) / (taperEnd - taperStart);
if (taper <= 0.02f) continue;
float hw = r.HalfW[k2] * (0.25f + 0.75f * taper);
int x0 = (int)MathF.Floor(r.Dense[k2].x - hw), x1 = (int)MathF.Ceiling(r.Dense[k2].x + hw);
int y0 = (int)MathF.Floor(r.Dense[k2].y - hw), y1 = (int)MathF.Ceiling(r.Dense[k2].y + hw);
for (int x = x0; x <= x1; x++)
{
if (x < 0 || x >= n) continue;
for (int y = y0; y <= y1; y++)
{
if (y < 0 || y >= n) continue;
if (wbid[x, y] != 0) continue; // never repaint existing water
float rx = x - r.Dense[k2].x, ry = y - r.Dense[k2].y;
if (rx * rx + ry * ry > hw * hw) continue;
float ddx = x - craterCx, ddy = y - craterCy;
if (ddx * ddx + ddy * ddy < coreSq) continue;
float h = height[x, y];
float sea = SeaAt(x, y);
if (h < sea) continue; // ocean/seam territory
// The tapered surface sits between the bed and the reach
// level, so the wet strip shallows as it narrows.
float localLevel = taper >= 1f
? level
: r.Bed[k2] + (level - r.Bed[k2]) * taper;
if (h >= localLevel) continue; // bank above the water line
wbid[x, y] = nextId;
reach.PixelCount++;
reach.Cx += x; reach.Cy += y;
}
}
}
if (reach.PixelCount > 0)
{
reach.Cx /= reach.PixelCount; reach.Cy /= reach.PixelCount;
reaches.Add(reach);
nextId++;
lastLevel = level;
}
i = j;
}
}
return reaches;
}
}