rustial-engine 0.0.1

//! Elevation grid data type and bilinear interpolation.
//!
//! # Overview
//!
//! An [`ElevationGrid`] is a regular rectangular grid of `f32` height
//! samples aligned to a single slippy-map [`TileId`].  Heights are in
//! **meters above the WGS-84 ellipsoid**.
//!
//! # Memory layout
//!
//! Samples are stored **row-major** (X varies fastest):
//!
//! ```text
//! index = y * width + x
//! ```
//!
//! The grid origin `(0, 0)` corresponds to the **north-west** corner of
//! the tile, matching the slippy-map convention where Y increases
//! southward.
//!
//! # UV convention
//!
//! The [`sample`](ElevationGrid::sample) method takes `(u, v)` in
//! `[0, 1]`:
//!
//! - `(0, 0)` = north-west corner of the tile
//! - `(1, 0)` = north-east corner
//! - `(0, 1)` = south-west corner
//! - `(1, 1)` = south-east corner
//!
//! Values outside `[0, 1]` are clamped, so queries near tile edges
//! never panic.
//!
//! # Consumers
//!
//! - `build_terrain_mesh()` in `rustial-engine` reads the grid to
//!   displace terrain vertices along Z.
//! - `TerrainManager::elevation_at()` delegates to
//!   [`sample_geo`](ElevationGrid::sample_geo) for point elevation
//!   queries.
//! - `FlatElevationSource` produces all-zero grids via
//!   [`flat`](ElevationGrid::flat).

use crate::coord::GeoCoord;
use crate::tile::{tile_to_geo, tile_xy_to_geo, TileId};

/// A regular grid of elevation samples aligned to a slippy-map tile.
///
/// Heights are in meters above the WGS-84 ellipsoid.  The grid is
/// always rectangular (`width * height` samples) and aligned to a
/// single [`TileId`].
///
/// Construct with [`flat`](Self::flat) for an all-zero grid or
/// [`from_data`](Self::from_data) to supply pre-decoded heights.
#[derive(Debug, Clone, PartialEq)]
#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
pub struct ElevationGrid {
    /// Number of samples along the X (east-west) axis.
    pub width: u32,
    /// Number of samples along the Y (north-south) axis.
    pub height: u32,
    /// Minimum elevation in the grid (meters).
    pub min_elev: f32,
    /// Maximum elevation in the grid (meters).
    pub max_elev: f32,
    /// Row-major elevation data (`width * height` samples, meters).
    ///
    /// Index `y * width + x` gives the sample at column `x`, row `y`,
    /// where row 0 is the northern-most row.
    pub data: Vec<f32>,
    /// The tile this grid is aligned to.
    pub tile: TileId,
}

impl ElevationGrid {
    /// Create a flat (all-zero) elevation grid for a tile.
    ///
    /// Useful as a no-op placeholder when terrain is disabled or while
    /// a real elevation tile is still loading.
    pub fn flat(tile: TileId, width: u32, height: u32) -> Self {
        Self {
            width,
            height,
            min_elev: 0.0,
            max_elev: 0.0,
            data: vec![0.0; (width * height) as usize],
            tile,
        }
    }

    /// Create an elevation grid from raw height data, computing min/max
    /// automatically.
    ///
    /// Returns `None` if `data.len() != width * height` (size mismatch).
    pub fn from_data(tile: TileId, width: u32, height: u32, mut data: Vec<f32>) -> Option<Self> {
        if data.len() != (width * height) as usize {
            return None;
        }
        let mut min_elev = f32::MAX;
        let mut max_elev = f32::MIN;
        for v in data.iter_mut() {
            *v = v.clamp(-500.0, 10_000.0);
            if *v < min_elev {
                min_elev = *v;
            }
            if *v > max_elev {
                max_elev = *v;
            }
        }
        Some(Self {
            width,
            height,
            min_elev,
            max_elev,
            data,
            tile,
        })
    }

    /// The elevation range (max - min) in meters.
    ///
    /// Returns 0.0 for flat grids.  Useful for terrain mesh LOD
    /// decisions: tiles with small elevation range can use coarser
    /// mesh resolution.
    #[inline]
    pub fn elevation_range(&self) -> f32 {
        self.max_elev - self.min_elev
    }

    /// Bilinear interpolation at fractional coordinates `(u, v)` in
    /// `[0, 1]`.
    ///
    /// - `(0, 0)` is the north-west corner of the tile.
    /// - `(1, 1)` is the south-east corner.
    ///
    /// Values outside `[0, 1]` are clamped to the grid boundary.
    ///
    /// Returns `None` if the grid is empty (`width == 0` or
    /// `height == 0`).
    pub fn sample(&self, u: f64, v: f64) -> Option<f32> {
        if self.data.is_empty() || self.width == 0 || self.height == 0 {
            return None;
        }

        let u = u.clamp(0.0, 1.0);
        let v = v.clamp(0.0, 1.0);

        let fx = u * (self.width - 1) as f64;
        let fy = v * (self.height - 1) as f64;

        let x0 = (fx.floor() as u32).min(self.width - 1);
        let y0 = (fy.floor() as u32).min(self.height - 1);
        let x1 = (x0 + 1).min(self.width - 1);
        let y1 = (y0 + 1).min(self.height - 1);

        let sx = (fx - x0 as f64) as f32;
        let sy = (fy - y0 as f64) as f32;

        let v00 = self.data[(y0 * self.width + x0) as usize];
        let v10 = self.data[(y0 * self.width + x1) as usize];
        let v01 = self.data[(y1 * self.width + x0) as usize];
        let v11 = self.data[(y1 * self.width + x1) as usize];

        let top = v00 * (1.0 - sx) + v10 * sx;
        let bot = v01 * (1.0 - sx) + v11 * sx;
        Some(top * (1.0 - sy) + bot * sy)
    }

    /// Sample elevation at a geographic coordinate.
    ///
    /// Converts the `GeoCoord` to tile-relative `(u, v)` using the
    /// tile's geographic bounds, then delegates to [`sample`](Self::sample).
    ///
    /// Returns `None` if:
    /// - The tile has zero geographic extent (degenerate).
    /// - The grid is empty.
    ///
    /// Coordinates outside the tile's bounds are clamped to the
    /// nearest edge (via the `[0, 1]` clamp in `sample`).
    pub fn sample_geo(&self, coord: &GeoCoord) -> Option<f32> {
        // NW corner = tile origin, SE corner = next tile diagonally.
        let nw = tile_to_geo(&self.tile);
        let se = tile_xy_to_geo(
            self.tile.zoom,
            self.tile.x as f64 + 1.0,
            self.tile.y as f64 + 1.0,
        );

        let lon_range = se.lon - nw.lon;
        let lat_range = nw.lat - se.lat;

        if lon_range.abs() < 1e-12 || lat_range.abs() < 1e-12 {
            return None;
        }

        let u = (coord.lon - nw.lon) / lon_range;
        // v is inverted: lat decreases southward but v increases.
        let v = (nw.lat - coord.lat) / lat_range;

        self.sample(u, v)
    }
}

// ---------------------------------------------------------------------------
// Tests
// ---------------------------------------------------------------------------

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

    // -- Construction ------------------------------------------------------

    #[test]
    fn flat_grid() {
        let grid = ElevationGrid::flat(TileId::new(0, 0, 0), 3, 3);
        assert_eq!(grid.data.len(), 9);
        assert_eq!(grid.min_elev, 0.0);
        assert_eq!(grid.max_elev, 0.0);
        assert_eq!(grid.sample(0.5, 0.5), Some(0.0));
    }

    #[test]
    fn from_data_wrong_size() {
        assert!(ElevationGrid::from_data(TileId::new(0, 0, 0), 2, 2, vec![0.0]).is_none());
    }

    #[test]
    fn min_max_elev() {
        let data = vec![-100.0, 50.0, 200.0, 0.0];
        let grid = ElevationGrid::from_data(TileId::new(0, 0, 0), 2, 2, data).unwrap();
        assert!((grid.min_elev - (-100.0)).abs() < 1e-6);
        assert!((grid.max_elev - 200.0).abs() < 1e-6);
    }

    #[test]
    fn from_data_clamps_extreme_elevations() {
        let data = vec![-32_768.0, -600.0, 10_500.0, 25.0];
        let grid = ElevationGrid::from_data(TileId::new(0, 0, 0), 2, 2, data).unwrap();
        assert!((grid.min_elev - (-500.0)).abs() < 1e-6);
        assert!((grid.max_elev - 10_000.0).abs() < 1e-6);
        assert_eq!(grid.data, vec![-500.0, -500.0, 10_000.0, 25.0]);
    }

    #[test]
    fn elevation_range() {
        let data = vec![-100.0, 50.0, 200.0, 0.0];
        let grid = ElevationGrid::from_data(TileId::new(0, 0, 0), 2, 2, data).unwrap();
        assert!((grid.elevation_range() - 300.0).abs() < 1e-6);
    }

    #[test]
    fn elevation_range_flat() {
        let grid = ElevationGrid::flat(TileId::new(0, 0, 0), 4, 4);
        assert!((grid.elevation_range()).abs() < 1e-6);
    }

    // -- Bilinear interpolation (sample) -----------------------------------

    #[test]
    fn sample_corners() {
        let data = vec![0.0, 10.0, 20.0, 30.0];
        let grid = ElevationGrid::from_data(TileId::new(0, 0, 0), 2, 2, data).unwrap();
        assert!((grid.sample(0.0, 0.0).unwrap() - 0.0).abs() < 1e-6);
        assert!((grid.sample(1.0, 0.0).unwrap() - 10.0).abs() < 1e-6);
        assert!((grid.sample(0.0, 1.0).unwrap() - 20.0).abs() < 1e-6);
        assert!((grid.sample(1.0, 1.0).unwrap() - 30.0).abs() < 1e-6);
    }

    #[test]
    fn bilinear_midpoint() {
        let data = vec![0.0, 10.0, 20.0, 30.0];
        let grid = ElevationGrid::from_data(TileId::new(0, 0, 0), 2, 2, data).unwrap();
        let mid = grid.sample(0.5, 0.5).unwrap();
        // Average of all four corners: (0+10+20+30)/4 = 15.
        assert!((mid - 15.0).abs() < 1e-6);
    }

    #[test]
    fn sample_1x1_grid() {
        let grid = ElevationGrid::from_data(TileId::new(0, 0, 0), 1, 1, vec![42.0]).unwrap();
        // Every UV should return the single sample.
        assert!((grid.sample(0.0, 0.0).unwrap() - 42.0).abs() < 1e-6);
        assert!((grid.sample(0.5, 0.5).unwrap() - 42.0).abs() < 1e-6);
        assert!((grid.sample(1.0, 1.0).unwrap() - 42.0).abs() < 1e-6);
    }

    #[test]
    fn sample_clamps_out_of_range_uv() {
        let data = vec![0.0, 10.0, 20.0, 30.0];
        let grid = ElevationGrid::from_data(TileId::new(0, 0, 0), 2, 2, data).unwrap();
        // u < 0 clamps to u=0, v < 0 clamps to v=0.
        assert!((grid.sample(-1.0, -1.0).unwrap() - 0.0).abs() < 1e-6);
        // u > 1 clamps to u=1, v > 1 clamps to v=1.
        assert!((grid.sample(2.0, 2.0).unwrap() - 30.0).abs() < 1e-6);
    }

    #[test]
    fn sample_empty_grid_returns_none() {
        let grid = ElevationGrid {
            width: 0,
            height: 0,
            min_elev: 0.0,
            max_elev: 0.0,
            data: vec![],
            tile: TileId::new(0, 0, 0),
        };
        assert!(grid.sample(0.5, 0.5).is_none());
    }

    // -- Geographic sampling (sample_geo) ----------------------------------

    #[test]
    fn sample_geo_tile_center() {
        // Create a 2x2 grid on tile 0/0/0 (covers the whole world).
        let data = vec![100.0, 200.0, 300.0, 400.0];
        let grid = ElevationGrid::from_data(TileId::new(0, 0, 0), 2, 2, data).unwrap();

        // The center of tile 0/0/0 is approximately (0, 0) lat/lon.
        let center = GeoCoord::from_lat_lon(0.0, 0.0);
        let elev = grid.sample_geo(&center).unwrap();
        // At the midpoint of a 2x2 grid, expect the average: 250.
        assert!((elev - 250.0).abs() < 1.0);
    }

    #[test]
    fn sample_geo_nw_corner() {
        // Tile 1/0/0 covers the NW quadrant.
        let data = vec![10.0, 20.0, 30.0, 40.0];
        let grid = ElevationGrid::from_data(TileId::new(1, 0, 0), 2, 2, data).unwrap();

        // NW corner of tile 1/0/0 should return the (0,0) sample.
        let nw = tile_to_geo(&TileId::new(1, 0, 0));
        let elev = grid.sample_geo(&nw).unwrap();
        assert!((elev - 10.0).abs() < 1e-3);
    }

    // -- PartialEq ---------------------------------------------------------

    #[test]
    fn partial_eq() {
        let a = ElevationGrid::flat(TileId::new(5, 10, 10), 4, 4);
        let b = ElevationGrid::flat(TileId::new(5, 10, 10), 4, 4);
        assert_eq!(a, b);
    }
}