jpmschweitzerandClaude Fable 5
4ef79b6b45
Hoshe finding 1: dropping HydrologySample.elevation was PROVEN unsafe (TerrainAnalysis.elev_pct is a rank percentile, not raw elevation; HydrologyResult carries no elevation) — the copy stays, and the truth moves into the docs instead: gen_queue's TerrainAnalysisCache sizing comment corrected to real 512x256 working-grid numbers (~1.57 -> ~2.62 MB/entry, capacity-8 worst case ~21 MB), clone-on-hit cost documented, byte-safety proof recorded on HydrologySample itself. Arc<TerrainAnalysis> follow-up filed as T-1187. Hoshe finding 2: the endorheic-split test now asserts the bowl basin's BasinOutcome actually diverges (Endorheic at moisture 0, Overflow at 100, straddling ENDORHEIC_MOISTURE_CEILING=60) plus basin-count sanity — mutation-verified by stubbing is_endorheic and watching it fail. Doc/test-only round: goldens byte-unchanged, full suite 2114 green. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
890 lines
32 KiB
Rust
890 lines
32 KiB
Rust
//! Geographic feature tag extraction — Layer 1 (D-209).
|
||
//!
|
||
//! After D8 drainage analysis (D-208), this module extracts the 7
|
||
//! `AttractorType` tags from the heightmap + river network. Each attractor has
|
||
//! a pixel position and a `strength` (integer 0–100) derived from local terrain
|
||
//! quality — computed in floating point, then quantized to an integer at the
|
||
//! extraction boundary so all downstream decisions stay integer-deterministic
|
||
//! (D-010).
|
||
//!
|
||
//! **D-209 / D-223 reconciliation:** D-209 reads ocean/lake polygons from
|
||
//! `markers.json`, but D-223 reduced markers to a names-only pool — those
|
||
//! polygons no longer exist. Coast and lake cells are therefore derived from
|
||
//! the heightmap itself: ocean = the largest connected below-sea-level water
|
||
//! body; lakes = smaller enclosed below-sea-level bodies.
|
||
//!
|
||
//! **Determinism (D-010 #4):** all collections iterate in sorted/row-major
|
||
//! order; the final attractor list is sorted by `(attractor_type, row, col)`.
|
||
//! No `HashMap`/`HashSet` iteration. `strength` (integer 0–100) is not part of
|
||
//! that ordering here, so float quantization can't perturb the sort; consumers
|
||
//! that rank by strength (e.g. `layer1::attach_feature_names`) do so on the
|
||
//! integer value.
|
||
|
||
use std::collections::{BTreeMap, VecDeque};
|
||
|
||
use crate::atlas::drainage::DrainageResult;
|
||
use crate::atlas::heightmap::BodyHeightmap;
|
||
use crate::simulation::generator::AttractorType;
|
||
|
||
/// 8-neighbor offsets (dr, dc). Columns wrap horizontally (equirectangular
|
||
/// globe); rows are bounds-clamped at the poles. Matches `drainage::D8`.
|
||
const NB8: [(i32, i32); 8] = [
|
||
(-1, 0),
|
||
(1, 0),
|
||
(0, 1),
|
||
(0, -1),
|
||
(-1, 1),
|
||
(-1, -1),
|
||
(1, 1),
|
||
(1, -1),
|
||
];
|
||
|
||
/// Minimum spacing (cells) between attractors of an areal type, so coastlines /
|
||
/// valleys / plains yield a sparse, placement-friendly set rather than one
|
||
/// attractor per pixel.
|
||
const MIN_SPACING: i32 = 12;
|
||
|
||
/// Hard cap on attractors per body (keeps the #955 matching tractable).
|
||
const MAX_ATTRACTORS: usize = 256;
|
||
|
||
/// A feature before sub-biome classification: position, type, strength.
|
||
/// `layer1` enriches these into `GeographicAttractor`s.
|
||
#[derive(Debug, Clone, Copy, PartialEq)]
|
||
pub struct RawAttractor {
|
||
pub row: u16,
|
||
pub col: u16,
|
||
pub attractor_type: AttractorType,
|
||
/// Strength on a 0–100 integer scale (100 = strongest). Quantized here from
|
||
/// the f32 flow-accumulation ratio — the single f32→integer boundary, after
|
||
/// which every ranking/scoring decision is integer (D-010, #955).
|
||
pub strength: i32,
|
||
}
|
||
|
||
/// Precomputed per-cell terrain fields, shared by feature extraction (D-209)
|
||
/// and sub-biome classification (D-210) so neither recomputes them.
|
||
#[derive(Debug, Clone)]
|
||
pub struct TerrainAnalysis {
|
||
pub w: usize,
|
||
pub h: usize,
|
||
/// `elev < sea_level` (any submerged cell).
|
||
pub ocean_mask: Vec<bool>,
|
||
/// Submerged cells not part of the largest water body (enclosed lakes/seas).
|
||
pub lake_mask: Vec<bool>,
|
||
/// Chebyshev distance (cells) to the nearest ocean cell or river mouth,
|
||
/// capped at `WATER_DIST_CAP`. Moisture proxy for habitability/sub-biome.
|
||
pub water_dist: Vec<u16>,
|
||
/// Local slope proxy in degrees: `atan(max |Δelev| over 8 neighbors)`.
|
||
/// Elevation is normalized [0,1]; this is a relative steepness measure.
|
||
pub slope_deg: Vec<f32>,
|
||
/// Elevation percentile [0,1] among land cells (ocean cells = 0.0).
|
||
pub elev_pct: Vec<f32>,
|
||
/// Settled-equilibrium hydrology sourcing (T-1184, D-227 amendment (4) /
|
||
/// D-255(f) seed-chaining mechanism B). `None` when hydrology hasn't been
|
||
/// solved for this analysis (e.g. every pre-T-1184 call site still using
|
||
/// bare [`TerrainAnalysis::analyze`], and every unit test that constructs
|
||
/// a `TerrainAnalysis` directly without going through the hydrology-aware
|
||
/// entry point) — callers MUST treat `None` as "fall through to the
|
||
/// `ocean_fraction_q` heuristic", never as an error. `Some` when
|
||
/// [`TerrainAnalysis::with_hydrology`] populated it from a real
|
||
/// [`crate::atlas::hydrology_equilibrium::HydrologyResult`].
|
||
pub hydrology: Option<HydrologySample>,
|
||
}
|
||
|
||
/// The two continuous working-grid fields `derive_morphology_zone`'s lake
|
||
/// sourcing bilinearly samples (T-1184) — never a discrete basin-membership
|
||
/// lookup (that gives blocky, non-refining lake edges, the exact D-166
|
||
/// magnified-composite artifact this design avoids; see D-227 amendment (4)
|
||
/// / D-255(f) mechanism B). Both fields are row-major, `w × h`, in the SAME
|
||
/// `[0.0, 1.0]` normalized domain the raw heightmap and `sea_level` already
|
||
/// share — so a bilinear sample of one is directly comparable to a bilinear
|
||
/// sample of the other, no rescaling at the call site.
|
||
///
|
||
/// **Size + clone cost (PR #200 review, Hoshe finding 1):** two `Vec<f32>` at
|
||
/// the real 512×256 working grid = ~1.05 MB/entry, added on top of
|
||
/// `TerrainAnalysis`'s pre-existing ~1.57 MB of dense fields (~2.62 MB total,
|
||
/// ×1.67 growth, not quite a doubling) — see the corrected sizing comment on
|
||
/// `GenWorkItem::DeriveWindow` (`gen_queue.rs`) for the full accounting and
|
||
/// the `TerrainAnalysisCache` cache-HIT clone-cost note (every hit
|
||
/// deep-copies both these `Vec`s, not just the first miss/insert).
|
||
///
|
||
/// **`elevation` is a deliberate, provably-necessary redundant copy, not an
|
||
/// oversight.** `TerrainAnalysis` has no OTHER field that retains the raw
|
||
/// `[0,1]` heightmap: `elev_pct` is a RANK PERCENTILE (`rank(elev[i]) /
|
||
/// (land_cell_count - 1)`, `compute_elev_percentile`'s own doc/impl) —
|
||
/// mathematically a different quantity from absolute elevation, and NOT
|
||
/// safe to compare against `filled` (two cells at different true elevations
|
||
/// can share adjacent ranks; ocean cells are forced to `0.0` regardless of
|
||
/// their real depth). `HydrologyResult` itself carries no elevation field
|
||
/// either (`hydrology_equilibrium.rs`: `basins`, `filled_scaled`,
|
||
/// `channel_depth_scaled`, `cliff_edge` — no `original`/`elevation` member).
|
||
/// So there is no existing bit-identical grid this field could point at
|
||
/// instead — carrying its own copy is the only byte-safe option today.
|
||
#[derive(Debug, Clone)]
|
||
pub struct HydrologySample {
|
||
/// The original (unfilled) heightmap elevation, `[0.0, 1.0]`. Not stored
|
||
/// anywhere else on `TerrainAnalysis` (`elev_pct` is a land-cell RANK
|
||
/// percentile, a different quantity — see its own doc) — this is the
|
||
/// literal `hm.data` the solver's `original` array was built from,
|
||
/// carried alongside `filled` so both halves of the lake comparison
|
||
/// sample from the identical grid at the identical resolution.
|
||
pub elevation: Vec<f32>,
|
||
/// `HydrologyResult.filled_scaled`, rescaled back from the solver's
|
||
/// `i64`-scaled integer domain to `[0.0, 1.0]` (dividing by the same
|
||
/// `ELEV_SCALE` the solver used to go the other way) — the settled
|
||
/// water-surface height at every working-grid cell (equal to
|
||
/// `elevation` wherever no lake exists).
|
||
pub filled: Vec<f32>,
|
||
}
|
||
|
||
const WATER_DIST_CAP: u16 = 255;
|
||
|
||
#[inline]
|
||
fn idx(r: usize, c: usize, w: usize) -> usize {
|
||
r * w + c
|
||
}
|
||
|
||
#[inline]
|
||
fn wrap_col(c: i32, w: i32) -> usize {
|
||
c.rem_euclid(w) as usize
|
||
}
|
||
|
||
impl TerrainAnalysis {
|
||
/// Compute all shared terrain fields for a body. O(w·h).
|
||
pub fn analyze(hm: &BodyHeightmap, drainage: &DrainageResult) -> TerrainAnalysis {
|
||
let w = hm.width as usize;
|
||
let h = hm.height as usize;
|
||
let n = w * h;
|
||
let elev = &hm.data;
|
||
let sea = hm.sea_level;
|
||
|
||
let ocean_mask: Vec<bool> = (0..n).map(|i| elev[i] < sea).collect();
|
||
let lake_mask = compute_lake_mask(&ocean_mask, w, h);
|
||
let water_dist = compute_water_dist(&ocean_mask, &drainage.river_network.mouths, w, h);
|
||
let slope_deg = compute_slope(elev, w, h);
|
||
let elev_pct = compute_elev_percentile(elev, &ocean_mask, w, h);
|
||
|
||
TerrainAnalysis {
|
||
w,
|
||
h,
|
||
ocean_mask,
|
||
lake_mask,
|
||
water_dist,
|
||
slope_deg,
|
||
elev_pct,
|
||
hydrology: None,
|
||
}
|
||
}
|
||
|
||
/// Populate the settled-hydrology sourcing fields (T-1184, D-227
|
||
/// amendment (4) / D-255(f) mechanism B) from a solved
|
||
/// [`crate::atlas::hydrology_equilibrium::HydrologyResult`].
|
||
///
|
||
/// Builder-style (consumes and returns `self`) rather than a constructor
|
||
/// parameter on [`TerrainAnalysis::analyze`] — `analyze` has ~20 call
|
||
/// sites across production code and tests that have no hydrology input
|
||
/// (and, per D-227, don't need one: hydrology sourcing is a lake-specific
|
||
/// refinement, not a precondition for every other terrain field this
|
||
/// struct carries). Keeping `analyze`'s signature untouched means every
|
||
/// existing caller keeps working byte-identically; only the two
|
||
/// production sites that actually solve hydrology
|
||
/// (`layer1::run_layer1`, `gen_queue::TerrainAnalysisCache::get_or_derive`)
|
||
/// opt in by chaining this call.
|
||
///
|
||
/// Panics if `result`'s grids aren't `self.w * self.h` cells — a
|
||
/// programmer error (mismatched working-grid resolution between the
|
||
/// heightmap this `TerrainAnalysis` was built from and the elevation grid
|
||
/// `solve()` was called on), never a legitimate runtime state.
|
||
pub fn with_hydrology(
|
||
mut self,
|
||
elevation: &[f32],
|
||
result: &crate::atlas::hydrology_equilibrium::HydrologyResult,
|
||
) -> TerrainAnalysis {
|
||
let n = self.w * self.h;
|
||
assert_eq!(
|
||
elevation.len(),
|
||
n,
|
||
"with_hydrology: elevation grid size does not match TerrainAnalysis dims"
|
||
);
|
||
assert_eq!(
|
||
result.filled_scaled.len(),
|
||
n,
|
||
"with_hydrology: HydrologyResult grid size does not match TerrainAnalysis dims"
|
||
);
|
||
let filled: Vec<f32> = result
|
||
.filled_scaled
|
||
.iter()
|
||
.map(|&s| crate::atlas::hydrology_equilibrium::scaled_to_fraction(s))
|
||
.collect();
|
||
self.hydrology = Some(HydrologySample {
|
||
elevation: elevation.to_vec(),
|
||
filled,
|
||
});
|
||
self
|
||
}
|
||
|
||
#[inline]
|
||
pub fn is_ocean(&self, r: usize, c: usize) -> bool {
|
||
self.ocean_mask[idx(r, c, self.w)]
|
||
}
|
||
|
||
/// Compass bearing toward the nearest water from cell `(r, c)`, quantized to
|
||
/// 8 octants (0=N, 45=NE … 315=NW); `360` = "no water in range" (#957, D-234).
|
||
///
|
||
/// Reads the `water_dist` field's local gradient — the 8-neighbour with the
|
||
/// smallest distance-to-water points toward water. Integer-only (no `atan2`)
|
||
/// for D-010 determinism. Returns `360` when the cell is itself water or no
|
||
/// neighbour is closer to water (flat/inland).
|
||
pub fn water_bearing(&self, r: usize, c: usize) -> u16 {
|
||
let here = self.water_dist[idx(r, c, self.w)];
|
||
if here == 0 || here >= WATER_DIST_CAP {
|
||
return NO_WATER_BEARING; // on water, or no water within range
|
||
}
|
||
let mut best = here;
|
||
let mut bdir = (0i32, 0i32);
|
||
for &(dr, dc) in &NB8 {
|
||
let nr = r as i32 + dr;
|
||
if nr < 0 || nr >= self.h as i32 {
|
||
continue;
|
||
}
|
||
let nc = wrap_col(c as i32 + dc, self.w as i32);
|
||
let nd = self.water_dist[idx(nr as usize, nc, self.w)];
|
||
if nd < best {
|
||
best = nd;
|
||
bdir = (dr, dc);
|
||
}
|
||
}
|
||
if bdir == (0, 0) {
|
||
NO_WATER_BEARING
|
||
} else {
|
||
octant_bearing(bdir.0, bdir.1)
|
||
}
|
||
}
|
||
}
|
||
|
||
/// Sentinel for [`TerrainAnalysis::water_bearing`] meaning "no water direction".
|
||
pub const NO_WATER_BEARING: u16 = 360;
|
||
|
||
/// Quantize a (Δrow, Δcol) step to a compass octant bearing (0=N … 315=NW).
|
||
/// `Δrow < 0` is north (rows increase downward). Integer-only (D-010).
|
||
fn octant_bearing(drow: i32, dcol: i32) -> u16 {
|
||
let (ar, ac) = (drow.abs(), dcol.abs());
|
||
let north = drow < 0;
|
||
let east = dcol > 0;
|
||
if ar >= ac * 2 {
|
||
if north {
|
||
0
|
||
} else {
|
||
180
|
||
}
|
||
} else if ac >= ar * 2 {
|
||
if east {
|
||
90
|
||
} else {
|
||
270
|
||
}
|
||
} else {
|
||
match (north, east) {
|
||
(true, true) => 45,
|
||
(true, false) => 315,
|
||
(false, true) => 135,
|
||
(false, false) => 225,
|
||
}
|
||
}
|
||
}
|
||
|
||
/// Largest connected below-sea-level component = ocean; all others = lakes.
|
||
/// Deterministic: BFS seeds scanned row-major; ties broken by lowest cell index.
|
||
fn compute_lake_mask(ocean_mask: &[bool], w: usize, h: usize) -> Vec<bool> {
|
||
let n = w * h;
|
||
let mut comp = vec![-1i32; n];
|
||
let mut comp_sizes: Vec<usize> = Vec::new();
|
||
let mut next_comp = 0i32;
|
||
|
||
for start in 0..n {
|
||
if !ocean_mask[start] || comp[start] >= 0 {
|
||
continue;
|
||
}
|
||
// Flood fill this component (row-major BFS = deterministic).
|
||
let mut size = 0usize;
|
||
let mut q = VecDeque::new();
|
||
comp[start] = next_comp;
|
||
q.push_back(start);
|
||
while let Some(cur) = q.pop_front() {
|
||
size += 1;
|
||
let (r, c) = (cur / w, cur % w);
|
||
for &(dr, dc) in &NB8 {
|
||
let nr = r as i32 + dr;
|
||
if nr < 0 || nr >= h as i32 {
|
||
continue;
|
||
}
|
||
let nc = wrap_col(c as i32 + dc, w as i32);
|
||
let ni = idx(nr as usize, nc, w);
|
||
if ocean_mask[ni] && comp[ni] < 0 {
|
||
comp[ni] = next_comp;
|
||
q.push_back(ni);
|
||
}
|
||
}
|
||
}
|
||
comp_sizes.push(size);
|
||
next_comp += 1;
|
||
}
|
||
|
||
if comp_sizes.is_empty() {
|
||
return vec![false; n]; // no water at all
|
||
}
|
||
|
||
// Largest component (tie → lowest comp id, which is the earliest row-major).
|
||
let mut ocean_comp = 0i32;
|
||
let mut best = 0usize;
|
||
for (cid, &sz) in comp_sizes.iter().enumerate() {
|
||
if sz > best {
|
||
best = sz;
|
||
ocean_comp = cid as i32;
|
||
}
|
||
}
|
||
|
||
// Lakes = submerged cells in any non-ocean component.
|
||
(0..n)
|
||
.map(|i| comp[i] >= 0 && comp[i] != ocean_comp)
|
||
.collect()
|
||
}
|
||
|
||
/// Multi-source BFS Chebyshev distance to nearest ocean cell or river mouth.
|
||
fn compute_water_dist(ocean_mask: &[bool], mouths: &[(u16, u16)], w: usize, h: usize) -> Vec<u16> {
|
||
let n = w * h;
|
||
let mut dist = vec![WATER_DIST_CAP; n];
|
||
let mut q = VecDeque::new();
|
||
// Seeds in row-major order for determinism.
|
||
for i in 0..n {
|
||
if ocean_mask[i] {
|
||
dist[i] = 0;
|
||
q.push_back(i);
|
||
}
|
||
}
|
||
for &(mr, mc) in mouths {
|
||
let i = idx(mr as usize, mc as usize, w);
|
||
if dist[i] != 0 {
|
||
dist[i] = 0;
|
||
q.push_back(i);
|
||
}
|
||
}
|
||
while let Some(cur) = q.pop_front() {
|
||
let d = dist[cur];
|
||
if d >= WATER_DIST_CAP {
|
||
continue;
|
||
}
|
||
let (r, c) = (cur / w, cur % w);
|
||
for &(dr, dc) in &NB8 {
|
||
let nr = r as i32 + dr;
|
||
if nr < 0 || nr >= h as i32 {
|
||
continue;
|
||
}
|
||
let nc = wrap_col(c as i32 + dc, w as i32);
|
||
let ni = idx(nr as usize, nc, w);
|
||
if dist[ni] > d + 1 {
|
||
dist[ni] = d + 1;
|
||
q.push_back(ni);
|
||
}
|
||
}
|
||
}
|
||
dist
|
||
}
|
||
|
||
/// Local slope proxy: `atan(max |Δelev| to 8 neighbors)` in degrees.
|
||
fn compute_slope(elev: &[f32], w: usize, h: usize) -> Vec<f32> {
|
||
let n = w * h;
|
||
let mut slope = vec![0.0f32; n];
|
||
for r in 0..h {
|
||
for c in 0..w {
|
||
let i = idx(r, c, w);
|
||
let e = elev[i];
|
||
let mut max_grad = 0.0f32;
|
||
for &(dr, dc) in &NB8 {
|
||
let nr = r as i32 + dr;
|
||
if nr < 0 || nr >= h as i32 {
|
||
continue;
|
||
}
|
||
let nc = wrap_col(c as i32 + dc, w as i32);
|
||
let g = (e - elev[idx(nr as usize, nc, w)]).abs();
|
||
if g > max_grad {
|
||
max_grad = g;
|
||
}
|
||
}
|
||
slope[i] = max_grad.atan().to_degrees();
|
||
}
|
||
}
|
||
slope
|
||
}
|
||
|
||
/// Elevation percentile [0,1] among land cells; ocean cells get 0.0.
|
||
fn compute_elev_percentile(elev: &[f32], ocean_mask: &[bool], w: usize, h: usize) -> Vec<f32> {
|
||
let n = w * h;
|
||
// (scaled_elev, idx) for land cells; integer key for deterministic sort.
|
||
let mut land: Vec<(i64, usize)> = (0..n)
|
||
.filter(|&i| !ocean_mask[i])
|
||
.map(|i| ((elev[i] as f64 * 1_000_000.0) as i64, i))
|
||
.collect();
|
||
land.sort_unstable_by(|a, b| a.0.cmp(&b.0).then(a.1.cmp(&b.1)));
|
||
let mut pct = vec![0.0f32; n];
|
||
let m = land.len();
|
||
if m <= 1 {
|
||
for &(_, i) in &land {
|
||
pct[i] = 0.5;
|
||
}
|
||
return pct;
|
||
}
|
||
for (rank, &(_, i)) in land.iter().enumerate() {
|
||
pct[i] = rank as f32 / (m - 1) as f32;
|
||
}
|
||
pct
|
||
}
|
||
|
||
/// Habitability score [0,1] from elevation band, flatness, and moisture.
|
||
/// Used by `ValleyFloor`/`PlainCenter` strengths and the D-210 classifier.
|
||
pub fn habitability(elev_pct: f32, slope_deg: f32, water_dist: u16) -> f32 {
|
||
let elev_score = (1.0 - (elev_pct - 0.35).abs() / 0.65).clamp(0.0, 1.0);
|
||
let flat_score = (1.0 - slope_deg / 15.0).clamp(0.0, 1.0);
|
||
let water_score = (1.0 - water_dist as f32 / 40.0).clamp(0.0, 1.0);
|
||
(0.4 * elev_score + 0.3 * flat_score + 0.3 * water_score).clamp(0.0, 1.0)
|
||
}
|
||
|
||
/// Extract the 7 D-209 attractor tags. Returns raw attractors (no sub-biome),
|
||
/// sorted by `(attractor_type, row, col)` for determinism. Higher-priority
|
||
/// types claim their cells first so a cell is tagged at most once.
|
||
pub fn extract_attractors(
|
||
hm: &BodyHeightmap,
|
||
drainage: &DrainageResult,
|
||
ta: &TerrainAnalysis,
|
||
) -> Vec<RawAttractor> {
|
||
let w = hm.width as usize;
|
||
let h = hm.height as usize;
|
||
let elev = &hm.data;
|
||
let accum = &drainage.flow_accumulation;
|
||
let max_accum = drainage.max_accumulation.max(1) as f32;
|
||
|
||
let mut claimed = vec![false; w * h];
|
||
let mut out: Vec<RawAttractor> = Vec::new();
|
||
|
||
// O(1) river-cell membership (avoids a binary_search per valley candidate).
|
||
let mut river_mask = vec![false; w * h];
|
||
for &(r, c) in &drainage.river_network.river_cells {
|
||
river_mask[idx(r as usize, c as usize, w)] = true;
|
||
}
|
||
|
||
let claim = |out: &mut Vec<RawAttractor>,
|
||
claimed: &mut [bool],
|
||
r: usize,
|
||
c: usize,
|
||
at: AttractorType,
|
||
strength: f32| {
|
||
let i = idx(r, c, w);
|
||
if claimed[i] {
|
||
return;
|
||
}
|
||
claimed[i] = true;
|
||
out.push(RawAttractor {
|
||
row: r as u16,
|
||
col: c as u16,
|
||
attractor_type: at,
|
||
// The single f32→integer boundary: quantize the 0.0–1.0 ratio to 0–100.
|
||
strength: (strength.clamp(0.0, 1.0) * 100.0).round() as i32,
|
||
});
|
||
};
|
||
|
||
// 1. RiverMouth — strength = accum / max_accum.
|
||
for &(r, c) in &drainage.river_network.mouths {
|
||
let i = idx(r as usize, c as usize, w);
|
||
let s = accum[i] as f32 / max_accum;
|
||
claim(
|
||
&mut out,
|
||
&mut claimed,
|
||
r as usize,
|
||
c as usize,
|
||
AttractorType::RiverMouth,
|
||
s,
|
||
);
|
||
}
|
||
|
||
// 2. RiverCrossing — confluences, strength = accum / max_accum * 0.7.
|
||
for &(r, c) in &drainage.river_network.confluences {
|
||
let i = idx(r as usize, c as usize, w);
|
||
let s = accum[i] as f32 / max_accum * 0.7;
|
||
claim(
|
||
&mut out,
|
||
&mut claimed,
|
||
r as usize,
|
||
c as usize,
|
||
AttractorType::RiverCrossing,
|
||
s,
|
||
);
|
||
}
|
||
|
||
// 3. CoastalAccess — land within 3 cells of ocean, thinned by spacing.
|
||
// strength = 0.6 + coast-density bonus (capped).
|
||
let mut coastal: Vec<(usize, usize, f32)> = Vec::new();
|
||
for r in 0..h {
|
||
for c in 0..w {
|
||
let i = idx(r, c, w);
|
||
// water_dist (precomputed) is a cheap pre-filter: only cells within
|
||
// 3 of water can be coastal, so skip the 49-cell scan for inland.
|
||
if ta.ocean_mask[i] || claimed[i] || ta.water_dist[i] > 3 {
|
||
continue;
|
||
}
|
||
let near = ocean_cells_within(ta, r, c, 3);
|
||
if near > 0 {
|
||
let bonus = (near as f32 / 24.0).min(0.3);
|
||
coastal.push((r, c, 0.6 + bonus));
|
||
}
|
||
}
|
||
}
|
||
for (r, c, s) in thin_by_spacing(coastal, &claimed, w) {
|
||
claim(
|
||
&mut out,
|
||
&mut claimed,
|
||
r,
|
||
c,
|
||
AttractorType::CoastalAccess,
|
||
s,
|
||
);
|
||
}
|
||
|
||
// 4. ValleyFloor — gentle slope, mid elevation, positive habitability.
|
||
let mut valleys: Vec<(usize, usize, f32)> = Vec::new();
|
||
for r in 0..h {
|
||
for c in 0..w {
|
||
let i = idx(r, c, w);
|
||
if ta.ocean_mask[i] || claimed[i] {
|
||
continue;
|
||
}
|
||
if river_mask[i] || ta.slope_deg[i] >= 5.0 {
|
||
continue;
|
||
}
|
||
if ta.elev_pct[i] < 0.10 || ta.elev_pct[i] > 0.60 {
|
||
continue;
|
||
}
|
||
let hab = habitability(ta.elev_pct[i], ta.slope_deg[i], ta.water_dist[i]);
|
||
if hab > 0.0 {
|
||
valleys.push((r, c, hab));
|
||
}
|
||
}
|
||
}
|
||
for (r, c, s) in thin_by_spacing(valleys, &claimed, w) {
|
||
claim(&mut out, &mut claimed, r, c, AttractorType::ValleyFloor, s);
|
||
}
|
||
|
||
// 5. PassEntrance — morphological saddles in higher terrain.
|
||
// strength = 1 - elev_pct (lower passes score higher).
|
||
let mut passes: Vec<(usize, usize, f32)> = Vec::new();
|
||
for r in 1..h.saturating_sub(1) {
|
||
for c in 0..w {
|
||
let i = idx(r, c, w);
|
||
if ta.ocean_mask[i] || claimed[i] || ta.elev_pct[i] < 0.5 {
|
||
continue;
|
||
}
|
||
if is_saddle(elev, r, c, w, h) {
|
||
passes.push((r, c, 1.0 - ta.elev_pct[i]));
|
||
}
|
||
}
|
||
}
|
||
for (r, c, s) in thin_by_spacing(passes, &claimed, w) {
|
||
claim(&mut out, &mut claimed, r, c, AttractorType::PassEntrance, s);
|
||
}
|
||
|
||
// 6. LakeShore — land adjacent to an enclosed lake. strength = 0.5.
|
||
let mut shores: Vec<(usize, usize, f32)> = Vec::new();
|
||
for r in 0..h {
|
||
for c in 0..w {
|
||
let i = idx(r, c, w);
|
||
if ta.ocean_mask[i] || ta.lake_mask[i] || claimed[i] {
|
||
continue;
|
||
}
|
||
if adjacent_to_lake(ta, r, c) {
|
||
shores.push((r, c, 0.5));
|
||
}
|
||
}
|
||
}
|
||
for (r, c, s) in thin_by_spacing(shores, &claimed, w) {
|
||
claim(&mut out, &mut claimed, r, c, AttractorType::LakeShore, s);
|
||
}
|
||
|
||
// 7. PlainCenter — very flat, away from everything. strength = hab * 0.4.
|
||
let mut plains: Vec<(usize, usize, f32)> = Vec::new();
|
||
for r in 0..h {
|
||
for c in 0..w {
|
||
let i = idx(r, c, w);
|
||
if ta.ocean_mask[i] || claimed[i] || ta.slope_deg[i] >= 2.0 {
|
||
continue;
|
||
}
|
||
let hab = habitability(ta.elev_pct[i], ta.slope_deg[i], ta.water_dist[i]);
|
||
plains.push((r, c, hab * 0.4));
|
||
}
|
||
}
|
||
for (r, c, s) in thin_by_spacing(plains, &claimed, w) {
|
||
claim(&mut out, &mut claimed, r, c, AttractorType::PlainCenter, s);
|
||
}
|
||
|
||
// Cap to MAX_ATTRACTORS while preserving type diversity. D-209 calls
|
||
// RiverMouth "always high-value", but its *normalized* strength
|
||
// (accum / max_accum) is tiny for all but the largest river, so a pure
|
||
// global-strength cap lets abundant ValleyFloor/CoastalAccess crowd every
|
||
// RiverMouth out. Instead: group by type, sort each group strongest-first,
|
||
// then round-robin across types so every present type keeps representation.
|
||
// Deterministic (BTreeMap type order, integer strength key, fixed rotation).
|
||
if out.len() > MAX_ATTRACTORS {
|
||
let mut by_type: std::collections::BTreeMap<u8, Vec<RawAttractor>> =
|
||
std::collections::BTreeMap::new();
|
||
for a in out.drain(..) {
|
||
by_type.entry(a.attractor_type as u8).or_default().push(a);
|
||
}
|
||
for group in by_type.values_mut() {
|
||
group.sort_by(|a, b| {
|
||
// strength is integer now — rank by it directly (descending).
|
||
b.strength
|
||
.cmp(&a.strength)
|
||
.then(a.row.cmp(&b.row))
|
||
.then(a.col.cmp(&b.col))
|
||
});
|
||
}
|
||
let mut kept: Vec<RawAttractor> = Vec::with_capacity(MAX_ATTRACTORS);
|
||
let mut depth = 0usize;
|
||
'fill: loop {
|
||
let mut progressed = false;
|
||
for group in by_type.values() {
|
||
if let Some(a) = group.get(depth) {
|
||
kept.push(*a);
|
||
progressed = true;
|
||
if kept.len() >= MAX_ATTRACTORS {
|
||
break 'fill;
|
||
}
|
||
}
|
||
}
|
||
if !progressed {
|
||
break;
|
||
}
|
||
depth += 1;
|
||
}
|
||
out = kept;
|
||
}
|
||
out.sort_by(|a, b| {
|
||
(a.attractor_type as u8, a.row, a.col).cmp(&(b.attractor_type as u8, b.row, b.col))
|
||
});
|
||
out
|
||
}
|
||
|
||
fn ocean_cells_within(ta: &TerrainAnalysis, r: usize, c: usize, radius: i32) -> usize {
|
||
let mut count = 0;
|
||
for dr in -radius..=radius {
|
||
let nr = r as i32 + dr;
|
||
if nr < 0 || nr >= ta.h as i32 {
|
||
continue;
|
||
}
|
||
for dc in -radius..=radius {
|
||
let nc = wrap_col(c as i32 + dc, ta.w as i32);
|
||
if ta.ocean_mask[idx(nr as usize, nc, ta.w)] {
|
||
count += 1;
|
||
}
|
||
}
|
||
}
|
||
count
|
||
}
|
||
|
||
fn adjacent_to_lake(ta: &TerrainAnalysis, r: usize, c: usize) -> bool {
|
||
for &(dr, dc) in &NB8 {
|
||
let nr = r as i32 + dr;
|
||
if nr < 0 || nr >= ta.h as i32 {
|
||
continue;
|
||
}
|
||
let nc = wrap_col(c as i32 + dc, ta.w as i32);
|
||
if ta.lake_mask[idx(nr as usize, nc, ta.w)] {
|
||
return true;
|
||
}
|
||
}
|
||
false
|
||
}
|
||
|
||
/// Morphological saddle: walking the 8-neighbor ring, the sign of
|
||
/// `(neighbor - cell)` alternates at least 4 times (≥2 higher sectors
|
||
/// separated by ≥2 lower sectors).
|
||
fn is_saddle(elev: &[f32], r: usize, c: usize, w: usize, h: usize) -> bool {
|
||
// Ring order (clockwise) so transitions are meaningful.
|
||
const RING: [(i32, i32); 8] = [
|
||
(-1, 0),
|
||
(-1, 1),
|
||
(0, 1),
|
||
(1, 1),
|
||
(1, 0),
|
||
(1, -1),
|
||
(0, -1),
|
||
(-1, -1),
|
||
];
|
||
let e = elev[idx(r, c, w)];
|
||
let mut signs = [0i8; 8];
|
||
for (k, &(dr, dc)) in RING.iter().enumerate() {
|
||
let nr = r as i32 + dr;
|
||
if nr < 0 || nr >= h as i32 {
|
||
return false; // poles can't be saddles in this scheme
|
||
}
|
||
let nc = wrap_col(c as i32 + dc, w as i32);
|
||
signs[k] = if elev[idx(nr as usize, nc, w)] > e {
|
||
1
|
||
} else {
|
||
-1
|
||
};
|
||
}
|
||
let mut transitions = 0;
|
||
for k in 0..8 {
|
||
if signs[k] != signs[(k + 1) % 8] {
|
||
transitions += 1;
|
||
}
|
||
}
|
||
transitions >= 4
|
||
}
|
||
|
||
/// Greedy spatial thinning: sort candidates by descending strength (ties by
|
||
/// row, col), keep one per `MIN_SPACING` Chebyshev neighborhood.
|
||
///
|
||
/// Uses a bucket grid (cell size = `MIN_SPACING`) so each candidate only checks
|
||
/// the 3×3 neighboring buckets — O(k) amortized rather than O(k²). The kept
|
||
/// order is fully determined by the sorted candidate iteration; the bucket map
|
||
/// is `BTreeMap` (the project bans `HashMap` for determinism) and is lookup-only
|
||
/// regardless.
|
||
fn thin_by_spacing(
|
||
mut cands: Vec<(usize, usize, f32)>,
|
||
_claimed: &[bool],
|
||
_w: usize,
|
||
) -> Vec<(usize, usize, f32)> {
|
||
// Deterministic order: strength desc, then row, col asc.
|
||
cands.sort_by(|a, b| {
|
||
let sa = (a.2 * 1e6) as i64;
|
||
let sb = (b.2 * 1e6) as i64;
|
||
sb.cmp(&sa).then(a.0.cmp(&b.0)).then(a.1.cmp(&b.1))
|
||
});
|
||
let sp = MIN_SPACING.max(1) as usize;
|
||
let mut buckets: BTreeMap<(usize, usize), Vec<(usize, usize)>> = BTreeMap::new();
|
||
let mut kept: Vec<(usize, usize, f32)> = Vec::new();
|
||
for (r, c, s) in cands {
|
||
let (br, bc) = (r / sp, c / sp);
|
||
let mut ok = true;
|
||
'scan: for nbr in br.saturating_sub(1)..=br + 1 {
|
||
for nbc in bc.saturating_sub(1)..=bc + 1 {
|
||
if let Some(pts) = buckets.get(&(nbr, nbc)) {
|
||
for &(kr, kc) in pts {
|
||
let dr = (kr as i32 - r as i32).abs();
|
||
let dc = (kc as i32 - c as i32).abs();
|
||
if dr.max(dc) < MIN_SPACING {
|
||
ok = false;
|
||
break 'scan;
|
||
}
|
||
}
|
||
}
|
||
}
|
||
}
|
||
if ok {
|
||
buckets.entry((br, bc)).or_default().push((r, c));
|
||
kept.push((r, c, s));
|
||
}
|
||
}
|
||
kept
|
||
}
|
||
|
||
#[cfg(test)]
|
||
mod tests {
|
||
use super::*;
|
||
use crate::atlas::drainage;
|
||
|
||
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;
|
||
1.0 - (r as f32 / h as f32 * 0.5 + c as f32 / w as f32 * 0.5)
|
||
})
|
||
.collect()
|
||
}
|
||
|
||
fn hm(data: Vec<f32>, w: u32, h: u32, sea: f32) -> BodyHeightmap {
|
||
BodyHeightmap {
|
||
body_id: "T".into(),
|
||
width: w,
|
||
height: h,
|
||
data,
|
||
sea_level: sea,
|
||
}
|
||
}
|
||
|
||
#[test]
|
||
fn deterministic_extraction() {
|
||
let h = hm(slope_grid(64, 32), 64, 32, 0.3);
|
||
let dr = drainage::analyze(&h.data, 64, 32, 0.3);
|
||
let ta = TerrainAnalysis::analyze(&h, &dr);
|
||
let a1 = extract_attractors(&h, &dr, &ta);
|
||
let a2 = extract_attractors(&h, &dr, &ta);
|
||
assert_eq!(a1, a2, "attractor extraction must be deterministic");
|
||
}
|
||
|
||
#[test]
|
||
fn attractors_sorted_and_bounded() {
|
||
let h = hm(slope_grid(128, 64), 128, 64, 0.3);
|
||
let dr = drainage::analyze(&h.data, 128, 64, 0.3);
|
||
let ta = TerrainAnalysis::analyze(&h, &dr);
|
||
let a = extract_attractors(&h, &dr, &ta);
|
||
assert!(a.len() <= MAX_ATTRACTORS);
|
||
// Sorted by (type as u8, row, col).
|
||
for win in a.windows(2) {
|
||
let ka = (win[0].attractor_type as u8, win[0].row, win[0].col);
|
||
let kb = (win[1].attractor_type as u8, win[1].row, win[1].col);
|
||
assert!(ka <= kb, "attractors must be sorted");
|
||
}
|
||
}
|
||
|
||
#[test]
|
||
fn percentile_in_range() {
|
||
let h = hm(slope_grid(32, 16), 32, 16, 0.3);
|
||
let dr = drainage::analyze(&h.data, 32, 16, 0.3);
|
||
let ta = TerrainAnalysis::analyze(&h, &dr);
|
||
assert!(ta.elev_pct.iter().all(|&p| (0.0..=1.0).contains(&p)));
|
||
assert_eq!(ta.slope_deg.len(), 32 * 16);
|
||
}
|
||
|
||
/// Multi-octave sine terrain (continents + many small coastal streams) —
|
||
/// produces > MAX_ATTRACTORS candidates with plenty of river mouths.
|
||
fn sine_grid(w: u32, h: u32) -> Vec<f32> {
|
||
use std::f32::consts::{PI, TAU};
|
||
(0..(w * h))
|
||
.map(|i| {
|
||
let r = (i / w) as f32;
|
||
let c = (i % w) as f32;
|
||
let x = c / w as f32 * TAU;
|
||
let y = r / h as f32 * PI;
|
||
(0.5 + 0.25 * (x * 3.0).sin() * (y * 2.0).sin()
|
||
+ 0.15 * (x * 7.0).cos() * (y * 5.0).sin()
|
||
+ 0.08 * (x * 13.0).sin() * (y * 11.0).cos()
|
||
+ 0.05 * (x * 23.0).cos() * (y * 19.0).sin())
|
||
.clamp(0.0, 1.0)
|
||
})
|
||
.collect()
|
||
}
|
||
|
||
#[test]
|
||
fn river_mouths_survive_cap() {
|
||
// D-209 + the type-aware cap: even though RiverMouth normalized strength
|
||
// is tiny, a body full of mouths must still keep RiverMouth attractors
|
||
// (a global-strength cap would drop all of them — the bug Hoshe caught).
|
||
let h = hm(sine_grid(512, 256), 512, 256, 0.40);
|
||
let dr = drainage::analyze(&h.data, 512, 256, 0.40);
|
||
assert!(
|
||
!dr.river_network.mouths.is_empty(),
|
||
"fixture must have mouths"
|
||
);
|
||
let ta = TerrainAnalysis::analyze(&h, &dr);
|
||
let a = extract_attractors(&h, &dr, &ta);
|
||
assert!(a.len() <= MAX_ATTRACTORS);
|
||
assert!(
|
||
a.iter()
|
||
.any(|x| x.attractor_type == AttractorType::RiverMouth),
|
||
"RiverMouth attractors must survive the cap when mouths exist"
|
||
);
|
||
}
|
||
}
|