settled-reach/server/src/atlas/drainage.rs

//! 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, RIVER_DOWNSTREAM_EDGE_DRAIN, RIVER_DOWNSTREAM_MOUTH,
};
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
];

/// `(dr, dc)` for D8 direction index `k` (0-7) — the same fixed priority-order
/// table [`extract_river_network`]/`flow_direction` use internally, exposed
/// `pub(crate)` so downstream consumers of [`crate::atlas::body_world_state::
/// RiverNetwork::river_downstream`] (T-1170's river course inventor) can walk
/// a river cell's D8 pointer without duplicating the table. Panics on an
/// out-of-range index — callers must check against the
/// `RIVER_DOWNSTREAM_MOUTH`/`RIVER_DOWNSTREAM_EDGE_DRAIN`/
/// `RIVER_DOWNSTREAM_TERMINAL` sentinels (values ≥ 8) before calling this.
pub(crate) fn d8_offset(k: u8) -> (i32, i32) {
    D8[k as usize]
}

/// 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. River-class banding (T-1156) anchors on its own
    // river-restricted max internally — see `extract_river_network` — not on
    // the grid-wide max computed below, so no dependency ordering between
    // the two is needed.
    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). This is
    // the grid-wide max (includes below-sea-level cells) — distinct from the
    // river-restricted max `extract_river_network` uses for its own T-1156
    // river-class banding.
    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, &deg) 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();

    // River-restricted max accumulation — the ceiling for the T-1156 log-band
    // classifier below. Deliberately NOT the grid-wide `max_accumulation`
    // (DrainageResult's D-209 normalization denominator, which includes
    // below-sea-level ocean cells where accumulation typically peaks, just
    // past a river's mouth): anchoring on that grid-wide value would classify
    // a wet, large-ocean body's actual wettest *river* cell short of trunk,
    // producing an entirely riverless District rung (Araminta's per-rung
    // table shows trunk only at District) on exactly the bodies with the
    // most river to show. Anchoring on the max among cells that passed the
    // `is_river` filter guarantees every body with any river cells has its
    // wettest one classified trunk, by construction — see
    // `classify_river_cell`'s doc comment.
    let river_max_accumulation = (0..n)
        .filter(|&i| is_river[i])
        .map(|i| accum[i])
        .max()
        .unwrap_or(RIVER_THRESHOLD + 1); // unused when river_cells is empty

    // River class per entry of `river_cells`, same order (T-1156 wave 1).
    let river_class: Vec<u8> = (0..n)
        .filter(|&i| is_river[i])
        .map(|i| classify_river_cell(accum[i], river_max_accumulation))
        .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_downstream (T-1170 Ruling 2a/2c/3f): a single pass over
    // `river_cells`, in the SAME order, computing both the D8-downstream
    // sentinel/pointer AND whether this cell is a mouth. `fdir[i]` is already
    // in scope here — capturing it as `river_downstream` costs one extra push
    // per river cell, no new grid pass (Ruling 2a's binding requirement).
    //
    // **Pole-edge exits are NOT mouths (Ruling 3f, binding).** A river cell
    // with no outflow because its D8 walk ran off the grid's top/bottom row
    // is a grid artifact — the equirectangular projection simply stops there,
    // there is no sea. The former code classified this the same as a real
    // sea-adjacent mouth, which rendered double-ring mouth markers in polar
    // ice with no sea in sight (Jeroen's capture question, T-1170 ticket).
    // `RIVER_DOWNSTREAM_EDGE_DRAIN` cells are excluded from `mouths` here;
    // T-1170's course inventor (Ruling 3f) ends their course geometry at the
    // last in-grid station with no mouth flag.
    //
    // A flat-peak interior cell (no outflow, but NOT at a pole row) is D8's
    // other `k < 0` case — vanishingly rare for a cell that also cleared
    // `RIVER_THRESHOLD`, but handled the same as `EDGE_DRAIN` (no downstream
    // neighbor to point at, not a sea mouth) rather than crashing the
    // pointer's "always points somewhere real" contract.
    let mut mouths: Vec<(u16, u16)> = Vec::new();
    let river_downstream: Vec<u8> = (0..n)
        .filter(|&i| is_river[i])
        .map(|i| {
            let r = i / w;
            let c = i % w;
            let k = fdir[i];
            if k < 0 {
                // No outflow at all — edge/flat-peak. Not a mouth (Ruling 3f).
                return RIVER_DOWNSTREAM_EDGE_DRAIN;
            }
            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 {
                // Flow direction points off the polar edge — a grid artifact,
                // not a mouth (Ruling 3f, the pole-edge-drain fix).
                return RIVER_DOWNSTREAM_EDGE_DRAIN;
            }
            if elevation[nr as usize * w + nc] < sea_level {
                // Flows into a sub-sea-level cell — a real mouth.
                mouths.push((r as u16, c as u16));
                return RIVER_DOWNSTREAM_MOUTH;
            }
            k as u8
        })
        .collect();

    RiverNetwork {
        river_cells,
        confluences,
        mouths,
        river_class,
        river_downstream,
    }
}

/// Bin a river cell's flow accumulation into a quantized class (T-1156 wave 1):
/// 0=stream, 1=tributary, 2=trunk. The client filters the ladder rung's river
/// draw by this class (Araminta's per-rung table: Region shows trunk only,
/// District adds tributary, Quarter shows everything) — no new wire field,
/// this is the sole carrier (Tyre's ruling).
///
/// **Binning: log-scaled fraction of the log-range between `RIVER_THRESHOLD`
/// (the accumulation floor below which a cell isn't a river cell at all) and
/// `river_max_accumulation` (the highest flow accumulation among this body's
/// own river cells), split into equal thirds.** Rationale for log rather than
/// linear: flow accumulation grows combinatorially downstream (each
/// confluence roughly sums its tributaries), so a linear split over-populates
/// the trunk band with anything past the halfway point and starves it on
/// modest bodies. Log-scaling spreads the bands evenly across orders of
/// magnitude instead, so a river's headwaters (streams), mid-course
/// tributaries, and lower trunk read as three roughly even bands on both a
/// wet, many-confluence body and a dry, single-channel one.
///
/// **The ceiling must be `river_max_accumulation` (max over cells that pass
/// the `is_river` filter — `accum > RIVER_THRESHOLD && elevation >=
/// sea_level`), never `DrainageResult::max_accumulation` (the grid-wide max
/// used elsewhere for D-209 strength normalization).** Flow accumulation
/// peaks right at a river's mouth, typically on the ocean-side cell just past
/// the coastline — a cell that is *never* a river cell by definition
/// (`is_river` requires `elevation >= sea_level`). Anchoring on the grid-wide
/// max therefore admits a ceiling no river cell can ever reach: on a wet body
/// with a large ocean, where accumulation piles up hardest past the
/// coastline, every actual river cell would land short of trunk and the
/// District rung (trunk-only per Araminta's table) would render riverless —
/// exactly backwards, since that is the body with the most river to show.
/// Anchoring on `river_max_accumulation` instead guarantees, by construction,
/// that a body's own wettest *river* cell — not its wettest cell overall —
/// always lands in the trunk band. Every body with any river cells gets a
/// trunk, scaled to its own wet/dry character, which is what "this body's
/// main river" should mean, and it holds unconditionally (not merely "if the
/// wettest water happens to be fluvial").
///
/// Using a per-body-relative ceiling at all (rather than an absolute multiple
/// of `RIVER_THRESHOLD`, e.g. trunk = accum ≥ 800) is itself deliberate: a
/// body whose single river barely clears the threshold would classify every
/// cell as `stream` under an absolute scheme, reading as "no real river"
/// even though it has exactly one.
///
/// Determinism (D-010/D-208): pure integer/float arithmetic on
/// `(accum, river_max_accumulation)`, no RNG, same body+seed → same class
/// every run. Monotonic by construction: `log` and the linear division into
/// thirds are both non-decreasing in `accum`, so a strictly higher
/// accumulation never produces a strictly lower class.
fn classify_river_cell(accum: i32, river_max_accumulation: i32) -> u8 {
    // Callers only invoke this for cells that passed `is_river` (accum >
    // RIVER_THRESHOLD == 200), and `river_max_accumulation` is the max over
    // that same cell set, so both logs below are well-defined (positive
    // arguments) and `river_max_accumulation > RIVER_THRESHOLD` always holds
    // when there is at least one river cell.
    let floor = (RIVER_THRESHOLD as f64).ln();
    let ceil = (river_max_accumulation as f64)
        .max(RIVER_THRESHOLD as f64 + 1.0)
        .ln();
    let span = (ceil - floor).max(f64::EPSILON);
    let frac = ((accum as f64).ln() - floor) / span;
    let frac = frac.clamp(0.0, 1.0);
    if frac >= 2.0 / 3.0 {
        2 // trunk
    } else if frac >= 1.0 / 3.0 {
        1 // tributary
    } else {
        0 // stream
    }
}

// ---------------------------------------------------------------------------
// 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");
    }

    // -----------------------------------------------------------------------
    // River class (T-1156 wave 1)
    // -----------------------------------------------------------------------

    #[test]
    fn every_river_cell_has_a_class() {
        let elev = slope_grid(512, 256);
        let result = analyze(&elev, 512, 256, 0.3);
        assert_eq!(
            result.river_network.river_cells.len(),
            result.river_network.river_class.len(),
            "river_class must be parallel/aligned with river_cells"
        );
        assert!(
            !result.river_network.river_cells.is_empty(),
            "test grid should produce river cells"
        );
    }

    #[test]
    fn river_class_monotonic_with_accumulation() {
        // A cell with higher accumulation must never have a lower class than
        // a cell with lower accumulation — the core binning contract.
        let elev = slope_grid(512, 256);
        let result = analyze(&elev, 512, 256, 0.3);
        let rn = &result.river_network;
        assert!(!rn.river_cells.is_empty());

        // Recover each river cell's accumulation and pair it with its class.
        let w = 512usize;
        let mut pairs: Vec<(i32, u8)> = rn
            .river_cells
            .iter()
            .zip(rn.river_class.iter())
            .map(|(&(r, c), &class)| {
                let idx = r as usize * w + c as usize;
                (result.flow_accumulation[idx], class)
            })
            .collect();
        pairs.sort_by_key(|&(accum, _)| accum);

        let mut max_class_seen = 0u8;
        for (_, class) in pairs {
            assert!(
                class >= max_class_seen,
                "monotonicity violated: saw class {class} after class {max_class_seen} \
                 in ascending-accumulation order"
            );
            max_class_seen = max_class_seen.max(class);
        }
    }

    #[test]
    fn at_least_one_trunk_cell_when_rivers_exist() {
        let elev = slope_grid(512, 256);
        let result = analyze(&elev, 512, 256, 0.3);
        assert!(!result.river_network.river_cells.is_empty());
        assert!(
            result.river_network.river_class.contains(&2),
            "a body with any rivers must have at least one trunk (class 2) cell — \
             this is the classify_river_cell river_max_accumulation-anchoring guarantee"
        );
    }

    #[test]
    fn at_least_one_trunk_cell_on_a_real_body_with_a_large_ocean() {
        // Regression for the grid-wide-max anchoring bug: GJ1c is exactly the
        // "wet body with a large ocean" shape where flow accumulation peaks
        // past the coastline (a non-river cell), which starved the trunk band
        // when the ceiling was anchored on the grid-wide max instead of the
        // river-restricted max. Same body + downsample as the cascade golden
        // (tests/golden/cascade_layer1.json) — 93 river cells there, so this
        // is a real, non-synthetic exercise of the guarantee.
        use crate::atlas::heightmap::load_heightmap_png;
        let src = std::path::PathBuf::from(env!("CARGO_MANIFEST_DIR"))
            .join("../wiki/star-systems/GJ-1/bodies/GJ1c/heightmap.png");
        let heightmap =
            load_heightmap_png(&src, "GJ1c", 0.3).expect("decode committed GJ1c heightmap");
        let small = heightmap.downsample(256, 128);
        let result = analyze(&small.data, small.width, small.height, small.sea_level);
        assert!(
            !result.river_network.river_cells.is_empty(),
            "GJ1c should have river cells at this downsample"
        );
        assert!(
            result.river_network.river_class.contains(&2),
            "GJ1c's own wettest river cell must classify as trunk — river-restricted \
             anchoring must not be starved by ocean-cell accumulation past the coastline"
        );
    }

    #[test]
    fn river_class_deterministic() {
        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_class, r2.river_network.river_class,
            "river_class must be deterministic (D-010/D-208)"
        );
    }

    #[test]
    fn classify_river_cell_barely_above_threshold_still_gets_a_trunk() {
        // A body whose single river barely clears RIVER_THRESHOLD must still
        // classify its own maximum as trunk — the whole point of anchoring
        // the log-range ceiling at river_max_accumulation instead of an
        // absolute multiple of RIVER_THRESHOLD.
        let river_max_accumulation = RIVER_THRESHOLD + 5;
        assert_eq!(
            classify_river_cell(river_max_accumulation, river_max_accumulation),
            2,
            "the body's own max river-cell accumulation must always classify as trunk"
        );
    }

    #[test]
    fn classify_river_cell_spans_all_three_classes_on_wide_range() {
        // Sanity check on the log-binning: a body with a wide dynamic range
        // (headwater trickles up to a major trunk) should exercise all three
        // classes, not collapse to two.
        let river_max_accumulation = 131_000;
        let low = classify_river_cell(RIVER_THRESHOLD + 1, river_max_accumulation);
        let mid = classify_river_cell(5_000, river_max_accumulation);
        let high = classify_river_cell(river_max_accumulation, river_max_accumulation);
        assert_eq!(low, 0, "just above threshold should be a stream");
        assert_eq!(mid, 1, "mid-range accumulation should be a tributary");
        assert_eq!(
            high, 2,
            "the body's max river-cell accumulation should be trunk"
        );
    }

    // -----------------------------------------------------------------------
    // river_downstream (T-1170 Ruling 2a/2c/3f)
    // -----------------------------------------------------------------------

    #[test]
    fn river_downstream_parallel_to_river_cells() {
        let elev = slope_grid(512, 256);
        let result = analyze(&elev, 512, 256, 0.3);
        let rn = &result.river_network;
        assert_eq!(
            rn.river_cells.len(),
            rn.river_downstream.len(),
            "river_downstream must be parallel/aligned with river_cells"
        );
        assert!(!rn.river_cells.is_empty());
    }

    #[test]
    fn river_downstream_values_are_direction_or_sentinel() {
        // Every entry is either a real D8 index (0-7) or one of the T-1170
        // sentinels (MOUTH=8, EDGE_DRAIN=9); TERMINAL=10 is reserved/unused.
        let elev = slope_grid(512, 256);
        let result = analyze(&elev, 512, 256, 0.3);
        for &v in &result.river_network.river_downstream {
            assert!(
                v <= RIVER_DOWNSTREAM_EDGE_DRAIN,
                "unexpected river_downstream value {v} (round-1 only emits 0-7, MOUTH=8, \
                 EDGE_DRAIN=9 — TERMINAL=10 is reserved and unused)"
            );
        }
    }

    #[test]
    fn river_downstream_mouth_sentinel_matches_mouths_list() {
        // Every river cell whose river_downstream is MOUTH must appear in
        // `mouths`, and every entry of `mouths` must have MOUTH as its
        // river_downstream — the two are the same underlying classification,
        // captured in the same pass (Ruling 2a/3f).
        use crate::atlas::heightmap::load_heightmap_png;
        let src = std::path::PathBuf::from(env!("CARGO_MANIFEST_DIR"))
            .join("../wiki/star-systems/GJ-1/bodies/GJ1c/heightmap.png");
        let heightmap =
            load_heightmap_png(&src, "GJ1c", 0.3).expect("decode committed GJ1c heightmap");
        let small = heightmap.downsample(256, 128);
        let result = analyze(&small.data, small.width, small.height, small.sea_level);
        let rn = &result.river_network;
        assert!(!rn.mouths.is_empty(), "GJ1c should have real sea mouths");

        let mouth_set: std::collections::BTreeSet<(u16, u16)> =
            rn.mouths.iter().copied().collect();
        for (i, &pos) in rn.river_cells.iter().enumerate() {
            let is_mouth_sentinel = rn.river_downstream[i] == RIVER_DOWNSTREAM_MOUTH;
            let is_in_mouths_list = mouth_set.contains(&pos);
            assert_eq!(
                is_mouth_sentinel, is_in_mouths_list,
                "cell {pos:?}: MOUTH sentinel ({is_mouth_sentinel}) must agree with \
                 mouths-list membership ({is_in_mouths_list})"
            );
        }
    }

    #[test]
    fn pole_edge_drains_are_not_mouths() {
        // Ruling 3f, binding: a river cell with no valid D8 outflow is a grid
        // artifact, not a river-meets-sea event.
        //
        // **Why this test targets the flat/no-outflow (`k < 0`) case, not a
        // literal "flow direction points past row 0" case:** `flow_direction`
        // (this file, `fn flow_direction`) bounds-checks every D8 candidate
        // BEFORE comparing drops (`if nr < 0 || nr >= h { continue; }`) — a
        // row-0 cell can therefore never even be ASSIGNED a north-pointing
        // `fdir` in the first place; the off-grid-direction branch in
        // `extract_river_network`'s `river_downstream` computation exists as
        // correct defensive code but is structurally unreachable given this
        // invariant. The real, reachable "pole-edge-drain" case (confirmed
        // against the committed GJ1c golden fixture, which has EDGE_DRAIN
        // cells at several rows including row 0) is `k < 0`: a cell with NO
        // neighbor at a strictly lower elevation — most commonly a flat
        // plateau at the grid's fringe, which the depression-fill/flow
        // algorithm cannot route off of. This fixture constructs exactly
        // that: a perfectly flat plateau at row 0 (identical elevation
        // across the whole top row, so no cell in it has a positive-drop
        // neighbor and `flow_direction` assigns `k=-1` — verified this
        // reproduces before ever reasoning about mouths) that river cells
        // from a converging valley drain into, with no ocean anywhere.
        let (w, h) = (64usize, 64usize);
        let n = w * h;
        let center_col = (w / 2) as f32;
        let elev: Vec<f32> = (0..n)
            .map(|i| {
                let r = i / w;
                let c = (i % w) as f32;
                if r == 0 {
                    return 0.15; // flat plateau — no cell here has a strictly lower neighbor
                }
                // V-shaped valley converging on center_col, sloping down
                // toward row 0 (but never reaching the plateau's own
                // elevation until row 1, so row-1 cells drain INTO the flat
                // row-0 plateau and then have nowhere further to go).
                let dist_from_center = (c - center_col).abs() / center_col;
                let valley = 0.15 + dist_from_center * 0.6;
                let pole_gradient = (r as f32 / h as f32) * 0.25;
                (valley + pole_gradient).clamp(0.0, 1.0)
            })
            .collect();
        let result = analyze(&elev, w as u32, h as u32, 0.0);
        let rn = &result.river_network;
        assert!(
            !rn.river_cells.is_empty(),
            "the converging-valley fixture must clear RIVER_THRESHOLD — if this starts \
             failing, the fixture (not the production code) needs retuning, since an \
             empty river_cells silently no-ops every assertion below"
        );
        assert!(
            rn.river_downstream
                .iter()
                .any(|&v| v == RIVER_DOWNSTREAM_EDGE_DRAIN),
            "expected at least one EDGE_DRAIN-sentinel river cell (the flat-plateau case) \
             on this fixture — downstream values were {:?}",
            rn.river_downstream
        );
        assert!(
            rn.mouths.is_empty(),
            "an all-land, pole-draining world must have zero mouths — got {:?}",
            rn.mouths
        );
        assert!(
            !rn.river_downstream
                .iter()
                .any(|&v| v == RIVER_DOWNSTREAM_MOUTH),
            "an all-land world must never emit a MOUTH sentinel"
        );
    }

    #[test]
    fn river_downstream_deterministic() {
        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_downstream, r2.river_network.river_downstream,
            "river_downstream must be deterministic (D-010/D-208)"
        );
    }
}