rustial-engine 0.0.1

//! Near-sided vertical perspective projection.
//!
//! This projection places the viewer above a configurable center point and
//! projects the visible hemisphere onto the tangent plane at that center.
//! It is a useful math primitive for future MapLibre-style vertical
//! perspective support.

use crate::bounds::GeoBounds;
use crate::coord::{GeoCoord, WorldCoord};
use crate::ellipsoid::Ellipsoid;
use crate::projection::Projection;
use glam::DVec3;

/// Near-sided vertical perspective projection onto a tangent plane.
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct VerticalPerspective {
    center: GeoCoord,
    camera_height: f64,
}

impl VerticalPerspective {
    /// Create a vertical perspective projection centered at `center`.
    ///
    /// `camera_height` is the viewer height above the ellipsoid surface, in
    /// meters, and must be positive.
    pub fn new(center: GeoCoord, camera_height: f64) -> Self {
        debug_assert!(camera_height.is_finite() && camera_height > 0.0);
        Self {
            center: GeoCoord::from_lat_lon(center.lat, center.lon),
            camera_height: if camera_height.is_finite() {
                camera_height.max(1.0)
            } else {
                1.0
            },
        }
    }

    /// Projection center on the globe.
    #[inline]
    pub fn center(&self) -> GeoCoord {
        self.center
    }

    /// Viewer height above the surface in meters.
    #[inline]
    pub fn camera_height(&self) -> f64 {
        self.camera_height
    }

    /// Earth radius used by the projection.
    #[inline]
    pub fn radius(&self) -> f64 {
        Ellipsoid::WGS84.a
    }

    /// Maximum visible central angle from the projection center, in radians.
    #[inline]
    pub fn horizon_central_angle(&self) -> f64 {
        (self.radius() / (self.radius() + self.camera_height)).acos()
    }

    /// Checked projection that rejects points beyond the visible horizon.
    pub fn project_checked(&self, geo: &GeoCoord) -> Option<WorldCoord> {
        let r = self.radius();
        let h = self.camera_height;
        let (east, north, up) = self.local_basis();
        let p = self.surface_unit(geo);
        let cos_c = up.dot(p);
        let visible_limit = r / (r + h);
        if cos_c < visible_limit {
            return None;
        }

        let scale = h / (r + h - r * cos_c);
        let x = scale * r * east.dot(p);
        let y = scale * r * north.dot(p);
        Some(WorldCoord::new(x, y, geo.alt))
    }

    fn local_basis(&self) -> (DVec3, DVec3, DVec3) {
        let lat0 = self.center.lat.to_radians();
        let lon0 = self.center.lon.to_radians();
        let (sin_lat0, cos_lat0) = lat0.sin_cos();
        let (sin_lon0, cos_lon0) = lon0.sin_cos();

        let east = DVec3::new(-sin_lon0, cos_lon0, 0.0);
        let north = DVec3::new(-sin_lat0 * cos_lon0, -sin_lat0 * sin_lon0, cos_lat0);
        let up = DVec3::new(cos_lat0 * cos_lon0, cos_lat0 * sin_lon0, sin_lat0);
        (east, north, up)
    }

    fn surface_unit(&self, geo: &GeoCoord) -> DVec3 {
        let lat = geo.lat.to_radians();
        let lon = geo.lon.to_radians();
        let (sin_lat, cos_lat) = lat.sin_cos();
        let (sin_lon, cos_lon) = lon.sin_cos();
        DVec3::new(cos_lat * cos_lon, cos_lat * sin_lon, sin_lat)
    }
}

impl Projection for VerticalPerspective {
    fn project(&self, geo: &GeoCoord) -> WorldCoord {
        self.project_checked(geo)
            .unwrap_or_else(|| WorldCoord::new(f64::NAN, f64::NAN, geo.alt))
    }

    fn unproject(&self, world: &WorldCoord) -> GeoCoord {
        let r = self.radius();
        let h = self.camera_height;
        let (east, north, up) = self.local_basis();
        let x = world.position.x;
        let y = world.position.y;
        let alt = world.position.z;

        let a = x * x + y * y + h * h;
        let b = -2.0 * h * (r + h);
        let c = h * (2.0 * r + h);
        let disc = b * b - 4.0 * a * c;
        if disc < 0.0 {
            return GeoCoord::new(f64::NAN, f64::NAN, alt);
        }

        let s = (-b - disc.sqrt()) / (2.0 * a);
        let local = DVec3::new(s * x, s * y, r + h - s * h);
        let ecef = east * local.x + north * local.y + up * local.z;
        let unit = ecef.normalize_or_zero();
        if unit.length_squared() <= 0.0 {
            return GeoCoord::new(f64::NAN, f64::NAN, alt);
        }

        let lat = unit.z.asin().to_degrees();
        let lon = unit.y.atan2(unit.x).to_degrees();
        GeoCoord::new(lat, lon, alt)
    }

    fn scale_factor(&self, geo: &GeoCoord) -> f64 {
        let r = self.radius();
        let h = self.camera_height;
        let p = self.surface_unit(geo);
        let (_, _, up) = self.local_basis();
        let cos_c = up.dot(p);
        let denom = r + h - r * cos_c;
        ((h * (r + h)) / (denom * denom)).max(0.0)
    }

    fn projection_bounds(&self) -> GeoBounds {
        let angle = self.horizon_central_angle().to_degrees();
        GeoBounds::new(
            GeoCoord::new((self.center.lat - angle).max(-90.0), -180.0, 0.0),
            GeoCoord::new((self.center.lat + angle).min(90.0), 180.0, 0.0),
        )
    }
}

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

    #[test]
    fn center_projects_to_origin() {
        let proj = VerticalPerspective::new(GeoCoord::from_lat_lon(0.0, 0.0), 5_000_000.0);
        let world = proj.project(&GeoCoord::from_lat_lon(0.0, 0.0));
        assert!(world.position.x.abs() < 1e-8);
        assert!(world.position.y.abs() < 1e-8);
    }

    #[test]
    fn roundtrip_visible_point() {
        let proj = VerticalPerspective::new(GeoCoord::from_lat_lon(0.0, 0.0), 8_000_000.0);
        let geo = GeoCoord::new(10.0, 15.0, 250.0);
        let world = proj.project_checked(&geo).expect("visible point");
        let back = proj.unproject(&world);
        assert!((back.lat - geo.lat).abs() < 1e-6);
        assert!((back.lon - geo.lon).abs() < 1e-6);
        assert!((back.alt - geo.alt).abs() < 1e-10);
    }

    #[test]
    fn checked_projection_rejects_far_side_point() {
        let proj = VerticalPerspective::new(GeoCoord::from_lat_lon(0.0, 0.0), 1_000_000.0);
        assert!(proj
            .project_checked(&GeoCoord::from_lat_lon(0.0, 120.0))
            .is_none());
    }

    #[test]
    fn projection_returns_nan_for_invisible_point() {
        let proj = VerticalPerspective::new(GeoCoord::from_lat_lon(0.0, 0.0), 1_000_000.0);
        let world = proj.project(&GeoCoord::from_lat_lon(0.0, 120.0));
        assert!(world.position.x.is_nan());
        assert!(world.position.y.is_nan());
    }

    #[test]
    fn horizon_angle_grows_with_camera_height() {
        let low = VerticalPerspective::new(GeoCoord::from_lat_lon(0.0, 0.0), 1_000_000.0);
        let high = VerticalPerspective::new(GeoCoord::from_lat_lon(0.0, 0.0), 20_000_000.0);
        assert!(high.horizon_central_angle() > low.horizon_central_angle());
    }
}