Files
settled-reach/server/src/atlas/features.rs
T
jpmschweitzerandClaude Opus 4.8 b73b4e6d92 feat(simulation): Layer-1 water bearing + D-234 waterfront rule (#957)
Close the last D-234 piece — terrain water-direction extraction wired through
to founding orientation and the quarter waterfront rule:

- Layer 1: TerrainAnalysis::water_bearing — 8-octant integer bearing toward the
  nearest water from the water_dist gradient (D-010, no atan2). Stored on
  GeographicAttractor.water_bearing (360 = none).
- #956 founding orientation: coastal/river settlements now get a real
  water-facing bearing (the anchoring attractor's), replacing the 0 stub.
- #957 waterfront rule (D-234b): the water-facing quarter edge (from the
  settlement's Coastal founding orientation) drops its block setback to 0 so
  buildings present flush to the quay (dock-orthogonal). Typed Edge + coastal_edge
  + per-block gating.

Golden + atlas_response fixture rebaked (additive water_bearing field only).
8 new tests. All integer-deterministic (D-010).

Pending: the waterfront rule reads context.founding_orientation, which
city_context_reader still stubs to Cardinal — real per-settlement orientation
reaches quarter generation once the Layer-3 placement -> Layer-4 GenerateSkeleton
dispatch is wired (the remaining cross-layer integration).

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
2026-06-06 10:22:14 +02:00

786 lines
27 KiB
Rust
Raw Blame History

This file contains ambiguous Unicode characters
This file contains Unicode characters that might be confused with other characters. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.
//! 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>,
}
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,
}
}
#[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"
);
}
}