summaryrefslogtreecommitdiff
path: root/crates/gipfel-core/src/panorama.rs
diff options
context:
space:
mode:
Diffstat (limited to 'crates/gipfel-core/src/panorama.rs')
-rw-r--r--crates/gipfel-core/src/panorama.rs502
1 files changed, 502 insertions, 0 deletions
diff --git a/crates/gipfel-core/src/panorama.rs b/crates/gipfel-core/src/panorama.rs
new file mode 100644
index 0000000..9e27251
--- /dev/null
+++ b/crates/gipfel-core/src/panorama.rs
@@ -0,0 +1,502 @@
+//! Panorama: viewpoint, hill visibility (occlusion, refraction) and the
+//! link between geographic coordinates and image coordinates.
+//!
+//! Faithful port of gipfel's `Panorama.cxx`.
+
+use crate::hill::{Hill, DEG2RAD, PI};
+use crate::projection::{projection_from_type, Projection, ProjectionType, ViewParams};
+
+const EARTH_RADIUS: f64 = 6371000.785;
+
+pub struct Panorama {
+ pub hills: Vec<Hill>,
+ /// indices into `hills`, sorted by alpha descending (rebuilt by update_angles)
+ pub alpha_order: Vec<usize>,
+ /// indices into `hills`, sorted by alpha descending
+ pub close: Vec<usize>,
+ /// indices into `hills`
+ pub visible: Vec<usize>,
+
+ pub view_name: Option<String>,
+ /// latitude, radians
+ pub view_phi: f64,
+ /// longitude, radians
+ pub view_lam: f64,
+ /// meters
+ pub view_height: f64,
+
+ pub height_dist_ratio: f64,
+ pub hide_value: f64,
+ /// Additional visibility limit: maximum angular distance (radians) from
+ /// the image center. The app sets this to the angle the projection maps
+ /// to the image edge (the lens' actual field of view, including fitted
+ /// distortion), so that directions outside the photo — which a distortion
+ /// model can fold back into the image — neither get markers nor compass
+ /// ticks. 0 = no additional limit (only the projection's own limit).
+ pub max_view_angle: f64,
+
+ pub parms: ViewParams,
+ pub projection_type: ProjectionType,
+ proj: Box<dyn Projection>,
+}
+
+impl Default for Panorama {
+ fn default() -> Self {
+ Self::new()
+ }
+}
+
+impl Panorama {
+ pub fn new() -> Self {
+ Panorama {
+ hills: Vec::new(),
+ alpha_order: Vec::new(),
+ close: Vec::new(),
+ visible: Vec::new(),
+ view_name: None,
+ view_phi: 0.0,
+ view_lam: 0.0,
+ view_height: 0.0,
+ height_dist_ratio: 0.07,
+ hide_value: 1.2,
+ max_view_angle: 0.0,
+ parms: ViewParams::new(),
+ projection_type: ProjectionType::Rectilinear,
+ proj: projection_from_type(ProjectionType::Rectilinear),
+ }
+ }
+
+ pub fn load_data(&mut self, path: &str) -> std::io::Result<()> {
+ let mut hills = crate::hill::load(path)?;
+ crate::hill::mark_duplicates(&mut hills, 0.00001);
+ self.hills = hills;
+ self.update_angles();
+ Ok(())
+ }
+
+ pub fn add_hills(&mut self, extra: &[Hill]) {
+ self.hills.extend_from_slice(extra);
+ crate::hill::mark_duplicates(&mut self.hills, 0.00001);
+ self.update_angles();
+ }
+
+ pub fn remove_hills(&mut self, flags: u32) {
+ self.hills.retain(|m| m.flags & flags == 0);
+ // indices shifted; rebuild all derived lists
+ self.alpha_order.clear();
+ self.close.clear();
+ self.visible.clear();
+ }
+
+ pub fn set_viewpoint(&mut self, idx: Option<usize>) {
+ let Some(idx) = idx else { return };
+ let m = &self.hills[idx];
+ self.view_phi = m.phi;
+ self.view_lam = m.lam;
+ self.view_height = m.height;
+ self.view_name = Some(m.name.clone());
+ self.update_angles();
+ }
+
+ /// Fit view parameters to the known hills (indices into `self.hills`).
+ /// `fit_distortion` enables the 7-parameter fit (k0/k1/x0) for 4+ hills;
+ /// otherwise the distortion params stay locked at their current values.
+ pub fn comp_params(&mut self, known: &[usize], fit_distortion: bool) -> Result<(), String> {
+ let ret = crate::lsq::comp_params(&*self.proj, &self.hills, known, &mut self.parms, fit_distortion);
+ if ret.is_ok() {
+ self.update_visible_mountains(Some(known));
+ }
+ ret
+ }
+
+ /// Like [`Panorama::comp_params`] but the observed hills are an external
+ /// slice (used by the golden fit tests, mirroring the C harness).
+ pub fn comp_params_with(&mut self, hills: &[Hill], known: &[usize], fit_distortion: bool) -> Result<(), String> {
+ let ret = crate::lsq::comp_params(&*self.proj, hills, known, &mut self.parms, fit_distortion);
+ if ret.is_ok() {
+ self.update_visible_mountains(None);
+ }
+ ret
+ }
+
+ pub fn set_center_angle(&mut self, deg: f64) {
+ self.parms.a_center = deg * DEG2RAD;
+ self.update_visible_mountains(None);
+ }
+
+ pub fn set_pitch_angle(&mut self, deg: f64) {
+ self.parms.a_pitch = deg * DEG2RAD;
+ self.update_coordinates(None);
+ }
+
+ pub fn set_roll_angle(&mut self, deg: f64) {
+ self.parms.a_roll = deg * DEG2RAD;
+ self.update_coordinates(None);
+ }
+
+ pub fn set_scale(&mut self, s: f64) {
+ self.parms.scale = s;
+ self.update_coordinates(None);
+ }
+
+ pub fn set_distortion_params(&mut self, k0: f64, k1: f64, x0: f64) {
+ self.parms.k0 = k0;
+ self.parms.k1 = k1;
+ self.parms.x0 = x0;
+ self.update_coordinates(None);
+ }
+
+ pub fn set_height_dist_ratio(&mut self, r: f64) {
+ self.height_dist_ratio = r;
+ self.update_close_mountains();
+ }
+
+ pub fn set_view_long(&mut self, deg: f64) {
+ self.view_lam = deg * DEG2RAD;
+ self.update_angles();
+ }
+
+ pub fn set_view_lat(&mut self, deg: f64) {
+ self.view_phi = deg * DEG2RAD;
+ self.update_angles();
+ }
+
+ pub fn set_view_height(&mut self, v: f64) {
+ self.view_height = v;
+ self.update_angles();
+ }
+
+ pub fn set_projection(&mut self, t: ProjectionType) {
+ self.projection_type = t;
+ self.proj = projection_from_type(t);
+ self.update_angles();
+ }
+
+ pub fn get_center_angle(&self) -> f64 {
+ self.parms.a_center / DEG2RAD
+ }
+ pub fn get_pitch_angle(&self) -> f64 {
+ self.parms.a_pitch / DEG2RAD
+ }
+ pub fn get_roll_angle(&self) -> f64 {
+ self.parms.a_roll / DEG2RAD
+ }
+ pub fn get_scale(&self) -> f64 {
+ self.parms.scale
+ }
+ pub fn get_view_long(&self) -> f64 {
+ self.view_lam / DEG2RAD
+ }
+ pub fn get_view_lat(&self) -> f64 {
+ self.view_phi / DEG2RAD
+ }
+ pub fn get_view_height(&self) -> f64 {
+ self.view_height
+ }
+
+ /// Last hill whose name matches exactly (C keeps the last match).
+ pub fn get_pos(&self, name: &str) -> Option<usize> {
+ let mut ret = None;
+ for (i, m) in self.hills.iter().enumerate() {
+ if m.name == name {
+ ret = Some(i);
+ }
+ }
+ ret
+ }
+
+ pub fn update_angles(&mut self) {
+ for i in 0..self.hills.len() {
+ let (phi, lam) = (self.hills[i].phi, self.hills[i].lam);
+ let d = self.distance(phi, lam);
+ self.hills[i].dist = d;
+ if phi != self.view_phi || lam != self.view_lam {
+ let a = self.alpha(i);
+ self.hills[i].alph = a;
+ }
+ }
+
+ // keep `hills` in insertion order (indices stay stable) and maintain a
+ // separate alpha-descending order, mirroring the C pointer-array sort
+ let mut idx: Vec<usize> = (0..self.hills.len()).collect();
+ idx.sort_by(|&a, &b| {
+ self.hills[b]
+ .alph
+ .partial_cmp(&self.hills[a].alph)
+ .unwrap_or(std::cmp::Ordering::Equal)
+ });
+ self.alpha_order = idx;
+
+ self.update_close_mountains();
+ }
+
+ /// Mark hills that are hidden behind a nearer, taller neighbour.
+ /// Faithful port including the original's `fabs(m->alph - n->alph > pi/2)`
+ /// oddity (comparison inside fabs, i.e. skip when m.alph - n.alph > pi/2).
+ pub fn mark_hidden(&mut self, hills: &[usize]) {
+ for i in 0..hills.len() {
+ let m = hills[i];
+ self.hills[m].flags &= !crate::hill::FLAG_HIDDEN;
+
+ if self.hills[m].flags & crate::hill::FLAG_DUPLIC != 0 {
+ continue;
+ }
+
+ for j in 0..hills.len() {
+ let n = hills[j];
+ if self.hills[n].flags & crate::hill::FLAG_DUPLIC != 0
+ || self.hills[n].flags & crate::hill::FLAG_TRACK_POINT != 0
+ {
+ continue;
+ }
+ if n == m || self.hills[m].alph - self.hills[n].alph > PI / 2.0 {
+ continue;
+ }
+ if self.hills[m].dist < self.hills[n].dist
+ || self.hills[m].a_elev > self.hills[n].a_elev
+ {
+ continue;
+ }
+
+ let h = (self.hills[n].a_elev - self.hills[m].a_elev)
+ / (self.hills[m].alph - self.hills[n].alph).abs();
+ if h.is_infinite() || h > self.hide_value {
+ self.hills[m].flags |= crate::hill::FLAG_HIDDEN;
+ }
+ }
+ }
+ }
+
+ pub fn update_close_mountains(&mut self) {
+ self.close.clear();
+
+ let order = self.alpha_order.clone();
+ for &i in &order {
+ let m = &self.hills[i];
+ if m.flags & crate::hill::FLAG_TRACK_POINT != 0
+ || ((m.phi != self.view_phi || m.lam != self.view_lam)
+ && (m.height / (m.dist * EARTH_RADIUS) > self.height_dist_ratio))
+ {
+ let elev = self.elevation(i);
+ self.hills[i].a_elev = elev;
+ self.close.push(i);
+ }
+ }
+
+ let close = self.close.clone();
+ self.mark_hidden(&close);
+ self.update_visible_mountains(None);
+ }
+
+ pub fn update_visible_mountains(&mut self, excluded: Option<&[usize]>) {
+ self.visible.clear();
+
+ for k in 0..self.close.len() {
+ let i = self.close[k];
+ if self.is_visible(self.hills[i].alph) {
+ self.visible.push(i);
+ self.hills[i].flags |= crate::hill::FLAG_VISIBLE;
+ } else {
+ self.hills[i].flags &= !crate::hill::FLAG_VISIBLE;
+ }
+ }
+
+ let owned;
+ let excl: Option<&[usize]> = match excluded {
+ Some(e) => {
+ owned = e.to_vec();
+ Some(&owned)
+ }
+ None => None,
+ };
+ self.update_coordinates(excl);
+ }
+
+ pub fn update_coordinates(&mut self, excluded: Option<&[usize]>) {
+ for k in 0..self.visible.len() {
+ let i = self.visible[k];
+ if let Some(e) = excluded {
+ if e.contains(&i) {
+ continue;
+ }
+ }
+ let (alph, elev) = (self.hills[i].alph, self.hills[i].a_elev);
+ let (x, y) = self.proj.get_coordinates(alph, elev, &self.parms);
+ self.hills[i].x = x;
+ self.hills[i].y = y;
+ }
+ }
+
+ /// angular distance (fraction of earth radius) viewpoint -> (phi, lam)
+ pub fn distance(&self, phi: f64, lam: f64) -> f64 {
+ let d_lam = self.view_lam - lam;
+
+ let t1 = (phi.cos() * d_lam.sin()).powi(2);
+ let t2 = (self.view_phi.cos() * phi.sin()
+ - self.view_phi.sin() * phi.cos() * d_lam.cos())
+ .powi(2);
+
+ (t1 + t2).sqrt().atan2(
+ self.view_phi.sin() * phi.sin() + self.view_phi.cos() * phi.cos() * d_lam.cos(),
+ )
+ }
+
+ pub fn alpha(&self, idx: usize) -> f64 {
+ let m = &self.hills[idx];
+ let sin_alph = (m.lam - self.view_lam).sin() * m.phi.cos() / m.dist.sin();
+ let cos_alph =
+ (m.phi.sin() - self.view_phi.sin() * m.dist.cos()) / (self.view_phi.cos() * m.dist.sin());
+
+ (sin_alph.atan2(cos_alph) + 2.0 * PI) % (2.0 * PI)
+ }
+
+ /// approximation of refraction effect as described by Tom Chester at
+ /// <http://tchester.org/sgm/analysis/peaks/refraction.html>
+ pub fn refraction(&self, idx: usize) -> f64 {
+ let (alpha, t0) = (6.5, 10.0);
+
+ let a = 2.9e-4 * (-self.view_height / 10000.0).exp() / (1.0 + 2.9 * t0 / 760.0);
+ let b = 2.9 * alpha / (760.0 * (1.0 + 2.9 * t0 / 760.0));
+ let c = a * (b - 1.0 / 10.0);
+
+ c * self.get_real_distance(idx) / (2000.0 * (1.0 + a))
+ }
+
+ pub fn elevation(&self, idx: usize) -> f64 {
+ let m = &self.hills[idx];
+ let theta = self.refraction(idx);
+
+ let b = m.height + self.get_earth_radius(m.phi);
+ let c = self.view_height + self.get_earth_radius(self.view_phi);
+
+ ((m.dist.cos() * b - c) / (m.dist.sin() * b)).atan() - theta
+ }
+
+ /// local distance to center of WGS84 ellipsoid
+ pub fn get_earth_radius(&self, phi: f64) -> f64 {
+ let a = 6378137.000f64;
+ let b = 6356752.315f64;
+ let ata = phi.tan();
+
+ a * (ata.powi(2) + 1.0).powf(1.0 / 2.0)
+ * b.abs()
+ * (b.powi(2) + a.powi(2) * ata.powi(2)).powf(-1.0 / 2.0)
+ }
+
+ /// real (line-of-sight) distance in meters
+ pub fn get_real_distance(&self, idx: usize) -> f64 {
+ let m = &self.hills[idx];
+ let a = self.view_height + self.get_earth_radius(self.view_phi);
+ let b = m.height + self.get_earth_radius(m.phi);
+ let gam = m.dist;
+
+ (a.powi(2) + b.powi(2) - 2.0 * a * b * gam.cos()).sqrt()
+ }
+
+ pub fn is_visible(&self, alph: f64) -> bool {
+ let mut center_dist = (alph - self.parms.a_center).abs();
+ if center_dist > PI {
+ center_dist = 2.0 * PI - center_dist;
+ }
+ if center_dist >= self.proj.get_view_angle() {
+ return false;
+ }
+ self.max_view_angle <= 0.0 || center_dist < self.max_view_angle
+ }
+
+ /// Forward-project (azimuth, elevation) to image-centered pixel coords.
+ /// `None` when the direction is not visible or the projection is
+ /// undefined there (mirrors C `Panorama::get_coordinates` returning 1,
+ /// plus the NaN guard of C `ScanImage::get_pixel_nearest`).
+ pub fn project_point(&self, alph: f64, elev: f64) -> Option<(f64, f64)> {
+ if !self.is_visible(alph) {
+ return None;
+ }
+ let (x, y) = self.proj.get_coordinates(alph, elev, &self.parms);
+ if x.is_nan() || y.is_nan() {
+ return None;
+ }
+ Some((x, y))
+ }
+
+ /// Parameter-only snapshot (no hill data) for resampling: same viewpoint
+ /// and projection parameters, but without the UI's lens-FOV visibility
+ /// limit (the stitch output range defines the FOV itself).
+ pub fn params_snapshot(&self) -> Panorama {
+ let mut p = Panorama::new();
+ p.projection_type = self.projection_type;
+ p.proj = projection_from_type(self.projection_type);
+ p.parms = self.parms;
+ p.view_phi = self.view_phi;
+ p.view_lam = self.view_lam;
+ p.view_height = self.view_height;
+ p
+ }
+
+ pub fn get_coordinates(&self, alph: f64, elev: f64) -> Option<(f64, f64)> {
+ if self.is_visible(alph) {
+ Some(self.proj.get_coordinates(alph, elev, &self.parms))
+ } else {
+ None
+ }
+ }
+
+ /// Export rows for the visible hills, mirroring the C `export_hills`
+ /// output loop: `name\theight\tx\ty\tdistance` with image-center origin.
+ /// When `require_export` is set, only hills with `FLAG_EXPORT` are listed.
+ pub fn export_hill_rows(&self, img_w: i64, img_h: i64, require_export: bool) -> Vec<String> {
+ let mut rows = Vec::new();
+ for &i in &self.visible {
+ let m = &self.hills[i];
+ if m.flags & crate::hill::FLAG_DUPLIC != 0
+ || m.flags & crate::hill::FLAG_HIDDEN != 0
+ || (require_export && m.flags & crate::hill::FLAG_EXPORT == 0)
+ {
+ continue;
+ }
+
+ let x = m.x.round_ties_even() as i64 + img_w / 2;
+ let y = m.y.round_ties_even() as i64 + img_h / 2;
+ if x < 0 || x > img_w || y < 0 || y > img_h {
+ continue;
+ }
+
+ rows.push(format!(
+ "{}\t{}\t{}\t{}\t{}",
+ m.name,
+ m.height.round_ties_even() as i64,
+ x,
+ y,
+ self.get_real_distance(i).round_ties_even() as i64
+ ));
+ }
+ rows
+ }
+}
+
+#[cfg(test)]
+mod tests {
+ use super::*;
+
+ #[test]
+ fn max_view_angle_limits_visibility() {
+ let mut pan = Panorama::new();
+ pan.set_center_angle(0.0);
+ let a30 = 30f64.to_radians();
+ let a70 = 70f64.to_radians();
+
+ // rectilinear default limit is +-60 deg
+ assert!(pan.is_visible(a30));
+ assert!(!pan.is_visible(a70));
+
+ // a tighter lens FOV (say +-35 deg) hides the 40 deg direction that
+ // the projection's own +-60 deg limit would allow
+ pan.max_view_angle = 35f64.to_radians();
+ assert!(pan.is_visible(a30));
+ assert!(!pan.is_visible(40f64.to_radians()));
+
+ // a wider limit never overrides the projection's own bound
+ pan.max_view_angle = 80f64.to_radians();
+ assert!(!pan.is_visible(a70));
+ }
+}
generated by cgit v1.2.3 (git 2.46.0) at 2026-10-30 21:02:20 +0000