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(
150                (self.center.lat - angle).max(-90.0),
151                -180.0,
152                0.0,
153            ),
154            GeoCoord::new(
155                (self.center.lat + angle).min(90.0),
156                180.0,
157                0.0,
158            ),
159        )
160    }
161}
162
163#[cfg(test)]
164mod tests {
165    use super::*;
166
167    #[test]
168    fn center_projects_to_origin() {
169        let proj = VerticalPerspective::new(GeoCoord::from_lat_lon(0.0, 0.0), 5_000_000.0);
170        let world = proj.project(&GeoCoord::from_lat_lon(0.0, 0.0));
171        assert!(world.position.x.abs() < 1e-8);
172        assert!(world.position.y.abs() < 1e-8);
173    }
174
175    #[test]
176    fn roundtrip_visible_point() {
177        let proj = VerticalPerspective::new(GeoCoord::from_lat_lon(0.0, 0.0), 8_000_000.0);
178        let geo = GeoCoord::new(10.0, 15.0, 250.0);
179        let world = proj.project_checked(&geo).expect("visible point");
180        let back = proj.unproject(&world);
181        assert!((back.lat - geo.lat).abs() < 1e-6);
182        assert!((back.lon - geo.lon).abs() < 1e-6);
183        assert!((back.alt - geo.alt).abs() < 1e-10);
184    }
185
186    #[test]
187    fn checked_projection_rejects_far_side_point() {
188        let proj = VerticalPerspective::new(GeoCoord::from_lat_lon(0.0, 0.0), 1_000_000.0);
189        assert!(proj.project_checked(&GeoCoord::from_lat_lon(0.0, 120.0)).is_none());
190    }
191
192    #[test]
193    fn projection_returns_nan_for_invisible_point() {
194        let proj = VerticalPerspective::new(GeoCoord::from_lat_lon(0.0, 0.0), 1_000_000.0);
195        let world = proj.project(&GeoCoord::from_lat_lon(0.0, 120.0));
196        assert!(world.position.x.is_nan());
197        assert!(world.position.y.is_nan());
198    }
199
200    #[test]
201    fn horizon_angle_grows_with_camera_height() {
202        let low = VerticalPerspective::new(GeoCoord::from_lat_lon(0.0, 0.0), 1_000_000.0);
203        let high = VerticalPerspective::new(GeoCoord::from_lat_lon(0.0, 0.0), 20_000_000.0);
204        assert!(high.horizon_central_angle() > low.horizon_central_angle());
205    }
206}
rustial_engine/math/vertical_perspective.rs

rustial_engine/math/
vertical_perspective.rs
