nordisk/facett

facett-core/src/engine/rte.rs

//! **Relative-to-eye (RTE) precision** — one fix for two jitter bugs.
//!
//! WGSL has no `f64`. A Web-Mercator coordinate at zoom 22 needs ~24 bits of
//! fraction *below* a value near 1.0, and an `f32` has 24 bits of mantissa
//! total — so the vertex shader quantises the world to a visible lattice and
//! roads shimmer. A graph drilled six levels deep has the identical failure with
//! different units.
//!
//! The fix is to never hand the shader an absolute coordinate. Split each `f64`
//! on the CPU into a **high** `f32` (the eye's neighbourhood) and a **low** `f32`
//! (the residual), and evaluate *relative to the eye* in the shader:
//!
//! ```text
//! P_eye = (P_high − C_eye_high) + (P_low − C_eye_low)
//! ```
//!
//! `P_high − C_eye_high` cancels the large magnitude exactly (both are the same
//! rounded lattice), leaving a small number the `f32` mantissa can carry in
//! full. This module is the CPU half; the WGSL half is two subtractions.

use super::source::WorldPos;

/// A `f64` split into two `f32`s whose sum reproduces it to `f32`-of-the-residual
/// precision: `value ≈ high + low`.
#[derive(Clone, Copy, Debug, Default, PartialEq)]
pub struct Split {
    pub high: f32,
    pub low: f32,
}

/// Split one `f64`. `high` is the `f32`-rounded value; `low` carries the part
/// that rounding threw away.
#[inline]
#[must_use]
pub fn split(v: f64) -> Split {
    let high = v as f32;
    let low = (v - high as f64) as f32;
    Split { high, low }
}

/// Reassemble — exact inverse of [`split`] up to the residual's own rounding.
#[inline]
#[must_use]
pub fn join(s: Split) -> f64 {
    s.high as f64 + s.low as f64
}

/// A world position in the form the vertex shader consumes: three `vec2<f32>`
/// worth of `(high, low)` pairs, `bytemuck`-able as two `[f32; 3]`.
#[derive(Clone, Copy, Debug, Default, PartialEq)]
pub struct RtePos {
    pub high: [f32; 3],
    pub low: [f32; 3],
}

impl RtePos {
    #[must_use]
    pub fn of(p: WorldPos) -> Self {
        let (sx, sy, sz) = (split(p.x), split(p.y), split(p.z));
        Self { high: [sx.high, sy.high, sz.high], low: [sx.low, sy.low, sz.low] }
    }

    /// The CPU mirror of the shader's two subtractions — **the parity oracle**.
    /// The CPU renderer must compute eye-relative positions exactly this way or
    /// the two lanes drift apart at high zoom (and only at high zoom, where no
    /// smoke test looks).
    #[must_use]
    pub fn relative_to(self, eye: RtePos) -> [f32; 3] {
        [
            (self.high[0] - eye.high[0]) + (self.low[0] - eye.low[0]),
            (self.high[1] - eye.high[1]) + (self.low[1] - eye.low[1]),
            (self.high[2] - eye.high[2]) + (self.low[2] - eye.low[2]),
        ]
    }
}

/// The **naive** `f32` evaluation — subtract after collapsing to `f32`. Kept
/// public **only** so a guard can red-prove itself: a test asserts this one is
/// visibly worse than [`RtePos::relative_to`] at high zoom. Never use it to draw.
#[must_use]
pub fn naive_relative(p: WorldPos, eye: WorldPos) -> [f32; 3] {
    [(p.x as f32) - (eye.x as f32), (p.y as f32) - (eye.y as f32), (p.z as f32) - (eye.z as f32)]
}

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn split_join_round_trips_to_f64_precision() {
        for v in [0.0_f64, 1.0, -1.0, 0.123_456_789_012_345, 1e12, -3.75e-9] {
            let back = join(split(v));
            let err = (back - v).abs();
            // high carries ~24 bits, low the next ~24 -> ~2^-48 relative.
            assert!(err <= v.abs() * 1e-14 + 1e-300, "v={v} back={back} err={err}");
        }
    }

    /// RED-PROVEN: this test fails if `relative_to` is replaced by
    /// `naive_relative` — the whole point of RTE, asserted rather than asserted-about.
    #[test]
    fn rte_beats_naive_f32_at_mercator_zoom_22() {
        // Two points 1e-9 of Mercator unit apart near x = 0.55 — about 4 cm at
        // the equator, i.e. two consecutive vertices of a footpath at z22.
        let eye = WorldPos::flat(0.554_321_098_765_432_1, 0.321_098_765_432_109_8);
        let step = 1e-9_f64;
        let p = WorldPos::flat(eye.x + step, eye.y);

        let rte = RtePos::of(p).relative_to(RtePos::of(eye));
        let naive = naive_relative(p, eye);

        let rte_err = (rte[0] as f64 - step).abs();
        let naive_err = (naive[0] as f64 - step).abs();

        // The naive lane collapses the step to exactly zero: both coordinates
        // round to the same f32. That IS the jitter.
        assert_eq!(naive[0], 0.0, "naive f32 was expected to lose the step entirely");
        assert!(rte_err < step * 1e-3, "RTE error {rte_err} should be <0.1% of the {step} step");
        assert!(
            naive_err > rte_err * 1000.0,
            "RTE must be >1000x better: rte_err={rte_err} naive_err={naive_err}"
        );
    }

    #[test]
    fn rte_is_exact_at_the_eye() {
        let eye = WorldPos::new(1e9, -2e9, 3.5);
        let r = RtePos::of(eye).relative_to(RtePos::of(eye));
        assert_eq!(r, [0.0, 0.0, 0.0]);
    }
}