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

//! 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.
///
/// `pub(crate)`: also the quantization band set for `layer_proxy`'s
/// `window_min_wl_m` (T-1150, zoom ladder design doc §5 — "quantize
/// `min_wl_m` to a small fixed set of bands per rung, **matching the rung's
/// own octave bands**"). District/quarter windows both derive via
/// `terrain_detail`, so this IS "the rung's own octave bands" for both rungs
/// today — one array, no duplicated magic numbers that could drift out of
/// sync with the actual cutoff behavior.
pub(crate) const OCTAVE_WAVELENGTHS_M: [f64; 4] = [32_768.0, 16_384.0, 8_192.0, 4_096.0];

/// Voxel-tier octave wavelengths in metres — the ≈0.13–1 km **sub-district** band
/// (all finer than the 2 km district planning unit) that the district-tier
/// [`terrain_detail`] (4–40 km) is too coarse to carry and the per-voxel scatter
/// (<64 m) is too fine to reach. This is the [`voxel_relief`] band: the rolling/ridged
/// hills a *walking character* navigates by (T-1081). The coarsest octave stays below
/// the district size so the relief never competes with `elev_q`'s district-scale role.
///
/// `pub(crate)` (T-1162): also the fine-tier band
/// [`crate::atlas::vegetation_invention`]'s texture field truncates against,
/// alongside [`OCTAVE_WAVELENGTHS_M`] — the same reuse-not-reinvent posture
/// that array already documents for `layer_proxy`.
pub(crate) const VOXEL_OCTAVE_WAVELENGTHS_M: [f64; 4] = [1_024.0, 512.0, 256.0, 128.0];

/// Sub-chunk mosaic octave wavelengths in metres — the ≈8–64 m band, finer than the
/// [`VOXEL_OCTAVE_WAVELENGTHS_M`] sub-district band, so a micro-habitat patch reads as
/// a distinct blob *inside* a single 64 m chunk (T-1084, D-246). Coarsest first.
const MOSAIC_OCTAVE_WAVELENGTHS_M: [f64; 4] = [64.0, 32.0, 16.0, 8.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.
/// `pub(crate)`: also the noise substrate of the T-1125 coastline invention
/// ([`crate::atlas::coast_invention`]), which salts its own seed stream.
pub(crate) 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.
///
/// `min_wavelength_m` (T-1149, zoom ladder §2): octaves whose wavelength is
/// below this cutoff are skipped entirely — the Nyquist truncation a coarse
/// sample density needs (no point paying for detail finer than the sample
/// spacing can resolve). `0.0` = no cutoff = every octave, byte-identical to
/// pre-T-1149 behavior.
pub fn terrain_detail(
    seed: u64,
    wx: f64,
    wy: f64,
    envelope: f64,
    ruggedness: f64,
    min_wavelength_m: f64,
) -> f64 {
    enveloped_fbm(
        seed,
        wx,
        wy,
        envelope,
        ruggedness,
        &OCTAVE_WAVELENGTHS_M,
        min_wavelength_m,
    )
}

/// The voxel-tier mid-scale relief perturbation in `[0,1]`-normalized units
/// (multiply by the body's relief span for metres) — the ≈0.25–2 km rolling/ridged
/// hills a walking character navigates by (T-1081). Identical envelope/ruggedness
/// contract as [`terrain_detail`]; only the octave band differs
/// ([`VOXEL_OCTAVE_WAVELENGTHS_M`] vs the district [`OCTAVE_WAVELENGTHS_M`]).
///
/// Seed must be **body-global** (constant across the body) — the position
/// `(wx, wy)` carries the variation. A per-voxel seed would make every voxel a
/// fresh lattice (white noise, not smooth hills).
///
/// `min_wavelength_m` — see [`terrain_detail`]; `0.0` = no cutoff.
pub fn voxel_relief(
    seed: u64,
    wx: f64,
    wy: f64,
    envelope: f64,
    ruggedness: f64,
    min_wavelength_m: f64,
) -> f64 {
    enveloped_fbm(
        seed,
        wx,
        wy,
        envelope,
        ruggedness,
        &VOXEL_OCTAVE_WAVELENGTHS_M,
        min_wavelength_m,
    )
}

/// The sub-chunk micro-habitat **selection field** in `[0, 1]` — a smooth,
/// body-global, position-keyed multi-octave value-noise scalar the mosaic pass
/// thresholds against a per-class palette's cumulative weights (T-1084, D-246).
///
/// Unlike [`voxel_relief`] this is not an elevation perturbation: it carries no
/// envelope and no ridge folding — it is a plain fBm remapped to `[0, 1]` whose only
/// job is to pick *which* micro-habitat a voxel lands in, as spatially-coherent blobs
/// (low-frequency → patches, never per-tile salt-and-pepper). The seed must be
/// **body-global** (constant across the body); the position `(wx, wy)` carries the
/// variation, exactly as for [`voxel_relief`].
pub fn voxel_mosaic(seed: u64, wx: f64, wy: f64) -> f64 {
    let mut sum = 0.0;
    let mut amp = 1.0;
    let mut norm = 0.0;
    for (i, &wl) in MOSAIC_OCTAVE_WAVELENGTHS_M.iter().enumerate() {
        let n = value_noise(
            seed.wrapping_add((i as u64).wrapping_mul(0x1000)),
            wx,
            wy,
            wl,
        );
        sum += n * amp;
        norm += amp;
        amp *= 0.5;
    }
    let fbm = sum / norm; // ≈ [-1, 1]
    ((fbm + 1.0) * 0.5).clamp(0.0, 1.0)
}

/// Shared adaptive-fBm core for [`terrain_detail`] and [`voxel_relief`] — the only
/// difference between the two tiers is the octave wavelength band. Returns the
/// enveloped, ruggedness-modulated perturbation (roughly `[-envelope, +envelope]`).
///
/// `min_wavelength_m` (T-1149): octaves with `wl < min_wavelength_m` are
/// skipped — HARD TRUNCATE, not amplitude-faded. `amp`/`norm` still advance
/// through the skipped octave's weight step (`i` keeps its position in the
/// wavelength array for the seed-salt term), so the surviving octaves keep
/// their same relative weighting as if the truncated tail were simply cut
/// off the sum, not renormalized against a smaller octave count. `0.0` = no
/// octave is ever skipped = today's behavior, byte-for-byte.
///
/// **Seam (not built, see T-1149/design doc §2 + §9 R1):** a hard truncate can
/// pop when the sample density crosses an octave boundary between two
/// requests (an octave present at one zoom step vanishes at the next,
/// discontinuously). The documented fix is fading the highest surviving
/// octave's amplitude toward zero as `wl` approaches `min_wavelength_m` from
/// above, rather than an on/off cut. Not implemented here — it needs a
/// visual A/B against a real client zoom ladder, which does not exist yet
/// (T-1150/T-1153); building it speculatively risks tuning against nothing.
fn enveloped_fbm(
    seed: u64,
    wx: f64,
    wy: f64,
    envelope: f64,
    ruggedness: f64,
    wavelengths: &[f64],
    min_wavelength_m: 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 wavelengths.iter().enumerate() {
        if wl < min_wavelength_m {
            // Below the per-sample Nyquist cutoff — skip the term entirely
            // (hard truncate), but still advance amp so later octaves (there
            // are none finer in these const arrays, but the rule is general)
            // keep their intended relative weight.
            amp *= 0.5 + 0.35 * rug;
            continue;
        }
        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;
    }
    if norm == 0.0 {
        return 0.0; // every octave cut by the cutoff → no invented relief left
    }
    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, 0.0);
        let b = terrain_detail(42, 12_345.0, -6_789.0, 0.6, 0.5, 0.0);
        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), 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, 0.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, 0.0).abs())
                .sum::<f64>()
                / 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, 0.0);
        let near = terrain_detail(11, 5_000.5, 5_000.0, 0.8, 0.6, 0.0);
        assert!(
            (base - near).abs() < 0.05,
            "0.5 m step jumped by {}",
            (base - near).abs()
        );
    }

    // ── min_wavelength_m cutoff (T-1149, zoom ladder §2) ──────────────────────

    #[test]
    fn cutoff_zero_matches_pre_t1149_behavior() {
        // 0.0 = no cutoff = every octave — this is the compatibility contract
        // every existing caller (derive_district's default) relies on.
        for i in 0..200 {
            let wx = i as f64 * 91.0;
            let wy = i as f64 * -53.0;
            let with_explicit_zero = terrain_detail(21, wx, wy, 0.6, 0.5, 0.0);
            // The finest OCTAVE_WAVELENGTHS_M entry is 4_096.0 — a cutoff below
            // that admits every octave too, and must agree exactly.
            let with_below_finest = terrain_detail(21, wx, wy, 0.6, 0.5, 1.0);
            assert_eq!(with_explicit_zero, with_below_finest);
        }
    }

    #[test]
    fn cutoff_truncates_octaves_below_it() {
        // A cutoff placed above the coarsest OCTAVE_WAVELENGTHS_M entry
        // (32_768.0) must skip every octave and fall back to 0.0 (the
        // norm==0.0 empty-sum guard), same as the flat_envelope_invents_nothing
        // envelope==0 case but reached via the cutoff instead.
        assert_eq!(
            terrain_detail(5, 1_000.0, 2_000.0, 0.6, 0.5, 100_000.0),
            0.0
        );
    }

    #[test]
    fn cutoff_changes_output_relative_to_uncut() {
        // A mid-band cutoff (drops the two finest octaves: 8_192.0, 4_096.0)
        // must produce DIFFERENT output than the uncut derive at the same
        // position — otherwise the cutoff parameter would be a no-op.
        let uncut = terrain_detail(17, 12_000.0, 9_000.0, 0.7, 0.6, 0.0);
        let cut = terrain_detail(17, 12_000.0, 9_000.0, 0.7, 0.6, 8_193.0);
        assert_ne!(
            uncut, cut,
            "a mid-band cutoff must change the derived output"
        );
    }

    // ── voxel_relief (T-1081): same contract, sub-district band ──────────────

    #[test]
    fn voxel_relief_deterministic_and_bounded() {
        let a = voxel_relief(42, 12_345.0, -6_789.0, 0.6, 0.5, 0.0);
        let b = voxel_relief(42, 12_345.0, -6_789.0, 0.6, 0.5, 0.0);
        assert_eq!(a, b);
        // Same envelope rule + amplitude ceiling as terrain_detail.
        assert_eq!(voxel_relief(7, 1_000.0, 2_000.0, 0.0, 1.0, 0.0), 0.0);
        for i in 0..400 {
            let v = voxel_relief(99, i as f64 * 137.0, i as f64 * -91.0, 0.5, 1.0, 0.0);
            assert!(v.abs() <= 0.5 + 1e-9, "v={v} exceeded envelope at {i}");
        }
    }

    #[test]
    fn voxel_relief_varies_at_sub_district_scale() {
        // The voxel band must produce DISTINCT relief within a single 2 km district —
        // exactly the variation terrain_detail's 4–40 km band cannot carry, and the
        // flatness T-1081 fixes. Body-global seed; position carries the variation. Steps
        // are deliberately non-aligned with the octave wavelengths to avoid aliasing.
        let seed = 1234;
        let vals: Vec<f64> = (0..16)
            .map(|i| voxel_relief(seed, i as f64 * 137.0, i as f64 * 89.0, 0.7, 0.6, 0.0))
            .collect();
        let min = vals.iter().cloned().fold(f64::INFINITY, f64::min);
        let max = vals.iter().cloned().fold(f64::NEG_INFINITY, f64::max);
        assert!(
            max - min > 0.1,
            "voxel relief is near-flat across a district: range {}",
            max - min
        );
    }

    // ── voxel_mosaic (T-1084): sub-chunk selection field ─────────────────────

    #[test]
    fn voxel_mosaic_deterministic_and_in_range() {
        let a = voxel_mosaic(42, 12_345.0, -6_789.0);
        let b = voxel_mosaic(42, 12_345.0, -6_789.0);
        assert_eq!(a, b);
        for i in 0..400 {
            let v = voxel_mosaic(99, i as f64 * 137.0, i as f64 * -91.0);
            assert!((0.0..=1.0).contains(&v), "v={v} out of [0,1] at {i}");
        }
    }

    #[test]
    fn voxel_mosaic_varies_at_sub_chunk_scale() {
        // The mosaic must produce DISTINCT values within a single 64 m chunk — the
        // blob scale that reads as patches, not the district-scale voxel_relief band.
        // Steps stay well under 64 m and non-aligned to the octave wavelengths.
        let seed = 7;
        let vals: Vec<f64> = (0..16)
            .map(|i| voxel_mosaic(seed, i as f64 * 3.7, i as f64 * 2.9))
            .collect();
        let min = vals.iter().cloned().fold(f64::INFINITY, f64::min);
        let max = vals.iter().cloned().fold(f64::NEG_INFINITY, f64::max);
        assert!(
            max - min > 0.1,
            "mosaic field is near-flat inside a chunk: range {}",
            max - min
        );
    }

    #[test]
    fn voxel_mosaic_is_smooth_not_salt_and_pepper() {
        // C¹ value noise → a small position step is a small field change, so patches
        // are coherent blobs, never the per-tile scatter the mosaic replaces.
        let base = voxel_mosaic(11, 5_000.0, 5_000.0);
        let near = voxel_mosaic(11, 5_000.5, 5_000.0);
        assert!(
            (base - near).abs() < 0.05,
            "0.5 m step jumped by {}",
            (base - near).abs()
        );
    }
}