rustial_engine/terrain/

mesh.rs

//! CPU-side terrain mesh generation.

use crate::camera_projection::CameraProjection;
use crate::tile_manager::TileTextureRegion;
use rustial_math::{ElevationGrid, TileId};
use std::sync::Arc;

/// Compute zoom-dependent skirt height in meters.
///
/// This follows the Mapbox GL JS heuristic: the skirt should be large
/// enough to hide seams between tiles at different LODs but small
/// enough to avoid visible vertical faces near the camera.  The value
/// decreases exponentially with zoom:
///
/// - zoom 0:  ~20 000 m
/// - zoom 11: ~520 m
/// - zoom 22: ~6 m
///
/// When `exaggeration` is greater than 1 the skirt is scaled
/// proportionally so that amplified elevation differences are still
/// fully hidden.
pub fn skirt_height(zoom: u8, exaggeration: f64) -> f64 {
    6.0 * 1.5_f64.powi(22 - zoom as i32) * exaggeration.max(1.0)
}

fn tile_vertex_geo(tile: &TileId, u: f64, v: f64) -> rustial_math::GeoCoord {
    let nw = rustial_math::tile_to_geo(tile);
    let se = rustial_math::tile_xy_to_geo(tile.zoom, tile.x as f64 + 1.0, tile.y as f64 + 1.0);
    let lat = nw.lat + (se.lat - nw.lat) * v;
    let lon = nw.lon + (se.lon - nw.lon) * u;
    rustial_math::GeoCoord::from_lat_lon(lat, lon)
}

fn project_tile_vertex(
    projection: CameraProjection,
    tile: &TileId,
    u: f64,
    v: f64,
    altitude: f64,
) -> [f64; 3] {
    let mut geo = tile_vertex_geo(tile, u, v);
    geo.alt = altitude;
    projection.project_position(&geo)
}

/// GPU-sampleable elevation texture payload for a terrain tile.
#[derive(Debug, Clone)]
pub struct TerrainElevationTexture {
    /// Width of the elevation texture in texels.
    pub width: u32,
    /// Height of the elevation texture in texels.
    pub height: u32,
    /// Minimum elevation in meters.
    pub min_elev: f32,
    /// Maximum elevation in meters.
    pub max_elev: f32,
    /// Elevation samples in meters, row-major, matching [`ElevationGrid`].
    ///
    /// Wrapped in [`Arc`] so that cloning terrain descriptors (e.g. for
    /// the frame-key cache in `TerrainManager`) is a cheap reference-count
    /// bump instead of a ~260 KB memcpy per tile.
    pub data: Arc<Vec<f32>>,
}

/// Remap a logical DEM region from interior tile UV space into the
/// full texture UV space of a border-expanded elevation grid.
///
/// Terrain DEM tiles are cached as `(W+2) x (H+2)` textures with a
/// 1-sample border used for seam patching.  `TileTextureRegion`
/// values, however, are defined over the original interior tile.  This
/// helper shifts and scales the region so CPU and GPU terrain sampling
/// address the correct sub-rectangle inside the expanded texture.
pub fn elevation_region_in_texture_space(
    region: TileTextureRegion,
    width: u32,
    height: u32,
) -> TileTextureRegion {
    if width < 4 || height < 4 {
        return region;
    }

    let map_axis = |min: f32, max: f32, extent: u32| {
        let denom = (extent - 1) as f32;
        let interior_span = (extent - 3) as f32;
        (
            (1.0 + min * interior_span) / denom,
            (1.0 + max * interior_span) / denom,
        )
    };

    let (u_min, u_max) = map_axis(region.u_min, region.u_max, width);
    let (v_min, v_max) = map_axis(region.v_min, region.v_max, height);
    TileTextureRegion {
        u_min,
        v_min,
        u_max,
        v_max,
    }
}

/// CPU-side terrain mesh data ready for GPU upload.
#[derive(Debug, Clone)]
pub struct TerrainMeshData {
    /// Tile this mesh represents.
    pub tile: TileId,
    /// Elevation source tile that backs this mesh.
    ///
    /// When terrain display zoom exceeds the DEM source max zoom, multiple
    /// child meshes reuse a parent DEM tile. This field identifies that parent.
    pub elevation_source_tile: TileId,
    /// Normalized sub-region of the elevation texture to sample for this mesh.
    ///
    /// `FULL` when the elevation texture matches the display tile. For
    /// overzoomed terrain meshes this is the child sub-rect inside the parent
    /// DEM tile.
    pub elevation_region: TileTextureRegion,
    /// Vertex positions in world-space (f64 for precision, cast to f32 at upload).
    ///
    /// For the reusable-grid path these positions represent the planar base
    /// surface before elevation displacement. Elevation is supplied separately
    /// via [`elevation_texture`] for GPU sampling.
    pub positions: Vec<[f64; 3]>,
    /// Texture coordinates.
    pub uvs: Vec<[f32; 2]>,
    /// Surface normals.
    ///
    /// For reusable-grid terrain these are base-grid normals. Final terrain
    /// shading may derive finer normals in the renderer from sampled heights.
    pub normals: Vec<[f32; 3]>,
    /// Triangle indices.
    pub indices: Vec<u32>,
    /// Generation counter from the terrain manager's elevation cache.
    ///
    /// Increments whenever the elevation source delivers new data.
    /// Renderers compare this against the value stored on their entity
    /// to detect when a mesh must be rebuilt (e.g. when a flat
    /// placeholder is replaced by real elevation data).
    pub generation: u64,
    /// Grid resolution used to build the reusable terrain mesh.
    pub grid_resolution: u16,
    /// Vertical exaggeration already selected for this tile.
    pub vertical_exaggeration: f32,
    /// Optional GPU-sampleable elevation payload for this tile.
    pub elevation_texture: Option<TerrainElevationTexture>,
}

/// Build a lightweight terrain descriptor that carries reusable-grid metadata
/// and GPU-sampleable elevation payload without eagerly generating displaced
/// CPU geometry.
pub fn build_terrain_descriptor(
    tile: &TileId,
    elevation: &ElevationGrid,
    resolution: u16,
    exaggeration: f64,
    generation: u64,
) -> TerrainMeshData {
    build_terrain_descriptor_with_source(
        tile,
        *tile,
        TileTextureRegion::FULL,
        elevation,
        resolution,
        exaggeration,
        generation,
    )
}

/// Build a lightweight terrain descriptor with explicit elevation-source mapping.
pub fn build_terrain_descriptor_with_source(
    tile: &TileId,
    elevation_source_tile: TileId,
    elevation_region: TileTextureRegion,
    elevation: &ElevationGrid,
    resolution: u16,
    exaggeration: f64,
    generation: u64,
) -> TerrainMeshData {
    // Compute min/max from the interior of the elevation grid, excluding
    // the 1-pixel border that may contain extreme values from neighboring
    // ocean tiles patched via backfill.  For overzoomed tiles this also
    // restricts to the subregion actually used by this child tile.
    let (min_elev, max_elev) = subregion_min_max(elevation, &elevation_region);

    TerrainMeshData {
        tile: *tile,
        elevation_source_tile,
        elevation_region,
        positions: Vec::new(),
        uvs: Vec::new(),
        normals: Vec::new(),
        indices: Vec::new(),
        generation,
        grid_resolution: resolution,
        vertical_exaggeration: exaggeration as f32,
        elevation_texture: Some(TerrainElevationTexture {
            width: elevation.width,
            height: elevation.height,
            min_elev,
            max_elev,
            data: Arc::new(elevation.data.clone()),
        }),
    }
}

/// Compute min/max elevation from a subregion of the elevation grid.
///
/// The grid may be border-expanded (W+2 x H+2) where the interior starts
/// at pixel (1,1).  Texture region coordinates (0..1) refer to the
/// original interior, so we must offset by 1 if a border is present.
fn subregion_min_max(elevation: &ElevationGrid, region: &TileTextureRegion) -> (f32, f32) {
    let w = elevation.width as usize;
    let h = elevation.height as usize;
    if w <= 2 || h <= 2 || elevation.data.is_empty() {
        return (elevation.min_elev, elevation.max_elev);
    }

    // Detect border-expanded grids: original interior is (w-2) x (h-2),
    // interior pixels start at (1, 1).
    let (interior_w, interior_h, offset) = if w >= 4 && h >= 4 {
        // Assume 1-pixel border expansion is always present when
        // the grid is used with overzoom texture regions.
        (w - 2, h - 2, 1usize)
    } else {
        (w, h, 0)
    };

    // Map normalized texture coordinates to interior pixel ranges,
    // then shift into the full (possibly expanded) grid.
    let x0 = (region.u_min as f64 * interior_w as f64).floor() as usize + offset;
    let x1 = ((region.u_max as f64 * interior_w as f64).ceil() as usize + offset).min(w);
    let y0 = (region.v_min as f64 * interior_h as f64).floor() as usize + offset;
    let y1 = ((region.v_max as f64 * interior_h as f64).ceil() as usize + offset).min(h);

    let mut lo = f32::MAX;
    let mut hi = f32::MIN;
    for y in y0..y1 {
        let row_start = y * w;
        for x in x0..x1 {
            let v = elevation.data[row_start + x];
            if v < lo {
                lo = v;
            }
            if v > hi {
                hi = v;
            }
        }
    }

    if lo > hi {
        (elevation.min_elev, elevation.max_elev)
    } else {
        (lo, hi)
    }
}

/// Materialize displaced CPU geometry from a lightweight terrain descriptor.
///
/// If geometry is already present it is returned unchanged.
pub fn materialize_terrain_mesh(
    mesh: &TerrainMeshData,
    projection: CameraProjection,
    skirt_depth: f64,
) -> TerrainMeshData {
    if !mesh.positions.is_empty() {
        return mesh.clone();
    }

    let Some(elevation_texture) = mesh.elevation_texture.as_ref() else {
        return mesh.clone();
    };
    let Some(elevation) = ElevationGrid::from_data(
        mesh.tile,
        elevation_texture.width,
        elevation_texture.height,
        elevation_texture.data.to_vec(),
    ) else {
        return mesh.clone();
    };

    build_terrain_mesh_with_source(
        &mesh.tile,
        mesh.elevation_source_tile,
        mesh.elevation_region,
        &elevation,
        projection,
        mesh.grid_resolution,
        mesh.vertical_exaggeration as f64,
        skirt_depth,
        mesh.generation,
    )
}

/// Build a terrain mesh for a tile, displacing vertex Z by elevation values.
///
/// `resolution` is the number of vertices per edge (e.g. 64 means a 64x64 grid).
/// `exaggeration` scales elevation values vertically.
/// `skirt_depth` adds a vertical skirt around the tile edges.
/// `generation` is stamped into the output so renderers can detect changes.
pub fn build_terrain_mesh(
    tile: &TileId,
    elevation: &ElevationGrid,
    projection: CameraProjection,
    resolution: u16,
    exaggeration: f64,
    skirt_depth: f64,
    generation: u64,
) -> TerrainMeshData {
    build_terrain_mesh_with_source(
        tile,
        *tile,
        TileTextureRegion::FULL,
        elevation,
        projection,
        resolution,
        exaggeration,
        skirt_depth,
        generation,
    )
}

/// Build a terrain mesh with explicit elevation-source mapping metadata.
#[allow(clippy::too_many_arguments)]
pub fn build_terrain_mesh_with_source(
    tile: &TileId,
    elevation_source_tile: TileId,
    elevation_region: TileTextureRegion,
    elevation: &ElevationGrid,
    projection: CameraProjection,
    resolution: u16,
    exaggeration: f64,
    skirt_depth: f64,
    generation: u64,
) -> TerrainMeshData {
    let res = resolution as usize;
    let vertex_count = res * res;
    let mut positions = Vec::with_capacity(vertex_count);
    let mut uvs = Vec::with_capacity(vertex_count);
    let sample_region = if *tile != elevation_source_tile {
        elevation_region_in_texture_space(elevation_region, elevation.width, elevation.height)
    } else {
        elevation_region
    };

    // Generate grid vertices.
    for row in 0..res {
        for col in 0..res {
            let u = col as f64 / (res - 1).max(1) as f64;
            let v = row as f64 / (res - 1).max(1) as f64;
            let sample_u =
                sample_region.u_min as f64 + (sample_region.u_max - sample_region.u_min) as f64 * u;
            let sample_v =
                sample_region.v_min as f64 + (sample_region.v_max - sample_region.v_min) as f64 * v;

            let raw_elev = elevation
                .sample(sample_u, sample_v)
                .unwrap_or(0.0)
                .clamp(-500.0, 10_000.0);
            let elev = raw_elev as f64 * exaggeration;

            positions.push(project_tile_vertex(projection, tile, u, v, elev));
            uvs.push([u as f32, v as f32]);
        }
    }

    // Generate indices.
    let quad_count = (res - 1) * (res - 1);
    let mut indices = Vec::with_capacity(quad_count * 6);
    for row in 0..(res - 1) {
        for col in 0..(res - 1) {
            let tl = (row * res + col) as u32;
            let tr = tl + 1;
            let bl = ((row + 1) * res + col) as u32;
            let br = bl + 1;
            indices.push(tl);
            indices.push(bl);
            indices.push(tr);
            indices.push(tr);
            indices.push(bl);
            indices.push(br);
        }
    }

    // Compute normals via central differences.
    let mut normals = vec![[0.0f32; 3]; vertex_count];
    for row in 0..res {
        for col in 0..res {
            let idx = row * res + col;
            let left = positions[if col > 0 { idx - 1 } else { idx }];
            let right = positions[if col < res - 1 { idx + 1 } else { idx }];
            let down = positions[if row < res - 1 { idx + res } else { idx }];
            let up = positions[if row > 0 { idx - res } else { idx }];

            let tangent_x = [
                (right[0] - left[0]) as f32,
                (right[1] - left[1]) as f32,
                (right[2] - left[2]) as f32,
            ];
            let tangent_y = [
                (up[0] - down[0]) as f32,
                (up[1] - down[1]) as f32,
                (up[2] - down[2]) as f32,
            ];

            let nx = tangent_y[1] * tangent_x[2] - tangent_y[2] * tangent_x[1];
            let ny = tangent_y[2] * tangent_x[0] - tangent_y[0] * tangent_x[2];
            let nz = tangent_y[0] * tangent_x[1] - tangent_y[1] * tangent_x[0];
            let len = (nx * nx + ny * ny + nz * nz).sqrt();
            normals[idx] = if len > 1e-6 {
                let mut normal = [nx / len, ny / len, nz / len];
                if normal[2] < 0.0 {
                    normal = [-normal[0], -normal[1], -normal[2]];
                }
                normal
            } else {
                [0.0, 0.0, 1.0]
            };
        }
    }

    if skirt_depth > 0.0 {
        add_skirt(
            &mut positions,
            &mut uvs,
            &mut normals,
            &mut indices,
            res,
            skirt_depth,
        );
    }

    TerrainMeshData {
        tile: *tile,
        elevation_source_tile,
        elevation_region,
        positions,
        uvs,
        normals,
        indices,
        generation,
        grid_resolution: resolution,
        vertical_exaggeration: exaggeration as f32,
        elevation_texture: Some(TerrainElevationTexture {
            width: elevation.width,
            height: elevation.height,
            min_elev: elevation.min_elev,
            max_elev: elevation.max_elev,
            data: Arc::new(elevation.data.clone()),
        }),
    }
}

fn add_skirt(
    positions: &mut Vec<[f64; 3]>,
    uvs: &mut Vec<[f32; 2]>,
    normals: &mut Vec<[f32; 3]>,
    indices: &mut Vec<u32>,
    res: usize,
    skirt_depth: f64,
) {
    // Compute a single base Z that ALL skirt vertices drop to.
    let min_z = positions[..res * res]
        .iter()
        .map(|p| p[2])
        .fold(f64::INFINITY, f64::min);
    let skirt_z = min_z - skirt_depth;

    // For each edge, duplicate the edge vertices at skirt_z
    // and add triangles connecting them.
    //
    // Outward-facing normals per edge so the skirt receives proper
    // hillshade lighting instead of appearing dark.
    let edges: [(&[f32; 3], Vec<usize>); 4] = [
        // Top edge (row 0) -- faces north (+Y).
        (&[0.0, 1.0, 0.0], (0..res).collect()),
        // Bottom edge (row res-1) -- faces south (-Y).
        (&[0.0, -1.0, 0.0], ((res - 1) * res..res * res).collect()),
        // Left edge (col 0) -- faces west (-X).
        (&[-1.0, 0.0, 0.0], (0..res).map(|r| r * res).collect()),
        // Right edge (col res-1) -- faces east (+X).
        (
            &[1.0, 0.0, 0.0],
            (0..res).map(|r| r * res + res - 1).collect(),
        ),
    ];

    for (normal, edge) in &edges {
        for i in 0..edge.len() - 1 {
            let a = edge[i] as u32;
            let b = edge[i + 1] as u32;

            let base_a = positions.len() as u32;
            let base_b = base_a + 1;

            // Skirt vertex for a: same XY, dropped to skirt_z.
            let pa = positions[edge[i]];
            positions.push([pa[0], pa[1], skirt_z]);
            uvs.push(uvs[edge[i]]);
            normals.push(**normal);

            // Skirt vertex for b: same XY, dropped to skirt_z.
            let pb = positions[edge[i + 1]];
            positions.push([pb[0], pb[1], skirt_z]);
            uvs.push(uvs[edge[i + 1]]);
            normals.push(**normal);

            // Two triangles for the skirt quad.
            indices.push(a);
            indices.push(base_a);
            indices.push(b);
            indices.push(b);
            indices.push(base_a);
            indices.push(base_b);
        }
    }
}

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

    #[test]
    fn flat_mesh_z_zero() {
        let tile = TileId::new(10, 100, 100);
        let elev = ElevationGrid::flat(tile, 4, 4);
        let mesh = build_terrain_mesh(&tile, &elev, CameraProjection::WebMercator, 4, 1.0, 0.0, 0);

        // 4x4 grid = 16 vertices, 3*3*2 = 18 triangles = 54 indices.
        assert_eq!(mesh.positions.len(), 16);
        assert_eq!(mesh.indices.len(), 54);

        // All Z should be 0 for flat grid.
        for pos in &mesh.positions {
            assert!((pos[2] - 0.0).abs() < 1e-6);
        }
    }

    #[test]
    fn sloped_mesh_z_nonzero() {
        let tile = TileId::new(10, 100, 100);
        let data = vec![0.0, 100.0, 200.0, 300.0];
        let elev = ElevationGrid::from_data(tile, 2, 2, data).unwrap();
        let mesh = build_terrain_mesh(&tile, &elev, CameraProjection::WebMercator, 2, 1.0, 0.0, 0);

        // Check that the Z values are displaced.
        assert!((mesh.positions[0][2] - 0.0).abs() < 1e-6);
        assert!((mesh.positions[1][2] - 100.0).abs() < 1e-6);
    }

    #[test]
    fn exaggeration() {
        let tile = TileId::new(10, 100, 100);
        let data = vec![100.0; 4];
        let elev = ElevationGrid::from_data(tile, 2, 2, data).unwrap();
        let mesh = build_terrain_mesh(&tile, &elev, CameraProjection::WebMercator, 2, 3.0, 0.0, 0);

        for pos in &mesh.positions {
            assert!((pos[2] - 300.0).abs() < 1e-6);
        }
    }

    #[test]
    fn skirt_adds_vertices() {
        let tile = TileId::new(10, 100, 100);
        let elev = ElevationGrid::flat(tile, 4, 4);
        let mesh_no_skirt =
            build_terrain_mesh(&tile, &elev, CameraProjection::WebMercator, 4, 1.0, 0.0, 0);
        let mesh_with_skirt = build_terrain_mesh(
            &tile,
            &elev,
            CameraProjection::WebMercator,
            4,
            1.0,
            100.0,
            0,
        );

        assert!(mesh_with_skirt.positions.len() > mesh_no_skirt.positions.len());
        assert!(mesh_with_skirt.indices.len() > mesh_no_skirt.indices.len());
    }

    #[test]
    fn normals_flat_point_up() {
        let tile = TileId::new(10, 100, 100);
        let elev = ElevationGrid::flat(tile, 4, 4);
        let mesh = build_terrain_mesh(&tile, &elev, CameraProjection::WebMercator, 4, 1.0, 0.0, 0);

        // Flat surface: all normals should point up (Z ~ 1).
        for n in &mesh.normals {
            assert!(n[2] > 0.99);
        }
    }

    #[test]
    fn adjacent_tiles_share_edge_positions() {
        // Two horizontally adjacent tiles at zoom 10 should have
        // identical vertex positions along their shared edge.
        let left = TileId::new(10, 100, 100);
        let right = TileId::new(10, 101, 100);

        let elev_l = ElevationGrid::flat(left, 4, 4);
        let elev_r = ElevationGrid::flat(right, 4, 4);

        let mesh_l = build_terrain_mesh(
            &left,
            &elev_l,
            CameraProjection::WebMercator,
            4,
            1.0,
            0.0,
            0,
        );
        let mesh_r = build_terrain_mesh(
            &right,
            &elev_r,
            CameraProjection::WebMercator,
            4,
            1.0,
            0.0,
            0,
        );

        // Right edge of left tile: col = 3 (last col), all rows.
        // Left edge of right tile: col = 0 (first col), all rows.
        let res = 4;
        for row in 0..res {
            let l_idx = row * res + (res - 1); // right edge of left tile
            let r_idx = row * res; // left edge of right tile

            let l_pos = mesh_l.positions[l_idx];
            let r_pos = mesh_r.positions[r_idx];

            assert!(
                (l_pos[0] - r_pos[0]).abs() < 1e-6
                    && (l_pos[1] - r_pos[1]).abs() < 1e-6
                    && (l_pos[2] - r_pos[2]).abs() < 1e-6,
                "row {row}: left right-edge {l_pos:?} != right left-edge {r_pos:?}",
            );
        }
    }

    #[test]
    fn adjacent_tiles_share_vertical_edge() {
        // Two vertically adjacent tiles should share positions along
        // bottom of upper / top of lower.
        let upper = TileId::new(10, 100, 100);
        let lower = TileId::new(10, 100, 101);

        let elev_u = ElevationGrid::flat(upper, 4, 4);
        let elev_d = ElevationGrid::flat(lower, 4, 4);

        let mesh_u = build_terrain_mesh(
            &upper,
            &elev_u,
            CameraProjection::WebMercator,
            4,
            1.0,
            0.0,
            0,
        );
        let mesh_d = build_terrain_mesh(
            &lower,
            &elev_d,
            CameraProjection::WebMercator,
            4,
            1.0,
            0.0,
            0,
        );

        let res = 4;
        for col in 0..res {
            let u_idx = (res - 1) * res + col; // bottom edge of upper tile
            let d_idx = col; // top edge of lower tile

            let u_pos = mesh_u.positions[u_idx];
            let d_pos = mesh_d.positions[d_idx];

            assert!(
                (u_pos[0] - d_pos[0]).abs() < 1e-6
                    && (u_pos[1] - d_pos[1]).abs() < 1e-6
                    && (u_pos[2] - d_pos[2]).abs() < 1e-6,
                "col {col}: upper bottom-edge {u_pos:?} != lower top-edge {d_pos:?}",
            );
        }
    }

    #[test]
    fn adjacent_tiles_share_edge_with_elevation() {
        // Verify that two horizontally adjacent tiles with non-flat
        // elevation still share exact positions along their common edge.
        // Uses a simple slope: elev = u * 100 (increasing east).
        let left = TileId::new(10, 100, 100);
        let right = TileId::new(10, 101, 100);

        let res: u16 = 4;
        let w = res as u32;
        let h = res as u32;

        // Left tile: elev rises from 0 (west) to 100 (east).
        let data_l: Vec<f32> = (0..h)
            .flat_map(|_row| (0..w).map(|col| col as f32 / (w - 1) as f32 * 100.0))
            .collect();
        let elev_l = ElevationGrid::from_data(left, w, h, data_l).unwrap();

        // Right tile starts at 100 (matching left's right edge).
        let data_r: Vec<f32> = (0..h)
            .flat_map(|_row| (0..w).map(|col| 100.0 + col as f32 / (w - 1) as f32 * 100.0))
            .collect();
        let elev_r = ElevationGrid::from_data(right, w, h, data_r).unwrap();

        let mesh_l = build_terrain_mesh(
            &left,
            &elev_l,
            CameraProjection::WebMercator,
            res,
            1.0,
            0.0,
            0,
        );
        let mesh_r = build_terrain_mesh(
            &right,
            &elev_r,
            CameraProjection::WebMercator,
            res,
            1.0,
            0.0,
            0,
        );

        let r = res as usize;
        for row in 0..r {
            let l_idx = row * r + (r - 1);
            let r_idx = row * r;

            let l_z = mesh_l.positions[l_idx][2];
            let r_z = mesh_r.positions[r_idx][2];

            assert!(
                (l_z - r_z).abs() < 1e-3,
                "row {row}: left right-edge Z={l_z:.4} != right left-edge Z={r_z:.4}",
            );
        }
    }

    #[test]
    fn skirt_drops_to_absolute_base() {
        // Verify that all skirt vertices share the same Z value
        // (min_surface_z - skirt_depth) regardless of the edge vertex
        // elevation.  This ensures adjacent tile skirts overlap.
        let tile = TileId::new(10, 100, 100);
        let data = vec![0.0, 100.0, 200.0, 300.0]; // sloped 2x2 grid
        let elev = ElevationGrid::from_data(tile, 2, 2, data).unwrap();
        let skirt_depth = 50.0;
        let mesh = build_terrain_mesh(
            &tile,
            &elev,
            CameraProjection::WebMercator,
            2,
            1.0,
            skirt_depth,
            0,
        );

        // Surface vertices: indices 0..4 with Z = 0, 100, 200, 300.
        // Min surface Z = 0. Skirt base = 0 - 50 = -50.
        let surface_count = 4; // 2x2 grid
        let skirt_vertices = &mesh.positions[surface_count..];
        assert!(!skirt_vertices.is_empty(), "should have skirt vertices");

        let expected_z = -50.0;
        for (i, sv) in skirt_vertices.iter().enumerate() {
            assert!(
                (sv[2] - expected_z).abs() < 1e-6,
                "skirt vertex {i}: Z={:.4}, expected {expected_z:.4}",
                sv[2],
            );
        }
    }

    #[test]
    fn adjacent_tile_skirts_overlap() {
        let right = TileId::new(10, 101, 100);

        let res: u16 = 4;
        let w = res as u32;
        let h = res as u32;

        // Right tile: flat at 1000m.
        let data_r = vec![1000.0f32; (w * h) as usize];
        let elev_r = ElevationGrid::from_data(right, w, h, data_r).unwrap();

        // With skirt_depth = 1200 (> 1000m height above Z=0), the right
        // tile's skirt base = 1000 - 1200 = -200, extending below any
        // neighbouring tile at Z=0.
        let mesh_r_deep = build_terrain_mesh(
            &right,
            &elev_r,
            CameraProjection::WebMercator,
            res,
            1.0,
            1200.0,
            0,
        );

        let surface_count = (res as usize) * (res as usize);
        let skirt_z_r: Vec<f64> = mesh_r_deep.positions[surface_count..]
            .iter()
            .map(|p| p[2])
            .collect();
        assert!(
            skirt_z_r.iter().all(|&z| z < 0.0),
            "right tile skirt should extend below left tile surface (Z=0)"
        );
    }

    #[test]
    fn equirectangular_projection_changes_xy_positions() {
        let tile = TileId::new(3, 4, 2);
        let elev = ElevationGrid::flat(tile, 4, 4);
        let merc = build_terrain_mesh(&tile, &elev, CameraProjection::WebMercator, 4, 1.0, 0.0, 0);
        let eq = build_terrain_mesh(
            &tile,
            &elev,
            CameraProjection::Equirectangular,
            4,
            1.0,
            0.0,
            0,
        );

        assert_eq!(merc.positions.len(), eq.positions.len());
        let different_xy = merc
            .positions
            .iter()
            .zip(eq.positions.iter())
            .any(|(a, b)| (a[0] - b[0]).abs() > 1.0 || (a[1] - b[1]).abs() > 1.0);
        assert!(different_xy);
    }

    #[test]
    fn elevation_region_maps_to_bordered_texture_space() {
        let full = elevation_region_in_texture_space(TileTextureRegion::FULL, 6, 6);
        assert!((full.u_min - 0.2).abs() < 1e-6);
        assert!((full.v_min - 0.2).abs() < 1e-6);
        assert!((full.u_max - 0.8).abs() < 1e-6);
        assert!((full.v_max - 0.8).abs() < 1e-6);

        let quarter = elevation_region_in_texture_space(
            TileTextureRegion {
                u_min: 0.0,
                v_min: 0.0,
                u_max: 0.5,
                v_max: 0.5,
            },
            6,
            6,
        );
        assert!((quarter.u_min - 0.2).abs() < 1e-6);
        assert!((quarter.v_min - 0.2).abs() < 1e-6);
        assert!((quarter.u_max - 0.5).abs() < 1e-6);
        assert!((quarter.v_max - 0.5).abs() < 1e-6);
    }

    #[test]
    fn overzoom_child_mesh_samples_only_child_region() {
        let child = TileId::new(1, 0, 0);
        let source = TileId::new(0, 0, 0);
        let data = vec![
            0.0, 0.0, 0.0, 100.0, 100.0, 100.0, 0.0, 0.0, 0.0, 100.0, 100.0, 100.0, 0.0, 0.0, 0.0,
            100.0, 100.0, 100.0, 200.0, 200.0, 200.0, 300.0, 300.0, 300.0, 200.0, 200.0, 200.0,
            300.0, 300.0, 300.0, 200.0, 200.0, 200.0, 300.0, 300.0, 300.0,
        ];
        let elev = ElevationGrid::from_data(source, 6, 6, data).unwrap();
        let mesh = build_terrain_mesh_with_source(
            &child,
            source,
            TileTextureRegion::from_tiles(&child, &source),
            &elev,
            CameraProjection::WebMercator,
            2,
            1.0,
            0.0,
            0,
        );

        let z_values: Vec<f64> = mesh.positions.iter().map(|p| p[2]).collect();
        let expected = [0.0, 50.0, 100.0, 150.0];
        for (actual, expected) in z_values.iter().zip(expected.iter()) {
            assert!(
                (actual - expected).abs() < 1e-3,
                "child mesh should sample the top-left parent subregion, got {z_values:?}"
            );
        }
    }
}