diff --git a/Core/Scripts/BlueprintFormat.cs b/Core/Scripts/BlueprintFormat.cs index 99d4597..5aad3c7 100644 --- a/Core/Scripts/BlueprintFormat.cs +++ b/Core/Scripts/BlueprintFormat.cs @@ -33,6 +33,7 @@ namespace IslaApocalypse.Core public const uint TAG_WATER_BODY_TABLE = 0x42544257; // "WBTB" public const uint TAG_WATER_SURFACE = 0x46525357; // "WSRF" public const uint TAG_TERRAIN_CURVE = 0x56524354; // "TCRV" + public const uint TAG_TERRAIN_DETAIL = 0x4C544454; // "TDTL" // WSRF quantization: u16, 0 reserved as the no-water sentinel. A real level L // (raw blueprint height units) encodes as 1 + round(L × 32768), so a genuine diff --git a/Core/Scripts/BlueprintWriter.cs b/Core/Scripts/BlueprintWriter.cs index ed2de85..27fe859 100644 --- a/Core/Scripts/BlueprintWriter.cs +++ b/Core/Scripts/BlueprintWriter.cs @@ -33,6 +33,8 @@ namespace IslaApocalypse.Core WriteSection(writer, BlueprintFormat.TAG_PARAMS, w => WriteParams(w, bp)); if (bp.TerrainCurve != null) WriteSection(writer, BlueprintFormat.TAG_TERRAIN_CURVE, w => WriteTerrainCurve(w, bp.TerrainCurve)); + if (bp.TerrainDetail != null) + WriteSection(writer, BlueprintFormat.TAG_TERRAIN_DETAIL, w => WriteTerrainDetail(w, bp.TerrainDetail)); WriteSection(writer, BlueprintFormat.TAG_HEIGHTS, w => WriteHeights(w, bp)); WriteSection(writer, BlueprintFormat.TAG_BIOMES, w => WriteBiomes(w, bp)); @@ -151,6 +153,16 @@ namespace IslaApocalypse.Core } } + private static void WriteTerrainDetail(BinaryWriter writer, TerrainDetailInfo d) + { + writer.Write(d.Version); + writer.Write(d.ReliefAmpM); writer.Write(d.ReliefFreqIslands); + writer.Write(d.IncK); writer.Write(d.IncP); writer.Write(d.IncCapM); + writer.Write(d.SeaClampRaw); writer.Write(d.CraterExclFactor); + writer.Write(d.ShelfIncWeight); + writer.Write(d.ReliefSeedOffset); + } + private static void WriteWaterBodyIds(BinaryWriter writer, WorldBlueprint bp) { int n = bp.MapSize; diff --git a/Core/Scripts/ConfigManager.cs b/Core/Scripts/ConfigManager.cs index c9a0033..14e233e 100644 --- a/Core/Scripts/ConfigManager.cs +++ b/Core/Scripts/ConfigManager.cs @@ -33,6 +33,13 @@ namespace IslaApocalypse.Core // Change this if your namespace is different // are retired. Default: v5. public static string TerrainCurve = "v5"; + // Terrain detail passes (task 10): "v1" = shelf micro-relief + drainage + // incision as one judged unit (requires the curve; no-op when it is off); + // "off" disables both. ShelfReliefAmp is the micro-relief amplitude in + // metres. Defaults: v1, 3 m. + public static string TerrainDetail = "v1"; + public static float ShelfReliefAmp = 3.0f; + public static void LoadConfig() { string path = "res://ServerConfig.json"; @@ -118,6 +125,20 @@ namespace IslaApocalypse.Core // Change this if your namespace is different GD.PrintErr($"[ConfigManager] Unknown TerrainCurve '{curve}'. Keeping '{TerrainCurve}'."); } + // Extract the terrain-detail gate + relief amplitude + if (data.ContainsKey("TerrainDetail")) + { + string detail = (string)data["TerrainDetail"]; + if (detail == "off" || detail == "v1") + TerrainDetail = detail; + else + GD.PrintErr($"[ConfigManager] Unknown TerrainDetail '{detail}'. Keeping '{TerrainDetail}'."); + } + if (data.ContainsKey("ShelfReliefAmp")) + { + ShelfReliefAmp = (float)data["ShelfReliefAmp"]; + } + switch (profile) { case "4K": diff --git a/Core/Scripts/MapDataParser.cs b/Core/Scripts/MapDataParser.cs index aadc5df..35e4695 100644 --- a/Core/Scripts/MapDataParser.cs +++ b/Core/Scripts/MapDataParser.cs @@ -87,6 +87,20 @@ namespace IslaApocalypse.Core public float K5, K6; } + /// + /// The terrain detail passes that shaped this blueprint's HGTS (v2 TDTL section, + /// terrain-water task 10): shelf micro-relief + drainage incision parameters. + /// Null when detail was off. Metadata only — heights are already detailed. + /// + public class TerrainDetailInfo + { + public ushort Version; + public float ReliefAmpM, ReliefFreqIslands; + public float IncK, IncP, IncCapM; + public float SeaClampRaw, CraterExclFactor, ShelfIncWeight; + public int ReliefSeedOffset; + } + public class WorldBlueprint { public int MapSize; @@ -114,6 +128,9 @@ namespace IslaApocalypse.Core // The curve that shaped HeightMap (v2 TCRV section); null = raw legacy profile. public TerrainCurveInfo TerrainCurve; + + // The detail passes that shaped HeightMap (TDTL section); null = no detail. + public TerrainDetailInfo TerrainDetail; } // 2. The Parser Utility @@ -237,6 +254,7 @@ namespace IslaApocalypse.Core else if (tag == BlueprintFormat.TAG_WATER_SURFACE) sectionOk = ParseWaterGrid(reader, blueprint, payloadLength, isSurface: true); else if (tag == BlueprintFormat.TAG_WATER_BODY_TABLE) sectionOk = ParseWaterBodyTable(reader, blueprint); else if (tag == BlueprintFormat.TAG_TERRAIN_CURVE) sectionOk = ParseTerrainCurve(reader, blueprint); + else if (tag == BlueprintFormat.TAG_TERRAIN_DETAIL) sectionOk = ParseTerrainDetail(reader, blueprint); else { // The property the redesign exists to buy: future sections (water, @@ -430,6 +448,19 @@ namespace IslaApocalypse.Core return true; } + private static bool ParseTerrainDetail(BinaryReader reader, WorldBlueprint blueprint) + { + var d = new TerrainDetailInfo(); + d.Version = reader.ReadUInt16(); + d.ReliefAmpM = reader.ReadSingle(); d.ReliefFreqIslands = reader.ReadSingle(); + d.IncK = reader.ReadSingle(); d.IncP = reader.ReadSingle(); d.IncCapM = reader.ReadSingle(); + d.SeaClampRaw = reader.ReadSingle(); d.CraterExclFactor = reader.ReadSingle(); + d.ShelfIncWeight = reader.ReadSingle(); + d.ReliefSeedOffset = reader.ReadInt32(); + blueprint.TerrainDetail = d; + return true; + } + private static bool ParseRoadTier(BinaryReader reader, List into) { int pathCount = reader.ReadInt32(); diff --git a/Tools/Scripts/MapGenerator.cs b/Tools/Scripts/MapGenerator.cs index 961b38a..ccbc7fb 100644 --- a/Tools/Scripts/MapGenerator.cs +++ b/Tools/Scripts/MapGenerator.cs @@ -46,10 +46,15 @@ public partial class MapGenerator : TextureRect private FastNoiseLite _plateauNoise; private FastNoiseLite _strengthNoise; - // v5: the selected knot preset (task-09 taste batch: COMPACT or BALANCED); - // null = curve off. + // v5: the selected knot preset; null = curve off. private CurveKnots _curveKnots; + // Task-10 detail passes (shelf micro-relief + drainage incision): gated by + // TerrainDetail, active only with the curve on (the masks are curve-band + // defined). The relief noise seeds from resolvedSeed + 7409. + private bool _detailOn; + private FastNoiseLite _reliefNoise; + // The seed's raw pre-curve height maximum (post noise/falloff/Trench/spine, // pre-carve) — the v2 curve's per-seed spike normalizer. Computed in // GenerateTopography pass 1; recorded in TCRV (effective, guard applied). @@ -112,6 +117,11 @@ public partial class MapGenerator : TextureRect _plateauNoise = MakeModulationNoise(HeightCurve.PLATEAU_SEED_OFFSET, HeightCurve.ELEV_FREQ_ISLANDS); _strengthNoise = MakeModulationNoise(HeightCurve.STRENGTH_SEED_OFFSET, HeightCurve.STRENGTH_FREQ_ISLANDS); } + _detailOn = _curveOn && ConfigManager.TerrainDetail == "v1"; + if (_detailOn) + _reliefNoise = MakeModulationNoise(TerrainDetailPass.RELIEF_SEED_OFFSET, TerrainDetailPass.RELIEF_FREQ_ISLANDS); + else if (ConfigManager.TerrainDetail == "v1" && !_curveOn) + GD.Print("[MapGenerator] TerrainDetail v1 requires the curve — no-op with TerrainCurve off."); // THE CRATER FIX: Push it into the ocean (scales via percentage of MapSize!) // Supposedly! We will have to test this manually on other map sizes to confirm the crater is properly scaled and submerged on the north coast! @@ -299,6 +309,18 @@ public partial class MapGenerator : TextureRect StrengthSeedOffset = HeightCurve.STRENGTH_SEED_OFFSET, PresetId = _curveKnots.PresetId, K5 = _curveKnots.K5, K6 = _curveKnots.K6 } : null, + TerrainDetail = _detailOn ? new TerrainDetailInfo + { + Version = TerrainDetailPass.VERSION, + ReliefAmpM = ConfigManager.ShelfReliefAmp, + ReliefFreqIslands = TerrainDetailPass.RELIEF_FREQ_ISLANDS, + IncK = TerrainDetailPass.INC_K, IncP = TerrainDetailPass.INC_P, + IncCapM = TerrainDetailPass.INC_CAP_M, + SeaClampRaw = TerrainDetailPass.SEA_CLAMP, + CraterExclFactor = TerrainDetailPass.CRATER_EXCL_FACTOR, + ShelfIncWeight = TerrainDetailPass.SHELF_INC_WEIGHT, + ReliefSeedOffset = TerrainDetailPass.RELIEF_SEED_OFFSET + } : null, FormatVersion = 2, Params = new BlueprintParams { @@ -462,30 +484,54 @@ public partial class MapGenerator : TextureRect // the curve. Identity at and below sea + this ordering preserve the // Trench/ocean-border guarantee and the crater by construction. classifyH // stays uncurved — see _heightMapClassify; hMaxSeed never touches it. - float physicalCraterRadius = _impactRadius * 0.80f; + // --- PASS 2a: curve + shelf micro-relief (task 10 pass A) --- + // classify stays RAW; curved gets the v5 curve plus, when TerrainDetail is on, + // the shelf-ness-weighted noise skin (risers and peaks untouched). + float reliefAmpRaw = ConfigManager.ShelfReliefAmp / 251f; for (int x = 0; x < MapSize; x++) { for (int y = 0; y < MapSize; y++) { float raw = _heightMap[x, y]; - float classifyH = raw; float curvedH; if (_curveOn) { - // v4: per-column shelf modulation — anchors and strength from the - // low-frequency fields; ordering safety by construction (amplitudes - // bounded; asserted at all 8 field-extreme corners per generation). + // per-column shelf modulation (v4) — anchors and strength from the + // low-frequency fields; ordering safety by construction. float benchLo = HeightCurve.BENCH_BASE + _benchNoise.GetNoise2D(x, y) * HeightCurve.BENCH_AMP; float plateauLo = HeightCurve.PLATEAU_BASE + _plateauNoise.GetNoise2D(x, y) * HeightCurve.PLATEAU_AMP; float shelfSpan = HeightCurve.ShelfSpan((_strengthNoise.GetNoise2D(x, y) + 1f) * 0.5f); curvedH = HeightCurve.Apply(raw, _hMaxSeed, benchLo, shelfSpan, plateauLo, shelfSpan, _curveKnots); + + if (_detailOn) + { + float wShelf = TerrainDetailPass.ShelfWeight(raw, _curveKnots); + if (wShelf > 0f) + curvedH += _reliefNoise.GetNoise2D(x, y) * reliefAmpRaw * wShelf; + } } else { curvedH = raw; } - // --- 5. CARVE THE CRATER (The Flooded Bay & Landbridge Fix!) --- + _heightMapClassify[x, y] = raw; // uncurved; carve joins in pass 2c + _heightMap[x, y] = curvedH; + } + } + + // --- PASS 2b: drainage incision (task 10 pass B) --- + if (_curveOn && _detailOn) + RunIncisionPass(); + + // --- PASS 2c: THE CRATER CARVE (The Flooded Bay & Landbridge Fix!) --- + // The carve remains the FINAL authority on its own terrain: applied after + // curve/relief/incision, to both maps, with the original expression. + float physicalCraterRadius = _impactRadius * 0.80f; + for (int x = 0; x < MapSize; x++) + { + for (int y = 0; y < MapSize; y++) + { float distToCrater = new Vector2(x, y).DistanceTo(_impactCenter); // We only carve the physical hole at 80% of the radius to guarantee a landbridge! @@ -494,16 +540,77 @@ public partial class MapGenerator : TextureRect float craterDepth = 1.0f - (distToCrater / physicalCraterRadius); // Dialed back to -0.15f as per your excellent instinct! float carveTarget = GetSeaLevel(_tempMap[x, y]) - 0.15f; - classifyH = Mathf.Lerp(classifyH, carveTarget, craterDepth * 0.9f); - curvedH = Mathf.Lerp(curvedH, carveTarget, craterDepth * 0.9f); + _heightMapClassify[x, y] = Mathf.Lerp(_heightMapClassify[x, y], carveTarget, craterDepth * 0.9f); + _heightMap[x, y] = Mathf.Lerp(_heightMap[x, y], carveTarget, craterDepth * 0.9f); } - - _heightMapClassify[x, y] = classifyH; - _heightMap[x, y] = curvedH; } } } + /// + /// Task-10 pass B: D8 flow accumulation over the curved+relieved land, then + /// depth = K · accum^p · slope, masked to the risers (shelves feathered to 30 %, + /// toe and peaks zero, crater excluded), capped, and clamped to sea + 1 m. + /// Prints its own MEASURED depth distribution — the tuning/report source. + /// + private void RunIncisionPass() + { + ulong t0 = Time.GetTicksMsec(); + int n = MapSize; + int total = n * n; + + float[] flat = new float[total]; + for (int x = 0; x < n; x++) + for (int y = 0; y < n; y++) + flat[x * n + y] = _heightMap[x, y]; + + int[] accum = TerrainDetailPass.FlowAccumulation(flat, n, out float[] drop); + double accumSeconds = (Time.GetTicksMsec() - t0) / 1000.0; + + float capRaw = TerrainDetailPass.INC_CAP_M / 251f; + float exclR = _impactRadius * TerrainDetailPass.CRATER_EXCL_FACTOR; + float featherR = _impactRadius * TerrainDetailPass.CRATER_FEATHER_FACTOR; + long incised = 0, clampHits = 0; + var depthsM = new System.Collections.Generic.List(1 << 20); + + for (int x = 0; x < n; x++) + { + for (int y = 0; y < n; y++) + { + float raw = _heightMapClassify[x, y]; + float w = TerrainDetailPass.IncisionWeight(raw, _curveKnots); + if (w <= 0f) continue; + + float distToCrater = new Vector2(x, y).DistanceTo(_impactCenter); + if (distToCrater < exclR) continue; + if (distToCrater < featherR) + w *= (distToCrater - exclR) / (featherR - exclR); + + int i = x * n + y; + float depth = TerrainDetailPass.INC_K + * Mathf.Pow(accum[i], TerrainDetailPass.INC_P) * drop[i]; + depth = Mathf.Min(depth, capRaw) * w; + if (depth <= 0f) continue; + + float nh = _heightMap[x, y] - depth; + if (nh < TerrainDetailPass.SEA_CLAMP) + { + nh = TerrainDetailPass.SEA_CLAMP; + clampHits++; + } + float realized = _heightMap[x, y] - nh; + if (realized * 251f >= 0.5f) { incised++; depthsM.Add(realized * 251f); } + _heightMap[x, y] = nh; + } + } + + depthsM.Sort(); + float P(double q) => depthsM.Count == 0 ? 0 : depthsM[Mathf.Clamp((int)(q * depthsM.Count), 0, depthsM.Count - 1)]; + double totalSeconds = (Time.GetTicksMsec() - t0) / 1000.0; + GD.Print($"{T()} [Incision] accumulation {accumSeconds:F1}s, total {totalSeconds:F1}s."); + GD.Print($"{T()} [Incision] incised cells (≥0.5 m): {incised}; depth m: p50 {P(0.5):F1}, p90 {P(0.9):F1}, p99 {P(0.99):F1}, max {(depthsM.Count > 0 ? depthsM[depthsM.Count - 1] : 0):F1}; sea-clamp hits {clampHits}."); + } + private void CalculateTrueOcean() { Queue queue = new Queue(); diff --git a/Tools/Scripts/RoundTripHarness.cs b/Tools/Scripts/RoundTripHarness.cs index 9c0a5cc..10128cb 100644 --- a/Tools/Scripts/RoundTripHarness.cs +++ b/Tools/Scripts/RoundTripHarness.cs @@ -121,6 +121,7 @@ public partial class RoundTripHarness : Node ok &= CompareRoads("Trail", a.TrailRoads, b.TrailRoads); ok &= CompareWater(a, b); ok &= CompareTerrainCurve(a, b); + ok &= CompareTerrainDetail(a, b); if (ok) GD.Print($"[Harness] Semantic equality holds: {a.MapSize}x{a.MapSize} grid, " + @@ -215,6 +216,31 @@ public partial class RoundTripHarness : Node return true; } + private bool CompareTerrainDetail(WorldBlueprint a, WorldBlueprint b) + { + if (a.TerrainDetail == null && b.TerrainDetail == null) + { + GD.Print("[Harness] TDTL: absent in source — nothing to compare (and none reappeared)."); + return true; + } + if (a.TerrainDetail == null || b.TerrainDetail == null) + { + GD.PrintErr("[Harness] TDTL presence mismatch between source and reread."); + return false; + } + var da = a.TerrainDetail; var db = b.TerrainDetail; + float[] fa = { da.Version, da.ReliefAmpM, da.ReliefFreqIslands, da.IncK, da.IncP, da.IncCapM, da.SeaClampRaw, da.CraterExclFactor, da.ShelfIncWeight, da.ReliefSeedOffset }; + float[] fb = { db.Version, db.ReliefAmpM, db.ReliefFreqIslands, db.IncK, db.IncP, db.IncCapM, db.SeaClampRaw, db.CraterExclFactor, db.ShelfIncWeight, db.ReliefSeedOffset }; + for (int i = 0; i < fa.Length; i++) + if (System.BitConverter.SingleToInt32Bits(fa[i]) != System.BitConverter.SingleToInt32Bits(fb[i])) + { + GD.PrintErr("[Harness] TDTL fields differ."); + return false; + } + GD.Print($"[Harness] TDTL equal (detail v{da.Version})."); + return true; + } + private bool CompareRoads(string tier, List a, List b) { if (a.Count != b.Count) diff --git a/Tools/Scripts/TerrainDetailPass.cs b/Tools/Scripts/TerrainDetailPass.cs new file mode 100644 index 0000000..756160d --- /dev/null +++ b/Tools/Scripts/TerrainDetailPass.cs @@ -0,0 +1,143 @@ +using Godot; +using System; + +/// +/// The terrain DETAIL passes (terrain-water task 10) — pure numeric array machinery +/// (D-035; a named future C++ candidate, kept standalone): +/// +/// PASS A — shelf micro-relief: a medium-frequency noise skin (±ShelfReliefAmp, +/// default 3 m) weighted by shelf-ness, so the compressed shelves get their +/// rolling texture back while risers and peaks stay untouched. +/// +/// PASS B — drainage incision: D8 steepest-descent flow routing + accumulation +/// over the curved terrain; depth = K · accum^p · localSlope (capped), masked to +/// the risers (feathered ~30 % onto shelves, zero on the toe and above the +/// plateau top, zero near the crater), clamped so carved terrain never drops +/// below sea + 1 m. The channels double as the future river routes (Phase C). +/// +/// Ordering (enforced by the caller): curve → micro-relief → incision → crater +/// carve. The classify map never sees any of it. +/// +public static class TerrainDetailPass +{ + public const ushort VERSION = 1; + + // Pass A — micro-relief. + public const float RELIEF_AMP_DEFAULT_M = 3f; // config dial: ShelfReliefAmp (metres) + public const float RELIEF_FREQ_ISLANDS = 40f; // ~40 undulations per island width (~200 m features) + public const int RELIEF_SEED_OFFSET = 7409; + + // Pass B — incision. K/p tuned against the depth targets (gullies 8–15 m, + // trunks ~25 m, cap 30 m); the tuning run's achieved distribution is in the + // task-10 report. + public const float INC_K = 0.78f; + public const float INC_P = 0.45f; // concave: many fingers, few deep trunks + public const float INC_CAP_M = 30f; // IncisionMax + public const float SEA_CLAMP = 0.15f + 1f / 251f; // carved height ≥ sea + 1 m + public const float SHELF_INC_WEIGHT = 0.3f; // shelves get washes, not gorges + public const float CRATER_EXCL_FACTOR = 1.2f; // zero incision inside this × CraterRadius + public const float CRATER_FEATHER_FACTOR = 1.4f; // ...feathering to full by this × CraterRadius + + /// + /// Shelf-ness weight from the RAW input height: 1 mid-shelf, feathering to 0 + /// through the risers (feather extends 30 % of the band half-width past each + /// shelf edge). Covers both shelves. + /// + public static float ShelfWeight(float raw, CurveKnots k) + { + return Mathf.Max(BandBump(raw, k.K3, k.K4), BandBump(raw, k.K5, k.K6)); + } + + private static float BandBump(float h, float lo, float hi) + { + float half = (hi - lo) * 0.5f; + float t = Mathf.Abs(h - (lo + half)) / half; // 0 centre, 1 at band edge + // full inside 60 % of the band, linear feather to zero at 130 % + return Mathf.Clamp(1f - (t - 0.6f) / 0.7f, 0f, 1f); + } + + /// + /// Incision mask from the RAW input height: 0 below the red-ceiling input (K2) + /// and above the plateau top (K6); 1 on the riser bands; SHELF_INC_WEIGHT on the + /// shelf bands; smooth feathers (15 % of the local band width) at every boundary. + /// + public static float IncisionWeight(float raw, CurveKnots k) + { + if (raw <= k.K2 || raw >= k.K6) return 0f; + if (raw < k.K3) // foothill riser: feather in from K2, feather toward shelf weight at K3 + return EdgeBlend(raw, k.K2, k.K3, 0f, 1f, SHELF_INC_WEIGHT); + if (raw < k.K4) // bench + return SHELF_INC_WEIGHT; + if (raw < k.K5) // mid riser + return EdgeBlend(raw, k.K4, k.K5, SHELF_INC_WEIGHT, 1f, SHELF_INC_WEIGHT); + // plateau band: shelf weight, feathering to zero at K6 + float w = (k.K6 - raw) / ((k.K6 - k.K5) * 0.15f); + return Mathf.Min(SHELF_INC_WEIGHT, Mathf.Clamp(w, 0f, 1f) * SHELF_INC_WEIGHT); + } + + private static float EdgeBlend(float h, float lo, float hi, float wIn, float wMid, float wOut) + { + float f = (hi - lo) * 0.15f; + if (h < lo + f) return Mathf.Lerp(wIn, wMid, (h - lo) / f); + if (h > hi - f) return Mathf.Lerp(wMid, wOut, (h - (hi - f)) / f); + return wMid; + } + + /// + /// D8 flow accumulation over a height field (row-major idx = x·n + y). + /// Steepest-descent routing (drop / distance, diagonals ÷√2), deterministic + /// tie-break (fixed neighbour order, first winner). Cells with no lower + /// neighbour are pits/outlets (no outflow). accum = upslope contributing cells + /// including self; steepestDrop = drop per pixel toward the chosen neighbour. + /// + public static int[] FlowAccumulation(float[] h, int n, out float[] steepestDrop) + { + int total = n * n; + int[] downstream = new int[total]; + steepestDrop = new float[total]; + int[] dx = { 1, -1, 0, 0, 1, 1, -1, -1 }; + int[] dy = { 0, 0, 1, -1, 1, -1, 1, -1 }; + float[] invDist = { 1f, 1f, 1f, 1f, 0.7071068f, 0.7071068f, 0.7071068f, 0.7071068f }; + + for (int x = 0; x < n; x++) + { + for (int y = 0; y < n; y++) + { + int i = x * n + y; + float hc = h[i]; + float best = 0f; + int bestIdx = -1; + for (int d = 0; d < 8; d++) + { + int nx = x + dx[d], ny = y + dy[d]; + if (nx < 0 || nx >= n || ny < 0 || ny >= n) continue; + int ni = nx * n + ny; + float grade = (hc - h[ni]) * invDist[d]; + if (grade > best) + { + best = grade; + bestIdx = ni; + } + } + downstream[i] = bestIdx; + steepestDrop[i] = best; + } + } + + // Height-descending order: each cell pushes its accumulated count downstream. + float[] keys = (float[])h.Clone(); + int[] order = new int[total]; + for (int i = 0; i < total; i++) order[i] = i; + Array.Sort(keys, order); // ascending + + int[] accum = new int[total]; + for (int i = 0; i < total; i++) accum[i] = 1; + for (int i = total - 1; i >= 0; i--) + { + int c = order[i]; + int d = downstream[c]; + if (d >= 0) accum[d] += accum[c]; + } + return accum; + } +}