beezm
148602b4c5
LOWGROUND is the gate verdict and the locked default ('short' stays available
behind the dial for the record). Every river polyline (upland stem + lowland
reach) is now decimated with RDP(4.0) and rounded with 4 Chaikin passes —
numerically mirroring MapGenerator.SmoothPath, since this pass is Godot-free —
before the bed profile is built, so the carved centreline carries none of the
8-connected Dijkstra 45° kinks (task-22 §4 caveat a). The monotone-descent and
sea-clamp constraints are applied to the smoothed course, unchanged.
Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
411 lines
16 KiB
C#
411 lines
16 KiB
C#
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;
|
||
}
|
||
|
||
// ---- 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[] 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 —
|
||
// then SMOOTHED (task 23) so the carved centreline carries no routing kinks.
|
||
var pts = new List<(float x, float y)>(upland);
|
||
if (lowlandRoute != null && lowlandRoute.Count > 1)
|
||
pts.AddRange(lowlandRoute.GetRange(1, lowlandRoute.Count - 1));
|
||
pts = SmoothCourse(pts);
|
||
|
||
// 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;
|
||
}
|
||
}
|