diff --git a/server/src/atlas/detail_scatter.rs b/server/src/atlas/detail_scatter.rs new file mode 100644 index 000000000..479361f39 --- /dev/null +++ b/server/src/atlas/detail_scatter.rs @@ -0,0 +1,162 @@ +//! Detail-scatter synthesis (D-243 §2) — the mid-scale terrain (≈2–40 km) that +//! the coarse heightmap (~40–78 km/pixel) is too coarse to carry and the voxel +//! generators (<2 km, inside a district) are too fine to reach. This is the +//! "stretch magic" that fills the gap so an interpolated heightmap reads as real +//! terrain instead of smooth ramps between samples. +//! +//! **Character (locked, T-1077):** +//! - **Adaptive** — ruggedness follows the local slope/morphology: gentle flats +//! stay soft (low amplitude, rolling), rugged highs gain ridge/valley energy. +//! - **Envelope + modulation** — the heightmap's local relief is the *ceiling* on +//! amplitude (smooth where the authored surface is smooth, so a lore body's +//! plain never sprouts a mountain), and morphology *modulates the character* +//! within that ceiling. +//! +//! **Determinism (D-227 / D-010):** a pure, seed+position-keyed value-noise +//! function — no stored state, no RNG, no platform-variant lookup. The output is +//! an `f64` *positional* perturbation; structural decisions downstream quantise it +//! to integer `slope_q`/`elev_q` (the truncate-before-decision convention). + +use crate::seed::splitmix64; + +/// Mid-scale octave wavelengths in metres — the 2–40 km band. Coarsest first. +/// Below the finest (~4 km) the district→voxel layers own the detail; above the +/// coarsest (~33 km) the heightmap itself carries the shape. +const OCTAVE_WAVELENGTHS_M: [f64; 4] = [32_768.0, 16_384.0, 8_192.0, 4_096.0]; + +/// Deterministic lattice value in `[-1, 1)` for an integer noise cell. +#[inline] +fn lattice(seed: u64, ix: i64, iy: i64) -> f64 { + let mixed = (ix as u64).wrapping_mul(0x9E37_79B9_7F4A_7C15) + ^ (iy as u64).wrapping_mul(0xBF58_476D_1CE4_E5B9); + let h = splitmix64(seed ^ splitmix64(mixed)); + // Top 53 bits → [0, 1) → [-1, 1). + ((h >> 11) as f64 / (1u64 << 53) as f64) * 2.0 - 1.0 +} + +/// Smooth (C¹) value-noise sample in `[-1, 1]` at world `(wx, wy)` for one +/// wavelength. Smoothstep interpolation keeps lattice cell boundaries crease-free. +fn value_noise(seed: u64, wx: f64, wy: f64, wavelength_m: f64) -> f64 { + let fx = wx / wavelength_m; + let fy = wy / wavelength_m; + let (x0, y0) = (fx.floor(), fy.floor()); + let (tx, ty) = (fx - x0, fy - y0); + let sx = tx * tx * (3.0 - 2.0 * tx); + let sy = ty * ty * (3.0 - 2.0 * ty); + let (ix, iy) = (x0 as i64, y0 as i64); + let c00 = lattice(seed, ix, iy); + let c10 = lattice(seed, ix + 1, iy); + let c01 = lattice(seed, ix, iy + 1); + let c11 = lattice(seed, ix + 1, iy + 1); + let a = c00 + (c10 - c00) * sx; + let b = c01 + (c11 - c01) * sx; + a + (b - a) * sy +} + +/// The adaptive mid-scale elevation perturbation in `[0,1]`-normalized units (the +/// same space as `elev_pct`); multiply by the body's elevation span for metres. +/// +/// - `wx`, `wy`: world-metre coordinates (absolute; keyed so features are +/// position-stable, never folded into a chunk frame). +/// - `envelope` (0–1): the heightmap's local relief — the **amplitude ceiling**. +/// - `ruggedness` (0–1): from local slope/morphology — **adaptive** ridge energy. +/// +/// Returns roughly `[-envelope, +envelope]`, smaller and smoother as `ruggedness` +/// drops toward 0 (gentle flats), larger and ridged as it rises toward 1. +pub fn terrain_detail(seed: u64, wx: f64, wy: f64, envelope: f64, ruggedness: f64) -> f64 { + let env = envelope.clamp(0.0, 1.0); + let rug = ruggedness.clamp(0.0, 1.0); + if env == 0.0 { + return 0.0; // perfectly flat heightmap → no invented relief (envelope rule) + } + + let mut sum = 0.0; + let mut amp = 1.0; + let mut norm = 0.0; + for (i, &wl) in OCTAVE_WAVELENGTHS_M.iter().enumerate() { + let mut n = value_noise( + seed.wrapping_add((i as u64).wrapping_mul(0x1000)), + wx, + wy, + wl, + ); + // Ridged folding blends in with ruggedness: flats stay smooth (rolling + // fBm), highs fold to sharp ridges/incised valleys (modulation). + if rug > 0.0 { + let ridged = 1.0 - 2.0 * n.abs(); // ridge crest at zero-crossings + n = n * (1.0 - rug) + ridged * rug; + } + sum += n * amp; + norm += amp; + // Higher octaves gain weight in rugged zones (rough detail); flats damp + // them toward a single gentle swell. + amp *= 0.5 + 0.35 * rug; + } + let fbm = sum / norm; // ≈ [-1, 1] + + // Envelope caps the amplitude; within the cap, ruggedness scales how much of + // the headroom is used (flats use ~40%, peaks ~100%). + fbm * env * (0.4 + 0.6 * rug) +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn deterministic() { + let a = terrain_detail(42, 12_345.0, -6_789.0, 0.6, 0.5); + let b = terrain_detail(42, 12_345.0, -6_789.0, 0.6, 0.5); + assert_eq!(a, b); + } + + #[test] + fn flat_envelope_invents_nothing() { + // envelope = 0 → the authored heightmap is flat here → no relief (the + // envelope rule: never sprout terrain on an authored plain). + for &rug in &[0.0, 0.5, 1.0] { + assert_eq!(terrain_detail(7, 1000.0, 2000.0, 0.0, rug), 0.0); + } + } + + #[test] + fn bounded_by_envelope() { + // Output magnitude never exceeds the envelope ceiling, anywhere. + for i in 0..400 { + let wx = (i as f64) * 137.0; + let wy = (i as f64) * -91.0; + let v = terrain_detail(99, wx, wy, 0.5, 1.0); + assert!(v.abs() <= 0.5 + 1e-9, "v={v} exceeded envelope at {i}"); + } + } + + #[test] + fn rugged_is_higher_amplitude_than_gentle() { + // Adaptive: averaged |relief| under high ruggedness exceeds that under low + // ruggedness for the same envelope (gentle flats vs rugged highs). + let mean_abs = |rug: f64| -> f64 { + let n = 500; + (0..n) + .map(|i| terrain_detail(3, i as f64 * 53.0, i as f64 * 71.0, 0.7, rug).abs()) + .sum::() + / n as f64 + }; + assert!( + mean_abs(0.9) > mean_abs(0.1), + "rugged terrain must carry more relief than gentle terrain" + ); + } + + #[test] + fn continuous_no_creases() { + // Small position steps produce small output changes (C¹ value noise) — no + // lattice creases that would read as grid artifacts. + let base = terrain_detail(11, 5_000.0, 5_000.0, 0.8, 0.6); + let near = terrain_detail(11, 5_000.5, 5_000.0, 0.8, 0.6); + assert!( + (base - near).abs() < 0.05, + "0.5 m step jumped by {}", + (base - near).abs() + ); + } +} diff --git a/server/src/atlas/district_profile.rs b/server/src/atlas/district_profile.rs index d0dd1e4d2..e5fbd8312 100644 --- a/server/src/atlas/district_profile.rs +++ b/server/src/atlas/district_profile.rs @@ -23,6 +23,7 @@ use std::collections::BTreeMap; use serde::{Deserialize, Serialize}; use crate::atlas::features::TerrainAnalysis; +use crate::atlas::scale; use crate::seed::SeedChain; use crate::simulation::generator::MorphologyZone; @@ -950,6 +951,28 @@ pub fn derive_district_profile( (0, 0, 0) }; + build_district_profile( + seed, + body_params, + climate, + slope_q, + elev_q, + ocean_fraction_q, + ) +} + +/// Build the climate + morphology fields of a `DistrictProfile` from its three +/// terrain primitives (`slope_q`, `elev_q`, `ocean_fraction_q`) — the shared tail +/// of every derivation path (cell-aggregate [`derive_district_profile`] and the +/// interpolation+scatter [`derive_district`]). Pure (T-1024, D-239 §2 / D-240). +fn build_district_profile( + seed: SeedChain, + body_params: &BodyParams, + climate: &ClimateConstants, + slope_q: i32, + elev_q: i32, + ocean_fraction_q: i32, +) -> DistrictProfile { let tectonic_class = derive_tectonic_class(body_params); // Climate derivation (T-1024, D-239 §2). Temperature lapse must vary by THIS @@ -959,13 +982,13 @@ pub fn derive_district_profile( // district's own normalized elevation (elev_q, 0–100) scaled to the body's // elevation span. district_latitude_deg is already per-district (set by the caller // / derive_all_districts). Per-cell refinement happens later at ChunkContext. - let region_climate_params = BodyParams { + let district_climate_params = BodyParams { elevation_km: (elev_q as f64 / 100.0) * MAX_REGION_ELEVATION_KM, ..body_params.clone() }; // D-240: body-scoped seed for the deterministic per-body temperature nudge. let body_seed = seed.seed(); - let temperature_c = derive_temperature_c(®ion_climate_params, climate, body_seed); + let temperature_c = derive_temperature_c(&district_climate_params, climate, body_seed); let moisture_q = derive_moisture_q(body_params); // Climate-derived fields: computed from temperature + moisture primitives @@ -1004,7 +1027,140 @@ pub fn derive_district_profile( } } -/// Derive district profiles for all districts covering the body. +/// On-demand 2 km district profile (D-243 §2, T-1077) — the corrected carrier. +/// +/// Maps the district to a fractional heightmap position via the body radius (the +/// elastic seam, D-204), bilinearly interpolates the Layer-1 terrain (the +/// continental *envelope*), then composes the adaptive detail-scatter +/// ([`crate::atlas::detail_scatter`]) for the mid-scale relief the heightmap is +/// too coarse to carry. `body_radius_km = None` (e.g. tiny test bodies) falls back +/// to direct heightmap indexing. +/// +/// Pure and deterministic (D-227/D-010): a function of +/// `(seed, body_params, terrain, district_pos)`; the f64 scatter is quantised to +/// integer `slope_q`/`elev_q` at the decision boundary. +pub fn derive_district( + seed: SeedChain, + body_params: &BodyParams, + ta: &TerrainAnalysis, + district_pos: DistrictPos, + climate: &ClimateConstants, +) -> DistrictProfile { + let (dx, dy) = district_pos; + + // District → fractional heightmap pixel + world-metre coordinate + latitude. + let (px, py, world_x_m, world_y_m, lat_deg) = match body_params.body_radius_km { + Some(r_km) if r_km > 0.0 => { + let circumference_m = std::f64::consts::TAU * r_km * 1000.0; + let meridian_m = std::f64::consts::PI * r_km * 1000.0; + let wx = dx as f64 * scale::DISTRICT_M as f64; + let wy = dy as f64 * scale::DISTRICT_M as f64; + // Longitude wraps; district (0,0) sits at lon 0 / the equator. + let px = (wx / circumference_m).rem_euclid(1.0) * ta.w as f64; + // Latitude: equator at py = h/2, clamped at the poles. + let lat_frac = (wy / meridian_m).clamp(-0.5, 0.5); // −0.5 = N pole, +0.5 = S + let py = (0.5 + lat_frac) * ta.h.saturating_sub(1) as f64; + (px, py, wx, wy, -lat_frac * 180.0) + } + _ => { + // No radius: the district grid IS the heightmap grid (tiny test bodies). + let px = (dx as f64).clamp(0.0, ta.w.saturating_sub(1) as f64); + let py = (dy as f64).clamp(0.0, ta.h.saturating_sub(1) as f64); + let lat_deg = if ta.h > 1 { + 90.0 - (py / (ta.h - 1) as f64) * 180.0 + } else { + 0.0 + }; + let dm = scale::DISTRICT_M as f64; + (px, py, px * dm, py * dm, lat_deg) + } + }; + + // Bilinear-interpolated Layer-1 envelope at the district position. + let elev_pct = bilinear(&ta.elev_pct, ta.w, ta.h, px, py) as f64; + let slope_deg = bilinear(&ta.slope_deg, ta.w, ta.h, px, py) as f64; + let ocean_frac = bilinear_bool(&ta.ocean_mask, ta.w, ta.h, px, py) as f64; + + // Envelope + adaptive ruggedness from the local heightmap slope (steep + // heightmap ⇒ relief headroom ⇒ rugged). Both 0–1. + let local_slope = (slope_deg / 45.0).clamp(0.0, 1.0); + let envelope = local_slope; + let ruggedness = local_slope; + let scatter = crate::atlas::detail_scatter::terrain_detail( + seed.seed(), + world_x_m, + world_y_m, + envelope, + ruggedness, + ); + + // Compose the primitives. Scatter perturbs elevation (mid-scale relief); slope + // gains a ruggedness-weighted bump so the morphology classifier responds + // (adaptive: rugged highs → mountain/pass, flats → plains). + let elev_q = (((elev_pct + scatter) * 100.0).round() as i32).clamp(0, 100); + let slope_q = + (((local_slope + ruggedness * scatter.abs()) * 100.0).round() as i32).clamp(0, 100); + let ocean_fraction_q = ((ocean_frac * 100.0).round() as i32).clamp(0, 100); + + let params = BodyParams { + district_latitude_deg: lat_deg, + ..body_params.clone() + }; + build_district_profile(seed, ¶ms, climate, slope_q, elev_q, ocean_fraction_q) +} + +/// Bilinear interpolation of a row-major `f32` field at fractional `(px, py)`. +/// Columns wrap (equirectangular); rows clamp at the poles. +fn bilinear(field: &[f32], w: usize, h: usize, px: f64, py: f64) -> f32 { + if w == 0 || h == 0 { + return 0.0; + } + let x0 = px.floor(); + let y0 = py.floor().clamp(0.0, (h - 1) as f64); + let tx = (px - x0) as f32; + let ty = (py - y0) as f32; + let ix0 = (x0 as i64).rem_euclid(w as i64) as usize; + let ix1 = (ix0 + 1) % w; + let iy0 = (y0 as usize).min(h - 1); + let iy1 = (iy0 + 1).min(h - 1); + let v00 = field[iy0 * w + ix0]; + let v10 = field[iy0 * w + ix1]; + let v01 = field[iy1 * w + ix0]; + let v11 = field[iy1 * w + ix1]; + let a = v00 + (v10 - v00) * tx; + let b = v01 + (v11 - v01) * tx; + a + (b - a) * ty +} + +/// Bilinear interpolation of a boolean mask as a 0–1 fraction (for ocean coverage). +fn bilinear_bool(mask: &[bool], w: usize, h: usize, px: f64, py: f64) -> f32 { + if w == 0 || h == 0 { + return 0.0; + } + let x0 = px.floor(); + let y0 = py.floor().clamp(0.0, (h - 1) as f64); + let tx = (px - x0) as f32; + let ty = (py - y0) as f32; + let ix0 = (x0 as i64).rem_euclid(w as i64) as usize; + let ix1 = (ix0 + 1) % w; + let iy0 = (y0 as usize).min(h - 1); + let iy1 = (iy0 + 1).min(h - 1); + let f = |r: usize, c: usize| mask[r * w + c] as i32 as f32; + let a = f(iy0, ix0) + (f(iy0, ix1) - f(iy0, ix0)) * tx; + let b = f(iy1, ix0) + (f(iy1, ix1) - f(iy1, ix0)) * tx; + a + (b - a) * ty +} + +/// Eagerly derive a coarse profile grid covering the body, by direct heightmap +/// tiling (`grid_cells_per_district` cells per cell). +/// +/// **Scale note (D-243, T-1077):** this is the *coarse* eager grid — one cell per +/// `gcpr` heightmap pixels (tens-to-hundreds of km) — kept as the Atlas zone +/// overlay source. The **corrected 2 km carrier** the voxel chain consumes is the +/// on-demand [`derive_district`] (heightmap interpolation + detail-scatter via the +/// elastic seam). Swapping production from this eager grid to on-demand +/// district/region surfacing is the chartered job of T-1046; T-1077 provides the +/// correct on-demand derivation, not the eager/Atlas restructure. /// /// Returns a `BTreeMap` covering the full /// heightmap at the given district-grid resolution. @@ -1101,6 +1257,96 @@ mod tests { assert_eq!(districts.len(), 32, "district count mismatch"); } + // --- derive_district (on-demand 2 km, interpolation + detail-scatter) ----- + + fn earth_params() -> BodyParams { + BodyParams { + hydrosphere: Some("ocean".into()), + atmosphere: Some("breathable".into()), + planet_class: Some("temperate".into()), + body_radius_km: Some(6371.0), + ..Default::default() + } + } + + #[test] + fn derive_district_is_deterministic() { + let hm = test_hm(); + let ta = test_ta(&hm); + let climate = ClimateConstants::default(); + let p = earth_params(); + let a = derive_district(test_seed(), &p, &ta, (1234, -567), &climate); + let b = derive_district(test_seed(), &p, &ta, (1234, -567), &climate); + assert_eq!(a.elev_q, b.elev_q); + assert_eq!(a.slope_q, b.slope_q); + assert_eq!(a.morphology_zone, b.morphology_zone); + assert_eq!(a.temperature_c, b.temperature_c); + } + + #[test] + fn derive_district_radius_maps_to_latitude_climate() { + // With a body radius, equatorial vs near-polar districts get different + // temperature (the seam maps district_y → latitude). Pole = colder. + let hm = test_hm(); + let ta = test_ta(&hm); + let climate = ClimateConstants::default(); + let p = earth_params(); + // meridian ≈ π·6371·1000 m; a district near the pole is ~quarter-meridian away. + let merid_districts = + (std::f64::consts::PI * 6371.0 * 1000.0 / scale::DISTRICT_M as f64) as i32; + let equator = derive_district(test_seed(), &p, &ta, (0, 0), &climate); + let high_lat = + derive_district(test_seed(), &p, &ta, (0, merid_districts / 2 - 2), &climate); + match (equator.temperature_c, high_lat.temperature_c) { + (Some(eq), Some(hi)) => { + assert!(hi < eq, "near-pole district must be colder ({hi} !< {eq})") + } + _ => panic!("breathable body must have a temperature"), + } + } + + #[test] + fn derive_district_flat_envelope_invents_no_relief() { + // A perfectly flat heightmap (no slope) → envelope 0 → scatter invents + // nothing → elev_q equals the interpolated base (the envelope rule). + let (w, h) = (64u32, 32u32); + let flat = BodyHeightmap { + body_id: "flat".into(), + width: w, + height: h, + data: vec![0.6; (w * h) as usize], + sea_level: 0.3, + }; + let ta = test_ta(&flat); + let climate = ClimateConstants::default(); + let p = earth_params(); + let d = derive_district(test_seed(), &p, &ta, (500, 100), &climate); + // Flat land everywhere → slope 0 → no invented relief, no ocean. + assert_eq!( + d.slope_q, 0, + "flat heightmap must yield zero district slope" + ); + assert_eq!(d.ocean_fraction_q, 0, "above sea level → no ocean"); + } + + #[test] + fn derive_district_no_radius_falls_back_to_direct_indexing() { + // body_radius_km = None (tiny test bodies): the district grid is the + // heightmap grid; derivation still succeeds and is deterministic. + let hm = test_hm(); + let ta = test_ta(&hm); + let climate = ClimateConstants::default(); + let p = BodyParams { + planet_class: Some("temperate".into()), + atmosphere: Some("breathable".into()), + ..Default::default() // body_radius_km: None + }; + let a = derive_district(test_seed(), &p, &ta, (20, 10), &climate); + let b = derive_district(test_seed(), &p, &ta, (20, 10), &climate); + assert_eq!(a.elev_q, b.elev_q); + assert!((0..=100).contains(&a.elev_q) && (0..=100).contains(&a.slope_q)); + } + #[test] fn derive_district_profile_is_deterministic() { let hm = test_hm(); diff --git a/server/src/atlas/mod.rs b/server/src/atlas/mod.rs index b34c7f340..9e028dd91 100644 --- a/server/src/atlas/mod.rs +++ b/server/src/atlas/mod.rs @@ -10,6 +10,7 @@ pub mod body_world_state; pub mod cascade; pub mod chunk_context; pub mod city_context_reader; +pub mod detail_scatter; pub mod district_mix; pub mod district_profile; pub mod domain_warp;