facett-map3d-flight/src/origin.rs
//! **Floating origin** — the foundational fix for *far-from-origin vertex jitter*
//! at planetary scale (FLIGHT-SIM #1).
//!
//! ## The problem
//! The world is WGS-84. Projected to Web-Mercator **metres** (EPSG:3857) a tile
//! near the antimeridian / high latitude sits ~20 000 000 m from the projection
//! origin. An `f32` has ~24 bits of mantissa → at 2×10⁷ m its resolution is
//! ≈ **2 m**. Every vertex of that tile therefore *snaps* to a ~2 m grid: the
//! building edges shimmer (Z-fight) and crawl (jitter) as the camera moves, and a
//! sub-metre DEM is impossible. This is the classic "the planet is too big for
//! `f32`" failure that every AAA planetary engine solves with a **floating
//! origin**.
//!
//! ## The fix (what this module is)
//! Keep the entire world simulation in **`f64` metres** ([`WorldPos`]). Each frame,
//! pick a **camera-centred origin** ([`FloatingOrigin::rebased`] from the camera's
//! `f64` world position) and upload every vertex as `f32` **relative to that
//! origin**. Because the subtraction happens in `f64` *before* the cast, the
//! `f32` only ever has to represent the *small* camera-relative offset (metres, not
//! ten-millions of metres), so it keeps **sub-millimetre** precision regardless of
//! where on Earth the camera is. The jitter is gone.
//!
//! Pure `f64`/`f32` math — no GPU, no egui — so the precision win is a directly
//! inject-assertable unit test ([`tests`]).
/// Web-Mercator (EPSG:3857) **equatorial radius**, metres. The projection scale
/// constant: `x = R·λ`, `y = R·ln(tan(π/4 + φ/2))`.
pub const EARTH_RADIUS_M: f64 = 6_378_137.0;
/// Web-Mercator latitude clamp (the projection diverges at the poles). Matches
/// [`facett_map::mercator`]'s `LAT_LIMIT`.
pub const LAT_LIMIT_DEG: f64 = 85.051_128_78;
/// A position in the **`f64` world frame**: Web-Mercator metres (x east, y north,
/// z up = altitude). This is the canonical sim-space — never truncated to `f32`
/// until it has been rebased relative to the camera origin.
#[derive(Clone, Copy, Debug, PartialEq)]
pub struct WorldPos {
/// Web-Mercator easting, metres from the prime meridian.
pub x: f64,
/// Web-Mercator northing, metres (north +). NB: standard EPSG:3857 has y north
/// +; the design-plane mapping flips to its own convention downstream.
pub y: f64,
/// Altitude / height, metres above the ellipsoid (the extrusion axis).
pub z: f64,
}
impl WorldPos {
#[must_use]
pub const fn new(x: f64, y: f64, z: f64) -> Self {
Self { x, y, z }
}
/// Project WGS-84 `(lon, lat)` degrees + an altitude (metres) into the `f64`
/// Web-Mercator world frame. This is the **spherical** Web-Mercator used by
/// every slippy-map tile server, so `(x, y)` here line up with OSM/GeoParquet
/// tile coordinates. Latitude is clamped to the Mercator limit.
#[must_use]
pub fn from_lonlat(lon_deg: f64, lat_deg: f64, altitude_m: f64) -> Self {
let lat = lat_deg.clamp(-LAT_LIMIT_DEG, LAT_LIMIT_DEG);
let lon_r = lon_deg.to_radians();
let lat_r = lat.to_radians();
let x = EARTH_RADIUS_M * lon_r;
// Mercator northing: R·ln(tan(π/4 + φ/2)) = R·asinh(tan φ).
let y = EARTH_RADIUS_M * (lat_r.tan() + 1.0 / lat_r.cos()).ln();
Self { x, y, z: altitude_m }
}
/// Inverse of [`WorldPos::from_lonlat`]: recover WGS-84 `(lon, lat)` degrees
/// from the Web-Mercator easting/northing. Altitude (`z`) is ignored. This is the
/// exact analytic inverse of the spherical Mercator forward map
/// (`lon = x/R`, `lat = atan(sinh(y/R))`), so `from_lonlat(lon, lat, _).to_lonlat()`
/// round-trips to full `f64` precision for any in-range latitude.
#[must_use]
pub fn to_lonlat(self) -> (f64, f64) {
let lon = (self.x / EARTH_RADIUS_M).to_degrees();
let lat = (self.y / EARTH_RADIUS_M).sinh().atan().to_degrees();
(lon, lat)
}
/// Component-wise `self - other` in `f64` (stays in the world frame).
#[must_use]
#[allow(clippy::should_implement_trait)] // public by-value world-frame op; a `Sub` impl would not replace the published method
pub fn sub(self, o: WorldPos) -> WorldPos {
WorldPos::new(self.x - o.x, self.y - o.y, self.z - o.z)
}
/// Horizontal (x,y) distance to another world position, metres (`f64`).
#[must_use]
pub fn ground_distance(self, o: WorldPos) -> f64 {
let dx = self.x - o.x;
let dy = self.y - o.y;
(dx * dx + dy * dy).sqrt()
}
}
/// The **per-frame floating origin**: the `f64` world point that maps to the
/// renderer's `(0,0,0)`. Rebuilt every frame from the camera so the camera (and
/// everything near it) sits at small coordinates the `f32` GPU pipeline renders
/// without jitter.
#[derive(Clone, Copy, Debug, PartialEq)]
pub struct FloatingOrigin {
/// The world point that is the renderer origin (usually the camera, snapped).
pub anchor: WorldPos,
}
impl FloatingOrigin {
/// A floating origin anchored exactly at a world position (typically the camera
/// eye). Everything is uploaded relative to this.
#[must_use]
pub const fn at(anchor: WorldPos) -> Self {
Self { anchor }
}
/// **Snapped** floating origin: the anchor is quantised to a `grid_m`-metre
/// lattice so the origin only jumps in discrete steps as the camera moves. A
/// *continuously* moving origin would re-translate *every* vertex every frame
/// (defeating any vertex caching and adding its own sub-pixel crawl); snapping
/// to a coarse grid keeps the camera-relative offsets bounded (≤ `grid_m`) and
/// the rebase stable frame-to-frame. `grid_m ≤ 0` ⇒ no snapping (anchor == eye).
#[must_use]
pub fn snapped(eye: WorldPos, grid_m: f64) -> Self {
if grid_m <= 0.0 {
return Self::at(eye);
}
let q = |v: f64| (v / grid_m).round() * grid_m;
Self::at(WorldPos::new(q(eye.x), q(eye.y), q(eye.z)))
}
/// **The rebase**: a world position expressed as `f32` **relative to the
/// anchor**. The subtraction is done in `f64` (full precision) and only the
/// *small* residual is cast to `f32`, so a vertex 20 000 km from the projection
/// origin but 50 m from the camera lands at `≈[…, 50.0, …]` with sub-mm error —
/// the whole point of the floating origin.
#[must_use]
pub fn rebase(self, world: WorldPos) -> [f32; 3] {
let r = world.sub(self.anchor);
[r.x as f32, r.y as f32, r.z as f32]
}
/// Rebase a whole slice of world positions into camera-relative `f32` — the
/// per-frame "upload geometry relative to the camera before the wgpu buffers"
/// step. Allocation-light helper for the streaming tile path.
#[must_use]
pub fn rebase_all(self, world: &[WorldPos]) -> Vec<[f32; 3]> {
world.iter().map(|&w| self.rebase(w)).collect()
}
}
#[cfg(test)]
mod tests {
use super::*;
/// The naive "cast world metres straight to `f32`" path — the bug the floating
/// origin fixes. Used by the precision test as the negative control.
fn naive_relative_f32(world: WorldPos, camera: WorldPos) -> [f32; 3] {
// Cast FIRST (lossy at planetary magnitude), subtract in f32.
let wx = world.x as f32;
let wy = world.y as f32;
let cx = camera.x as f32;
let cy = camera.y as f32;
[wx - cx, wy - cy, (world.z - camera.z) as f32]
}
/// **INJECT-ASSERT (the headline precision win):** a tile ~20 000 km from the
/// Mercator origin, with the camera 50 m away, must rebase to `f32` coords that
/// are still accurate to **well under a millimetre** — while the naive
/// cast-then-subtract path loses *metres* of precision (the jitter). We compute
/// the exact `f64` offset as ground truth and compare both paths to it.
#[test]
fn far_tile_stays_precise_after_camera_relative_rebase() {
// A far-from-origin spot: high latitude near the antimeridian.
let cam_world = WorldPos::from_lonlat(179.5, 70.0, 1000.0);
// The far magnitude is real (tens of millions of metres out).
assert!(cam_world.x.abs() > 1.9e7, "camera is far from origin (x={})", cam_world.x);
// A cluster of vertices a few metres from the camera, with *fractional*
// offsets — the f32 grid step at this magnitude is ~2 m, so a fractional
// offset is exactly what the naive cast quantises away (the jitter). (A
// whole-number offset can land on the grid and look fine by luck; the
// sub-metre detail is what crawls.)
let origin = FloatingOrigin::at(cam_world);
let offsets = [(50.4f64, 30.7, 12.3), (0.9, -0.9, 0.4), (1.7, 3.3, -0.6), (120.6, -77.1, 5.5)];
let mut worst_good = 0.0f32;
let mut worst_naive = 0.0f32;
for &(dx, dy, dz) in &offsets {
let vtx_world = WorldPos::new(cam_world.x + dx, cam_world.y + dy, cam_world.z + dz);
let truth = [dx as f32, dy as f32, dz as f32];
// Floating origin: subtract in f64, THEN cast → sub-mm.
let good = origin.rebase(vtx_world);
worst_good = worst_good.max((0..3).map(|i| (good[i] - truth[i]).abs()).fold(0.0f32, f32::max));
// Naive: cast the 2×10⁷ m magnitudes to f32 first → metre-scale error.
let bad = naive_relative_f32(vtx_world, cam_world);
worst_naive = worst_naive.max((0..2).map(|i| (bad[i] - truth[i]).abs()).fold(0.0f32, f32::max));
}
assert!(worst_good < 1.0e-3, "floating-origin rebase keeps sub-mm precision, worst err={worst_good} m");
assert!(
worst_naive > 0.5,
"the naive cast-then-subtract path loses >0.5 m (the jitter) — worst err={worst_naive} m"
);
// The floating origin is orders of magnitude tighter than the naive cast.
assert!(worst_good * 100.0 < worst_naive, "floating origin is orders-of-magnitude tighter");
}
/// **INJECT-ASSERT (precision holds across the globe):** sweep a spread of
/// far-flung camera anchors; at each, a cluster of nearby vertices must rebase
/// to within a millimetre of their true camera-relative offset. Proves the win
/// is not specific to one lucky coordinate.
#[test]
fn rebase_precision_holds_across_a_global_sweep() {
let spots = [
(-179.9, -84.0),
(179.9, 84.0),
(139.69, 35.68), // Tokyo
(-122.42, 37.77), // SF
(18.07, 59.33), // Stockholm
(0.0, 0.0),
];
for &(lon, lat) in &spots {
let eye = WorldPos::from_lonlat(lon, lat, 500.0);
let origin = FloatingOrigin::snapped(eye, 64.0);
// The snapped origin is within one grid cell of the eye (so camera-local
// offsets stay small + bounded).
let cam_rel = origin.rebase(eye);
let cam_off = cam_rel.iter().map(|v| v.abs()).fold(0.0, f32::max);
assert!(cam_off <= 64.0 + 1e-3, "snapped origin keeps the camera within a grid cell ({cam_off} m)");
for (k, &(dx, dy, dz)) in [(7.5, -3.2, 4.0), (250.0, 110.0, 0.0), (-90.0, 600.0, 25.0)].iter().enumerate() {
let v = WorldPos::new(eye.x + dx, eye.y + dy, eye.z + dz);
let r = origin.rebase(v);
// True offset relative to the SNAPPED anchor.
let truth = v.sub(origin.anchor);
let err = ((r[0] as f64 - truth.x).abs())
.max((r[1] as f64 - truth.y).abs())
.max((r[2] as f64 - truth.z).abs());
assert!(err < 1.0e-3, "vertex {k} at {lon},{lat} rebases sub-mm (err={err} m)");
}
}
}
/// **INJECT-ASSERT (the world transform is the real Web-Mercator):** projecting
/// `(lon, lat)` agrees with [`facett_map::mercator`]'s normalised `[0,1]²` form
/// scaled by the Earth circumference — i.e. the `f64` metric world frame and the
/// existing map projection describe the SAME globe (so OSM tiles line up).
#[test]
fn world_transform_matches_the_map_mercator() {
let circ = 2.0 * std::f64::consts::PI * EARTH_RADIUS_M; // full Mercator world span
for &(lon, lat) in &[(9.5, 47.1), (-122.42, 37.77), (139.69, 35.68)] {
let w = WorldPos::from_lonlat(lon, lat, 0.0);
let (mx, my) = facett_map::mercator(lon, lat); // normalised [0,1], y south
// Our x is metres from the prime meridian; the normalised x is from the
// antimeridian (lon=-180), so shift by half a world.
let want_x = (mx - 0.5) * circ;
// Our y is north-up metres; the normalised y is south-down from the top.
let want_y = (0.5 - my) * circ;
assert!((w.x - want_x).abs() < 1.0, "x metres match map mercator at {lon},{lat}");
assert!((w.y - want_y).abs() < 1.0, "y metres match map mercator at {lon},{lat}");
}
}
/// **INJECT-ASSERT (the forward projection hits hand-computed Mercator metres):**
/// pin the spherical Web-Mercator forward map to closed-form reference values so a
/// regression in the constants/formula is caught, not just relative agreement.
/// * the prime-meridian/equator point is the exact origin `(0,0)`;
/// * a full ±180° of longitude is ±π·R metres (the world half-span, `WORLD_HALF`);
/// * easting is exactly linear in longitude (`x = R·λ`);
/// * the northing at the Mercator latitude limit is the same π·R half-span (the
/// projection is square), and latitude beyond the limit is clamped, not NaN.
#[test]
fn forward_projection_hits_known_mercator_values() {
const WORLD_HALF: f64 = std::f64::consts::PI * EARTH_RADIUS_M; // π·R ≈ 20 037 508 m
// Equator × prime meridian is the exact projection origin.
let o = WorldPos::from_lonlat(0.0, 0.0, 0.0);
assert!(o.x == 0.0 && o.y == 0.0, "equator/prime-meridian is (0,0), got ({},{})", o.x, o.y);
// ±180° longitude → ±π·R (the world half-span). Antimeridian is symmetric.
let east = WorldPos::from_lonlat(180.0, 0.0, 0.0);
let west = WorldPos::from_lonlat(-180.0, 0.0, 0.0);
assert!((east.x - WORLD_HALF).abs() < 1e-6, "lon +180 → +π·R, got {}", east.x);
assert!((west.x + WORLD_HALF).abs() < 1e-6, "lon −180 → −π·R, got {}", west.x);
// x is exactly linear in longitude: x = R·λ (metres per degree = R·π/180).
let per_deg = EARTH_RADIUS_M * std::f64::consts::PI / 180.0;
for lon in [-123.4, -17.0, 1.0, 45.0, 137.9] {
let x = WorldPos::from_lonlat(lon, 12.3, 0.0).x;
assert!((x - lon * per_deg).abs() < 1e-6, "x is linear in lon at {lon}° (x={x})");
}
// At the Mercator latitude limit the northing equals the world half-span (the
// projected world is the square [-π·R, π·R]²).
let top = WorldPos::from_lonlat(0.0, LAT_LIMIT_DEG, 0.0);
assert!((top.y - WORLD_HALF).abs() < 1.0, "lat limit → +π·R northing, got {}", top.y);
// Latitude past the limit is clamped (finite), not run to the divergent pole.
let polar = WorldPos::from_lonlat(0.0, 89.9, 0.0);
assert!(polar.y.is_finite(), "over-limit latitude is clamped, not NaN/inf (y={})", polar.y);
assert!((polar.y - top.y).abs() < 1e-6, "89.9° clamps to the same northing as the limit");
}
/// **INJECT-ASSERT (lon/lat → Web-Mercator → lon/lat round-trips):** the analytic
/// inverse [`WorldPos::to_lonlat`] recovers the input to full `f64` precision across
/// a global sweep — equator/prime-meridian zero point, the antimeridian, and high
/// latitudes — and an over-limit latitude round-trips to the *clamped* limit (the
/// only lossy case, and it is lossy by design, not by error).
#[test]
fn lonlat_round_trips_through_web_mercator() {
for &(lon, lat) in &[
(0.0, 0.0), // the zero point
(18.0686, 59.3293), // Stockholm
(-179.999, 84.9), // near the antimeridian, high north
(179.999, -84.9), // near the antimeridian, high south
(-122.42, 37.77), // SF
(139.69, 35.68), // Tokyo
] {
let (rlon, rlat) = WorldPos::from_lonlat(lon, lat, 250.0).to_lonlat();
assert!((rlon - lon).abs() < 1e-9, "lon round-trips at {lon},{lat} (got {rlon})");
assert!((rlat - lat).abs() < 1e-9, "lat round-trips at {lon},{lat} (got {rlat})");
}
// An over-limit latitude is clamped on the way out, so it round-trips to the
// Mercator limit — exactly, not approximately.
let (rlon, rlat) = WorldPos::from_lonlat(10.0, 88.0, 0.0).to_lonlat();
assert!((rlon - 10.0).abs() < 1e-9, "longitude still round-trips past the lat clamp");
assert!((rlat - LAT_LIMIT_DEG).abs() < 1e-6, "over-limit lat clamps to the Mercator limit (got {rlat})");
}
/// **INJECT-ASSERT (the "~2 m f32 grid" claim in the module docs, as data):** at the
/// far edge of Web-Mercator (~2×10⁷ m) the spacing between representable `f32`
/// values — one ULP — really is ~2 m, so world coordinates cast straight to `f32`
/// snap to a ~2 m lattice. This is the concrete magnitude the floating origin exists
/// to defeat; the rebase tests prove it does.
#[test]
fn f32_resolution_at_planetary_scale_is_metres() {
// A coordinate out near the antimeridian world edge.
let far = WorldPos::from_lonlat(179.9, 0.0, 0.0).x;
assert!(far > 1.9e7, "the sample really is far from the origin (x={far})");
// The next representable f32 above `far as f32` is one ULP away.
let a = far as f32;
let ulp = f32::from_bits(a.to_bits() + 1) - a;
assert!((1.0..=4.0).contains(&ulp), "f32 ULP at ~2×10⁷ m is metre-scale (~2 m), got {ulp} m");
// Concretely: a sub-metre detail (+0.4 m) is *invisible* to a straight f32 cast
// at this magnitude — it lands back on the same grid point (the jitter source).
let snapped = (far + 0.4) as f32;
assert_eq!(snapped, a, "a 0.4 m detail vanishes in a raw f32 cast at planetary scale");
}
/// INJECT-ASSERT: origin snapping is deterministic + idempotent — two cameras a
/// few metres apart inside the same grid cell pick the SAME anchor (so vertex
/// offsets don't crawl every frame), and re-snapping an anchor is a no-op.
#[test]
fn snapped_origin_is_stable_within_a_cell() {
let a = WorldPos::from_lonlat(10.0, 50.0, 100.0);
let b = WorldPos::new(a.x + 10.0, a.y - 12.0, a.z); // < 64 m away
let oa = FloatingOrigin::snapped(a, 64.0);
let ob = FloatingOrigin::snapped(b, 64.0);
assert_eq!(oa.anchor, ob.anchor, "nearby cameras share the snapped origin");
// Re-snapping the (already grid-aligned) anchor is idempotent.
let re = FloatingOrigin::snapped(oa.anchor, 64.0);
assert_eq!(re.anchor, oa.anchor, "snapping is idempotent on a grid point");
}
}