beezm/islaApocalypse
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

feat: shelf micro-relief + drainage incision (terrain-water task 10)

Browse source 
Two paired detail passes on the locked v5 skeleton, output-heights
only (classify path untouched):

PASS A — shelf micro-relief: +-ShelfReliefAmp (default 3 m) Simplex
skin (seed resolved+7409, ~40 undulations/island ~ 200 m features),
weighted by shelf-ness (full mid-shelf, feathered to zero 30% of a
band half-width into the risers) — both shelves get their rolling
texture back; risers and peaks untouched.

PASS B — drainage incision: D8 steepest-descent routing + height-
ordered flow accumulation over the curved+relieved land
(TerrainDetailPass — pure array machinery, named C++ candidate per
D-035); depth = K*accum^p*slope (K=0.78, p=0.45 concave, cap 30 m),
masked full on risers / 30% on shelves / zero on the toe, above the
plateau top, and within 1.2x CraterRadius (feathered to 1.4x); hard
clamp: carved height >= sea + 1 m. Ordering: curve -> relief ->
incision -> crater carve (the carve stays the final authority; carve
moved to its own pass 2c, bit-identical expression on both maps).
The pass prints its own MEASURED depth distribution.

Gate TerrainDetail "off"|"v1" (default v1; both passes one judged
unit; no-op without the curve) + ShelfReliefAmp dial. New TDTL section
(38 B: version + relief/incision params + seed offset) across
writer/parser/harness — blueprints stay self-describing.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
This commit is contained in:
Stewart Howe 2026-08-08 18:08:28 -04:00
parent c71026fd94
commit 45fe2ce034
7 changed files with 354 additions and 13 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

1
Core/Scripts/BlueprintFormat.cs
View file

@ -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

12
Core/Scripts/BlueprintWriter.cs
View file

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

21
Core/Scripts/ConfigManager.cs
View file

@ -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":

31
Core/Scripts/MapDataParser.cs
View file

@ -87,6 +87,20 @@ namespace IslaApocalypse.Core
public float K5, K6;
}
/// <summary>
/// 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.
/// </summary>
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<Vector2[]> into)
{
int pathCount = reader.ReadInt32();

133
Tools/Scripts/MapGenerator.cs
View file

@ -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;
}
}
}
/// <summary>
/// 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.
/// </summary>
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<float>(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<Vector2I> queue = new Queue<Vector2I>();

26
Tools/Scripts/RoundTripHarness.cs
View file

@ -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<Vector2[]> a, List<Vector2[]> b)
{
if (a.Count != b.Count)

143
Tools/Scripts/TerrainDetailPass.cs Normal file
View file

@ -0,0 +1,143 @@
using Godot;
using System;
/// <summary>
/// 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.
/// </summary>
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
/// <summary>
/// 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.
/// </summary>
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);
}
/// <summary>
/// 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.
/// </summary>
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;
}
/// <summary>
/// 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.
/// </summary>
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;
}
}