rustial_engine/math/

tile.rs

//! Tile addressing for slippy-map tile pyramids.
//!
//! # Tile grid convention
//!
//! This module uses the [OpenStreetMap / Slippy Map][osm] tile numbering
//! scheme, which is shared by virtually all web map services (Mapbox, Google
//! Maps, Leaflet, etc.):
//!
//! - **Origin** is the top-left (north-west) corner of the world.
//! - **X** increases eastward (column index).
//! - **Y** increases southward (row index).
//! - At zoom level `z`, the world is divided into `2^z * 2^z` tiles.
//!
//! Internally the tile grid is backed by the Web Mercator projection
//! (EPSG:3857), so latitudes beyond approximately +/-85.06 degrees are not
//! representable in the tile grid.
//!
//! # Coordinate precision
//!
//! All conversions between geographic coordinates and tile coordinates are
//! performed in f64.  The fractional [`TileCoord`] type preserves sub-tile
//! precision, which is important for the renderer's camera-relative mesh
//! placement (avoids f32 jitter when the camera is far from tile corners).
//!
//! [osm]: https://wiki.openstreetmap.org/wiki/Slippy_map_tilenames

use crate::bounds::WorldBounds;
use crate::coord::{GeoCoord, WorldCoord};
use crate::mercator::WebMercator;
use glam::{DMat4, DVec3, DVec4};
use std::f64::consts::PI;
use std::fmt;

/// Maximum zoom level for slippy-map tiles.
///
/// At zoom 22 the tile grid is 4 194 304 x 4 194 304, giving roughly 2 cm
/// ground resolution at the equator.  Higher levels are impractical for
/// raster tiles and would risk `u32` overflow in `axis_tiles()` at zoom 32.
pub const MAX_ZOOM: u8 = 22;

// ---------------------------------------------------------------------------
// Flat-view tile selection
// ---------------------------------------------------------------------------

/// Policy knobs for footprint-aware flat raster tile selection.
///
/// These values control when the selector switches from conservative
/// [`visible_tiles`] AABB coverage to the more precise sampled ground
/// footprint, and how aggressively that footprint is sampled/clamped.
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct FlatTileSelectionConfig {
    /// Minimum pitch (radians) before footprint filtering is applied.
    pub footprint_pitch_threshold_rad: f64,
    /// Minimum raw tile count before footprint filtering is applied.
    pub footprint_min_tiles: usize,
    /// Number of sample segments per screen edge when building the
    /// ground-footprint polygon.
    pub footprint_edge_steps: usize,
    /// Upper bound in meters for ground-intersection ray length.
    pub max_ground_distance: f64,
    /// Upper bound in meters for sky-pointing ray XY projection.
    pub max_sky_distance: f64,
}

impl Default for FlatTileSelectionConfig {
    fn default() -> Self {
        Self {
            footprint_pitch_threshold_rad: 0.3,
            footprint_min_tiles: 64,
            footprint_edge_steps: 24,
            max_ground_distance: 500_000.0,
            max_sky_distance: 400_000.0,
        }
    }
}

/// Camera parameters needed to select flat raster tiles for a pitched
/// perspective map view.
///
/// The selection algorithm starts from the conservative axis-aligned
/// [`visible_tiles`] result computed from the world-space viewport bounds,
/// then filters that set against the sampled ground footprint of the actual
/// screen frustum. This keeps coverage correct while avoiding severe tile
/// over-selection at steep pitch angles where the visible ground forms a thin,
/// rotated trapezoid rather than a broad rectangle.
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct FlatTileView {
    /// Camera target in Web Mercator world coordinates (meters).
    pub target_world: WorldCoord,
    /// Camera orbit distance in meters.
    pub distance: f64,
    /// Camera pitch in radians.
    pub pitch: f64,
    /// Camera yaw in radians.
    pub yaw: f64,
    /// Vertical field of view in radians.
    pub fov_y: f64,
    /// Viewport width in logical pixels.
    pub viewport_width: u32,
    /// Viewport height in logical pixels.
    pub viewport_height: u32,
}

impl FlatTileView {
    /// Construct flat-view selection parameters.
    #[inline]
    pub fn new(
        target_world: WorldCoord,
        distance: f64,
        pitch: f64,
        yaw: f64,
        fov_y: f64,
        viewport_width: u32,
        viewport_height: u32,
    ) -> Self {
        Self {
            target_world,
            distance,
            pitch,
            yaw,
            fov_y,
            viewport_width,
            viewport_height,
        }
    }
}

/// Return visible flat raster tiles for a pitched perspective view.
///
/// This is the production tile-selection path for **flat tile quads**.
/// It preserves complete coverage while reducing the extreme overfetch that
/// occurs when a steeply pitched view uses only an axis-aligned Mercator AABB.
///
/// # Algorithm
///
/// 1. Compute the conservative AABB tile set with [`visible_tiles`].
/// 2. For sufficiently pitched views, sample the full screen border
///    against the ground plane to build a world-space footprint polygon.
/// 3. Keep only tiles whose bounds intersect that sampled footprint.
///
/// Low-pitch and small tile sets bypass the extra filtering step.
pub fn visible_tiles_flat_view(
    bounds: &WorldBounds,
    zoom: u8,
    view: &FlatTileView,
) -> Vec<TileId> {
    visible_tiles_flat_view_with_config(bounds, zoom, view, &FlatTileSelectionConfig::default())
}

/// Return visible flat raster tiles for a pitched perspective view using
/// an explicit flat-tile selection policy.
///
/// This is the configurable form of [`visible_tiles_flat_view`].
pub fn visible_tiles_flat_view_with_config(
    bounds: &WorldBounds,
    zoom: u8,
    view: &FlatTileView,
    config: &FlatTileSelectionConfig,
) -> Vec<TileId> {
    let mut tiles = visible_tiles(bounds, zoom);

    if view.pitch <= config.footprint_pitch_threshold_rad
        || tiles.len() <= config.footprint_min_tiles
    {
        return tiles;
    }

    let footprint = sampled_ground_footprint(view, config);
    if footprint.len() < 3 {
        return tiles;
    }

    tiles.retain(|tile| tile_intersects_ground_footprint(*tile, &footprint));
    tiles
}

fn refine_nearby_flat_tiles(
    tiles: Vec<TileId>,
    zoom: u8,
    view: &FlatTileView,
    _footprint: &[(f64, f64)],
) -> Vec<TileId> {
    if zoom >= MAX_ZOOM || tiles.is_empty() {
        return tiles;
    }

    let cam_x = view.target_world.position.x;
    let cam_y = view.target_world.position.y;
    let refine_radius = (view.distance * 2.5).max(60_000.0);
    let refine_radius_sq = refine_radius * refine_radius;

    let mut refined = Vec::with_capacity(tiles.len() * 3);
    for tile in tiles {
        let bounds = tile_bounds_world(&tile);
        let cx = (bounds.min.position.x + bounds.max.position.x) * 0.5;
        let cy = (bounds.min.position.y + bounds.max.position.y) * 0.5;
        let dx = cx - cam_x;
        let dy = cy - cam_y;
        let dist_sq = dx * dx + dy * dy;

        if dist_sq <= refine_radius_sq {
            refined.extend_from_slice(&tile.children());
        } else {
            refined.push(tile);
        }
    }

    refined.sort();
    refined.dedup();
    refined
}

fn sampled_ground_footprint(
    view: &FlatTileView,
    config: &FlatTileSelectionConfig,
) -> Vec<(f64, f64)> {
    let w = view.viewport_width as f64;
    let h = view.viewport_height as f64;
    if w <= 0.0 || h <= 0.0 {
        return Vec::new();
    }

    let edge_steps = config.footprint_edge_steps.max(1);

    let half_fov = view.fov_y / 2.0;
    let max_angle = (view.pitch + half_fov).min(std::f64::consts::FRAC_PI_2 - 0.01);
    let ground_far = view.distance * max_angle.tan().abs().max(1.0);

    // Scale the distance caps with camera altitude so the footprint
    // polygon covers the full viewport during high-altitude zoom-out.
    // The config values act as a floor for close-up views; at high
    // altitude the camera-derived limit takes over.  This matches
    // MapLibre's behaviour where tile selection adapts to altitude.
    let altitude_ground_cap = view.distance * 6.0;
    let max_ground_dist = (ground_far * 2.0)
        .min(config.max_ground_distance.max(altitude_ground_cap));
    let altitude_sky_cap = view.distance * 4.0;
    let max_sky_dist = (view.distance * view.pitch.tan().abs().max(1.0) * 4.0)
        .min(config.max_sky_distance.max(altitude_sky_cap));

    let mut samples = Vec::with_capacity((edge_steps + 1) * 4);
    let mut push_hit = |px: f64, py: f64| {
        let (origin, dir) = screen_to_ray(view, px, py);

        if dir.z.abs() < 1e-12 {
            let xy_len = (dir.x * dir.x + dir.y * dir.y).sqrt();
            if xy_len > 1e-12 {
                samples.push((
                    origin.x + (dir.x / xy_len) * max_sky_dist,
                    origin.y + (dir.y / xy_len) * max_sky_dist,
                ));
            }
            return;
        }

        let t = -origin.z / dir.z;
        if t < 0.0 {
            let xy_len = (dir.x * dir.x + dir.y * dir.y).sqrt();
            if xy_len > 1e-12 {
                samples.push((
                    origin.x + (dir.x / xy_len) * max_sky_dist,
                    origin.y + (dir.y / xy_len) * max_sky_dist,
                ));
            }
            return;
        }

        let t_clamped = t.min(max_ground_dist);
        samples.push((origin.x + dir.x * t_clamped, origin.y + dir.y * t_clamped));
    };

    for i in 0..=edge_steps {
        let t = i as f64 / edge_steps as f64;
        push_hit(t * w, 0.0);
    }
    for i in 1..=edge_steps {
        let t = i as f64 / edge_steps as f64;
        push_hit(w, t * h);
    }
    for i in (0..edge_steps).rev() {
        let t = i as f64 / edge_steps as f64;
        push_hit(t * w, h);
    }
    for i in (1..edge_steps).rev() {
        let t = i as f64 / edge_steps as f64;
        push_hit(0.0, t * h);
    }

    dedupe_nearby_points(samples, 1.0)
}

fn screen_to_ray(view: &FlatTileView, px: f64, py: f64) -> (DVec3, DVec3) {
    let w = view.viewport_width.max(1) as f64;
    let h = view.viewport_height.max(1) as f64;

    let target_world = DVec3::new(
        view.target_world.position.x,
        view.target_world.position.y,
        view.target_world.position.z,
    );
    let view_m = view_matrix(view, target_world);
    let proj_m = perspective_matrix(view);
    let vp_inv = (proj_m * view_m).inverse();

    let ndc_x = (2.0 * px / w) - 1.0;
    let ndc_y = 1.0 - (2.0 * py / h);

    let near_ndc = DVec4::new(ndc_x, ndc_y, -1.0, 1.0);
    let far_ndc = DVec4::new(ndc_x, ndc_y, 1.0, 1.0);

    let near_world = vp_inv * near_ndc;
    let far_world = vp_inv * far_ndc;

    if near_world.w.abs() < 1e-12 || far_world.w.abs() < 1e-12 {
        return (DVec3::ZERO, -DVec3::Z);
    }

    let near = DVec3::new(
        near_world.x / near_world.w,
        near_world.y / near_world.w,
        near_world.z / near_world.w,
    );
    let far = DVec3::new(
        far_world.x / far_world.w,
        far_world.y / far_world.w,
        far_world.z / far_world.w,
    );

    let dir = (far - near).normalize();
    if dir.is_nan() {
        return (DVec3::ZERO, -DVec3::Z);
    }
    (near, dir)
}

fn eye_offset(view: &FlatTileView) -> DVec3 {
    let (sp, cp) = view.pitch.sin_cos();
    let (sy, cy) = view.yaw.sin_cos();
    DVec3::new(
        -view.distance * sp * sy,
        -view.distance * sp * cy,
        view.distance * cp,
    )
}

fn view_matrix(view: &FlatTileView, target_world: DVec3) -> DMat4 {
    let eye = target_world + eye_offset(view);
    const BLEND_RAD: f64 = 0.15;

    let (sy, cy) = view.yaw.sin_cos();
    let yaw_up = DVec3::new(sy, cy, 0.0);
    let right = DVec3::new(cy, -sy, 0.0);
    let look = (target_world - eye).normalize_or_zero();
    let pitched_up = right.cross(look).normalize_or_zero();

    let t = (view.pitch / BLEND_RAD).clamp(0.0, 1.0);
    let up = (pitched_up * t + yaw_up * (1.0 - t)).normalize_or_zero();
    let up = if up.length_squared() < 0.5 { DVec3::Z } else { up };

    DMat4::look_at_rh(eye, target_world, up)
}

fn perspective_matrix(view: &FlatTileView) -> DMat4 {
    let aspect = view.viewport_width as f64 / view.viewport_height.max(1) as f64;
    let near = (view.distance * 0.001).max(0.01);
    let pitch_far_scale = if view.pitch > 0.01 {
        (1.0 / view.pitch.cos().abs().max(0.05)).min(100.0)
    } else {
        1.0
    };
    let far = view.distance * 10.0 * pitch_far_scale;
    DMat4::perspective_rh(view.fov_y, aspect, near, far)
}

fn dedupe_nearby_points(points: Vec<(f64, f64)>, epsilon: f64) -> Vec<(f64, f64)> {
    let mut out = Vec::with_capacity(points.len());
    let eps2 = epsilon * epsilon;
    for p in points {
        let keep = out
            .last()
            .map(|q: &(f64, f64)| {
                let dx = p.0 - q.0;
                let dy = p.1 - q.1;
                dx * dx + dy * dy > eps2
            })
            .unwrap_or(true);
        if keep {
            out.push(p);
        }
    }
    if out.len() >= 2 {
        let first = out[0];
        if let Some(&last) = out.last() {
            let dx = first.0 - last.0;
            let dy = first.1 - last.1;
            if dx * dx + dy * dy <= eps2 {
                out.pop();
            }
        }
    }
    out
}

fn tile_intersects_ground_footprint(tile: TileId, footprint: &[(f64, f64)]) -> bool {
    let bounds = tile_bounds_world(&tile);
    let min_x = bounds.min.position.x;
    let min_y = bounds.min.position.y;
    let max_x = bounds.max.position.x;
    let max_y = bounds.max.position.y;

    if footprint
        .iter()
        .any(|&(x, y)| x >= min_x && x <= max_x && y >= min_y && y <= max_y)
    {
        return true;
    }

    let corners = [(min_x, min_y), (max_x, min_y), (max_x, max_y), (min_x, max_y)];
    if corners.iter().any(|&p| point_in_polygon(p, footprint)) {
        return true;
    }

    for i in 0..footprint.len() {
        let a = footprint[i];
        let b = footprint[(i + 1) % footprint.len()];
        if segment_intersects_aabb(a, b, min_x, min_y, max_x, max_y) {
            return true;
        }
    }

    false
}

fn point_in_polygon(point: (f64, f64), polygon: &[(f64, f64)]) -> bool {
    let (px, py) = point;
    let mut inside = false;
    let mut j = polygon.len() - 1;
    for i in 0..polygon.len() {
        let (xi, yi) = polygon[i];
        let (xj, yj) = polygon[j];
        if (yi > py) != (yj > py) {
            let denom = yj - yi;
            if denom.abs() > 1e-12 {
                let x_at_py = (xj - xi) * (py - yi) / denom + xi;
                if px < x_at_py {
                    inside = !inside;
                }
            }
        }
        j = i;
    }
    inside
}

fn segment_intersects_aabb(
    a: (f64, f64),
    b: (f64, f64),
    min_x: f64,
    min_y: f64,
    max_x: f64,
    max_y: f64,
) -> bool {
    if (a.0 >= min_x && a.0 <= max_x && a.1 >= min_y && a.1 <= max_y)
        || (b.0 >= min_x && b.0 <= max_x && b.1 >= min_y && b.1 <= max_y)
    {
        return true;
    }

    let rect = [
        ((min_x, min_y), (max_x, min_y)),
        ((max_x, min_y), (max_x, max_y)),
        ((max_x, max_y), (min_x, max_y)),
        ((min_x, max_y), (min_x, min_y)),
    ];

    rect.iter().any(|&(r0, r1)| segments_intersect(a, b, r0, r1))
}

fn segments_intersect(a1: (f64, f64), a2: (f64, f64), b1: (f64, f64), b2: (f64, f64)) -> bool {
    fn orient(a: (f64, f64), b: (f64, f64), c: (f64, f64)) -> f64 {
        (b.0 - a.0) * (c.1 - a.1) - (b.1 - a.1) * (c.0 - a.0)
    }
    fn on_segment(a: (f64, f64), b: (f64, f64), p: (f64, f64)) -> bool {
        p.0 >= a.0.min(b.0) - 1e-9
            && p.0 <= a.0.max(b.0) + 1e-9
            && p.1 >= a.1.min(b.1) - 1e-9
            && p.1 <= a.1.max(b.1) + 1e-9
    }

    let o1 = orient(a1, a2, b1);
    let o2 = orient(a1, a2, b2);
    let o3 = orient(b1, b2, a1);
    let o4 = orient(b1, b2, a2);

    if (o1 > 0.0 && o2 < 0.0 || o1 < 0.0 && o2 > 0.0)
        && (o3 > 0.0 && o4 < 0.0 || o3 < 0.0 && o4 > 0.0)
    {
        return true;
    }

    (o1.abs() <= 1e-9 && on_segment(a1, a2, b1))
        || (o2.abs() <= 1e-9 && on_segment(a1, a2, b2))
        || (o3.abs() <= 1e-9 && on_segment(b1, b2, a1))
        || (o4.abs() <= 1e-9 && on_segment(b1, b2, a2))
}

// ---------------------------------------------------------------------------
// TileId
// ---------------------------------------------------------------------------

/// A tile identifier in a slippy-map tile grid (zoom / x / y).
///
/// Zoom level 0 has a single tile covering the world.  At zoom `z` there
/// are `2^z` tiles along each axis, giving `4^z` tiles total.
///
/// `TileId` is `Copy`, `Hash`, and `Ord`, so it can be used directly as a
/// key in `HashMap`, `BTreeMap`, or sorted `Vec`.
#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash, PartialOrd, Ord)]
#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
pub struct TileId {
    /// Zoom level (0-22).
    pub zoom: u8,
    /// Column index (0 is the western-most column).
    pub x: u32,
    /// Row index (0 is the northern-most row).
    pub y: u32,
}

impl TileId {
    /// Create a new tile identifier.
    ///
    /// # Panics (debug only)
    ///
    /// Debug-asserts that `zoom <= 22` and that `x`, `y` are within range
    /// for the zoom level.  In release builds the values are unchecked --
    /// use [`new_checked`](Self::new_checked) when the inputs come from
    /// untrusted sources.
    #[inline]
    pub fn new(zoom: u8, x: u32, y: u32) -> Self {
        debug_assert!(zoom <= MAX_ZOOM, "zoom {zoom} exceeds maximum {MAX_ZOOM}");
        debug_assert!(
            x < Self::axis_tiles(zoom),
            "x={x} out of range for zoom {zoom}"
        );
        debug_assert!(
            y < Self::axis_tiles(zoom),
            "y={y} out of range for zoom {zoom}"
        );
        Self { zoom, x, y }
    }

    /// Checked constructor that returns `None` if parameters are out of range.
    ///
    /// Validates `zoom <= 22` and `x, y < 2^zoom`.
    #[inline]
    pub fn new_checked(zoom: u8, x: u32, y: u32) -> Option<Self> {
        if zoom > MAX_ZOOM {
            return None;
        }
        let n = Self::axis_tiles(zoom);
        if x >= n || y >= n {
            return None;
        }
        Some(Self { zoom, x, y })
    }

    /// Total number of tiles along one axis at this zoom level (`2^zoom`).
    ///
    /// # Contract
    ///
    /// Callers should pass zooms in `0..=MAX_ZOOM`. Higher zooms are not
    /// meaningful for this crate and may overflow shift semantics on some
    /// targets/toolchains.
    #[inline]
    pub fn axis_tiles(zoom: u8) -> u32 {
        1u32 << zoom
    }

    /// Return the parent tile (one zoom level up).  Returns `None` at zoom 0.
    ///
    /// The parent is the tile that fully contains this tile at the next
    /// coarser zoom level.  Used by `TileManager` for fallback rendering
    /// while a tile is still loading.
    #[inline]
    pub fn parent(&self) -> Option<TileId> {
        if self.zoom == 0 {
            None
        } else {
            Some(TileId {
                zoom: self.zoom - 1,
                x: self.x / 2,
                y: self.y / 2,
            })
        }
    }

    /// Return the four children (one zoom level down).
    ///
    /// The children are ordered: top-left, top-right, bottom-left,
    /// bottom-right.
    ///
    /// # Note
    ///
    /// Calling this at `zoom == MAX_ZOOM` produces tiles at zoom 23,
    /// which is beyond the validated range.  The engine never does this,
    /// but callers at the boundary should be aware.
    #[inline]
    pub fn children(&self) -> [TileId; 4] {
        let z = self.zoom + 1;
        let x = self.x * 2;
        let y = self.y * 2;
        [
            TileId { zoom: z, x, y },
            TileId {
                zoom: z,
                x: x + 1,
                y,
            },
            TileId {
                zoom: z,
                x,
                y: y + 1,
            },
            TileId {
                zoom: z,
                x: x + 1,
                y: y + 1,
            },
        ]
    }

    /// Encode this tile as a [Bing Maps-style quadkey][qk] string.
    ///
    /// Returns an empty string at zoom 0.  Each character is `'0'`-`'3'`,
    /// encoding two bits per level: bit 0 from X, bit 1 from Y.
    ///
    /// [qk]: https://learn.microsoft.com/en-us/bingmaps/articles/bing-maps-tile-system
    pub fn quadkey(&self) -> String {
        let mut key = String::with_capacity(self.zoom as usize);
        for i in (1..=self.zoom).rev() {
            let mut digit: u8 = b'0';
            let mask = 1u32 << (i - 1);
            if (self.x & mask) != 0 {
                digit += 1;
            }
            if (self.y & mask) != 0 {
                digit += 2;
            }
            key.push(digit as char);
        }
        key
    }

    /// Decode a [Bing Maps-style quadkey][qk] string into a tile ID.
    ///
    /// Returns `None` if the string contains invalid characters or if the
    /// resulting zoom exceeds [`MAX_ZOOM`].
    ///
    /// An empty quadkey decodes to the root tile `0/0/0`.
    ///
    /// [qk]: https://learn.microsoft.com/en-us/bingmaps/articles/bing-maps-tile-system
    pub fn from_quadkey(key: &str) -> Option<Self> {
        let zoom = key.len() as u8;
        if zoom > MAX_ZOOM {
            return None;
        }

        let mut x = 0u32;
        let mut y = 0u32;

        for (i, ch) in key.bytes().enumerate() {
            let bit = 1u32 << (zoom as usize - 1 - i);
            match ch {
                b'0' => {}
                b'1' => x |= bit,
                b'2' => y |= bit,
                b'3' => {
                    x |= bit;
                    y |= bit;
                }
                _ => return None,
            }
        }

        Self::new_checked(zoom, x, y)
    }

    /// Return the neighbouring tile in a cardinal direction at the same
    /// zoom level.
    ///
    /// `dx` and `dy` are offsets: `(-1, 0)` = west, `(1, 0)` = east,
    /// `(0, -1)` = north, `(0, 1)` = south.  Diagonal offsets such as
    /// `(-1, -1)` (north-west) are also supported.
    ///
    /// The X axis wraps around (longitude wrap): moving west from column 0
    /// yields column `2^zoom - 1`, and vice versa.  The Y axis does **not**
    /// wrap: moving north from row 0 or south from the last row returns
    /// `None` (there is no tile beyond the poles).
    ///
    /// Returns `None` if the resulting tile would be out of the valid
    /// Y range.
    #[inline]
    pub fn neighbor(&self, dx: i32, dy: i32) -> Option<TileId> {
        let n = Self::axis_tiles(self.zoom) as i64;
        // X wraps around (longitude wrapping).
        let nx = ((self.x as i64 + dx as i64) % n + n) % n;
        // Y does not wrap (no tiles beyond the poles).
        let ny = self.y as i64 + dy as i64;
        if ny < 0 || ny >= n {
            return None;
        }
        Some(TileId {
            zoom: self.zoom,
            x: nx as u32,
            y: ny as u32,
        })
    }
}

impl fmt::Display for TileId {
    /// Formats as `"zoom/x/y"` (the standard OSM URL path component).
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        write!(f, "{}/{}/{}", self.zoom, self.x, self.y)
    }
}

// ---------------------------------------------------------------------------
// TileCoord
// ---------------------------------------------------------------------------

/// A fractional position within the tile grid at a given zoom level.
///
/// Where [`TileId`] addresses an integer tile, `TileCoord` addresses an
/// exact point within the grid -- for example `(zoom=10, x=560.7, y=342.3)`
/// means "70% across tile column 560, 30% down tile row 342".
///
/// Useful for sub-tile precision when mapping world coordinates to tiles.
#[derive(Debug, Clone, Copy, PartialEq)]
#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
pub struct TileCoord {
    /// Zoom level.
    pub zoom: u8,
    /// Fractional column position (0.0 = west edge, `2^zoom` = east edge).
    pub x: f64,
    /// Fractional row position (0.0 = north edge, `2^zoom` = south edge).
    pub y: f64,
}

impl TileCoord {
    /// Create a new fractional tile coordinate.
    #[inline]
    pub fn new(zoom: u8, x: f64, y: f64) -> Self {
        Self { zoom, x, y }
    }

    /// Snap to the integer tile that contains this coordinate.
    ///
    /// Values are clamped to the valid tile range `[0, 2^zoom - 1]`, so
    /// this always returns a valid `TileId` -- even if `x` or `y` are
    /// slightly out of range due to floating-point rounding.
    #[inline]
    pub fn tile_id(&self) -> TileId {
        let n = TileId::axis_tiles(self.zoom);
        // Clamp to [0, n-1].  The lower clamp guards against negative
        // values (which `as u32` would saturate to 0 on stable Rust,
        // but explicitly clamping makes intent clear and safe).
        let x = (self.x.max(0.0) as u32).min(n.saturating_sub(1));
        let y = (self.y.max(0.0) as u32).min(n.saturating_sub(1));
        TileId {
            zoom: self.zoom,
            x,
            y,
        }
    }
}

impl fmt::Display for TileCoord {
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        write!(f, "{}/{:.3}/{:.3}", self.zoom, self.x, self.y)
    }
}

// ---------------------------------------------------------------------------
// Conversion functions
// ---------------------------------------------------------------------------

/// Convert a geographic coordinate to a fractional tile coordinate at the
/// given zoom.
///
/// Uses the standard [Mercator tile formula][osm]:
///
/// ```text
/// x = (lon + 180) / 360 * 2^zoom
/// y = (1 - ln(tan(lat) + sec(lat)) / pi) / 2 * 2^zoom
/// ```
///
/// The input must be within the Web Mercator valid range (latitude
/// approximately +/-85.06 degrees).  Latitudes at or beyond the poles
/// produce infinite or NaN `y` values.
///
/// [osm]: https://wiki.openstreetmap.org/wiki/Slippy_map_tilenames#Mathematics
pub fn geo_to_tile(geo: &GeoCoord, zoom: u8) -> TileCoord {
    let n = TileId::axis_tiles(zoom) as f64;
    let lat_rad = geo.lat.to_radians();
    let x = (geo.lon + 180.0) / 360.0 * n;
    let y = (1.0 - (lat_rad.tan() + 1.0 / lat_rad.cos()).ln() / PI) / 2.0 * n;
    TileCoord::new(zoom, x, y)
}

/// Checked variant of [`geo_to_tile`].
///
/// Returns `None` if `zoom > MAX_ZOOM` or if `geo` is outside Web Mercator
/// latitude range.
pub fn geo_to_tile_checked(geo: &GeoCoord, zoom: u8) -> Option<TileCoord> {
    if zoom > MAX_ZOOM || !geo.is_web_mercator_valid() {
        return None;
    }
    Some(geo_to_tile(geo, zoom))
}

/// Convert a tile ID (top-left corner) to the geographic coordinate of its
/// north-west corner.
///
/// This is the inverse of snapping `geo_to_tile(...).tile_id()` at integer
/// boundaries.
pub fn tile_to_geo(tile: &TileId) -> GeoCoord {
    tile_xy_to_geo(tile.zoom, tile.x as f64, tile.y as f64)
}

/// Convert fractional tile coordinates to geographic.
///
/// Accepts fractional `(x, y)` values that may sit on or beyond the tile
/// grid boundary (for example `x = n` to get the east edge of the last
/// column). This is useful for deriving geographic tile corners before
/// reprojecting them into other planar world spaces.
pub fn tile_xy_to_geo(zoom: u8, x: f64, y: f64) -> GeoCoord {
    let n = TileId::axis_tiles(zoom) as f64;
    let lon = x / n * 360.0 - 180.0;
    let lat_rad = (PI * (1.0 - 2.0 * y / n)).sinh().atan();
    GeoCoord::from_lat_lon(lat_rad.to_degrees(), lon)
}

/// Return the world-space bounding box (Web Mercator meters) for a tile.
///
/// The returned `WorldBounds` has `min` at the south-west corner and `max`
/// at the north-east corner, matching the convention used by therenderer's
/// quad mesh builder.
pub fn tile_bounds_world(tile: &TileId) -> WorldBounds {
    // NW corner of this tile, SE corner is the NW of the next tile (+1,+1).
    let nw = tile_to_geo(tile);
    let se = tile_xy_to_geo(tile.zoom, tile.x as f64 + 1.0, tile.y as f64 + 1.0);
    let min_world = WebMercator::project(&GeoCoord::from_lat_lon(se.lat, nw.lon));
    let max_world = WebMercator::project(&GeoCoord::from_lat_lon(nw.lat, se.lon));
    WorldBounds::new(min_world, max_world)
}

/// Return all tiles that intersect the given world-space bounding box at a
/// zoom level.
///
/// # Antimeridian handling
///
/// When the bounding box spans the antimeridian (min.x > max.x in world
/// space after unprojection), the X iteration wraps around, producing tiles
/// on both sides of the dateline.  The Y range never wraps because latitude
/// is bounded by the Mercator limit.
///
/// # Ordering and duplicates
///
/// Output is deterministic and row-major: `y` increases outermost, `x`
/// advances left-to-right within each row (with wrap when needed).
/// Each `(z, x, y)` appears at most once.
///
/// # Allocation
///
/// Returns an owned `Vec` with capacity pre-allocated from the tile count.
/// At zoom 14 a typical viewport might produce ~200 tiles; at zoom 18 with
/// a 4K display, up to ~4000.
pub fn visible_tiles(bounds: &WorldBounds, zoom: u8) -> Vec<TileId> {
    let extent = WebMercator::max_extent();
    let world_size = WebMercator::world_size();
    let n = TileId::axis_tiles(zoom);

    // If the viewport spans the full Mercator world width (or more), all
    // columns at this zoom are visible regardless of wrapped longitude.
    let spans_full_world_x = (bounds.max.position.x - bounds.min.position.x) >= world_size - 1.0;

    let x_min_raw: i64;
    let x_count: u32;
    if spans_full_world_x {
        x_min_raw = 0;
        x_count = n;
    } else {
        let tile_world_width = world_size / n as f64;
        x_min_raw = ((bounds.min.position.x + extent) / tile_world_width).floor() as i64;
        let x_max_raw = ((bounds.max.position.x + extent) / tile_world_width).floor() as i64;
        x_count = (x_max_raw - x_min_raw + 1).clamp(0, n as i64) as u32;
    }

    let world_y_to_tile_y = |world_y: f64| {
        (((extent - world_y.clamp(-extent, extent)) / world_size) * n as f64)
            .floor()
            .clamp(0.0, n.saturating_sub(1) as f64) as u32
    };

    let y_min = world_y_to_tile_y(bounds.max.position.y);
    let y_max = world_y_to_tile_y(bounds.min.position.y);

    let mut tiles = Vec::with_capacity((x_count * (y_max - y_min + 1)) as usize);
    for y in y_min..=y_max {
        for i in 0..x_count {
            let x = if spans_full_world_x {
                i
            } else {
                (x_min_raw + i as i64).rem_euclid(n as i64) as u32
            };
            tiles.push(TileId { zoom, x, y });
        }
    }
    tiles
}

/// Checked variant of [`visible_tiles`].
pub fn visible_tiles_checked(bounds: &WorldBounds, zoom: u8) -> Option<Vec<TileId>> {
    if zoom > MAX_ZOOM {
        return None;
    }
    Some(visible_tiles(bounds, zoom))
}

/// Return tiles at multiple zoom levels based on distance from the camera.
pub fn visible_tiles_lod(
    bounds: &WorldBounds,
    base_zoom: u8,
    camera_world: (f64, f64),
    near_threshold: f64,
    mid_threshold: f64,
    max_tiles: usize,
) -> Vec<TileId> {
    use std::collections::HashSet;

    let near_zoom = (base_zoom + 1).min(MAX_ZOOM);
    let far_zoom = base_zoom.saturating_sub(1);

    let near_tiles = visible_tiles(bounds, near_zoom);
    let mid_tiles = visible_tiles(bounds, base_zoom);
    let far_tiles = visible_tiles(bounds, far_zoom);

    let near_sq = near_threshold * near_threshold;
    let mid_sq = mid_threshold * mid_threshold;

    fn tile_dist_sq(tile: &TileId, cam: (f64, f64)) -> f64 {
        let b = tile_bounds_world(tile);
        let cx = (b.min.position.x + b.max.position.x) * 0.5;
        let cy = (b.min.position.y + b.max.position.y) * 0.5;
        let dx = cx - cam.0;
        let dy = cy - cam.1;
        dx * dx + dy * dy
    }

    let mut result: Vec<(TileId, f64)> = Vec::new();
    let mut seen = HashSet::new();

    for tile in &near_tiles {
        let d2 = tile_dist_sq(tile, camera_world);
        if d2 <= near_sq && seen.insert(*tile) {
            result.push((*tile, d2));
        }
    }

    let near_parent_set: HashSet<TileId> = if near_zoom > base_zoom {
        result.iter().filter_map(|(t, _)| t.parent()).collect()
    } else {
        HashSet::new()
    };

    for tile in &mid_tiles {
        let d2 = tile_dist_sq(tile, camera_world);
        if d2 <= mid_sq && seen.insert(*tile) {
            if near_parent_set.contains(tile) {
                let children = tile.children();
                if children.iter().all(|c| seen.contains(c)) {
                    continue;
                }
            }
            result.push((*tile, d2));
        }
    }

    let mid_parent_set: HashSet<TileId> = if base_zoom > far_zoom {
        result
            .iter()
            .filter(|(t, _)| t.zoom == base_zoom)
            .filter_map(|(t, _)| t.parent())
            .collect()
    } else {
        HashSet::new()
    };

    for tile in &far_tiles {
        let d2 = tile_dist_sq(tile, camera_world);
        if d2 > mid_sq && seen.insert(*tile) {
            if mid_parent_set.contains(tile) {
                let children = tile.children();
                if children.iter().all(|c| seen.contains(c)) {
                    continue;
                }
            }
            result.push((*tile, d2));
        }
    }

    result.sort_by(|a, b| a.1.partial_cmp(&b.1).unwrap_or(std::cmp::Ordering::Equal));
    result.truncate(max_tiles);

    result.into_iter().map(|(t, _)| t).collect()
}

/// Return visible tiles by depth-first quadtree traversal with frustum culling.
pub fn visible_tiles_frustum(
    frustum: &crate::frustum::Frustum,
    target_zoom: u8,
    max_tiles: usize,
    camera_world: (f64, f64),
) -> Vec<TileId> {
    let target_zoom = target_zoom.min(MAX_ZOOM);

    struct StackEntry {
        tile: TileId,
        fully_visible: bool,
    }

    let mut stack = Vec::with_capacity(64);
    stack.push(StackEntry {
        tile: TileId::new(0, 0, 0),
        fully_visible: false,
    });

    let mut result: Vec<(TileId, f64)> = Vec::new();

    while let Some(entry) = stack.pop() {
        let tile = entry.tile;
        let bounds = tile_bounds_world(&tile);

        if !entry.fully_visible {
            let mut all_inside = true;
            let mut any_outside = false;
            let min = bounds.min.position;
            let max = bounds.max.position;

            for plane in frustum.planes() {
                let px = if plane.normal()[0] >= 0.0 { max.x } else { min.x };
                let py = if plane.normal()[1] >= 0.0 { max.y } else { min.y };
                let pz = if plane.normal()[2] >= 0.0 { max.z } else { min.z };

                if plane.distance_to_point(px, py, pz) < 0.0 {
                    any_outside = true;
                    break;
                }

                let nx = if plane.normal()[0] >= 0.0 { min.x } else { max.x };
                let ny = if plane.normal()[1] >= 0.0 { min.y } else { max.y };
                let nz = if plane.normal()[2] >= 0.0 { min.z } else { max.z };

                if plane.distance_to_point(nx, ny, nz) < 0.0 {
                    all_inside = false;
                }
            }

            if any_outside {
                continue;
            }

            if tile.zoom >= target_zoom {
                let cx = (bounds.min.position.x + bounds.max.position.x) * 0.5;
                let cy = (bounds.min.position.y + bounds.max.position.y) * 0.5;
                let dx = cx - camera_world.0;
                let dy = cy - camera_world.1;
                result.push((tile, dx * dx + dy * dy));
                continue;
            }

            for child in tile.children().iter().rev() {
                stack.push(StackEntry {
                    tile: *child,
                    fully_visible: all_inside,
                });
            }
        } else {
            if tile.zoom >= target_zoom {
                let cx = (bounds.min.position.x + bounds.max.position.x) * 0.5;
                let cy = (bounds.min.position.y + bounds.max.position.y) * 0.5;
                let dx = cx - camera_world.0;
                let dy = cy - camera_world.1;
                result.push((tile, dx * dx + dy * dy));
                continue;
            }

            for child in tile.children().iter().rev() {
                stack.push(StackEntry {
                    tile: *child,
                    fully_visible: true,
                });
            }
        }
    }

    result.sort_by(|a, b| a.1.partial_cmp(&b.1).unwrap_or(std::cmp::Ordering::Equal));
    result.truncate(max_tiles);
    result.into_iter().map(|(t, _)| t).collect()
}

/// Return visible flat raster tiles for a pitched perspective view,
/// capped to a maximum count and prioritised by distance to the camera target.
pub fn visible_tiles_flat_view_capped(
    bounds: &WorldBounds,
    zoom: u8,
    view: &FlatTileView,
    max_tiles: usize,
) -> Vec<TileId> {
    visible_tiles_flat_view_capped_with_config(
        bounds,
        zoom,
        view,
        &FlatTileSelectionConfig::default(),
        max_tiles,
    )
}

/// Return visible flat raster tiles for a pitched perspective view,
/// capped to a maximum count and prioritised by distance to the camera target,
/// using an explicit flat-tile selection policy.
pub fn visible_tiles_flat_view_capped_with_config(
    bounds: &WorldBounds,
    zoom: u8,
    view: &FlatTileView,
    config: &FlatTileSelectionConfig,
    max_tiles: usize,
) -> Vec<TileId> {
    let mut tiles = visible_tiles_flat_view_with_config(bounds, zoom, view, config);
    if tiles.len() <= max_tiles {
        return tiles;
    }

    let cam_x = view.target_world.position.x;
    let cam_y = view.target_world.position.y;
    tiles.sort_by(|a, b| {
        let ad = tile_distance_sq(*a, cam_x, cam_y);
        let bd = tile_distance_sq(*b, cam_x, cam_y);
        ad.partial_cmp(&bd)
            .unwrap_or(std::cmp::Ordering::Equal)
            .then_with(|| a.zoom.cmp(&b.zoom))
            .then_with(|| a.y.cmp(&b.y))
            .then_with(|| a.x.cmp(&b.x))
    });
    tiles.truncate(max_tiles);
    tiles
}

fn tile_distance_sq(tile: TileId, cam_x: f64, cam_y: f64) -> f64 {
    let bounds = tile_bounds_world(&tile);
    let cx = (bounds.min.position.x + bounds.max.position.x) * 0.5;
    let cy = (bounds.min.position.y + bounds.max.position.y) * 0.5;
    let dx = cx - cam_x;
    let dy = cy - cam_y;
    dx * dx + dy * dy
}

/// Return visible flat raster tiles for a pitched perspective view with
/// an extra near-camera refinement pass, capped to a maximum count.
pub fn visible_tiles_flat_view_refined_capped(
    bounds: &WorldBounds,
    zoom: u8,
    view: &FlatTileView,
    max_tiles: usize,
) -> Vec<TileId> {
    visible_tiles_flat_view_refined_capped_with_config(
        bounds,
        zoom,
        view,
        &FlatTileSelectionConfig::default(),
        max_tiles,
    )
}

/// Configurable form of [`visible_tiles_flat_view_refined_capped`].
pub fn visible_tiles_flat_view_refined_capped_with_config(
    bounds: &WorldBounds,
    zoom: u8,
    view: &FlatTileView,
    config: &FlatTileSelectionConfig,
    max_tiles: usize,
) -> Vec<TileId> {
    let footprint_filtered = visible_tiles_flat_view_with_config(bounds, zoom, view, config);
    if footprint_filtered.is_empty() {
        return footprint_filtered;
    }

    let mut refined = footprint_filtered;
    if refined.len() < max_tiles
        && zoom < MAX_ZOOM
        && view.pitch > config.footprint_pitch_threshold_rad
    {
        let footprint = sampled_ground_footprint(view, config);
        if footprint.len() >= 3 {
            refined = refine_nearby_flat_tiles(refined, zoom, view, &footprint);
        }
    }

    if refined.len() <= max_tiles {
        return refined;
    }

    let cam_x = view.target_world.position.x;
    let cam_y = view.target_world.position.y;
    refined.sort_by(|a, b| {
        let ad = tile_distance_sq(*a, cam_x, cam_y);
        let bd = tile_distance_sq(*b, cam_x, cam_y);
        ad.partial_cmp(&bd)
            .unwrap_or(std::cmp::Ordering::Equal)
            .then_with(|| a.zoom.cmp(&b.zoom))
            .then_with(|| a.y.cmp(&b.y))
            .then_with(|| a.x.cmp(&b.x))
    });
    refined.truncate(max_tiles);
    refined
}

// ---------------------------------------------------------------------------
// Covering-tiles: MapLibre-equivalent quadtree traversal with variable zoom
// ---------------------------------------------------------------------------

/// Options controlling the covering-tiles quadtree traversal.
///
/// This is the Rustial equivalent of MapLibre's `CoveringTilesOptionsInternal`.
/// It couples frustum culling, per-tile variable zoom heuristics, wrap
/// handling, and tile-budget enforcement into a single covering-selection
/// call.
#[derive(Debug, Clone)]
pub struct CoveringTilesOptions {
    /// Smallest allowed tile zoom (inclusive).  Tiles at zooms below this
    /// are never emitted but may still be traversed.
    pub min_zoom: u8,
    /// Largest allowed tile zoom (inclusive).  The traversal never descends
    /// past this zoom regardless of the per-tile heuristic.
    pub max_zoom: u8,
    /// Whether to round (`true`) or floor (`false`) the computed tile zoom.
    pub round_zoom: bool,
    /// Source tile size in screen pixels (typically 256 or 512).
    pub tile_size: u32,
    /// Maximum number of result tiles.
    pub max_tiles: usize,
    /// Whether to allow per-tile variable zoom heuristics at steep pitch.
    /// When `false` all tiles use the single nominal zoom level.
    pub allow_variable_zoom: bool,
    /// Whether to emit world-copy tiles across the antimeridian.
    pub render_world_copies: bool,
}

impl Default for CoveringTilesOptions {
    fn default() -> Self {
        Self {
            min_zoom: 0,
            max_zoom: MAX_ZOOM,
            round_zoom: false,
            tile_size: 256,
            max_tiles: 512,
            allow_variable_zoom: true,
            render_world_copies: true,
        }
    }
}

/// Camera state needed by the covering-tiles traversal.
///
/// All distances / positions are in **Mercator normalised coordinates**
/// (range 0..1), matching the MapLibre convention where `worldSize` = 1.
/// The caller is responsible for converting meter-based engine coordinates
/// before calling [`visible_tiles_covering`].
#[derive(Debug, Clone, Copy)]
pub struct CoveringCamera {
    /// Camera position X in Mercator normalised coords.
    pub camera_x: f64,
    /// Camera position Y in Mercator normalised coords.
    pub camera_y: f64,
    /// Absolute vertical distance from camera to map center, in Mercator
    /// normalised coords.
    pub camera_to_center_z: f64,
    /// Center point X in Mercator normalised coords.
    pub center_x: f64,
    /// Center point Y in Mercator normalised coords.
    pub center_y: f64,
    /// Camera pitch in radians.
    pub pitch_rad: f64,
    /// Vertical field-of-view in **degrees** (matching MapLibre convention).
    pub fov_deg: f64,
    /// The logical zoom level requested at the center of the viewport
    /// (before tile-size adjustment).
    pub zoom: f64,
    /// Display tile size in screen pixels (typically 256).
    pub display_tile_size: u32,
}

/// Compute the per-tile desired zoom using a pitch-aware heuristic.
///
/// This is a direct Rust port of MapLibre's `createCalculateTileZoomFunction`
/// with `maxZoomLevelsOnScreen = 9.314` and `tileCountMaxMinRatio = 3.0`.
///
/// # Arguments
///
/// * `requested_center_zoom` - nominal zoom at the center.
/// * `dist_2d` - 2D distance from camera to tile center (Mercator units).
/// * `dist_z` - vertical distance from camera to center (Mercator units).
/// * `dist_center_3d` - 3D distance from camera to center (Mercator units).
/// * `fov_deg` - camera vertical FOV in degrees.
fn default_calculate_tile_zoom(
    requested_center_zoom: f64,
    dist_2d: f64,
    dist_z: f64,
    dist_center_3d: f64,
    fov_deg: f64,
) -> f64 {
    const MAX_ZOOM_LEVELS_ON_SCREEN: f64 = 9.314;
    const TILE_COUNT_MAX_MIN_RATIO: f64 = 3.0;
    const MAX_MERCATOR_HORIZON_ANGLE: f64 = 85.051129; // degrees

    fn scale_zoom(s: f64) -> f64 {
        s.log2()
    }

    fn integral_cos_x_by_p(p: f64, x1: f64, x2: f64) -> f64 {
        let num_points = 10usize;
        let dx = (x2 - x1) / num_points as f64;
        let mut sum = 0.0;
        for i in 0..num_points {
            let x = x1 + (i as f64 + 0.5) / num_points as f64 * (x2 - x1);
            sum += dx * x.cos().powf(p);
        }
        sum
    }

    let horizon_rad = (MAX_MERCATOR_HORIZON_ANGLE).to_radians();
    let fov_rad = fov_deg.to_radians();

    let pitch_tile_loading_behavior = 2.0
        * ((MAX_ZOOM_LEVELS_ON_SCREEN - 1.0)
            / scale_zoom(
                (horizon_rad - fov_rad).cos()
                    / horizon_rad.cos(),
            )
            - 1.0);

    let center_pitch = if dist_center_3d > 1e-15 {
        (dist_z / dist_center_3d).clamp(-1.0, 1.0).acos()
    } else {
        0.0
    };

    let half_fov_rad = fov_rad / 2.0;
    let _tile_count_pitch0 =
        2.0 * integral_cos_x_by_p(pitch_tile_loading_behavior - 1.0, 0.0, half_fov_rad);
    let highest_pitch = horizon_rad.min(center_pitch + half_fov_rad);
    let lowest_pitch = highest_pitch.min(center_pitch - half_fov_rad);
    let tile_count =
        integral_cos_x_by_p(pitch_tile_loading_behavior - 1.0, lowest_pitch, highest_pitch);
    let tile_count_pitch0 =
        2.0 * integral_cos_x_by_p(pitch_tile_loading_behavior - 1.0, 0.0, half_fov_rad);

    let this_tile_pitch = if dist_z.abs() > 1e-15 {
        (dist_2d / dist_z).atan()
    } else {
        std::f64::consts::FRAC_PI_2
    };
    let dist_tile_3d = (dist_2d * dist_2d + dist_z * dist_z).sqrt();

    let mut desired_z = requested_center_zoom;
    if dist_tile_3d > 1e-15 {
        desired_z += scale_zoom(
            dist_center_3d
                / dist_tile_3d
                / half_fov_rad.cos().max(0.5),
        );
    }
    desired_z += pitch_tile_loading_behavior * scale_zoom(this_tile_pitch.cos()) / 2.0;
    desired_z -=
        scale_zoom((tile_count / tile_count_pitch0 / TILE_COUNT_MAX_MIN_RATIO).max(1.0)) / 2.0;

    desired_z
}

/// Return the nominal covering zoom level for a source tile size.
///
/// Equivalent to MapLibre's `coveringZoomLevel`.
fn covering_zoom_level(cam: &CoveringCamera, opts: &CoveringTilesOptions, round: bool) -> f64 {
    fn scale_zoom(s: f64) -> f64 {
        s.log2()
    }
    let raw = cam.zoom + scale_zoom(cam.display_tile_size as f64 / opts.tile_size as f64);
    let z = if round { raw.round() } else { raw.floor() };
    z.max(0.0)
}

/// 2D distance from a point to the nearest point on a tile (Mercator normalised).
fn dist_to_tile_2d(cx: f64, cy: f64, tx: u32, ty: u32, z: u8) -> f64 {
    let n = (1u64 << z as u64) as f64;
    let inv = 1.0 / n;
    let tile_min_x = tx as f64 * inv;
    let tile_min_y = ty as f64 * inv;
    let tile_max_x = (tx + 1) as f64 * inv;
    let tile_max_y = (ty + 1) as f64 * inv;

    let dx = if cx < tile_min_x {
        tile_min_x - cx
    } else if cx > tile_max_x {
        cx - tile_max_x
    } else {
        0.0
    };
    let dy = if cy < tile_min_y {
        tile_min_y - cy
    } else if cy > tile_max_y {
        cy - tile_max_y
    } else {
        0.0
    };
    (dx * dx + dy * dy).sqrt()
}

/// Return visible tiles using MapLibre-equivalent depth-first quadtree
/// traversal with per-tile variable zoom heuristics and frustum culling.
///
/// This is the production covering-tiles path that matches MapLibre's
/// `coveringTiles()` algorithm:
///
/// 1. Derive a nominal target zoom from the camera zoom and source tile size.
/// 2. Perform a depth-first quadtree traversal from zoom 0.
/// 3. At each node, test visibility against the camera frustum.
/// 4. Optionally compute a per-tile desired zoom based on the tile's 3D
///    distance and pitch angle (variable zoom).
/// 5. When a tile reaches its desired zoom, emit it as a result.
/// 6. Optionally traverse world copies for antimeridian wrap.
///
/// Results are sorted by ascending distance from the center point and
/// capped to `opts.max_tiles`.
pub fn visible_tiles_covering(
    frustum: &crate::frustum::Frustum,
    cam: &CoveringCamera,
    opts: &CoveringTilesOptions,
) -> Vec<TileId> {
    let desired_z = covering_zoom_level(cam, opts, opts.round_zoom);
    let nominal_z = (desired_z as u8).min(opts.max_zoom).max(0);
    let num_tiles_f = (1u64 << nominal_z as u64) as f64;

    let center_tx = num_tiles_f * cam.center_x;
    let center_ty = num_tiles_f * cam.center_y;

    let dist_center_2d =
        ((cam.center_x - cam.camera_x).powi(2) + (cam.center_y - cam.camera_y).powi(2)).sqrt();
    let dist_z = cam.camera_to_center_z;
    let dist_center_3d = (dist_center_2d * dist_center_2d + dist_z * dist_z).sqrt();

    let requested_center_zoom = cam.zoom
        + if cam.display_tile_size > 0 && opts.tile_size > 0 {
            (cam.display_tile_size as f64 / opts.tile_size as f64).log2()
        } else {
            0.0
        };

    struct StackEntry {
        zoom: u8,
        x: u32,
        y: u32,
        wrap: i32,
        fully_visible: bool,
    }

    let mut stack: Vec<StackEntry> = Vec::with_capacity(64);

    // World copies for antimeridian wrap handling
    if opts.render_world_copies {
        for i in 1..=3i32 {
            stack.push(StackEntry {
                zoom: 0,
                x: 0,
                y: 0,
                wrap: -i,
                fully_visible: false,
            });
            stack.push(StackEntry {
                zoom: 0,
                x: 0,
                y: 0,
                wrap: i,
                fully_visible: false,
            });
        }
    }
    stack.push(StackEntry {
        zoom: 0,
        x: 0,
        y: 0,
        wrap: 0,
        fully_visible: false,
    });

    struct ResultEntry {
        tile: TileId,
        wrap: i32,
        dist_sq: f64,
    }

    let mut result: Vec<ResultEntry> = Vec::new();
    let world_size = WebMercator::world_size();

    while let Some(entry) = stack.pop() {
        let z = entry.zoom;
        let x = entry.x;
        let y = entry.y;
        let mut fully_visible = entry.fully_visible;

        // Compute the world-space AABB for this tile (with wrap offset)
        let tile_for_bounds = TileId {
            zoom: z,
            x: x % TileId::axis_tiles(z).max(1),
            y: y.min(TileId::axis_tiles(z).saturating_sub(1)),
        };
        let base_bounds = tile_bounds_world(&tile_for_bounds);
        let wrap_offset = entry.wrap as f64 * world_size;
        let bounds = WorldBounds::new(
            WorldCoord::new(
                base_bounds.min.position.x + wrap_offset,
                base_bounds.min.position.y,
                0.0,
            ),
            WorldCoord::new(
                base_bounds.max.position.x + wrap_offset,
                base_bounds.max.position.y,
                0.0,
            ),
        );

        if !fully_visible {
            if !frustum.intersects_aabb(&bounds) {
                continue;
            }
            // Check if fully inside (all corners inside all planes)
            let min = bounds.min.position;
            let max = bounds.max.position;
            let mut all_inside = true;
            for plane in frustum.planes() {
                let nx = if plane.normal()[0] >= 0.0 { min.x } else { max.x };
                let ny = if plane.normal()[1] >= 0.0 { min.y } else { max.y };
                let nz = if plane.normal()[2] >= 0.0 { min.z } else { max.z };
                if plane.distance_to_point(nx, ny, nz) < 0.0 {
                    all_inside = false;
                    break;
                }
            }
            fully_visible = all_inside;
        }

        // Compute the per-tile desired zoom
        let this_tile_desired_z = if opts.allow_variable_zoom && cam.pitch_rad > 0.05 {
            let d2d = dist_to_tile_2d(cam.camera_x, cam.camera_y, x, y, z);
            let z_val = default_calculate_tile_zoom(
                requested_center_zoom,
                d2d,
                dist_z,
                dist_center_3d,
                cam.fov_deg,
            );
            let z_rounded = if opts.round_zoom {
                z_val.round()
            } else {
                z_val.floor()
            };
            (z_rounded.max(0.0) as u8).min(opts.max_zoom)
        } else {
            nominal_z
        };

        // Have we reached the target depth for this tile?
        if z >= this_tile_desired_z {
            if z < opts.min_zoom {
                continue;
            }
            // Distance from center for sorting
            let dz_shift = nominal_z.saturating_sub(z);
            let tile_center_x = (x as f64 + 0.5) * (1u64 << dz_shift) as f64;
            let tile_center_y = (y as f64 + 0.5) * (1u64 << dz_shift) as f64;
            let dx = center_tx - tile_center_x;
            let dy = center_ty - tile_center_y;

            // Wrap x into valid tile range
            let n = TileId::axis_tiles(z);
            let wrapped_x = if n > 0 {
                ((x as i64).rem_euclid(n as i64)) as u32
            } else {
                0
            };

            result.push(ResultEntry {
                tile: TileId {
                    zoom: z,
                    x: wrapped_x,
                    y: y.min(n.saturating_sub(1)),
                },
                wrap: entry.wrap,
                dist_sq: dx * dx + dy * dy,
            });
            continue;
        }

        // Subdivide into 4 children
        let child_z = z + 1;
        if child_z > MAX_ZOOM {
            continue;
        }
        for i in (0..4u32).rev() {
            let cx = (x << 1) + (i % 2);
            let cy = (y << 1) + (i >> 1);
            stack.push(StackEntry {
                zoom: child_z,
                x: cx,
                y: cy,
                wrap: entry.wrap,
                fully_visible,
            });
        }
    }

    // Sort by distance and cap.
    //
    // Important: when `render_world_copies` is enabled, the same canonical
    // (z/x/y) can be visible in multiple wrapped worlds at once. Collapsing
    // them into a single entry (dedup by z/x/y) can cause tiles to disappear
    // as the camera crosses world boundaries, because the renderer needs
    // both wrapped copies to stay concurrently visible.
    result.sort_by(|a, b| {
        a.dist_sq
            .partial_cmp(&b.dist_sq)
            .unwrap_or(std::cmp::Ordering::Equal)
    });

    // Deduplicate exact duplicates (same wrap + canonical id) only.
    let mut seen = std::collections::HashSet::new();
    let mut final_tiles = Vec::with_capacity(result.len().min(opts.max_tiles));
    for entry in result {
        if seen.insert((entry.wrap, entry.tile.zoom, entry.tile.x, entry.tile.y)) {
            final_tiles.push(entry.tile);
            if final_tiles.len() >= opts.max_tiles {
                break;
            }
        }
    }

    final_tiles
}

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

    #[test]
    fn visible_tiles_flat_view_preserves_coverage() {
        let bounds = WorldBounds::new(
            WebMercator::project(&GeoCoord::from_lat_lon(37.7749, -122.4194)),
            WebMercator::project(&GeoCoord::from_lat_lon(37.8049, -122.3894)),
        );
        let zoom = 12;
        let view = FlatTileView::new(
            bounds.center(),
            1000.0,
            std::f64::consts::FRAC_PI_4,
            0.0,
            std::f64::consts::FRAC_PI_4,
            800,
            600,
        );

        let tiles = visible_tiles_flat_view(&bounds, zoom, &view);

        // Basic checks
        assert!(!tiles.is_empty());
        assert!(tiles.iter().all(|t| t.zoom == zoom));
    }

    #[test]
    fn visible_tiles_flat_view_capped_does_not_exceed_limit() {
        let bounds = WorldBounds::new(
            WebMercator::project(&GeoCoord::from_lat_lon(37.7749, -122.4194)),
            WebMercator::project(&GeoCoord::from_lat_lon(37.8049, -122.3894)),
        );
        let zoom = 12;
        let view = FlatTileView::new(
            bounds.center(),
            1000.0,
            std::f64::consts::FRAC_PI_4,
            0.0,
            std::f64::consts::FRAC_PI_4,
            800,
            600,
        );

        let max_tiles = 10;
        let tiles = visible_tiles_flat_view_capped(&bounds, zoom, &view, max_tiles);

        // Check that the number of tiles does not exceed the cap
        assert!(tiles.len() <= max_tiles);
    }

    #[test]
    fn visible_tiles_covering_top_down_returns_uniform_zoom() {
        use crate::frustum::Frustum;
        use glam::DMat4;

        // Top-down camera (no pitch) should give all tiles at the same zoom.
        let proj = DMat4::perspective_rh(std::f64::consts::FRAC_PI_4, 1.5, 1.0, 500_000.0);
        let eye = glam::DVec3::new(0.0, 0.0, 50_000.0);
        let target = glam::DVec3::ZERO;
        let view = DMat4::look_at_rh(eye, target, glam::DVec3::Y);
        let frustum = Frustum::from_view_projection(&(proj * view));

        let world_size = crate::mercator::WebMercator::world_size();
        let cam = CoveringCamera {
            camera_x: 0.5,
            camera_y: 0.5,
            camera_to_center_z: 50_000.0 / world_size,
            center_x: 0.5,
            center_y: 0.5,
            pitch_rad: 0.0,
            fov_deg: 45.0,
            zoom: 4.0,
            display_tile_size: 256,
        };
        let opts = CoveringTilesOptions {
            min_zoom: 0,
            max_zoom: MAX_ZOOM,
            round_zoom: false,
            tile_size: 256,
            max_tiles: 512,
            allow_variable_zoom: true,
            render_world_copies: false,
        };

        let tiles = visible_tiles_covering(&frustum, &cam, &opts);
        assert!(!tiles.is_empty());
        let first_zoom = tiles[0].zoom;
        // All tiles at same zoom for a top-down view
        assert!(tiles.iter().all(|t| t.zoom == first_zoom));
    }

    #[test]
    fn visible_tiles_covering_pitched_produces_variable_zoom() {
        use crate::frustum::Frustum;
        use glam::DMat4;

        // Steeply pitched camera should produce tiles at different zoom levels.
        let proj = DMat4::perspective_rh(std::f64::consts::FRAC_PI_4, 1.5, 1.0, 500_000.0);
        let eye = glam::DVec3::new(0.0, -40_000.0, 20_000.0);
        let target = glam::DVec3::ZERO;
        let view = DMat4::look_at_rh(eye, target, glam::DVec3::Z);
        let frustum = Frustum::from_view_projection(&(proj * view));

        let world_size = crate::mercator::WebMercator::world_size();
        let pitch_rad = 1.1; // ~63 degrees
        let cam = CoveringCamera {
            camera_x: 0.5,
            camera_y: 0.5 - 40_000.0 / world_size,
            camera_to_center_z: 20_000.0 / world_size,
            center_x: 0.5,
            center_y: 0.5,
            pitch_rad,
            fov_deg: 45.0,
            zoom: 8.0,
            display_tile_size: 256,
        };
        let opts = CoveringTilesOptions {
            min_zoom: 0,
            max_zoom: MAX_ZOOM,
            round_zoom: false,
            tile_size: 256,
            max_tiles: 512,
            allow_variable_zoom: true,
            render_world_copies: false,
        };

        let tiles = visible_tiles_covering(&frustum, &cam, &opts);
        assert!(!tiles.is_empty());
        let zooms: std::collections::HashSet<u8> = tiles.iter().map(|t| t.zoom).collect();
        // At steep pitch with variable zoom enabled, we expect multiple zoom levels
        assert!(
            zooms.len() > 1,
            "Expected multiple zoom levels at steep pitch, got: {:?}",
            zooms
        );
    }

    #[test]
    fn visible_tiles_covering_respects_max_tiles() {
        use crate::frustum::Frustum;
        use glam::DMat4;

        let vp = DMat4::orthographic_rh(-500_000.0, 500_000.0, -500_000.0, 500_000.0, -500_000.0, 500_000.0);
        let frustum = Frustum::from_view_projection(&vp);

        let cam = CoveringCamera {
            camera_x: 0.5,
            camera_y: 0.5,
            camera_to_center_z: 0.01,
            center_x: 0.5,
            center_y: 0.5,
            pitch_rad: 0.0,
            fov_deg: 45.0,
            zoom: 10.0,
            display_tile_size: 256,
        };
        let opts = CoveringTilesOptions {
            min_zoom: 0,
            max_zoom: MAX_ZOOM,
            round_zoom: false,
            tile_size: 256,
            max_tiles: 5,
            allow_variable_zoom: false,
            render_world_copies: false,
        };

        let tiles = visible_tiles_covering(&frustum, &cam, &opts);
        assert!(tiles.len() <= 5);
    }

    #[test]
    fn visible_tiles_covering_no_duplicates() {
        use crate::frustum::Frustum;
        use glam::DMat4;

        let proj = DMat4::perspective_rh(std::f64::consts::FRAC_PI_4, 1.5, 1.0, 500_000.0);
        let eye = glam::DVec3::new(0.0, 0.0, 50_000.0);
        let view = DMat4::look_at_rh(eye, glam::DVec3::ZERO, glam::DVec3::Y);
        let frustum = Frustum::from_view_projection(&(proj * view));

        let cam = CoveringCamera {
            camera_x: 0.5,
            camera_y: 0.5,
            camera_to_center_z: 0.001,
            center_x: 0.5,
            center_y: 0.5,
            pitch_rad: 0.0,
            fov_deg: 45.0,
            zoom: 4.0,
            display_tile_size: 256,
        };
        let opts = CoveringTilesOptions {
            render_world_copies: true,
            ..CoveringTilesOptions::default()
        };

        let tiles = visible_tiles_covering(&frustum, &cam, &opts);
        let unique: std::collections::HashSet<_> = tiles.iter().collect();
        assert_eq!(tiles.len(), unique.len(), "Covering tiles should have no duplicates");
    }

    #[test]
    fn high_altitude_pitched_flat_view_selects_enough_tiles() {
        // Reproduce the fly_to zoom-out scenario: zoom 7, distance ~1.4 M m,
        // pitch ~48 degrees.  Before the fix the footprint was clamped to
        // 500 km which left most of the viewport uncovered.
        let center = WebMercator::project(&GeoCoord::from_lat_lon(40.839, -72.9));
        let distance = 1_406_424.0; // ~zoom 7 altitude
        let pitch = 48.4_f64.to_radians();
        let yaw = -20.3_f64.to_radians();
        let fov_y = 45.0_f64.to_radians();
        let view = FlatTileView::new(center, distance, pitch, yaw, fov_y, 1280, 720);

        // Build viewport bounds that match the camera altitude.
        let half = distance * 4.0;
        let bounds = WorldBounds::new(
            WorldCoord::new(center.position.x - half, center.position.y - half, 0.0),
            WorldCoord::new(center.position.x + half, center.position.y + half, 0.0),
        );

        let tiles = visible_tiles_flat_view_with_config(
            &bounds,
            7,
            &view,
            &FlatTileSelectionConfig::default(),
        );

        // At zoom 7 over the Atlantic/US the viewport should contain many
        // tiles -- certainly more than the 5 that were produced when the
        // footprint was clamped to 500 km.
        assert!(
            tiles.len() >= 10,
            "Expected >= 10 tiles at zoom 7 high altitude, got {}",
            tiles.len()
        );
    }

    #[test]
    fn visible_tiles_covering_variable_zoom_disabled_gives_uniform_zoom() {
        use crate::frustum::Frustum;
        use glam::DMat4;

        let proj = DMat4::perspective_rh(std::f64::consts::FRAC_PI_4, 1.5, 1.0, 500_000.0);
        let eye = glam::DVec3::new(0.0, -40_000.0, 20_000.0);
        let target = glam::DVec3::ZERO;
        let view = DMat4::look_at_rh(eye, target, glam::DVec3::Z);
        let frustum = Frustum::from_view_projection(&(proj * view));

        let world_size = crate::mercator::WebMercator::world_size();
        let cam = CoveringCamera {
            camera_x: 0.5,
            camera_y: 0.5 - 40_000.0 / world_size,
            camera_to_center_z: 20_000.0 / world_size,
            center_x: 0.5,
            center_y: 0.5,
            pitch_rad: 1.1,
            fov_deg: 45.0,
            zoom: 6.0,
            display_tile_size: 256,
        };
        let opts = CoveringTilesOptions {
            allow_variable_zoom: false,
            render_world_copies: false,
            ..CoveringTilesOptions::default()
        };

        let tiles = visible_tiles_covering(&frustum, &cam, &opts);
        assert!(!tiles.is_empty());
        let first_zoom = tiles[0].zoom;
        assert!(
            tiles.iter().all(|t| t.zoom == first_zoom),
            "With variable zoom disabled, all tiles should be at the same zoom"
        );
    }
}