rustial_engine/math/
vertical_perspective.rs

1//! Near-sided vertical perspective projection.
2//!
3//! This projection places the viewer above a configurable center point and
4//! projects the visible hemisphere onto the tangent plane at that center.
5//! It is a useful math primitive for future MapLibre-style vertical
6//! perspective support.
7
8use crate::bounds::GeoBounds;
9use crate::coord::{GeoCoord, WorldCoord};
10use crate::ellipsoid::Ellipsoid;
11use crate::projection::Projection;
12use glam::DVec3;
13
14/// Near-sided vertical perspective projection onto a tangent plane.
15#[derive(Debug, Clone, Copy, PartialEq)]
16pub struct VerticalPerspective {
17    center: GeoCoord,
18    camera_height: f64,
19}
20
21impl VerticalPerspective {
22    /// Create a vertical perspective projection centered at `center`.
23    ///
24    /// `camera_height` is the viewer height above the ellipsoid surface, in
25    /// meters, and must be positive.
26    pub fn new(center: GeoCoord, camera_height: f64) -> Self {
27        debug_assert!(camera_height.is_finite() && camera_height > 0.0);
28        Self {
29            center: GeoCoord::from_lat_lon(center.lat, center.lon),
30            camera_height: if camera_height.is_finite() {
31                camera_height.max(1.0)
32            } else {
33                1.0
34            },
35        }
36    }
37
38    /// Projection center on the globe.
39    #[inline]
40    pub fn center(&self) -> GeoCoord {
41        self.center
42    }
43
44    /// Viewer height above the surface in meters.
45    #[inline]
46    pub fn camera_height(&self) -> f64 {
47        self.camera_height
48    }
49
50    /// Earth radius used by the projection.
51    #[inline]
52    pub fn radius(&self) -> f64 {
53        Ellipsoid::WGS84.a
54    }
55
56    /// Maximum visible central angle from the projection center, in radians.
57    #[inline]
58    pub fn horizon_central_angle(&self) -> f64 {
59        (self.radius() / (self.radius() + self.camera_height)).acos()
60    }
61
62    /// Checked projection that rejects points beyond the visible horizon.
63    pub fn project_checked(&self, geo: &GeoCoord) -> Option<WorldCoord> {
64        let r = self.radius();
65        let h = self.camera_height;
66        let (east, north, up) = self.local_basis();
67        let p = self.surface_unit(geo);
68        let cos_c = up.dot(p);
69        let visible_limit = r / (r + h);
70        if cos_c < visible_limit {
71            return None;
72        }
73
74        let scale = h / (r + h - r * cos_c);
75        let x = scale * r * east.dot(p);
76        let y = scale * r * north.dot(p);
77        Some(WorldCoord::new(x, y, geo.alt))
78    }
79
80    fn local_basis(&self) -> (DVec3, DVec3, DVec3) {
81        let lat0 = self.center.lat.to_radians();
82        let lon0 = self.center.lon.to_radians();
83        let (sin_lat0, cos_lat0) = lat0.sin_cos();
84        let (sin_lon0, cos_lon0) = lon0.sin_cos();
85
86        let east = DVec3::new(-sin_lon0, cos_lon0, 0.0);
87        let north = DVec3::new(-sin_lat0 * cos_lon0, -sin_lat0 * sin_lon0, cos_lat0);
88        let up = DVec3::new(cos_lat0 * cos_lon0, cos_lat0 * sin_lon0, sin_lat0);
89        (east, north, up)
90    }
91
92    fn surface_unit(&self, geo: &GeoCoord) -> DVec3 {
93        let lat = geo.lat.to_radians();
94        let lon = geo.lon.to_radians();
95        let (sin_lat, cos_lat) = lat.sin_cos();
96        let (sin_lon, cos_lon) = lon.sin_cos();
97        DVec3::new(cos_lat * cos_lon, cos_lat * sin_lon, sin_lat)
98    }
99}
100
101impl Projection for VerticalPerspective {
102    fn project(&self, geo: &GeoCoord) -> WorldCoord {
103        self.project_checked(geo)
104            .unwrap_or_else(|| WorldCoord::new(f64::NAN, f64::NAN, geo.alt))
105    }
106
107    fn unproject(&self, world: &WorldCoord) -> GeoCoord {
108        let r = self.radius();
109        let h = self.camera_height;
110        let (east, north, up) = self.local_basis();
111        let x = world.position.x;
112        let y = world.position.y;
113        let alt = world.position.z;
114
115        let a = x * x + y * y + h * h;
116        let b = -2.0 * h * (r + h);
117        let c = h * (2.0 * r + h);
118        let disc = b * b - 4.0 * a * c;
119        if disc < 0.0 {
120            return GeoCoord::new(f64::NAN, f64::NAN, alt);
121        }
122
123        let s = (-b - disc.sqrt()) / (2.0 * a);
124        let local = DVec3::new(s * x, s * y, r + h - s * h);
125        let ecef = east * local.x + north * local.y + up * local.z;
126        let unit = ecef.normalize_or_zero();
127        if unit.length_squared() <= 0.0 {
128            return GeoCoord::new(f64::NAN, f64::NAN, alt);
129        }
130
131        let lat = unit.z.asin().to_degrees();
132        let lon = unit.y.atan2(unit.x).to_degrees();
133        GeoCoord::new(lat, lon, alt)
134    }
135
136    fn scale_factor(&self, geo: &GeoCoord) -> f64 {
137        let r = self.radius();
138        let h = self.camera_height;
139        let p = self.surface_unit(geo);
140        let (_, _, up) = self.local_basis();
141        let cos_c = up.dot(p);
142        let denom = r + h - r * cos_c;
143        ((h * (r + h)) / (denom * denom)).max(0.0)
144    }
145
146    fn projection_bounds(&self) -> GeoBounds {
147        let angle = self.horizon_central_angle().to_degrees();
148        GeoBounds::new(
149            GeoCoord::new((self.center.lat - angle).max(-90.0), -180.0, 0.0),
150            GeoCoord::new((self.center.lat + angle).min(90.0), 180.0, 0.0),
151        )
152    }
153}
154
155#[cfg(test)]
156mod tests {
157    use super::*;
158
159    #[test]
160    fn center_projects_to_origin() {
161        let proj = VerticalPerspective::new(GeoCoord::from_lat_lon(0.0, 0.0), 5_000_000.0);
162        let world = proj.project(&GeoCoord::from_lat_lon(0.0, 0.0));
163        assert!(world.position.x.abs() < 1e-8);
164        assert!(world.position.y.abs() < 1e-8);
165    }
166
167    #[test]
168    fn roundtrip_visible_point() {
169        let proj = VerticalPerspective::new(GeoCoord::from_lat_lon(0.0, 0.0), 8_000_000.0);
170        let geo = GeoCoord::new(10.0, 15.0, 250.0);
171        let world = proj.project_checked(&geo).expect("visible point");
172        let back = proj.unproject(&world);
173        assert!((back.lat - geo.lat).abs() < 1e-6);
174        assert!((back.lon - geo.lon).abs() < 1e-6);
175        assert!((back.alt - geo.alt).abs() < 1e-10);
176    }
177
178    #[test]
179    fn checked_projection_rejects_far_side_point() {
180        let proj = VerticalPerspective::new(GeoCoord::from_lat_lon(0.0, 0.0), 1_000_000.0);
181        assert!(proj
182            .project_checked(&GeoCoord::from_lat_lon(0.0, 120.0))
183            .is_none());
184    }
185
186    #[test]
187    fn projection_returns_nan_for_invisible_point() {
188        let proj = VerticalPerspective::new(GeoCoord::from_lat_lon(0.0, 0.0), 1_000_000.0);
189        let world = proj.project(&GeoCoord::from_lat_lon(0.0, 120.0));
190        assert!(world.position.x.is_nan());
191        assert!(world.position.y.is_nan());
192    }
193
194    #[test]
195    fn horizon_angle_grows_with_camera_height() {
196        let low = VerticalPerspective::new(GeoCoord::from_lat_lon(0.0, 0.0), 1_000_000.0);
197        let high = VerticalPerspective::new(GeoCoord::from_lat_lon(0.0, 0.0), 20_000_000.0);
198        assert!(high.horizon_central_angle() > low.horizon_central_angle());
199    }
200}
rustial_engine/math/vertical_perspective.rs

rustial_engine/math/
vertical_perspective.rs
