nx_terrain_erode_lib.nx source

1// nx_terrain_erode_lib.nx -- TERRAIN-SCALE EROSION COMPOSITION + ITS REFEREE (PG20, procgen.plan). 2// 3// WHY THIS LIB EXISTS. nx_water_erosion.nx ships all seven Mei 2007 stages and is gate-proven 4// (nx_water_erosion_gate, PG18) -- and the ranker measured it LIB-UNIMPORTED: nothing in the estate 5// composed it. PG20's contract is that composition onto a terrain heightfield, plus the thermal-talus 6// companion Musgrave 1989 pairs with hydraulic erosion, plus the referee that decides whether the 7// result reads as terrain rather than as blobs. 8// 9// THREE THINGS LIVE HERE AND NOTHING ELSE: 10// 11// 1. te_thermal_pass -- THE ONE TALUS PASS IN THE ESTATE. 12// nx_pets_voxel3d.vp_erode_pass already implemented the voxel form (shed exactly one unit to every 13// 4-neighbour overhung by more than talus, mass-conserving by construction). That is a real 14// lower-layer kernel and it is NOT deleted: this function GENERALISES it by one parameter and 15// REPRODUCES it exactly at flux_q14 == 0, which is what makes retiring the private copy onto this 16// owner provable rather than asserted. flux_q14 > 0 selects the Musgrave proportional law that a 17// continuous heightfield needs, where a fixed one-unit shed would be a voxel quantum applied to a 18// quantity that has no voxels. 19// INDEX ORDER IS NOT A CONVENTION HERE, IT IS A PROVEN PROPERTY. vp_erode_pass indexes x*gn+z and 20// nx_water_erosion indexes y*w+x. On a square grid the four neighbours of cell k are k-1, k+1, 21// k-gn, k+gn under BOTH, the boundary guards map onto each other exactly, and the deltas are 22// accumulated additively from the pre-pass state -- so the two orders produce byte-identical 23// output. The gate asserts that against a reference written in the other order rather than 24// leaving it as a comment nobody can check. 25// 26// 2. te_hypsometric_permil -- THE REFEREE. Strahler's hypsometric integral, 27// (mean - min) / (max - min), in permille. It is the right referee for a procedural world for a 28// reason that is decidable before any measurement is taken: IT IS DIMENSIONLESS. The plan named 29// two published bands, hypsometric curve and DRAINAGE DENSITY -- and drainage density is 30// kilometres of channel per square kilometre, so quoting a published band for it against a world 31// with no metric scale binding would be a number about a subject we do not have. That half is 32// declared open by name in the gate, not half-measured. 33// The integral is invariant under h -> a*h + b for a > 0, which is exactly what makes it 34// non-gameable: an implementation that merely lowers the terrain, or flattens it, or adds noise, 35// moves it NOWHERE. Only moving mass from highs into lows moves it, and that is erosion. 36// A flat field has no hypsometry, so this ABSTAINS with -1 rather than returning a number. 37// 38// 3. te_erode_grid -- the composed loop: seed the Mei state from a caller-owned height grid, tick the 39// full seven-stage pipeline, interleave the talus pass, read the bed back. Every buffer is 40// caller-owned and sized by te_state_bytes / te_grid_bytes, so there is no cap to guess and no 41// allocation inside any loop. 42// 43// SCALE: this lib is SCALE-FREE. nx_water_erosion names its units metres; every quantity that decides 44// behaviour here is a RATIO (slope = height delta / cell size), so a caller passing height and cell 45// size in the same arbitrary world unit gets the same physics with no metre assumption anywhere. 46// 47// license_tier: ORIGINAL 48// genealogy_id: mei_2007_fast_hydraulic_erosion_gpu + musgrave_1989_eroded_fractal_terrains + 49// strahler_1952_hypsometric_analysis 50 51import "nx_syscalls.nx" 52import "nx_tier.nx" 53import "nx_water_erosion.nx" 54 55// Permille scale for the hypsometric integral. Named for its PURPOSE (the estate reports ratios in 56// permille) and used nowhere else in this file. 57const TE_PERMIL: i64 = 1000 58// The abstention value: "this field has no hypsometry", never a height and never a ratio. 59const TE_HYPSO_UNOBSERVABLE: i64 = 0 - 1 60 61// ===== buffer sizing, DERIVED from the composed organ's own stride ===== 62func te_state_bytes(gn: i64) -> i64 { return gn * gn * NX_WATER_STRIDE * NX_SIZEOF_NX_INT } 63func te_grid_bytes(gn: i64) -> i64 { return gn * gn * NX_SIZEOF_NX_INT } 64 65// ===== the talus law, one site ===== 66// Sheds material from src to dst when the overhang exceeds talus. Returns the amount moved (0 = none). 67// flux_q14 == 0 -> the VOXEL law: exactly one unit, reproducing nx_pets_voxel3d.vp_erode_pass. 68// flux_q14 > 0 -> the MUSGRAVE law: a fraction of the excess, halved because the pair of cells shares 69// the excess between them (moving all of it would overshoot past level). 70func te_shed(hg: *i64, delta: *i64, src: i64, dst: i64, h: i64, talus: i64, flux_q14: i64) -> i64 { 71 let dh: i64 = h - hg[dst] 72 if dh <= talus { return 0 } 73 var amt: i64 = 1 74 if flux_q14 > 0 { amt = (dh - talus) * flux_q14 / NX_WATER_Q / 2 } 75 if amt <= 0 { return 0 } 76 delta[src] = delta[src] - amt 77 delta[dst] = delta[dst] + amt 78 return amt 79} 80 81// One thermal-erosion pass over a square gn x gn grid, row-major (idx = y*gn + x). 82// Deltas are computed from the PRE-PASS state and applied afterwards, so the result does not depend 83// on iteration order. Every -amt has a matching +amt, so sum(hg) is invariant EXACTLY. 84// Returns the total material moved, so a caller can tell a real pass from a no-op. 85func te_thermal_pass(hg: *i64, gn: i64, talus: i64, flux_q14: i64, delta: *i64) -> i64 { 86 let n: i64 = gn * gn 87 var i: i64 = 0 88 while i < n { delta[i] = 0; i = i + 1 } 89 var moved: i64 = 0 90 var y: i64 = 0 91 while y < gn { 92 var x: i64 = 0 93 while x < gn { 94 let k: i64 = y * gn + x 95 let h: i64 = hg[k] 96 if x > 0 { moved = moved + te_shed(hg, delta, k, k - 1, h, talus, flux_q14) } 97 if x < gn - 1 { moved = moved + te_shed(hg, delta, k, k + 1, h, talus, flux_q14) } 98 if y > 0 { moved = moved + te_shed(hg, delta, k, k - gn, h, talus, flux_q14) } 99 if y < gn - 1 { moved = moved + te_shed(hg, delta, k, k + gn, h, talus, flux_q14) } 100 x = x + 1 101 } 102 y = y + 1 103 } 104 var j: i64 = 0 105 while j < n { hg[j] = hg[j] + delta[j]; j = j + 1 } 106 return moved 107} 108 109// ===== the referee ===== 110// Strahler's hypsometric integral in permille. Returns TE_HYPSO_UNOBSERVABLE on a flat or empty 111// field: a field with no relief has no hypsometry, and reporting a number there would acquit on a 112// measurement that could not be taken. 113func te_hypsometric_permil(hg: *i64, n: i64) -> i64 { 114 if n <= 0 { return TE_HYPSO_UNOBSERVABLE } 115 var mn: i64 = hg[0] 116 var mx: i64 = hg[0] 117 var sum: i64 = 0 118 var i: i64 = 0 119 while i < n { 120 let v: i64 = hg[i] 121 if v < mn { mn = v } 122 if v > mx { mx = v } 123 sum = sum + v 124 i = i + 1 125 } 126 if mx == mn { return TE_HYPSO_UNOBSERVABLE } 127 let mean: i64 = sum / n 128 return (mean - mn) * TE_PERMIL / (mx - mn) 129} 130 131// ===== state marshalling: only the BED field crosses, water and sediment persist across ticks ===== 132func te_seed_bed(state: *i64, hg: *i64, n: i64) -> i64 { 133 var i: i64 = 0 134 while i < n { state[i * NX_WATER_STRIDE + NX_WATER_OFF_H] = hg[i]; i = i + 1 } 135 return n 136} 137func te_read_bed(state: *i64, hg: *i64, n: i64) -> i64 { 138 var i: i64 = 0 139 while i < n { hg[i] = state[i * NX_WATER_STRIDE + NX_WATER_OFF_H]; i = i + 1 } 140 return n 141} 142 143// ===== the composed loop ===== 144// Every buffer is caller-owned: state (te_state_bytes), scratch and delta (te_grid_bytes each). 145// nx_water_erosion_init runs ONCE, so water, sediment and flux carry across ticks as the pipe model 146// requires; only the bed round-trips each tick so the talus pass can act on it. 147// Returns the total material the talus passes moved over the whole run. 148func te_erode_grid(hg: *i64, gn: i64, ticks: i64, rain_q14: i64, talus: i64, flux_q14: i64, 149 params: *i64, state: *i64, scratch: *i64, delta: *i64) -> i64 { 150 let n: i64 = gn * gn 151 te_seed_bed(state, hg, n) 152 nx_water_erosion_init(state, gn, gn) 153 var moved: i64 = 0 154 var t: i64 = 0 155 while t < ticks { 156 nx_water_erosion_tick(state, scratch, gn, gn, rain_q14, params) 157 te_read_bed(state, hg, n) 158 moved = moved + te_thermal_pass(hg, gn, talus, flux_q14, delta) 159 te_seed_bed(state, hg, n) 160 t = t + 1 161 } 162 te_read_bed(state, hg, n) 163 return moved 164}