facett-map3d-flight/src/terrain.rs
//! **Flight terrain mesh** — builds the flight engine's 3D terrain from real
//! elevations (roadmap §6): breaks the flat `y = 0` ground plane by sampling a
//! Digital Elevation Model at every terrain vertex and lifting it to true height.
//!
//! The **elevation source** ([`Dem`], the [`ProceduralTerrain`] fallback, the
//! [`TerrainField`] seam, the `.hgt` loader / [`parse_hgt_name`]) is **not forked
//! here** — it is the single implementation in [`facett_map3d::terrain`], which
//! flight already depends on; this module `pub use`s it so `flight::terrain::Dem`
//! (and the crate-root re-exports) keep resolving to the one type. What is genuinely
//! flight-specific stays here:
//! - [`GeoRegion`] — the WGS-84 rectangle the mesh is built over.
//! - [`build_terrain_mesh`] — samples elevation on an `nx × ny` grid over a
//! [`GeoRegion`], lifts each vertex in the **`f64` world frame**, **rebases**
//! through the [`FloatingOrigin`] (jitter-free planetary precision), and computes
//! per-vertex **normals** → a [`TerrainMesh`] (`#[cfg(feature = "egui")]` it
//! converts to the shared [`facett_map3d::Mesh3`], reusing the wgpu pipeline).
use crate::origin::{FloatingOrigin, WorldPos};
// The DEM elevation source is the single home in `facett_map3d::terrain`; re-exported
// so flight's own `terrain::{Dem,Elevation,TerrainField,…}` paths keep working with
// one implementation (audit cluster C — the ~340-LOC verbatim fork is deleted here).
pub use facett_map3d::terrain::{Dem, Elevation, HGT_VOID, ProceduralTerrain, TerrainField, parse_hgt_name};
// ───────────────────────────────────────────────────────────────────────────────
// GeoRegion + TerrainMesh — build the 3D terrain from real heights.
// ───────────────────────────────────────────────────────────────────────────────
/// A WGS-84 lon/lat rectangle the terrain mesh is built over.
#[derive(Clone, Copy, Debug, PartialEq)]
pub struct GeoRegion {
/// West edge longitude (°).
pub min_lon: f64,
/// South edge latitude (°).
pub min_lat: f64,
/// East edge longitude (°).
pub max_lon: f64,
/// North edge latitude (°).
pub max_lat: f64,
}
impl GeoRegion {
#[must_use]
pub const fn new(min_lon: f64, min_lat: f64, max_lon: f64, max_lat: f64) -> Self {
Self { min_lon, min_lat, max_lon, max_lat }
}
/// The whole Sweden extent (matches [`crate::sweden::SwedenExtent`]).
#[must_use]
pub const fn sweden() -> Self {
use crate::sweden::SwedenExtent;
Self::new(SwedenExtent::MIN_LON, SwedenExtent::MIN_LAT, SwedenExtent::MAX_LON, SwedenExtent::MAX_LAT)
}
/// The centre `(lon, lat)` of the region.
#[must_use]
pub fn center(&self) -> (f64, f64) {
((self.min_lon + self.max_lon) * 0.5, (self.min_lat + self.max_lat) * 0.5)
}
}
/// The **3D terrain mesh** built from real elevations: positions + per-vertex normals
/// (parallel) + triangle indices, plus the axis-aligned bounds. Positions are in the
/// renderer's `f32` frame **relative to a [`FloatingOrigin`]** (metres), so the terrain
/// slots straight into the jitter-free planetary pipeline. Colours/materials are the
/// render backend's job (terrain shades by height/slope in the shader), so this stays
/// a lean, egui-free, headless-testable geometry payload; [`TerrainMesh::to_mesh3`]
/// (feature `egui`) converts it to the shared [`facett_map3d::Mesh3`] for the reused
/// wgpu path.
#[derive(Clone, Debug, Default, PartialEq)]
pub struct TerrainMesh {
/// `(x, y, z)` per vertex, metres relative to the floating origin (`y` = up/height).
pub positions: Vec<[f32; 3]>,
/// Unit surface normal per vertex (parallel to `positions`), for lighting.
pub normals: Vec<[f32; 3]>,
/// Triangle indices (3 per face) into `positions`.
pub indices: Vec<u32>,
/// Grid columns (vertices west→east).
pub nx: usize,
/// Grid rows (vertices south→north).
pub ny: usize,
/// Axis-aligned minimum corner of `positions`.
pub min: [f32; 3],
/// Axis-aligned maximum corner of `positions`.
pub max: [f32; 3],
}
impl TerrainMesh {
/// Vertex count (`nx · ny`).
#[must_use]
pub fn vertex_count(&self) -> usize {
self.positions.len()
}
/// Triangle count (`(nx−1)·(ny−1)·2`).
#[must_use]
pub fn triangle_count(&self) -> usize {
self.indices.len() / 3
}
/// Convert to the shared [`facett_map3d::Mesh3`] (reused wgpu pipeline), painting
/// every triangle a single base `color` (real shading — height/slope tint — is a
/// shader concern the backend applies). Feature-gated on `egui` because `Mesh3`'s
/// per-triangle colour is an `egui::Color32`.
#[cfg(feature = "egui")]
#[must_use]
pub fn to_mesh3(&self, color: egui::Color32) -> facett_map3d::Mesh3 {
facett_map3d::Mesh3 {
positions: self.positions.clone(),
normals: self.normals.clone(),
colors: vec![color; self.triangle_count()],
indices: self.indices.clone(),
// Terrain is never a house — see `facett_map3d::extrude::Mesh3::building`.
building: vec![false; self.positions.len()],
}
}
}
/// **Build the 3D terrain mesh from real heights** (roadmap §6): sample `elev` on an
/// `nx × ny` vertex grid over `region`, lift each vertex to its elevation (× the
/// `vertical_exaggeration` scale) in the `f64` world frame, **rebase** through
/// `origin` to the renderer's `f32` frame, and compute per-vertex **normals** by
/// central differences of the neighbouring world positions (so lighting reads the real
/// slope). This is the replacement for the flat `y = 0` quad.
///
/// `nx, ny ≥ 2`. Vertices are laid out row 0 = **south**, row `ny−1` = **north**,
/// columns west→east. Triangles are wound CCW seen from above (front-facing up).
///
/// # Panics
/// If `nx < 2` or `ny < 2`.
#[must_use]
pub fn build_terrain_mesh(
elev: &impl Elevation,
region: GeoRegion,
nx: usize,
ny: usize,
origin: FloatingOrigin,
vertical_exaggeration: f64,
) -> TerrainMesh {
assert!(nx >= 2 && ny >= 2, "terrain grid needs ≥2×2 vertices");
let n = nx * ny;
let mut positions = Vec::with_capacity(n);
let idx = |i: usize, j: usize| j * nx + i;
// Vertices: row j maps south→north, column i west→east. The renderer's frame is
// **Y-up** (the flight camera's `up` basis is +Y and roadmap §6 displaces terrain
// along `y`), whereas the `f64` [`WorldPos`] frame is (x=easting, y=northing,
// z=altitude). So we rebase in the world frame, then map to render axes
// `[easting, elevation, northing]` — elevation is the vertical (Y) axis the camera
// and the shared `Mesh3`/wgpu pipeline expect.
for j in 0..ny {
let ty = j as f64 / (ny - 1) as f64;
let lat = region.min_lat + (region.max_lat - region.min_lat) * ty;
for i in 0..nx {
let tx = i as f64 / (nx - 1) as f64;
let lon = region.min_lon + (region.max_lon - region.min_lon) * tx;
let h = f64::from(elev.elevation_at(lon, lat)) * vertical_exaggeration;
let r = origin.rebase(WorldPos::from_lonlat(lon, lat, h)); // [easting, northing, elev]
positions.push([r[0], r[2], r[1]]); // → Y-up render frame [easting, elev, northing]
}
}
// Per-vertex normals from central differences of the neighbour positions. Using the
// already-rebased local metres keeps the differences small + precise; we force the
// normal to point up (+y) so a consistent winding isn't required for correct lighting.
let mut normals = vec![[0.0f32, 1.0, 0.0]; n];
for j in 0..ny {
for i in 0..nx {
let l = positions[idx(i.saturating_sub(1), j)];
let r = positions[idx((i + 1).min(nx - 1), j)];
let d = positions[idx(i, j.saturating_sub(1))];
let u = positions[idx(i, (j + 1).min(ny - 1))];
let du = [r[0] - l[0], r[1] - l[1], r[2] - l[2]]; // east tangent
let dv = [u[0] - d[0], u[1] - d[1], u[2] - d[2]]; // north tangent
// north × east points up for the south→north / west→east layout.
let mut nrm = [
dv[1] * du[2] - dv[2] * du[1],
dv[2] * du[0] - dv[0] * du[2],
dv[0] * du[1] - dv[1] * du[0],
];
let len = (nrm[0] * nrm[0] + nrm[1] * nrm[1] + nrm[2] * nrm[2]).sqrt();
if len > 1e-9 {
let s = if nrm[1] < 0.0 { -1.0 / len } else { 1.0 / len };
nrm = [nrm[0] * s, nrm[1] * s, nrm[2] * s];
normals[idx(i, j)] = nrm;
}
}
}
// Triangle indices (two per grid cell), CCW seen from above.
let mut indices = Vec::with_capacity((nx - 1) * (ny - 1) * 6);
for j in 0..ny - 1 {
for i in 0..nx - 1 {
let a = idx(i, j) as u32;
let b = idx(i + 1, j) as u32;
let c = idx(i, j + 1) as u32;
let d = idx(i + 1, j + 1) as u32;
indices.extend_from_slice(&[a, b, c, b, d, c]);
}
}
// Axis-aligned bounds.
let mut min = [f32::INFINITY; 3];
let mut max = [f32::NEG_INFINITY; 3];
for p in &positions {
for k in 0..3 {
min[k] = min[k].min(p[k]);
max[k] = max[k].max(p[k]);
}
}
TerrainMesh { positions, normals, indices, nx, ny, min, max }
}
#[cfg(test)]
mod tests {
use super::*;
/// **INJECT-ASSERT (`load_hgt_or_synthetic` falls back when no tile is present):**
/// pointing the loader at a nonexistent path yields the synthetic fallback (not an
/// error), so CI with no DEM tiles on disk is still green + deterministic.
#[test]
fn missing_tile_falls_back_to_synthetic() {
let (field, real) = TerrainField::load_hgt_or_synthetic("/no/such/tile/N59E017.hgt");
assert!(!real, "no real tile loaded");
assert!(!field.is_real(), "field is the synthetic fallback");
assert!(field.elevation_at(17.0, 59.0).is_finite());
}
/// **INJECT-ASSERT (mesh geometry from real heights):** building the terrain mesh
/// over a varied heightfield yields the exact vertex/normal/index counts for the
/// grid, unit normals that point **up**, finite bounds, and — the headline — a
/// **non-flat** surface (the `y` extent is non-zero), i.e. real topography, not the
/// old flat `y = 0` plane.
#[test]
fn mesh_generation_counts_normals_and_relief() {
let dem = Dem::synthetic(10.0, 60.0, 2.0, 2.0, 64, 64, |lon, lat| {
// A hill: metres, varying with position.
500.0 + 400.0 * (lon * 3.0).sin() as f32 * (lat * 2.0).cos() as f32
});
let region = GeoRegion::new(10.2, 58.2, 11.8, 59.8);
let (nx, ny) = (48usize, 40usize);
let origin = FloatingOrigin::snapped(WorldPos::from_lonlat(11.0, 59.0, 0.0), 256.0);
let mesh = build_terrain_mesh(&dem, region, nx, ny, origin, 1.0);
assert_eq!(mesh.vertex_count(), nx * ny, "vertex count = nx*ny");
assert_eq!(mesh.normals.len(), mesh.positions.len(), "one normal per vertex");
assert_eq!(mesh.triangle_count(), (nx - 1) * (ny - 1) * 2, "two triangles per cell");
assert_eq!(mesh.indices.len(), (nx - 1) * (ny - 1) * 6);
// Every index is in range.
assert!(mesh.indices.iter().all(|&i| (i as usize) < mesh.positions.len()));
// Normals are unit-length and point up.
for nrm in &mesh.normals {
let len = (nrm[0] * nrm[0] + nrm[1] * nrm[1] + nrm[2] * nrm[2]).sqrt();
assert!((len - 1.0).abs() < 1e-4, "unit normal (len={len})");
assert!(nrm[1] > 0.0, "normal points up (y={})", nrm[1]);
}
// Real relief: the terrain is NOT flat (y extent > 1 m).
assert!(mesh.max[1] - mesh.min[1] > 1.0, "terrain has vertical relief (Δy={})", mesh.max[1] - mesh.min[1]);
// Bounds are finite.
assert!(mesh.min.iter().chain(mesh.max.iter()).all(|v| v.is_finite()));
}
/// **INJECT-ASSERT (a flat DEM ⇒ the degenerate flat plane):** a constant-height
/// DEM builds a mesh with (near-)zero vertical relief and all normals straight up —
/// the old flat-ground behaviour falls out as the zero-relief special case, proving
/// the height is what drives the topography.
#[test]
fn flat_dem_reproduces_the_flat_plane() {
let dem = Dem::from_grid(10.0, 60.0, 1.0, 1.0, 2, 2, vec![100.0; 4]);
let origin = FloatingOrigin::at(WorldPos::from_lonlat(10.5, 59.5, 0.0));
let mesh = build_terrain_mesh(&dem, GeoRegion::new(10.1, 59.1, 10.9, 59.9), 16, 16, origin, 1.0);
assert!(mesh.max[1] - mesh.min[1] < 1e-3, "constant height ⇒ flat (Δy={})", mesh.max[1] - mesh.min[1]);
for nrm in &mesh.normals {
assert!(nrm[1] > 0.999, "flat ground normals point straight up (y={})", nrm[1]);
}
}
/// **INJECT-ASSERT (terrain rides the floating origin):** the same DEM region built
/// against a camera-snapped origin lands near the renderer's `(0,0,0)` (small
/// coordinates) — the terrain inherits the planetary jitter fix instead of carrying
/// 2×10⁷ m Mercator magnitudes into `f32`.
#[test]
fn terrain_is_rebased_onto_the_floating_origin() {
let region = GeoRegion::new(17.9, 59.1, 18.1, 59.3); // near Stockholm
let cam = WorldPos::from_lonlat(18.0, 59.2, 3000.0);
let origin = FloatingOrigin::snapped(cam, 512.0);
let p = ProceduralTerrain::default();
let mesh = build_terrain_mesh(&p, region, 32, 32, origin, 1.0);
// The region is ~10 km on a side, so rebased coords are ≲ tens of km, NOT the
// ~2×10⁶ m raw Mercator northing at 59°N.
let worst = mesh
.positions
.iter()
.flat_map(|p| p.iter().map(|v| v.abs()))
.fold(0.0f32, f32::max);
assert!(worst < 50_000.0, "rebased terrain sits near the origin ({worst} m), not at planetary magnitude");
}
/// **INJECT-ASSERT (`TerrainField` dispatches to the real DEM):** the enum forwards
/// `elevation_at` to a loaded DEM, so callers query one type whether the source is
/// real or synthetic.
#[test]
fn terrain_field_dispatches_to_the_real_dem() {
let dem = Dem::from_grid(0.0, 1.0, 1.0, 1.0, 2, 2, vec![10.0, 20.0, 30.0, 40.0]);
let field = TerrainField::Dem(Box::new(dem.clone()));
assert!(field.is_real());
assert_eq!(field.elevation_at(0.25, 0.75), dem.elevation_at(0.25, 0.75));
}
}