jpmschweitzerandClaude Opus 4.8
0263b68df9
T-1044: run_layer1 now returns TerrainAnalysis (carried transiently on CascadeSnapshot, dropped after the district + road-graph passes), eliminating the redundant per-body drainage::analyze + TerrainAnalysis::analyze re-run flagged by PERF/TODO(T-1044). Not persisted on the LRU-cached state (D-203/T-1048 size concern). T-1047: basin_direction is now derived from the real D8 thalweg. run_layer1 aggregates a per-district dominant D8 direction from the live fdir grid (carried transiently on DrainageResult), threaded via Layer1Output.district_basin_dirs -> derive_all_districts -> DistrictProfile.basin_direction; derive_chunk_context reads it directly. Removed the false derive_basin_direction (it branched on ocean_fraction_q then read seed bits despite a doc comment claiming an elev_q/slope_q D8 proxy) + corrected the module contract. D-239 §8 (D8 thalweg) now actually honoured. 1559 tests pass; golden byte-identical (district_basin_dirs is #[serde(skip)], transient). Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
810 lines
27 KiB
Rust
810 lines
27 KiB
Rust
//! D8 drainage routing — flow direction, flow accumulation, river network
|
||
//! extraction, and drainage basin delineation (D-208).
|
||
//!
|
||
//! **Determinism (D-010, D-208):** All flow-direction comparisons use integer
|
||
//! arithmetic on scaled elevation values (`(elev * 1_000_000.0) as i64`) to
|
||
//! avoid f32 comparison non-determinism. Tie-breaking uses a fixed D8 neighbor
|
||
//! priority order. The result is bit-identical across runs on the same inputs.
|
||
//!
|
||
//! **Algorithm:**
|
||
//! 1. Scale f32 elevation to i64 integers.
|
||
//! 2. Priority-flood depression fill (iterative, convergence in ≤10 passes).
|
||
//! 3. D8 flow direction: steepest descent, 8-neighbor, wraps horizontally.
|
||
//! 4. Flow accumulation via topological sort of the D8 DAG.
|
||
//! 5. River network extraction: cells with accumulation > RIVER_THRESHOLD.
|
||
//! 6. Basin labeling: flood-fill seeded at pour points.
|
||
//!
|
||
//! The grid is row-major. Row 0 is the north pole; row H-1 is the south pole.
|
||
//! Columns wrap horizontally (the globe is equirectangular).
|
||
|
||
use std::collections::VecDeque;
|
||
|
||
use crate::atlas::body_world_state::{DrainageBasin, RiverNetwork};
|
||
use crate::simulation::generator::TerritorialStatus;
|
||
|
||
/// A cell is a river cell when its flow accumulation exceeds this threshold (D-208).
|
||
pub const RIVER_THRESHOLD: i32 = 200;
|
||
|
||
/// Scale factor for converting f32 elevation to integer for deterministic comparison.
|
||
const ELEV_SCALE: f64 = 1_000_000.0;
|
||
|
||
// D8 neighbor offsets (dr, dc) in fixed priority order for deterministic tie-breaking.
|
||
// Priority: cardinal directions first (N, S, E, W), then diagonals (NE, NW, SE, SW).
|
||
const D8: [(i32, i32); 8] = [
|
||
(-1, 0), // N
|
||
(1, 0), // S
|
||
(0, 1), // E
|
||
(0, -1), // W
|
||
(-1, 1), // NE
|
||
(-1, -1), // NW
|
||
(1, 1), // SE
|
||
(1, -1), // SW
|
||
];
|
||
|
||
/// Result of the full D8 drainage analysis for one body.
|
||
#[derive(Debug, Clone)]
|
||
pub struct DrainageResult {
|
||
pub river_network: RiverNetwork,
|
||
pub drainage_basins: Vec<DrainageBasin>,
|
||
/// Per-cell flow accumulation (row-major, `w × h`): the upstream cell count
|
||
/// draining through each cell. Exposed for D-209 attractor-strength
|
||
/// normalization (`flow_accumulation[cell] / max_accumulation`).
|
||
pub flow_accumulation: Vec<i32>,
|
||
/// Maximum flow accumulation across the grid — the denominator for
|
||
/// normalized attractor strength (D-209). Always ≥ 1.
|
||
pub max_accumulation: i32,
|
||
/// Per-cell D8 flow-direction index into `D8` (0–7), or -1 for no outflow
|
||
/// (edge, flat peak, or ocean). Row-major, `w × h`.
|
||
///
|
||
/// **Transient — used within the Layer-1 pass only.** The caller aggregates a
|
||
/// per-district dominant direction from this grid (T-1047) and carries that
|
||
/// compact result on `Layer1Output.district_basin_dirs`; the full 131 KB
|
||
/// grid is NOT persisted on `BodyWorldState` or the LRU cache (D-203).
|
||
pub fdir: Vec<i8>,
|
||
}
|
||
|
||
// ---------------------------------------------------------------------------
|
||
// Public entry point
|
||
// ---------------------------------------------------------------------------
|
||
|
||
/// Run the full D8 drainage analysis on an elevation grid.
|
||
///
|
||
/// `elevation` is a row-major float32 grid of shape `height × width`, values
|
||
/// in [0.0, 1.0]. `sea_level` is the fraction below which terrain is ocean.
|
||
///
|
||
/// Returns `DrainageResult` with the river network and drainage basins.
|
||
pub fn analyze(elevation: &[f32], width: u32, height: u32, sea_level: f32) -> DrainageResult {
|
||
let w = width as usize;
|
||
let h = height as usize;
|
||
// 1. Scale to integers.
|
||
let scaled: Vec<i64> = elevation
|
||
.iter()
|
||
.map(|&e| (e as f64 * ELEV_SCALE) as i64)
|
||
.collect();
|
||
|
||
// 2. Depression fill.
|
||
let filled = depression_fill(&scaled, w, h);
|
||
|
||
// 3. D8 flow direction. -1 = no outflow (edge or flat peak).
|
||
let fdir = flow_direction(&filled, w, h);
|
||
|
||
// 4. Flow accumulation.
|
||
let accum = flow_accumulation(&fdir, w, h);
|
||
|
||
// 5. River network.
|
||
let river_network = extract_river_network(&accum, &fdir, w, h, sea_level, elevation);
|
||
|
||
// 6. Basin labeling.
|
||
let labels = label_basins(&fdir, &accum, w, h);
|
||
|
||
// 7. Merge small basins + clamp count to [4, 12].
|
||
let labels = merge_small_basins(labels, w, h, 4, 12);
|
||
|
||
// 8. Build DrainageBasin structs.
|
||
let drainage_basins = build_basins(&labels, w, h);
|
||
|
||
// Max accumulation for D-209 strength normalization (clamped ≥ 1 so the
|
||
// division is always well-defined, even on a flat/empty world).
|
||
let max_accumulation = accum.iter().copied().max().unwrap_or(1).max(1);
|
||
|
||
DrainageResult {
|
||
river_network,
|
||
drainage_basins,
|
||
flow_accumulation: accum,
|
||
max_accumulation,
|
||
fdir,
|
||
}
|
||
}
|
||
|
||
// ---------------------------------------------------------------------------
|
||
// Step 2: Depression fill
|
||
// ---------------------------------------------------------------------------
|
||
|
||
fn depression_fill(scaled: &[i64], w: usize, h: usize) -> Vec<i64> {
|
||
let mut filled = scaled.to_vec();
|
||
for _ in 0..10 {
|
||
let mut changed = false;
|
||
for r in 1..h.saturating_sub(1) {
|
||
for c in 0..w {
|
||
let mut nbr_min = i64::MAX;
|
||
for &(dr, dc) in &D8 {
|
||
let nr = r as i32 + dr;
|
||
let nc = (c as i32 + dc).rem_euclid(w as i32) as usize;
|
||
if nr >= 0 && nr < h as i32 {
|
||
let val = filled[nr as usize * w + nc];
|
||
if val < nbr_min {
|
||
nbr_min = val;
|
||
}
|
||
}
|
||
}
|
||
if filled[r * w + c] < nbr_min {
|
||
filled[r * w + c] = nbr_min + 1;
|
||
changed = true;
|
||
}
|
||
}
|
||
}
|
||
if !changed {
|
||
break;
|
||
}
|
||
}
|
||
filled
|
||
}
|
||
|
||
// ---------------------------------------------------------------------------
|
||
// Step 3: D8 flow direction
|
||
// ---------------------------------------------------------------------------
|
||
|
||
/// Returns per-cell flow direction index into D8 (0–7), or -1 for no outflow.
|
||
fn flow_direction(filled: &[i64], w: usize, h: usize) -> Vec<i8> {
|
||
let mut fdir = vec![-1i8; w * h];
|
||
for r in 0..h {
|
||
for c in 0..w {
|
||
let elev = filled[r * w + c];
|
||
let mut best_drop = 0i64;
|
||
let mut best_k: i8 = -1;
|
||
for (k, &(dr, dc)) in D8.iter().enumerate() {
|
||
let nr = r as i32 + dr;
|
||
let nc = (c as i32 + dc).rem_euclid(w as i32) as usize;
|
||
if nr < 0 || nr >= h as i32 {
|
||
continue;
|
||
}
|
||
let drop = elev - filled[nr as usize * w + nc];
|
||
if drop > best_drop {
|
||
best_drop = drop;
|
||
best_k = k as i8;
|
||
}
|
||
}
|
||
fdir[r * w + c] = best_k;
|
||
}
|
||
}
|
||
fdir
|
||
}
|
||
|
||
// ---------------------------------------------------------------------------
|
||
// Step 4: Flow accumulation
|
||
// ---------------------------------------------------------------------------
|
||
|
||
fn flow_accumulation(fdir: &[i8], w: usize, h: usize) -> Vec<i32> {
|
||
let n = w * h;
|
||
let mut in_degree = vec![0i32; n];
|
||
|
||
for r in 0..h {
|
||
for c in 0..w {
|
||
let k = fdir[r * w + c];
|
||
if k < 0 {
|
||
continue;
|
||
}
|
||
let (dr, dc) = D8[k as usize];
|
||
let nr = r as i32 + dr;
|
||
let nc = (c as i32 + dc).rem_euclid(w as i32) as usize;
|
||
if nr >= 0 && nr < h as i32 {
|
||
in_degree[nr as usize * w + nc] += 1;
|
||
}
|
||
}
|
||
}
|
||
|
||
let mut queue = VecDeque::new();
|
||
for (i, °) in in_degree.iter().enumerate().take(n) {
|
||
if deg == 0 {
|
||
queue.push_back(i);
|
||
}
|
||
}
|
||
|
||
let mut accum = vec![1i32; n];
|
||
while let Some(idx) = queue.pop_front() {
|
||
let r = idx / w;
|
||
let c = idx % w;
|
||
let k = fdir[idx];
|
||
if k < 0 {
|
||
continue;
|
||
}
|
||
let (dr, dc) = D8[k as usize];
|
||
let nr = r as i32 + dr;
|
||
let nc = (c as i32 + dc).rem_euclid(w as i32) as usize;
|
||
if nr >= 0 && nr < h as i32 {
|
||
let ni = nr as usize * w + nc;
|
||
accum[ni] += accum[idx];
|
||
in_degree[ni] -= 1;
|
||
if in_degree[ni] == 0 {
|
||
queue.push_back(ni);
|
||
}
|
||
}
|
||
}
|
||
|
||
accum
|
||
}
|
||
|
||
// ---------------------------------------------------------------------------
|
||
// Step 5: River network extraction
|
||
// ---------------------------------------------------------------------------
|
||
|
||
fn extract_river_network(
|
||
accum: &[i32],
|
||
fdir: &[i8],
|
||
w: usize,
|
||
h: usize,
|
||
sea_level: f32,
|
||
elevation: &[f32],
|
||
) -> RiverNetwork {
|
||
let n = w * h;
|
||
|
||
// River cells: above threshold AND above sea level.
|
||
let is_river: Vec<bool> = (0..n)
|
||
.map(|i| accum[i] > RIVER_THRESHOLD && elevation[i] >= sea_level)
|
||
.collect();
|
||
|
||
let river_cells: Vec<(u16, u16)> = (0..n)
|
||
.filter(|&i| is_river[i])
|
||
.map(|i| ((i / w) as u16, (i % w) as u16))
|
||
.collect();
|
||
|
||
// Confluences: river cells with 2+ river neighbors flowing into them.
|
||
let mut inflow_count = vec![0u8; n];
|
||
for r in 0..h {
|
||
for c in 0..w {
|
||
let i = r * w + c;
|
||
if !is_river[i] {
|
||
continue;
|
||
}
|
||
let k = fdir[i];
|
||
if k < 0 {
|
||
continue;
|
||
}
|
||
let (dr, dc) = D8[k as usize];
|
||
let nr = r as i32 + dr;
|
||
let nc = (c as i32 + dc).rem_euclid(w as i32) as usize;
|
||
if nr >= 0 && nr < h as i32 {
|
||
let ni = nr as usize * w + nc;
|
||
if is_river[ni] {
|
||
inflow_count[ni] = inflow_count[ni].saturating_add(1);
|
||
}
|
||
}
|
||
}
|
||
}
|
||
let confluences: Vec<(u16, u16)> = (0..n)
|
||
.filter(|&i| is_river[i] && inflow_count[i] >= 2)
|
||
.map(|i| ((i / w) as u16, (i % w) as u16))
|
||
.collect();
|
||
|
||
// Mouths: river cells that flow to a sea cell or to the polar edge.
|
||
let mouths: Vec<(u16, u16)> = (0..n)
|
||
.filter(|&i| {
|
||
if !is_river[i] {
|
||
return false;
|
||
}
|
||
let r = i / w;
|
||
let c = i % w;
|
||
let k = fdir[i];
|
||
if k < 0 {
|
||
return true; // no outflow — edge
|
||
}
|
||
let (dr, dc) = D8[k as usize];
|
||
let nr = r as i32 + dr;
|
||
let nc = (c as i32 + dc).rem_euclid(w as i32) as usize;
|
||
if nr < 0 || nr >= h as i32 {
|
||
return true; // polar edge
|
||
}
|
||
// Flows into a sub-sea-level cell = mouth
|
||
elevation[nr as usize * w + nc] < sea_level
|
||
})
|
||
.map(|i| ((i / w) as u16, (i % w) as u16))
|
||
.collect();
|
||
|
||
RiverNetwork {
|
||
river_cells,
|
||
confluences,
|
||
mouths,
|
||
}
|
||
}
|
||
|
||
// ---------------------------------------------------------------------------
|
||
// Step 6: Basin labeling
|
||
// ---------------------------------------------------------------------------
|
||
|
||
fn label_basins(fdir: &[i8], accum: &[i32], w: usize, h: usize) -> Vec<i32> {
|
||
let n = w * h;
|
||
let mut labels = vec![-1i32; n];
|
||
|
||
// Pour points: local accumulation maxima above river threshold.
|
||
let mut pour_pts: Vec<usize> = Vec::new();
|
||
for i in 0..n {
|
||
if accum[i] <= RIVER_THRESHOLD {
|
||
continue;
|
||
}
|
||
let r = i / w;
|
||
let c = i % w;
|
||
let mut is_max = true;
|
||
for &(dr, dc) in &D8 {
|
||
let nr = r as i32 + dr;
|
||
let nc = (c as i32 + dc).rem_euclid(w as i32) as usize;
|
||
if nr >= 0 && nr < h as i32 && accum[nr as usize * w + nc] > accum[i] {
|
||
is_max = false;
|
||
break;
|
||
}
|
||
}
|
||
if is_max {
|
||
pour_pts.push(i);
|
||
}
|
||
}
|
||
|
||
if pour_pts.is_empty() {
|
||
// Flat/ocean world — single basin.
|
||
labels.iter_mut().for_each(|l| *l = 0);
|
||
return labels;
|
||
}
|
||
|
||
for (basin_id, &idx) in pour_pts.iter().enumerate() {
|
||
labels[idx] = basin_id as i32;
|
||
}
|
||
|
||
// Trace remaining cells: follow fdir until a labeled cell is reached.
|
||
for start in 0..n {
|
||
if labels[start] >= 0 {
|
||
continue;
|
||
}
|
||
// Walk forward, accumulate path.
|
||
let mut path: Vec<usize> = Vec::new();
|
||
let mut cur = start;
|
||
let label = loop {
|
||
if labels[cur] >= 0 {
|
||
break labels[cur];
|
||
}
|
||
path.push(cur);
|
||
let k = fdir[cur];
|
||
if k < 0 {
|
||
break 0; // no outflow — assign to basin 0
|
||
}
|
||
let r = cur / w;
|
||
let c = cur % w;
|
||
let (dr, dc) = D8[k as usize];
|
||
let nr = r as i32 + dr;
|
||
let nc = (c as i32 + dc).rem_euclid(w as i32) as usize;
|
||
if nr < 0 || nr >= h as i32 {
|
||
break 0; // polar edge
|
||
}
|
||
let next = nr as usize * w + nc;
|
||
// Cycle guard: if we're visiting a cell already in path, stop.
|
||
if path.contains(&next) {
|
||
break 0;
|
||
}
|
||
cur = next;
|
||
};
|
||
for idx in path {
|
||
labels[idx] = label;
|
||
}
|
||
}
|
||
|
||
labels
|
||
}
|
||
|
||
// ---------------------------------------------------------------------------
|
||
// Step 7: Merge small basins
|
||
// ---------------------------------------------------------------------------
|
||
|
||
/// Union-find root with path compression.
|
||
fn uf_find(parent: &mut [i32], x: i32) -> i32 {
|
||
let mut root = x;
|
||
while parent[root as usize] != root {
|
||
root = parent[root as usize];
|
||
}
|
||
let mut cur = x;
|
||
while parent[cur as usize] != root {
|
||
let next = parent[cur as usize];
|
||
parent[cur as usize] = root;
|
||
cur = next;
|
||
}
|
||
root
|
||
}
|
||
|
||
/// Merge small basins into their largest neighbor until the count is in
|
||
/// `[min_count, max_count]` and every basin holds ≥ 2% of the surface.
|
||
///
|
||
/// Builds a basin adjacency graph + sizes in a single grid pass, then performs
|
||
/// all merges as union-find operations on that graph — the grid is rewritten
|
||
/// exactly once at the end. This replaces the former O(merges × n) loop (which
|
||
/// rescanned the whole grid per merge: ~250ms at 512×256) with O(n + merges).
|
||
/// Determinism: smallest basin chosen by `(size, id)`, largest neighbor by
|
||
/// `(size, then lowest id)` — both fixed orders.
|
||
fn merge_small_basins(
|
||
mut labels: Vec<i32>,
|
||
w: usize,
|
||
h: usize,
|
||
min_count: usize,
|
||
max_count: usize,
|
||
) -> Vec<i32> {
|
||
use std::collections::BTreeSet;
|
||
let n = w * h;
|
||
let min_frac = 0.02f64; // 2% minimum basin area
|
||
|
||
let max_label = labels.iter().copied().max().unwrap_or(0);
|
||
let nb = (max_label + 1) as usize;
|
||
if nb <= 1 {
|
||
return labels; // single basin — nothing to merge
|
||
}
|
||
|
||
// One pass: basin sizes + adjacency (neighbor labels per basin).
|
||
let mut size = vec![0usize; nb];
|
||
let mut adj: Vec<BTreeSet<i32>> = vec![BTreeSet::new(); nb];
|
||
for r in 0..h {
|
||
for c in 0..w {
|
||
let l = labels[r * w + c];
|
||
size[l as usize] += 1;
|
||
for &(dr, dc) in &D8 {
|
||
let nr = r as i32 + dr;
|
||
let nc = (c as i32 + dc).rem_euclid(w as i32) as usize;
|
||
if nr >= 0 && nr < h as i32 {
|
||
let nl = labels[nr as usize * w + nc];
|
||
if nl != l {
|
||
adj[l as usize].insert(nl);
|
||
}
|
||
}
|
||
}
|
||
}
|
||
}
|
||
|
||
let mut parent: Vec<i32> = (0..nb as i32).collect();
|
||
let mut active: BTreeSet<i32> = (0..nb as i32).collect();
|
||
|
||
while active.len() > min_count {
|
||
// Smallest active basin (tie → lowest id; BTreeSet iterates ascending).
|
||
let smallest = *active
|
||
.iter()
|
||
.min_by_key(|&&b| (size[b as usize], b))
|
||
.unwrap();
|
||
let smallest_size = size[smallest as usize];
|
||
if active.len() <= max_count && smallest_size as f64 / n as f64 >= min_frac {
|
||
break;
|
||
}
|
||
|
||
// Largest active neighbor (tie → lowest id).
|
||
let mut best: i32 = -1;
|
||
let mut best_size = 0usize;
|
||
for &nb_lbl in &adj[smallest as usize] {
|
||
let rep = uf_find(&mut parent, nb_lbl);
|
||
if rep == smallest {
|
||
continue;
|
||
}
|
||
let s = size[rep as usize];
|
||
if s > best_size || (s == best_size && (best < 0 || rep < best)) {
|
||
best_size = s;
|
||
best = rep;
|
||
}
|
||
}
|
||
// No neighbor (isolated basin) → merge into the next smallest active.
|
||
let merge_into = if best >= 0 {
|
||
best
|
||
} else {
|
||
match active.iter().find(|&&b| b != smallest) {
|
||
Some(&other) => other,
|
||
None => break,
|
||
}
|
||
};
|
||
|
||
// Union smallest → merge_into; fold size and adjacency.
|
||
parent[smallest as usize] = merge_into;
|
||
size[merge_into as usize] += smallest_size;
|
||
let small_adj = std::mem::take(&mut adj[smallest as usize]);
|
||
for nb_lbl in small_adj {
|
||
let rep = uf_find(&mut parent, nb_lbl);
|
||
if rep != merge_into {
|
||
adj[merge_into as usize].insert(rep);
|
||
}
|
||
}
|
||
active.remove(&smallest);
|
||
}
|
||
|
||
// Resolve every cell to its basin representative (single pass).
|
||
for l in labels.iter_mut() {
|
||
*l = uf_find(&mut parent, *l);
|
||
}
|
||
|
||
// Renumber contiguously from 0.
|
||
let unique: BTreeSet<i32> = labels.iter().copied().collect();
|
||
let remap: std::collections::BTreeMap<i32, i32> = unique
|
||
.iter()
|
||
.enumerate()
|
||
.map(|(new, &old)| (old, new as i32))
|
||
.collect();
|
||
for l in labels.iter_mut() {
|
||
*l = remap[l];
|
||
}
|
||
|
||
labels
|
||
}
|
||
|
||
// ---------------------------------------------------------------------------
|
||
// Step 8: Build DrainageBasin structs
|
||
// ---------------------------------------------------------------------------
|
||
|
||
fn build_basins(labels: &[i32], w: usize, h: usize) -> Vec<DrainageBasin> {
|
||
let n = w * h;
|
||
let mut basin_map: std::collections::BTreeMap<i32, Vec<usize>> =
|
||
std::collections::BTreeMap::new();
|
||
|
||
for (i, &l) in labels.iter().enumerate() {
|
||
basin_map.entry(l).or_default().push(i);
|
||
}
|
||
|
||
let mut basins: Vec<DrainageBasin> = Vec::with_capacity(basin_map.len());
|
||
let mut ids: Vec<i32> = basin_map.keys().copied().collect();
|
||
ids.sort();
|
||
|
||
for basin_id in ids {
|
||
let cells = &basin_map[&basin_id];
|
||
let area_pct = cells.len() as f32 / n as f32;
|
||
|
||
// Outer boundary as an ordered, non-self-crossing contour via Moore-
|
||
// neighbour tracing from the basin's first (row-major) cell. The previous
|
||
// angle-from-centroid sort produced star-shaped, self-crossing polygons for
|
||
// concave basins, which rendered as straight chords across the map (#960).
|
||
let start = cells.iter().copied().min().unwrap_or(0);
|
||
let mut boundary = trace_outer_boundary(labels, w, h, basin_id, start);
|
||
// Decimate to ≤500 points, preserving traversal order (and thus shape).
|
||
if boundary.len() > 500 {
|
||
let step = boundary.len() / 500;
|
||
boundary = boundary.into_iter().step_by(step).collect();
|
||
}
|
||
|
||
basins.push(DrainageBasin {
|
||
basin_id: basin_id as u32,
|
||
boundary,
|
||
area_pct,
|
||
// Default; the cascade sets the real status from dominant_faction
|
||
// after Layer 1 (D-212, #956).
|
||
territorial_status: TerritorialStatus::FrontierUnclaimed,
|
||
});
|
||
}
|
||
|
||
basins
|
||
}
|
||
|
||
/// Trace the outer boundary of the connected component of `basin_id` containing
|
||
/// `start` (a row-major cell index), clockwise, via Moore-neighbour tracing.
|
||
/// Produces an ordered, 8-connected, non-self-crossing perimeter. The grid edge is
|
||
/// treated as background (no x-wrap) — this is for atlas visualization, not flow.
|
||
fn trace_outer_boundary(
|
||
labels: &[i32],
|
||
w: usize,
|
||
h: usize,
|
||
basin_id: i32,
|
||
start: usize,
|
||
) -> Vec<(u16, u16)> {
|
||
// Moore-neighbourhood offsets in clockwise order: N, NE, E, SE, S, SW, W, NW.
|
||
const DIRS: [(i32, i32); 8] = [
|
||
(-1, 0),
|
||
(-1, 1),
|
||
(0, 1),
|
||
(1, 1),
|
||
(1, 0),
|
||
(1, -1),
|
||
(0, -1),
|
||
(-1, -1),
|
||
];
|
||
let is_fg = |r: i32, c: i32| -> bool {
|
||
r >= 0
|
||
&& r < h as i32
|
||
&& c >= 0
|
||
&& c < w as i32
|
||
&& labels[r as usize * w + c as usize] == basin_id
|
||
};
|
||
let dir_index =
|
||
|dr: i32, dc: i32| -> usize { DIRS.iter().position(|&o| o == (dr, dc)).unwrap_or(0) };
|
||
|
||
let sr = (start / w) as i32;
|
||
let sc = (start % w) as i32;
|
||
let s = (sr, sc);
|
||
let mut boundary: Vec<(u16, u16)> = vec![(sr as u16, sc as u16)];
|
||
let mut p = s;
|
||
// Backtrack starts west of `start`: it is the first cell in scan order, so its
|
||
// western neighbour is background. Consecutive Moore neighbours are 8-adjacent,
|
||
// so the new backtrack stays adjacent to the new boundary cell each step.
|
||
let mut b = (sr, sc - 1);
|
||
let max_steps = w * h * 8 + 16;
|
||
for _ in 0..max_steps {
|
||
let b_idx = dir_index(b.0 - p.0, b.1 - p.1);
|
||
let mut prev = b;
|
||
let mut advanced = false;
|
||
for k in 1..=8 {
|
||
let d = (b_idx + k) % 8;
|
||
let cand = (p.0 + DIRS[d].0, p.1 + DIRS[d].1);
|
||
if is_fg(cand.0, cand.1) {
|
||
if cand == s {
|
||
return boundary; // closed the loop (start already at index 0)
|
||
}
|
||
boundary.push((cand.0 as u16, cand.1 as u16));
|
||
b = prev;
|
||
p = cand;
|
||
advanced = true;
|
||
break;
|
||
}
|
||
prev = cand;
|
||
}
|
||
if !advanced {
|
||
break; // isolated cell — no foreground neighbour
|
||
}
|
||
}
|
||
boundary
|
||
}
|
||
|
||
// ---------------------------------------------------------------------------
|
||
// Tests
|
||
// ---------------------------------------------------------------------------
|
||
|
||
#[cfg(test)]
|
||
mod tests {
|
||
use super::*;
|
||
|
||
fn flat_grid(w: u32, h: u32, val: f32) -> Vec<f32> {
|
||
vec![val; (w * h) as usize]
|
||
}
|
||
|
||
fn slope_grid(w: u32, h: u32) -> Vec<f32> {
|
||
let n = (w * h) as usize;
|
||
(0..n)
|
||
.map(|i| {
|
||
let r = i / w as usize;
|
||
let c = i % w as usize;
|
||
// Slope: higher in top-left, drains toward bottom-right.
|
||
1.0 - (r as f32 / h as f32 * 0.5 + c as f32 / w as f32 * 0.5)
|
||
})
|
||
.collect()
|
||
}
|
||
|
||
#[test]
|
||
fn trace_outer_boundary_is_an_adjacent_contour() {
|
||
// Concave (L-shaped) basin (id=1) in a 4x4 grid; -1 is background.
|
||
let labels: Vec<i32> = vec![
|
||
1, 1, -1, -1, //
|
||
1, 1, -1, -1, //
|
||
1, 1, 1, 1, //
|
||
1, 1, 1, 1, //
|
||
];
|
||
let boundary = trace_outer_boundary(&labels, 4, 4, 1, 0);
|
||
assert!(
|
||
boundary.len() >= 8,
|
||
"expected a real perimeter, got {boundary:?}"
|
||
);
|
||
// The defining property the angle-sort violated: consecutive boundary
|
||
// points are 8-adjacent (a genuine contour walk, not crossing chords).
|
||
for w in boundary.windows(2) {
|
||
let dr = (w[0].0 as i32 - w[1].0 as i32).abs();
|
||
let dc = (w[0].1 as i32 - w[1].1 as i32).abs();
|
||
assert!(
|
||
dr <= 1 && dc <= 1 && dr + dc > 0,
|
||
"non-adjacent step {:?} -> {:?}",
|
||
w[0],
|
||
w[1]
|
||
);
|
||
}
|
||
// Deterministic (D-010): same input → same trace.
|
||
assert_eq!(boundary, trace_outer_boundary(&labels, 4, 4, 1, 0));
|
||
}
|
||
|
||
#[test]
|
||
fn flat_grid_produces_single_basin() {
|
||
let elev = flat_grid(16, 8, 0.5);
|
||
let result = analyze(&elev, 16, 8, 0.3);
|
||
// Flat world → no pour points → single basin
|
||
assert_eq!(result.drainage_basins.len(), 1);
|
||
assert!((result.drainage_basins[0].area_pct - 1.0).abs() < 0.01);
|
||
}
|
||
|
||
#[test]
|
||
fn slope_grid_has_no_river_cells_below_threshold_by_default() {
|
||
// Small 8×4 grid: max flow_accum ≤ 32, below RIVER_THRESHOLD (200).
|
||
let elev = slope_grid(8, 4);
|
||
let result = analyze(&elev, 8, 4, 0.3);
|
||
// River cells may be empty on this tiny grid — that is acceptable.
|
||
// What matters: no panic and basin count ≥ 1.
|
||
assert!(!result.drainage_basins.is_empty());
|
||
}
|
||
|
||
#[test]
|
||
fn large_grid_river_cells_nonempty() {
|
||
// 512×256: max flow accumulation ~131K >> RIVER_THRESHOLD.
|
||
let elev = slope_grid(512, 256);
|
||
let result = analyze(&elev, 512, 256, 0.3);
|
||
assert!(
|
||
!result.river_network.river_cells.is_empty(),
|
||
"Expected river cells on a large sloped grid"
|
||
);
|
||
}
|
||
|
||
#[test]
|
||
fn basin_area_pcts_sum_to_one() {
|
||
let elev = slope_grid(64, 32);
|
||
let result = analyze(&elev, 64, 32, 0.3);
|
||
let total: f32 = result.drainage_basins.iter().map(|b| b.area_pct).sum();
|
||
assert!(
|
||
(total - 1.0).abs() < 0.01,
|
||
"Basin area fractions must sum to 1, got {}",
|
||
total
|
||
);
|
||
}
|
||
|
||
#[test]
|
||
fn basin_count_within_target_range() {
|
||
let elev = slope_grid(128, 64);
|
||
let result = analyze(&elev, 128, 64, 0.3);
|
||
let n = result.drainage_basins.len();
|
||
assert!(
|
||
(1..=12).contains(&n),
|
||
"Basin count {} out of expected range [1, 12]",
|
||
n
|
||
);
|
||
}
|
||
|
||
#[test]
|
||
fn determinism() {
|
||
// Running analyze twice on the same input must produce identical results.
|
||
let elev = slope_grid(64, 32);
|
||
let r1 = analyze(&elev, 64, 32, 0.3);
|
||
let r2 = analyze(&elev, 64, 32, 0.3);
|
||
assert_eq!(
|
||
r1.river_network.river_cells, r2.river_network.river_cells,
|
||
"River cells must be deterministic"
|
||
);
|
||
assert_eq!(
|
||
r1.drainage_basins.len(),
|
||
r2.drainage_basins.len(),
|
||
"Basin count must be deterministic"
|
||
);
|
||
}
|
||
|
||
#[test]
|
||
fn flow_accumulation_deterministic_and_clamped() {
|
||
// flow_accumulation/max_accumulation are the D-209 strength denominator —
|
||
// a silent drift corrupts every attractor strength. Lock them down.
|
||
let elev = slope_grid(64, 32);
|
||
let r1 = analyze(&elev, 64, 32, 0.3);
|
||
let r2 = analyze(&elev, 64, 32, 0.3);
|
||
assert_eq!(r1.flow_accumulation, r2.flow_accumulation);
|
||
assert_eq!(r1.max_accumulation, r2.max_accumulation);
|
||
assert!(
|
||
r1.max_accumulation >= 1,
|
||
"max_accumulation must be clamped ≥ 1"
|
||
);
|
||
// Flat / all-ocean world: still well-defined (no division by zero).
|
||
let flat = flat_grid(16, 8, 0.5);
|
||
let rf = analyze(&flat, 16, 8, 0.9); // sea_level above all terrain
|
||
assert!(rf.max_accumulation >= 1);
|
||
}
|
||
|
||
#[test]
|
||
fn isolated_basins_no_panic() {
|
||
// Two land patches split by an ocean band (rows 3-4 below sea level):
|
||
// exercises basin labeling/merge on a disconnected world.
|
||
let (w, h) = (32usize, 8usize);
|
||
let mut elev = vec![0.1f32; w * h]; // ocean everywhere
|
||
for r in [0, 1, 2, 5, 6, 7] {
|
||
for c in 0..w {
|
||
// two raised land bands, sloped so they drain internally
|
||
elev[r * w + c] = 0.5 + (c as f32 / w as f32) * 0.3;
|
||
}
|
||
}
|
||
let res = analyze(&elev, w as u32, h as u32, 0.3);
|
||
let n = res.drainage_basins.len();
|
||
assert!((1..=12).contains(&n), "basin count {n} out of range");
|
||
}
|
||
}
|