Skip to main content

rustial_engine/math/
elevation.rs

1//! Elevation grid data type and bilinear interpolation.
2//!
3//! # Overview
4//!
5//! An [`ElevationGrid`] is a regular rectangular grid of `f32` height
6//! samples aligned to a single slippy-map [`TileId`].  Heights are in
7//! **meters above the WGS-84 ellipsoid**.
8//!
9//! # Memory layout
10//!
11//! Samples are stored **row-major** (X varies fastest):
12//!
13//! ```text
14//! index = y * width + x
15//! ```
16//!
17//! The grid origin `(0, 0)` corresponds to the **north-west** corner of
18//! the tile, matching the slippy-map convention where Y increases
19//! southward.
20//!
21//! # UV convention
22//!
23//! The [`sample`](ElevationGrid::sample) method takes `(u, v)` in
24//! `[0, 1]`:
25//!
26//! - `(0, 0)` = north-west corner of the tile
27//! - `(1, 0)` = north-east corner
28//! - `(0, 1)` = south-west corner
29//! - `(1, 1)` = south-east corner
30//!
31//! Values outside `[0, 1]` are clamped, so queries near tile edges
32//! never panic.
33//!
34//! # Consumers
35//!
36//! - `build_terrain_mesh()` in `rustial-engine` reads the grid to
37//!   displace terrain vertices along Z.
38//! - `TerrainManager::elevation_at()` delegates to
39//!   [`sample_geo`](ElevationGrid::sample_geo) for point elevation
40//!   queries.
41//! - `FlatElevationSource` produces all-zero grids via
42//!   [`flat`](ElevationGrid::flat).
43
44use crate::coord::GeoCoord;
45use crate::tile::{tile_to_geo, tile_xy_to_geo, TileId};
46
47/// A regular grid of elevation samples aligned to a slippy-map tile.
48///
49/// Heights are in meters above the WGS-84 ellipsoid.  The grid is
50/// always rectangular (`width * height` samples) and aligned to a
51/// single [`TileId`].
52///
53/// Construct with [`flat`](Self::flat) for an all-zero grid or
54/// [`from_data`](Self::from_data) to supply pre-decoded heights.
55#[derive(Debug, Clone, PartialEq)]
56#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
57pub struct ElevationGrid {
58    /// Number of samples along the X (east-west) axis.
59    pub width: u32,
60    /// Number of samples along the Y (north-south) axis.
61    pub height: u32,
62    /// Minimum elevation in the grid (meters).
63    pub min_elev: f32,
64    /// Maximum elevation in the grid (meters).
65    pub max_elev: f32,
66    /// Row-major elevation data (`width * height` samples, meters).
67    ///
68    /// Index `y * width + x` gives the sample at column `x`, row `y`,
69    /// where row 0 is the northern-most row.
70    pub data: Vec<f32>,
71    /// The tile this grid is aligned to.
72    pub tile: TileId,
73}
74
75impl ElevationGrid {
76    /// Create a flat (all-zero) elevation grid for a tile.
77    ///
78    /// Useful as a no-op placeholder when terrain is disabled or while
79    /// a real elevation tile is still loading.
80    pub fn flat(tile: TileId, width: u32, height: u32) -> Self {
81        Self {
82            width,
83            height,
84            min_elev: 0.0,
85            max_elev: 0.0,
86            data: vec![0.0; (width * height) as usize],
87            tile,
88        }
89    }
90
91    /// Create an elevation grid from raw height data, computing min/max
92    /// automatically.
93    ///
94    /// Returns `None` if `data.len() != width * height` (size mismatch).
95    pub fn from_data(tile: TileId, width: u32, height: u32, mut data: Vec<f32>) -> Option<Self> {
96        if data.len() != (width * height) as usize {
97            return None;
98        }
99        let mut min_elev = f32::MAX;
100        let mut max_elev = f32::MIN;
101        for v in data.iter_mut() {
102            *v = v.clamp(-500.0, 10_000.0);
103            if *v < min_elev {
104                min_elev = *v;
105            }
106            if *v > max_elev {
107                max_elev = *v;
108            }
109        }
110        Some(Self {
111            width,
112            height,
113            min_elev,
114            max_elev,
115            data,
116            tile,
117        })
118    }
119
120    /// The elevation range (max - min) in meters.
121    ///
122    /// Returns 0.0 for flat grids.  Useful for terrain mesh LOD
123    /// decisions: tiles with small elevation range can use coarser
124    /// mesh resolution.
125    #[inline]
126    pub fn elevation_range(&self) -> f32 {
127        self.max_elev - self.min_elev
128    }
129
130    /// Bilinear interpolation at fractional coordinates `(u, v)` in
131    /// `[0, 1]`.
132    ///
133    /// - `(0, 0)` is the north-west corner of the tile.
134    /// - `(1, 1)` is the south-east corner.
135    ///
136    /// Values outside `[0, 1]` are clamped to the grid boundary.
137    ///
138    /// Returns `None` if the grid is empty (`width == 0` or
139    /// `height == 0`).
140    pub fn sample(&self, u: f64, v: f64) -> Option<f32> {
141        if self.data.is_empty() || self.width == 0 || self.height == 0 {
142            return None;
143        }
144
145        let u = u.clamp(0.0, 1.0);
146        let v = v.clamp(0.0, 1.0);
147
148        let fx = u * (self.width - 1) as f64;
149        let fy = v * (self.height - 1) as f64;
150
151        let x0 = (fx.floor() as u32).min(self.width - 1);
152        let y0 = (fy.floor() as u32).min(self.height - 1);
153        let x1 = (x0 + 1).min(self.width - 1);
154        let y1 = (y0 + 1).min(self.height - 1);
155
156        let sx = (fx - x0 as f64) as f32;
157        let sy = (fy - y0 as f64) as f32;
158
159        let v00 = self.data[(y0 * self.width + x0) as usize];
160        let v10 = self.data[(y0 * self.width + x1) as usize];
161        let v01 = self.data[(y1 * self.width + x0) as usize];
162        let v11 = self.data[(y1 * self.width + x1) as usize];
163
164        let top = v00 * (1.0 - sx) + v10 * sx;
165        let bot = v01 * (1.0 - sx) + v11 * sx;
166        Some(top * (1.0 - sy) + bot * sy)
167    }
168
169    /// Sample elevation at a geographic coordinate.
170    ///
171    /// Converts the `GeoCoord` to tile-relative `(u, v)` using the
172    /// tile's geographic bounds, then delegates to [`sample`](Self::sample).
173    ///
174    /// Returns `None` if:
175    /// - The tile has zero geographic extent (degenerate).
176    /// - The grid is empty.
177    ///
178    /// Coordinates outside the tile's bounds are clamped to the
179    /// nearest edge (via the `[0, 1]` clamp in `sample`).
180    pub fn sample_geo(&self, coord: &GeoCoord) -> Option<f32> {
181        // NW corner = tile origin, SE corner = next tile diagonally.
182        let nw = tile_to_geo(&self.tile);
183        let se = tile_xy_to_geo(
184            self.tile.zoom,
185            self.tile.x as f64 + 1.0,
186            self.tile.y as f64 + 1.0,
187        );
188
189        let lon_range = se.lon - nw.lon;
190        let lat_range = nw.lat - se.lat;
191
192        if lon_range.abs() < 1e-12 || lat_range.abs() < 1e-12 {
193            return None;
194        }
195
196        let u = (coord.lon - nw.lon) / lon_range;
197        // v is inverted: lat decreases southward but v increases.
198        let v = (nw.lat - coord.lat) / lat_range;
199
200        self.sample(u, v)
201    }
202}
203
204// ---------------------------------------------------------------------------
205// Tests
206// ---------------------------------------------------------------------------
207
208#[cfg(test)]
209mod tests {
210    use super::*;
211
212    // -- Construction ------------------------------------------------------
213
214    #[test]
215    fn flat_grid() {
216        let grid = ElevationGrid::flat(TileId::new(0, 0, 0), 3, 3);
217        assert_eq!(grid.data.len(), 9);
218        assert_eq!(grid.min_elev, 0.0);
219        assert_eq!(grid.max_elev, 0.0);
220        assert_eq!(grid.sample(0.5, 0.5), Some(0.0));
221    }
222
223    #[test]
224    fn from_data_wrong_size() {
225        assert!(ElevationGrid::from_data(TileId::new(0, 0, 0), 2, 2, vec![0.0]).is_none());
226    }
227
228    #[test]
229    fn min_max_elev() {
230        let data = vec![-100.0, 50.0, 200.0, 0.0];
231        let grid = ElevationGrid::from_data(TileId::new(0, 0, 0), 2, 2, data).unwrap();
232        assert!((grid.min_elev - (-100.0)).abs() < 1e-6);
233        assert!((grid.max_elev - 200.0).abs() < 1e-6);
234    }
235
236        #[test]
237        fn from_data_clamps_extreme_elevations() {
238            let data = vec![-32_768.0, -600.0, 10_500.0, 25.0];
239            let grid = ElevationGrid::from_data(TileId::new(0, 0, 0), 2, 2, data).unwrap();
240            assert!((grid.min_elev - (-500.0)).abs() < 1e-6);
241            assert!((grid.max_elev - 10_000.0).abs() < 1e-6);
242            assert_eq!(grid.data, vec![-500.0, -500.0, 10_000.0, 25.0]);
243        }
244
245    #[test]
246    fn elevation_range() {
247        let data = vec![-100.0, 50.0, 200.0, 0.0];
248        let grid = ElevationGrid::from_data(TileId::new(0, 0, 0), 2, 2, data).unwrap();
249        assert!((grid.elevation_range() - 300.0).abs() < 1e-6);
250    }
251
252    #[test]
253    fn elevation_range_flat() {
254        let grid = ElevationGrid::flat(TileId::new(0, 0, 0), 4, 4);
255        assert!((grid.elevation_range()).abs() < 1e-6);
256    }
257
258    // -- Bilinear interpolation (sample) -----------------------------------
259
260    #[test]
261    fn sample_corners() {
262        let data = vec![0.0, 10.0, 20.0, 30.0];
263        let grid = ElevationGrid::from_data(TileId::new(0, 0, 0), 2, 2, data).unwrap();
264        assert!((grid.sample(0.0, 0.0).unwrap() - 0.0).abs() < 1e-6);
265        assert!((grid.sample(1.0, 0.0).unwrap() - 10.0).abs() < 1e-6);
266        assert!((grid.sample(0.0, 1.0).unwrap() - 20.0).abs() < 1e-6);
267        assert!((grid.sample(1.0, 1.0).unwrap() - 30.0).abs() < 1e-6);
268    }
269
270    #[test]
271    fn bilinear_midpoint() {
272        let data = vec![0.0, 10.0, 20.0, 30.0];
273        let grid = ElevationGrid::from_data(TileId::new(0, 0, 0), 2, 2, data).unwrap();
274        let mid = grid.sample(0.5, 0.5).unwrap();
275        // Average of all four corners: (0+10+20+30)/4 = 15.
276        assert!((mid - 15.0).abs() < 1e-6);
277    }
278
279    #[test]
280    fn sample_1x1_grid() {
281        let grid = ElevationGrid::from_data(TileId::new(0, 0, 0), 1, 1, vec![42.0]).unwrap();
282        // Every UV should return the single sample.
283        assert!((grid.sample(0.0, 0.0).unwrap() - 42.0).abs() < 1e-6);
284        assert!((grid.sample(0.5, 0.5).unwrap() - 42.0).abs() < 1e-6);
285        assert!((grid.sample(1.0, 1.0).unwrap() - 42.0).abs() < 1e-6);
286    }
287
288    #[test]
289    fn sample_clamps_out_of_range_uv() {
290        let data = vec![0.0, 10.0, 20.0, 30.0];
291        let grid = ElevationGrid::from_data(TileId::new(0, 0, 0), 2, 2, data).unwrap();
292        // u < 0 clamps to u=0, v < 0 clamps to v=0.
293        assert!((grid.sample(-1.0, -1.0).unwrap() - 0.0).abs() < 1e-6);
294        // u > 1 clamps to u=1, v > 1 clamps to v=1.
295        assert!((grid.sample(2.0, 2.0).unwrap() - 30.0).abs() < 1e-6);
296    }
297
298    #[test]
299    fn sample_empty_grid_returns_none() {
300        let grid = ElevationGrid {
301            width: 0,
302            height: 0,
303            min_elev: 0.0,
304            max_elev: 0.0,
305            data: vec![],
306            tile: TileId::new(0, 0, 0),
307        };
308        assert!(grid.sample(0.5, 0.5).is_none());
309    }
310
311    // -- Geographic sampling (sample_geo) ----------------------------------
312
313    #[test]
314    fn sample_geo_tile_center() {
315        // Create a 2x2 grid on tile 0/0/0 (covers the whole world).
316        let data = vec![100.0, 200.0, 300.0, 400.0];
317        let grid = ElevationGrid::from_data(TileId::new(0, 0, 0), 2, 2, data).unwrap();
318
319        // The center of tile 0/0/0 is approximately (0, 0) lat/lon.
320        let center = GeoCoord::from_lat_lon(0.0, 0.0);
321        let elev = grid.sample_geo(&center).unwrap();
322        // At the midpoint of a 2x2 grid, expect the average: 250.
323        assert!((elev - 250.0).abs() < 1.0);
324    }
325
326    #[test]
327    fn sample_geo_nw_corner() {
328        // Tile 1/0/0 covers the NW quadrant.
329        let data = vec![10.0, 20.0, 30.0, 40.0];
330        let grid = ElevationGrid::from_data(TileId::new(1, 0, 0), 2, 2, data).unwrap();
331
332        // NW corner of tile 1/0/0 should return the (0,0) sample.
333        let nw = tile_to_geo(&TileId::new(1, 0, 0));
334        let elev = grid.sample_geo(&nw).unwrap();
335        assert!((elev - 10.0).abs() < 1e-3);
336    }
337
338    // -- PartialEq ---------------------------------------------------------
339
340    #[test]
341    fn partial_eq() {
342        let a = ElevationGrid::flat(TileId::new(5, 10, 10), 4, 4);
343        let b = ElevationGrid::flat(TileId::new(5, 10, 10), 4, 4);
344        assert_eq!(a, b);
345    }
346}