From 375d87a9b0ec1cf2fffd705adb870c86144d2400 Mon Sep 17 00:00:00 2001 From: Johannes Hofmann Date: Mon, 26 Oct 2026 20:05:20 +0100 Subject: Replace the C++ implementation with the Rust port The Rust port (formerly gipfel-rs) is now the official successor to the FLTK/GSL C++ code base. This squashes the port into a single commit and removes the C++ sources (they remain in the git history below). Layout: crates/gipfel-core headless core: hills, panorama geometry, projections (maxima-generated derivatives), pure-Rust LM fit, gipfel: JPEG comment + EXIF I/O, profiles, stitching crates/gipfel-cli command-line tools (export, info, stitch, profiles) crates/gipfel-app Cosmic DE (libcosmic) GUI, binary name `gipfel` codegen/ maxima scripts + regen.sh for the projection code tests/ golden files + C++ reference harnesses (historical) The original README is kept for now; it will be updated to describe the Rust build. gipfel.dat, AUTHORS, ChangeLog and NEWS are unchanged. --- crates/gipfel-app/Cargo.toml | 19 + crates/gipfel-app/data/org.gipfel.App.desktop | 13 + crates/gipfel-app/src/compass.rs | 611 +++++++ crates/gipfel-app/src/dump.rs | 285 ++++ crates/gipfel-app/src/main.rs | 1805 ++++++++++++++++++++ crates/gipfel-app/src/viewport.rs | 659 +++++++ crates/gipfel-cli/Cargo.toml | 13 + crates/gipfel-cli/src/main.rs | 371 ++++ crates/gipfel-core/Cargo.toml | 13 + crates/gipfel-core/examples/genjpg.rs | 24 + crates/gipfel-core/src/hill.rs | 224 +++ crates/gipfel-core/src/lib.rs | 23 + crates/gipfel-core/src/lsq.rs | 392 +++++ crates/gipfel-core/src/metadata.rs | 697 ++++++++ crates/gipfel-core/src/panorama.rs | 502 ++++++ crates/gipfel-core/src/profile.rs | 362 ++++ crates/gipfel-core/src/projection/cylindrical.rs | 86 + .../src/projection/generated/cylindrical_funcs.rs | 10 + .../src/projection/generated/rectilinear_funcs.rs | 16 + crates/gipfel-core/src/projection/mod.rs | 174 ++ crates/gipfel-core/src/projection/rectilinear.rs | 104 ++ crates/gipfel-core/src/stitch.rs | 323 ++++ crates/gipfel-core/src/track.rs | 248 +++ crates/gipfel-core/tests/golden.rs | 511 ++++++ 24 files changed, 7485 insertions(+) create mode 100644 crates/gipfel-app/Cargo.toml create mode 100644 crates/gipfel-app/data/org.gipfel.App.desktop create mode 100644 crates/gipfel-app/src/compass.rs create mode 100644 crates/gipfel-app/src/dump.rs create mode 100644 crates/gipfel-app/src/main.rs create mode 100644 crates/gipfel-app/src/viewport.rs create mode 100644 crates/gipfel-cli/Cargo.toml create mode 100644 crates/gipfel-cli/src/main.rs create mode 100644 crates/gipfel-core/Cargo.toml create mode 100644 crates/gipfel-core/examples/genjpg.rs create mode 100644 crates/gipfel-core/src/hill.rs create mode 100644 crates/gipfel-core/src/lib.rs create mode 100644 crates/gipfel-core/src/lsq.rs create mode 100644 crates/gipfel-core/src/metadata.rs create mode 100644 crates/gipfel-core/src/panorama.rs create mode 100644 crates/gipfel-core/src/profile.rs create mode 100644 crates/gipfel-core/src/projection/cylindrical.rs create mode 100644 crates/gipfel-core/src/projection/generated/cylindrical_funcs.rs create mode 100644 crates/gipfel-core/src/projection/generated/rectilinear_funcs.rs create mode 100644 crates/gipfel-core/src/projection/mod.rs create mode 100644 crates/gipfel-core/src/projection/rectilinear.rs create mode 100644 crates/gipfel-core/src/stitch.rs create mode 100644 crates/gipfel-core/src/track.rs create mode 100644 crates/gipfel-core/tests/golden.rs (limited to 'crates') diff --git a/crates/gipfel-app/Cargo.toml b/crates/gipfel-app/Cargo.toml new file mode 100644 index 0000000..22589c1 --- /dev/null +++ b/crates/gipfel-app/Cargo.toml @@ -0,0 +1,19 @@ +[package] +name = "gipfel-app" +version = "0.1.0" +edition = "2021" +description = "gipfel port to Cosmic DE" +license = "GPL-2.0-or-later" + +[[bin]] +name = "gipfel" +path = "src/main.rs" + +[dependencies] +gipfel-core = { path = "../gipfel-core" } +libcosmic = { git = "https://github.com/pop-os/libcosmic.git", rev = "430f3d48a97d8779f9c1985473912e594ccca1ae", default-features = false, features = ["winit", "wgpu", "tokio", "a11y", "xdg-portal"] } +image.workspace = true +fontdue = "0.9" +cosmic-text = "0.19" +fontdb = "0.23" +tokio = { version = "1", features = ["time"] } diff --git a/crates/gipfel-app/data/org.gipfel.App.desktop b/crates/gipfel-app/data/org.gipfel.App.desktop new file mode 100644 index 0000000..b20c468 --- /dev/null +++ b/crates/gipfel-app/data/org.gipfel.App.desktop @@ -0,0 +1,13 @@ +[Desktop Entry] +Type=Application +Name=gipfel +GenericName=Panorama annotator +Comment=Annotate mountain panoramas +Exec=gipfel %U +TryExec=gipfel +Icon=org.gipfel.App +Terminal=false +Categories=Graphics;Photography; +Keywords=mountain;panorama;photo;annotation; +MimeType=image/jpeg;image/png; +StartupNotify=true diff --git a/crates/gipfel-app/src/compass.rs b/crates/gipfel-app/src/compass.rs new file mode 100644 index 0000000..a17cd5c --- /dev/null +++ b/crates/gipfel-app/src/compass.rs @@ -0,0 +1,611 @@ +//! Experimental compass strip: shows N/E/S/W at their projected positions +//! relative to the view direction; drag horizontally to change the direction. +//! +//! Two experimental styles (switchable from the View menu): +//! - `Aligned`: the strip shares the image's horizontal transform, so each +//! direction appears exactly above its position in the image. Grab-drag +//! keeps the direction under the cursor fixed. +//! - `Ruler`: a fixed 0..360° ruler across the widget width; press or drag +//! sets the direction to the value under the cursor. The shaded band marks +//! the current horizontal field of view. +//! +//! Layout (top to bottom): direction labels, center marker, tick scale on a +//! baseline. The whole strip sits on a secondary-container background so it +//! is visually distinct from the image (pan) area. + +use std::cell::RefCell; +use std::rc::Rc; + +use cosmic::iced::advanced::graphics::geometry::Renderer as _; +use cosmic::iced::advanced::renderer::Style as RendererStyle; +use cosmic::iced::advanced::text::{self as adv_text, paragraph::Plain, LineHeight}; +use cosmic::iced::advanced::widget::{tree, Tree, Widget as WidgetTrait}; +use cosmic::iced::advanced::{Clipboard, Layout, Shell}; +use cosmic::iced::mouse::{self, Button, Cursor, Interaction}; +use cosmic::iced::widget::canvas::{Frame, Path, Stroke}; +use cosmic::iced::advanced::graphics::geometry::Style as GeoStyle; +use cosmic::iced::{alignment, Color, Event, Length, Pixels, Point, Rectangle, Size}; +use cosmic::Renderer as R; +use cosmic::Theme; + +use gipfel_core::projection::{ + projection_from_type, Projection, ProjArgs, ProjectionType, ViewParams, +}; + +pub const STRIP_H: f32 = 56.0; + +#[derive(Clone, Copy, Debug, PartialEq, Eq)] +pub enum CompassStyle { + Off, + Aligned, + Ruler, +} + +/// inverse: azimuth whose image x is `x_img`, measured with `c_ref` as +/// center angle. Analytic base inverse + damped Newton refinement so lens +/// distortion is honoured. +/// Damped, clamped Newton iteration for a scalar inverse projection: find +/// `a` such that `f(a) ≈ target`. `f` returns `None` where the forward +/// projection is undefined (the iteration stops there). The step is clamped +/// to ±0.3 rad and `a` kept in `[lo, hi]` to stay off the projection poles. +pub(crate) fn newton_inverse( + target: f64, + mut f: impl FnMut(f64) -> Option, + seed: f64, + lo: f64, + hi: f64, +) -> f64 { + let mut a = seed.clamp(lo, hi); + for _ in 0..6 { + let Some(y) = f(a) else { break }; + let resid = y - target; + let h = 1e-5; + let (Some(y1), Some(y2)) = (f(a + h), f(a - h)) else { + break; + }; + let d = (y1 - y2) / (2.0 * h); + if !d.is_finite() || d.abs() < 1e-9 { + break; + } + let step = (resid / d).clamp(-0.3, 0.3); + a = (a - step).clamp(lo, hi); + if step.abs() < 1e-10 { + break; + } + } + a +} + +pub fn alph_at_img_x( + ptype: ProjectionType, + proj: &dyn Projection, + parms: &ViewParams, + x_img: f64, + c_ref: f64, +) -> f64 { + let mut p = *parms; + p.a_center = c_ref; + // analytic base inverse (ignores distortion); Newton refines it + let seed = match ptype { + ProjectionType::Rectilinear => c_ref + (x_img / p.scale).atan(), + ProjectionType::Cylindrical => c_ref + x_img / p.scale, + }; + // rectilinear has a pole at |alph - center| = 90°; keep Newton inside + let (lo, hi) = match ptype { + ProjectionType::Rectilinear => (c_ref - 1.55, c_ref + 1.55), + ProjectionType::Cylindrical => (c_ref - 3.1, c_ref + 3.1), + }; + newton_inverse( + x_img, + |a| Some(proj.mac_x(ProjArgs::from_params(&p, a, 0.0))), + seed, + lo, + hi, + ) +} + +pub struct Compass { + pub style: CompassStyle, + pub center_deg: f64, + pub zoom: f32, + pub pan_x: f32, + pub parms: ViewParams, + pub ptype: ProjectionType, + /// horizontal field of view of the image, degrees (Ruler shading) + pub fov_deg: f64, + /// lens field-of-view half-angle limit in radians (0 = none). Ticks and + /// direction labels beyond this are dropped: a distortion model can fold + /// out-of-view directions back onto the strip, showing bogus ticks. + pub max_view_angle: f64, + pub on_change: Rc M>, +} + +impl Clone for Compass { + fn clone(&self) -> Self { + Compass { + style: self.style, + center_deg: self.center_deg, + zoom: self.zoom, + pan_x: self.pan_x, + parms: self.parms, + ptype: self.ptype, + fov_deg: self.fov_deg, + max_view_angle: self.max_view_angle, + on_change: self.on_change.clone(), + } + } +} + +#[derive(Default)] +struct CompassState { + /// (press x, center angle at press) + drag: Option<(f32, f64)>, + /// strong refs: the text pipeline only keeps a Weak to each paragraph + paragraphs: RefCell::Paragraph>>>, +} + +const DIRS: [(f64, &str); 8] = [ + (0.0, "N"), + (45.0, "NE"), + (90.0, "E"), + (135.0, "SE"), + (180.0, "S"), + (225.0, "SW"), + (270.0, "W"), + (315.0, "NW"), +]; + +/// one text to draw: content, x, y, center around x?, color +struct Txt { + text: String, + x: f32, + y: f32, + centered: bool, + color: Color, +} + +fn stroke_tick(frame: &mut Frame, x: f32, y0: f32, h: f32, w: f32, color: Color) { + if !x.is_finite() { + return; + } + // snap thin lines to integer coords for crisp 1 px strokes + let x = if w < 1.5 { x.round() } else { x }; + let p = Path::new(|p| { + p.move_to(Point::new(x, y0)); + p.line_to(Point::new(x, y0 - h)); + }); + frame.stroke( + &p, + Stroke { + width: w, + style: GeoStyle::Solid(color), + ..Default::default() + }, + ); +} + +impl Compass { + /// image-space x of azimuth `alph` (on the horizon), current parms + fn x_img_at(&self, proj: &dyn Projection, alph: f64) -> f64 { + proj.get_coordinates(alph, 0.0, &self.parms).0 + } + + fn alph_at_img_x( + &self, + proj: &dyn Projection, + x_img: f64, + c_ref: f64, + ) -> f64 { + alph_at_img_x(self.ptype, proj, &self.parms, x_img, c_ref) + } + + /// screen x of the image center for the Aligned style + fn center_x(&self, b: Rectangle) -> f32 { + b.x + b.width / 2.0 + self.pan_x + } +} + +impl WidgetTrait for Compass { + fn size(&self) -> Size { + Size::new(Length::Fill, Length::Fixed(STRIP_H)) + } + + fn tag(&self) -> tree::Tag { + tree::Tag::of::() + } + + fn state(&self) -> tree::State { + tree::State::new(CompassState::default()) + } + + fn layout( + &mut self, + _tree: &mut Tree, + _renderer: &R, + limits: &cosmic::iced::advanced::layout::Limits, + ) -> cosmic::iced::advanced::layout::Node { + cosmic::iced::advanced::layout::Node::new(Size::new(limits.max().width, STRIP_H)) + } + + fn draw( + &self, + tree: &Tree, + renderer: &mut R, + theme: &Theme, + _style: &RendererStyle, + layout: Layout<'_>, + _cursor: Cursor, + _viewport: &Rectangle, + ) { + use cosmic::iced::advanced::text::Renderer as _; + + let b = layout.bounds(); + let _ = theme; + // dark instrument panel, clearly distinct from the image area + let on: Color = Color::from_rgb(0.85, 0.85, 0.85); + let bg: Color = Color::from_rgb(0.16, 0.16, 0.16); + let divider: Color = Color::from_rgb(0.32, 0.32, 0.32); + let red = Color::from_rgb(1.0, 0.4, 0.25); + + // geometry rows (snapped so lines stay crisp) + let label_y = (b.y + 5.0).round(); + let y0 = (b.y + STRIP_H - 10.0).round(); // baseline + let tri_y = (b.y + 25.0).round(); + let z = self.zoom; + let proj = projection_from_type(self.ptype); + + let mut texts: Vec = Vec::new(); + let mut frame = Frame::with_bounds(renderer, b); + + // panel background + top divider: clearly a separate, draggable area + let panel = Path::new(|p| p.rectangle(b.position(), b.size())); + frame.fill(&panel, bg); + let sep = Path::new(|p| { + p.move_to(Point::new(b.x, b.y + 0.5)); + p.line_to(Point::new(b.x + b.width, b.y + 0.5)); + }); + frame.stroke( + &sep, + Stroke { + width: 1.0, + style: GeoStyle::Solid(divider), + ..Default::default() + }, + ); + + match self.style { + CompassStyle::Off => return, + CompassStyle::Aligned => { + let cx = self.center_x(b); + let c = self.parms.a_center; + // pixels per degree at the center + let ppd = ((self.x_img_at(&*proj, c + 1f64.to_radians()) + - self.x_img_at(&*proj, c - 1f64.to_radians())) + * z as f64 + / 2.0) + .abs(); + let step = [1.0f64, 2.0, 5.0, 10.0, 15.0, 30.0, 45.0, 90.0] + .into_iter() + .find(|s| s * ppd >= 6.0) + .unwrap_or(90.0); + + let a_lo = self + .alph_at_img_x(&*proj, ((b.x - cx) / z) as f64, c) + .min(c - std::f64::consts::PI); + let a_hi = self + .alph_at_img_x(&*proj, ((b.x + b.width - cx) / z) as f64, c) + .max(c + std::f64::consts::PI); + // the rectilinear projection folds directions behind the + // camera onto the same pixels (tan has 180 deg period); + // never draw beyond +-88 deg of the center. Cylindrical is + // linear in azimuth and may legitimately span > 180 deg. + let fold = match self.ptype { + ProjectionType::Rectilinear => 1.535, + ProjectionType::Cylindrical => std::f64::consts::PI, + }; + // additionally clamp to the lens FOV so distortion-folded + // out-of-view directions do not get ticks/labels + let fold = if self.max_view_angle > 0.0 { + fold.min(self.max_view_angle) + } else { + fold + }; + let a_lo = a_lo.max(c - fold); + let a_hi = a_hi.min(c + fold); + + // baseline + let line = Path::new(|p| { + p.move_to(Point::new(b.x, y0 + 0.5)); + p.line_to(Point::new(b.x + b.width, y0 + 0.5)); + }); + frame.stroke( + &line, + Stroke { + width: 1.0, + style: GeoStyle::Solid(on), + ..Default::default() + }, + ); + + // ticks (rise from the baseline, stay below the labels); + // `a` is radians, `step` degrees -> convert once + let step_rad = step.to_radians(); + let mut a = (a_lo / step_rad).floor() * step_rad; + while a <= a_hi { + let x = cx + self.x_img_at(&*proj, a) as f32 * z; + let rem = a.to_degrees().rem_euclid(45.0); + let major = rem.min(45.0 - rem) < step / 2.0; + stroke_tick(&mut frame, x, y0, if major { 12.0 } else { 6.0 }, 1.0, on); + a += step_rad; + } + + // direction labels on the top row; skip directions outside + // the +-88 deg window (they would fold onto the same pixel) + for (deg, name) in DIRS { + let alph = deg.to_radians(); + let rel = (alph - c + std::f64::consts::PI).rem_euclid(std::f64::consts::TAU) + - std::f64::consts::PI; + if rel.abs() > fold { + continue; + } + let x = cx + self.x_img_at(&*proj, alph) as f32 * z; + if !x.is_finite() || x < b.x - 20.0 || x > b.x + b.width + 20.0 { + continue; + } + let color = if deg == 0.0 { red } else { on }; + texts.push(Txt { + text: name.to_string(), + x, + y: label_y, + centered: true, + color, + }); + } + + // center marker: triangle + hairline (tick area only) + let tri = Path::new(|p| { + p.move_to(Point::new(cx - 6.0, tri_y)); + p.line_to(Point::new(cx + 6.0, tri_y)); + p.line_to(Point::new(cx, tri_y + 8.0)); + p.close(); + }); + frame.fill(&tri, red); + let hair = Path::new(|p| { + p.move_to(Point::new(cx.round(), tri_y + 8.0)); + p.line_to(Point::new(cx.round(), y0)); + }); + frame.stroke( + &hair, + Stroke { + width: 1.0, + style: GeoStyle::Solid(red), + ..Default::default() + }, + ); + + let deg = self.center_deg.rem_euclid(360.0); + texts.push(Txt { + text: format!("{deg:05.1}°"), + x: cx + 10.0, + y: tri_y - 3.0, + centered: false, + color: red, + }); + } + CompassStyle::Ruler => { + let bx = b.x; + let f = |alph_deg: f64| { + bx + (alph_deg.rem_euclid(360.0) / 360.0) as f32 * b.width + }; + let center = self.center_deg.rem_euclid(360.0); + let mx = f(center); + + // FOV band (between label row and baseline) + let half = (self.fov_deg / 2.0 / 360.0) as f32 * b.width; + let x0 = (mx - half).max(bx); + let x1 = (mx + half).min(bx + b.width); + if x1 > x0 { + let band = Path::new(|p| { + p.rectangle( + Point::new(x0, tri_y), + Size::new(x1 - x0, y0 - tri_y), + ); + }); + frame.fill(&band, Color::from_rgba(0.4, 0.65, 1.0, 0.20)); + } + + // baseline + ticks every 5 deg + let line = Path::new(|p| { + p.move_to(Point::new(bx, y0 + 0.5)); + p.line_to(Point::new(bx + b.width, y0 + 0.5)); + }); + frame.stroke( + &line, + Stroke { + width: 1.0, + style: GeoStyle::Solid(on), + ..Default::default() + }, + ); + for deg in (0..360).step_by(5) { + let major = deg % 45 == 0; + stroke_tick( + &mut frame, + f(deg as f64), + y0, + if major { 12.0 } else { 6.0 }, + 1.0, + on, + ); + } + for (deg, name) in DIRS { + let color = if deg == 0.0 { red } else { on }; + texts.push(Txt { + text: name.to_string(), + x: f(deg), + y: label_y, + centered: true, + color, + }); + } + + // marker + let tri = Path::new(|p| { + p.move_to(Point::new(mx - 6.0, tri_y)); + p.line_to(Point::new(mx + 6.0, tri_y)); + p.line_to(Point::new(mx, tri_y + 8.0)); + p.close(); + }); + frame.fill(&tri, red); + texts.push(Txt { + text: format!("{center:05.1}°"), + x: mx + 10.0, + y: tri_y - 3.0, + centered: false, + color: red, + }); + } + } + + renderer.draw_geometry(frame.into_geometry()); + + // text (kept alive in state until next frame; positions rounded + // so glyphs land on pixel boundaries) + let font = renderer.default_font(); + let state = tree.state.downcast_ref::(); + { + let mut paras = state.paragraphs.borrow_mut(); + paras.clear(); + for t in texts.iter() { + paras.push(Plain::<::Paragraph>::new(adv_text::Text { + content: t.text.clone(), + bounds: Size::new(120.0, 20.0), + size: Pixels(13.0), + line_height: LineHeight::default(), + font, + align_x: alignment::Horizontal::Left.into(), + align_y: alignment::Vertical::Top.into(), + shaping: adv_text::Shaping::Basic, + wrapping: adv_text::Wrapping::None, + ellipsize: adv_text::Ellipsize::None, + })); + } + } + let paras = state.paragraphs.borrow(); + for (gi, t) in texts.iter().enumerate() { + let para = ¶s[gi]; + let w = para.min_bounds().width; + let x = if t.centered { t.x - w / 2.0 } else { t.x }; + renderer.fill_paragraph( + para.raw(), + Point::new(x.round(), t.y.round()), + t.color, + b, + ); + } + } + + fn update( + &mut self, + tree: &mut Tree, + event: &Event, + layout: Layout<'_>, + cursor: Cursor, + _renderer: &R, + _clipboard: &mut dyn Clipboard, + shell: &mut Shell<'_, M>, + _viewport: &Rectangle, + ) { + let state = tree.state.downcast_mut::(); + let b = layout.bounds(); + let proj = projection_from_type(self.ptype); + + match event { + Event::Mouse(mouse::Event::ButtonPressed(Button::Left)) => { + let Some(p) = cursor.position_over(b) else { return }; + state.drag = Some((p.x, self.center_deg)); + if self.style == CompassStyle::Ruler { + let alph = ((p.x - b.x) / b.width * 360.0) as f64; + shell.publish((self.on_change)(alph)); + } + } + Event::Mouse(mouse::Event::ButtonReleased(Button::Left)) => { + state.drag = None; + } + Event::Mouse(mouse::Event::CursorMoved { .. }) => { + let Some((x0, c0)) = state.drag else { return }; + let Some(p) = cursor.position() else { return }; + match self.style { + CompassStyle::Aligned => { + let cx = self.center_x(b); + let z = self.zoom; + let a0 = + self.alph_at_img_x(&*proj, ((x0 - cx) / z) as f64, c0.to_radians()); + let a1 = + self.alph_at_img_x(&*proj, ((p.x - cx) / z) as f64, c0.to_radians()); + let new_c = (c0.to_radians() - (a1 - a0)).to_degrees(); + shell.publish((self.on_change)(new_c)); + } + CompassStyle::Ruler => { + let alph = ((p.x - b.x) / b.width * 360.0) as f64; + shell.publish((self.on_change)(alph)); + } + CompassStyle::Off => {} + } + } + _ => {} + } + } + + fn mouse_interaction( + &self, + _tree: &Tree, + layout: Layout<'_>, + cursor: Cursor, + _viewport: &Rectangle, + _renderer: &R, + ) -> Interaction { + if cursor.position_over(layout.bounds()).is_some() { + Interaction::Grab + } else { + Interaction::None + } + } +} + +#[cfg(test)] +mod tests { + use super::*; + + fn roundtrip(ptype: ProjectionType, k0: f64, k1: f64, x0: f64) { + let proj = projection_from_type(ptype); + let mut parms = ViewParams::new(); + parms.a_center = 0.3; + parms.scale = 3000.0; + parms.k0 = k0; + parms.k1 = k1; + parms.x0 = x0; + for x in [-2500.0, -800.0, -50.0, 0.0, 120.0, 900.0, 2400.0] { + let a = alph_at_img_x(ptype, &*proj, &parms, x, parms.a_center); + let (x2, _) = proj.get_coordinates(a, 0.0, &parms); + assert!( + (x2 - x).abs() < 1e-6, + "{ptype:?} x={x}: round-trips to {x2} (alph {a})" + ); + } + } + + #[test] + fn inverse_rectilinear() { + roundtrip(ProjectionType::Rectilinear, 0.0, 0.0, 0.0); + } + #[test] + fn inverse_rectilinear_distorted() { + // realistic magnitudes: x0 is in normalized (pre-scale) units + roundtrip(ProjectionType::Rectilinear, 0.002, -0.0003, 0.01); + } + #[test] + fn inverse_cylindrical() { + roundtrip(ProjectionType::Cylindrical, 0.0, 0.0, 0.0); + } +} diff --git a/crates/gipfel-app/src/dump.rs b/crates/gipfel-app/src/dump.rs new file mode 100644 index 0000000..0bcfe56 --- /dev/null +++ b/crates/gipfel-app/src/dump.rs @@ -0,0 +1,285 @@ +//! Screen dump: render markers/labels onto the photo and save as PNG/JPEG. +//! Replaces FLTK's window capture from the original "Screen Dump" menu item. + +use image::RgbaImage; +use std::path::Path; + +pub struct DumpMarker { + pub x: f64, + pub y: f64, + pub name: String, + pub known: bool, + pub hidden: bool, +} + +pub struct DumpTrackPt { + pub x: f64, + pub y: f64, + pub hidden: bool, + pub width: f64, +} + +fn px(img: &mut RgbaImage, x: i64, y: i64, c: [u8; 4]) { + if x >= 0 && y >= 0 && (x as u32) < img.width() && (y as u32) < img.height() { + img.put_pixel(x as u32, y as u32, image::Rgba(c)); + } +} + +fn line(img: &mut RgbaImage, x0: f64, y0: f64, x1: f64, y1: f64, c: [u8; 4], w: i64) { + let dx = x1 - x0; + let dy = y1 - y0; + let n = (dx.hypot(dy) as i64).max(1); + for i in 0..=n { + let t = i as f64 / n as f64; + let (x, y) = (x0 + dx * t, y0 + dy * t); + for oy in -w..=w { + for ox in -w..=w { + px(img, x.round() as i64 + ox, y.round() as i64 + oy, c); + } + } + } +} + +fn circle(img: &mut RgbaImage, cx: f64, cy: f64, r: f64, c: [u8; 4]) { + let n = (2.0 * std::f64::consts::PI * r.max(1.0)) as i64; + for i in 0..=n { + let a = i as f64 / n as f64 * 2.0 * std::f64::consts::PI; + px(img, (cx + r * a.cos()).round() as i64, (cy + r * a.sin()).round() as i64, c); + px(img, (cx + (r + 1.0) * a.cos()).round() as i64, (cy + (r + 1.0) * a.sin()).round() as i64, c); + } +} + +fn tri(img: &mut RgbaImage, p: [(f64, f64); 3], c: [u8; 4]) { + let ymin = p.iter().map(|q| q.y() as i64).min().unwrap(); + let ymax = p.iter().map(|q| q.y() as i64).max().unwrap(); + for y in ymin..=ymax { + let mut xmin = i64::MAX; + let mut xmax = i64::MIN; + for i in 0..3 { + let a = p[i]; + let b = p[(i + 1) % 3]; + if (a.y() as i64) == (b.y() as i64) { + continue; + } + if y >= (a.y() as i64).min(b.y() as i64) && y <= (a.y() as i64).max(b.y() as i64) { + let t = (y as f64 - a.y()) / (b.y() - a.y()); + let x = a.x() + t * (b.x() - a.x()); + xmin = xmin.min(x.round() as i64); + xmax = xmax.max(x.round() as i64); + } + } + for x in xmin..=xmax { + px(img, x, y, c); + } + } +} + +trait XY { + fn x(self) -> f64; + fn y(self) -> f64; +} +impl XY for (f64, f64) { + fn x(self) -> f64 { + self.0 + } + fn y(self) -> f64 { + self.1 + } +} + +fn load_font(family: &str) -> Option { + // Resolve the COSMIC interface font (or system sans-serif fallback) through + // the same font discovery cosmic's text engine uses, then hand the bytes to + // fontdue for the bitmap dump. + static FS: std::sync::OnceLock> = + std::sync::OnceLock::new(); + let mut fs = FS + .get_or_init(|| std::sync::Mutex::new(cosmic_text::FontSystem::new())) + .lock() + .unwrap(); + let fam = [fontdb::Family::Name(family)]; + let id = fs + .db() + .query(&fontdb::Query { families: &fam, ..Default::default() }) + .or_else(|| fs.db().query(&fontdb::Query::default()))?; + let font = fs.get_font(id, fontdb::Weight::NORMAL)?; + fontdue::Font::from_bytes( + font.data().to_vec(), + fontdue::FontSettings { scale: 64.0, ..Default::default() }, + ) + .ok() +} + +struct GlyphBlit { + x: f32, + y: f32, + width: usize, + #[allow(dead_code)] + height: usize, + pixels: Vec, +} + +fn layout_text(font: &fontdue::Font, text: &str, size: f32) -> Vec { + use fontdue::layout::{CoordinateSystem, Layout as FontLayout, TextStyle}; + let mut l = FontLayout::new(CoordinateSystem::PositiveYDown); + l.append(&[font], &TextStyle::new(text, size, 0)); + l.glyphs() + .iter() + .map(|g| { + let (m, pixels) = font.rasterize_config(g.key); + GlyphBlit { x: g.x, y: g.y, width: m.width, height: m.height, pixels } + }) + .collect() +} + +fn draw_text(img: &mut RgbaImage, font: &fontdue::Font, x: f64, y: f64, text: &str, c: [u8; 4], size: f32) { + for glyph in layout_text(font, text, size) { + let gx0 = (x + glyph.x as f64).round() as i64; + let gy0 = (y + glyph.y as f64).round() as i64; + for (i, &alpha) in glyph.pixels.iter().enumerate() { + if alpha == 0 { + continue; + } + let gx = gx0 + (i % glyph.width) as i64; + let gy = gy0 + (i / glyph.width) as i64; + if gx < 0 || gy < 0 || gx >= img.width() as i64 || gy >= img.height() as i64 { + continue; + } + let a = alpha as u32; + let old = img.get_pixel(gx as u32, gy as u32).0; + let blend = |o: u8, n: u8| ((o as u32 * (255 - a) + n as u32 * a) / 255) as u8; + px(img, gx, gy, [blend(old[0], c[0]), blend(old[1], c[1]), blend(old[2], c[2]), 255]); + } + } +} + +fn text_width(font: &fontdue::Font, text: &str, size: f32) -> f64 { + layout_text(font, text, size) + .iter() + .map(|g| (g.x + g.width as f32) as f64) + .fold(0.0, f64::max) +} + +/// Compose the dump: photo + crosses/leaders/labels + track. +/// +/// `font_size` is the on-screen label font size in screen pixels, `scale` +/// converts screen pixels to image pixels (1 / zoom): the overlay geometry +/// is specified in screen pixels (as seen on the display) but the dump is +/// rendered at the image's native resolution, so every size must be scaled +/// accordingly. Label placement runs the exact same collision-avoiding +/// layout as the on-screen viewport (in screen space, then mapped back to +/// image space), so the dump shows the labels at the same positions. +pub fn render( + img: &RgbaImage, + markers: &[DumpMarker], + track: &[DumpTrackPt], + out: &Path, + font_size: f32, + scale: f64, + font_family: &str, +) -> Result<(), String> { + let mut canvas = img.clone(); + let (w, h) = (img.width() as f64, img.height() as f64); + let (cx, cy) = (w / 2.0, h / 2.0); + let s = scale.max(0.1); + let font_img = f64::from(font_size) * s; + let font = load_font(font_family); + + // track + for win in track.windows(2) { + let (a, b) = (&win[0], &win[1]); + let c = if a.hidden { [51, 102, 255, 255] } else { [255, 51, 51, 255] }; + line(&mut canvas, cx + a.x, cy + a.y, cx + b.x, cy + b.y, c, (a.width as i64).max(1)); + } + + // label layout: the exact same collision-avoiding code as the on-screen + // viewport (see viewport::layout_label_y). It runs in screen space + // (marker positions / scale, screen font size); the resulting offsets + // are mapped back to image space with s. Needed for the leader lines + // and the label drawing below. + let line_h = f64::from(font_size) * 1.3 + 2.0; + let in_img: Vec = markers + .iter() + .map(|m| m.x.abs() * 2.0 <= w && m.y.abs() * 2.0 <= h) + .collect(); + let layout: Option> = font.as_ref().map(|font| { + let sx: Vec = markers.iter().map(|m| (m.x / s) as f32).collect(); + let sy: Vec = markers.iter().map(|m| (m.y / s) as f32).collect(); + let widths: Vec = markers + .iter() + .map(|m| text_width(font, &m.name, font_size) as f32) + .collect(); + crate::viewport::layout_label_y(&sx, &sy, &widths, &in_img, line_h as f32) + }); + + let yellow = [255, 230, 51, 255]; + // ~1 screen pixel stroke at image resolution + let lw = ((s - 1.0) / 2.0).round().max(0.0) as i64; + // like the screen: known-hill flags in red + let flag_color = [255, 64, 64, 255]; + for (i, m) in markers.iter().enumerate() { + let (x, y) = (cx + m.x, cy + m.y); + if x < -10.0 || y < -10.0 || x > w + 10.0 || y > h + 10.0 || !in_img[i] { + continue; + } + // yellow cross + leader line up to the label (screen section 3) + let d = 4.0 * s; + line(&mut canvas, x - d, y, x + d, y, yellow, lw); + if let Some(lay) = &layout { + let (lx, ly) = lay[i]; + let label_line = y + (f64::from(ly) - line_h) * s; + line(&mut canvas, x, y, x, label_line, yellow, lw); + line(&mut canvas, x, label_line, x + f64::from(lx) * s, label_line, yellow, lw); + } + } + // known-hill flags (pole + pennant + circle), on top of all leaders + for (i, m) in markers.iter().enumerate() { + if !m.known || !in_img[i] { + continue; + } + let (x, y) = (cx + m.x, cy + m.y); + if x < -10.0 || y < -10.0 || x > w + 10.0 || y > h + 10.0 { + continue; + } + let pole = 16.0 * s; + line(&mut canvas, x, y, x, y - pole, flag_color, (s.round() as i64).max(1)); + tri( + &mut canvas, + [(x, y - pole), (x + 10.0 * s, y - pole * 0.8), (x, y - pole * 0.625)], + flag_color, + ); + circle(&mut canvas, x, y, 3.0 * s, flag_color); + } + + if let (Some(font), Some(layout)) = (&font, &layout) { + for (i, m) in markers.iter().enumerate() { + if !in_img[i] { + continue; + } + let (_, ly) = layout[i]; + let (x, y) = (cx + m.x + 2.0 * s, cy + m.y + (f64::from(ly) - line_h) * s); + if x < 0.0 || y < 0.0 || x > w || y > h { + continue; + } + let tw = text_width(font, &m.name, font_img as f32); + for yy in (y - 2.0 * s) as i64..=(y + font_img + 2.0 * s) as i64 { + for xx in (x - 1.0 * s) as i64..=(x + tw + 1.0 * s) as i64 { + px(&mut canvas, xx, yy, [0, 0, 0, 160]); + } + } + let c = if m.hidden { [102, 153, 255, 255] } else { [255, 255, 255, 255] }; + draw_text(&mut canvas, font, x, y, &m.name, c, font_img as f32); + } + } + + if out + .extension() + .is_some_and(|e| e.eq_ignore_ascii_case("jpg") || e.eq_ignore_ascii_case("jpeg")) + { + let rgb = image::DynamicImage::ImageRgba8(canvas).to_rgb8(); + rgb.save(out.with_extension("jpg")).map_err(|e| e.to_string())?; + } else { + canvas.save(out.with_extension("png")).map_err(|e| e.to_string())?; + } + Ok(()) +} diff --git a/crates/gipfel-app/src/main.rs b/crates/gipfel-app/src/main.rs new file mode 100644 index 0000000..4cf0c63 --- /dev/null +++ b/crates/gipfel-app/src/main.rs @@ -0,0 +1,1805 @@ +//! gipfel for the COSMIC desktop: annotate mountain panoramas. +//! +//! M2/M3: full app skeleton (menus, dialogs, settings sidebar) plus the +//! interactive viewport (marker placement, refit, context menu, track). + +mod compass; +mod dump; +mod viewport; + +use std::marker::PhantomData; +use std::path::{Path, PathBuf}; +use std::rc::Rc; + +use cosmic::app::{Core, Task}; +use cosmic::iced::keyboard::{self, Key}; +use cosmic::iced::window::Event as WinEvent; +use cosmic::iced::{event, Color, Event, Length, Subscription, Vector}; +use cosmic::widget::{self, menu, menu::action::MenuAction}; +use cosmic::{Application, ApplicationExt, Element, Theme}; + +use gipfel_core::hill::{self, Hill, FLAG_DUPLIC, FLAG_HIDDEN, FLAG_TRACK_POINT, FLAG_VISIBLE}; +use gipfel_core::metadata::ImageMetaData; +use gipfel_core::panorama::Panorama; +use gipfel_core::projection::ProjectionType; +use gipfel_core::stitch; +use gipfel_core::track; + +use compass::{Compass, CompassStyle, STRIP_H}; + +use viewport::{Marker, TrackPt, ViewportWidget, VpMsg}; + +/// actions of the right-click popover menu on the viewport. +/// The click position / hit index are baked into the action at menu-build +/// time (like cosmic-files bakes its entity), so the handler never depends +/// on state that may have changed while the menu was open. +#[derive(Clone, Copy, PartialEq)] +pub enum CtxAction { + /// image-space position of the right-click + FindPeak(f64, f64), + /// hill index of the clicked marker + ToggleHidden(usize), +} + +impl Eq for CtxAction {} + +impl MenuAction for CtxAction { + type Message = Message; + + fn message(&self) -> Message { + match self { + CtxAction::FindPeak(x, y) => Message::ContextFindPeak(*x, *y), + CtxAction::ToggleHidden(i) => Message::ContextToggleHidden(*i), + } + } +} + +pub struct Flags { + pub data: Option, + pub image: Option, +} + +pub struct Gipfel { + core: Core, + pan: Panorama, + project: Option, + /// hill indices of the known (flagged) markers + known: Vec, + /// hill indices of the GPS track points, in file order + track: Vec, + // viewport state + zoom: f32, + fit: bool, + panv: Vector, + hover: Option, + ctrl: bool, + version: u64, + win_w: f64, + win_h: f64, + vp_w: f64, + vp_h: f64, + // sidebar state + lat_in: String, + lon_in: String, + height_in: String, + k0_in: String, + k1_in: String, + x0_in: String, + profile_in: String, + profiles: Vec, + profile_overwrite: bool, + track_width: f64, + show_hidden: bool, + fit_distortion: bool, + label_font: f32, + compass_style: CompassStyle, + status: String, + // overlays + picker: Option, + context: Option, + /// peak info shown in the status bar while the mouse hovers a marker + hover_info: Option, + about: bool, +} + +struct Project { + path: PathBuf, + handle: cosmic::widget::image::Handle, + w: u32, + h: u32, + md: ImageMetaData, +} + +struct Picker { + title: String, + query: String, + mode: PickerMode, +} + +#[derive(Clone, Copy)] +enum PickerMode { + Viewpoint, + FindPeak { img_x: f64, img_y: f64 }, +} + +struct ContextState { + img_x: f64, + img_y: f64, + hit: Option, +} + +#[derive(Debug, Clone)] +pub enum Message { + Vp(VpMsg), + Event(Event), + WinResized(f64, f64), + FileDropped(Vec), + + OpenImage, + OpenPath(Option), + ImageLoaded(Box>), + SaveImage, + SavePath(Option), + ImageSaved(Option), + DumpImage, + ImageDumped(Box>), + StitchTiff, + StitchDone(Box>), + LoadTrack, + TrackLoaded(Box, String>>), + + ChooseViewpoint, + PickerQuery(String), + PickerPick(usize), + PickerClose, + ContextFindPeak(f64, f64), + ContextToggleHidden(usize), + Surface(cosmic::surface::Action), + Redraw, + + Center(f64), + Pitch(f64), + Roll(f64), + LabelFont(f32), + Focal(f64), + Visibility(f64), + TrackWidth(f64), + Lat(String), + Lon(String), + Height(String), + ViewApply, + K0(String), + K1(String), + X0(String), + DistortionApply, + ProfileName(String), + ProfileSave, + ProfileLoad, + ProfileSelect(usize), + Projection(ProjectionType), + ToggleShowHidden, + FitDistortion(bool), + ToggleNav, + CompassStyle(CompassStyle), + + Fit, + ZoomIn, + ZoomOut, + About(bool), + Quit, +} + +#[derive(Clone, Copy, Debug, Eq, PartialEq)] +enum Action { + Open, + Save, + Dump, + Stitch, + LoadTrack, + Quit, + ChooseViewpoint, + ShowHidden, + Fit, + ZoomIn, + ZoomOut, + ProjNormal, + ProjPanoramic, + CompassAligned, + CompassRuler, + CompassOff, + About, +} + +impl MenuAction for Action { + type Message = Message; + fn message(&self) -> Message { + match self { + Action::Open => Message::OpenImage, + Action::Save => Message::SaveImage, + Action::Dump => Message::DumpImage, + Action::Stitch => Message::StitchTiff, + Action::LoadTrack => Message::LoadTrack, + Action::Quit => Message::Quit, + Action::ChooseViewpoint => Message::ChooseViewpoint, + Action::ShowHidden => Message::ToggleShowHidden, + Action::Fit => Message::Fit, + Action::ZoomIn => Message::ZoomIn, + Action::ZoomOut => Message::ZoomOut, + Action::ProjNormal => Message::Projection(ProjectionType::Rectilinear), + Action::ProjPanoramic => Message::Projection(ProjectionType::Cylindrical), + Action::CompassAligned => Message::CompassStyle(CompassStyle::Aligned), + Action::CompassRuler => Message::CompassStyle(CompassStyle::Ruler), + Action::CompassOff => Message::CompassStyle(CompassStyle::Off), + Action::About => Message::About(true), + } + } +} + +impl From for Message { + fn from(m: VpMsg) -> Self { + Message::Vp(m) + } +} + +#[derive(Clone)] +pub struct LoadedImage { + path: PathBuf, + w: u32, + h: u32, + handle: cosmic::widget::image::Handle, + md: ImageMetaData, +} + +impl std::fmt::Debug for LoadedImage { + fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { + write!(f, "LoadedImage({})", self.path.display()) + } +} + +/// Load a GPS track from a `.gpx` (GPS Exchange) or `.csv`/`.dat` (gipfel +/// hill-CSV) file, returning hills flagged as track points. +fn load_track_file(path: &Path) -> Result, String> { + if path.extension().is_some_and(|e| e.eq_ignore_ascii_case("gpx")) { + track::load_gpx(&path.to_string_lossy(), 2000) + } else { + let mut hs = hill::load(&path.to_string_lossy()).map_err(|e| e.to_string())?; + for h in hs.iter_mut() { + h.flags |= FLAG_TRACK_POINT; + } + Ok(hs) + } +} + +fn data_candidates(flag: Option<&Path>) -> Vec { + let mut v = Vec::new(); + if let Some(p) = flag { + v.push(p.to_path_buf()); + } + if let Ok(e) = std::env::var("GIPFEL_DATA") { + v.push(PathBuf::from(e)); + } + v.push(PathBuf::from("gipfel.dat")); + if let Some(home) = std::env::var_os("HOME") { + v.push(PathBuf::from(home).join(".local/share/gipfel/gipfel.dat")); + } + v.push(PathBuf::from("/usr/local/share/gipfel/gipfel.dat")); + v.push(PathBuf::from("/usr/share/gipfel/gipfel.dat")); + v +} + +/// Inverse projection for the vertical axis: the elevation whose projected +/// pixel row (image-centered coords, center column) is `y_img`. Shares the +/// damped Newton iteration with the compass strip (`compass::newton_inverse`). +fn elev_at_img_y(pan: &Panorama, y_img: f64) -> f64 { + let c = pan.parms.a_center; + let p = pan.parms.a_pitch; + let seed = p - (y_img / pan.parms.scale).atan(); + compass::newton_inverse( + y_img, + |e| pan.project_point(c, e).map(|v| v.1), + seed, + p - 1.535, + p + 1.535, + ) +} + +impl Gipfel { + fn bump(&mut self) { + // keep the visibility limit in sync with the lens' actual field of + // view: the angle at which the (distorted) projection reaches the + // image edge. Directions beyond it cannot be inside the photo, even + // if the distortion model folds their projection back into it. + let limit = self.view_angle_limit(); + if (limit - self.pan.max_view_angle).abs() > 1e-9 { + self.pan.max_view_angle = limit; + // keep the user-placed marker positions (observations) intact + self.pan.update_visible_mountains(Some(&self.known)); + } + self.version += 1; + } + + /// Reload the stored profile names from the prefs file. + fn refresh_profiles(&mut self) { + self.profiles = gipfel_core::profile::list().into_iter().map(|(n, _)| n).collect(); + } + + /// Load a stored profile and apply its distortion parameters. + fn apply_profile(&mut self, name: &str) { + match gipfel_core::profile::load(name) { + Ok(p) => { + self.pan.set_distortion_params(p.k0, p.k1, p.x0); + self.sync_inputs(); + self.bump(); + self.status = format!("loaded profile '{name}'"); + } + Err(e) => self.status = format!("profile '{name}': {e}"), + } + } + + /// half-angle (radians) the projection maps the image's right edge to, + /// i.e. the horizontal FOV of the lens including fitted distortion. + /// 0 = no image loaded / not computable. + fn view_angle_limit(&self) -> f64 { + let Some(p) = &self.project else { + return 0.0; + }; + let proj = gipfel_core::projection::projection_from_type(self.pan.projection_type); + let c = self.pan.parms.a_center; + let a_edge = compass::alph_at_img_x( + self.pan.projection_type, + &*proj, + &self.pan.parms, + p.w as f64 / 2.0, + c, + ); + (a_edge - c).abs() + } + + fn img_max(&self) -> f64 { + self.project + .as_ref() + .map(|p| p.w.max(p.h) as f64) + .unwrap_or(1000.0) + } + + fn focal_35mm(&self) -> f64 { + self.pan.get_scale() * 35.0 / self.img_max() + } + + fn set_focal_35mm(&mut self, f: f64) { + self.pan.set_scale(f * self.img_max() / 35.0); + } + + fn sync_inputs(&mut self) { + self.lat_in = format!("{:.5}", self.pan.get_view_lat()); + self.lon_in = format!("{:.5}", self.pan.get_view_long()); + self.height_in = format!("{:.0}", self.pan.get_view_height()); + self.k0_in = format!("{:.5}", self.pan.parms.k0); + self.k1_in = format!("{:.5}", self.pan.parms.k1); + self.x0_in = format!("{:.1}", self.pan.parms.x0); + } + + fn fit_zoom(&self) -> f32 { + let p = match &self.project { + Some(p) => p, + None => return 1.0, + }; + if self.vp_w < 10.0 || self.vp_h < 10.0 { + return 1.0; + } + (self.vp_w / p.w as f64).min(self.vp_h / p.h as f64) as f32 + } + + fn markers(&self) -> Rc> { + let mut v = Vec::new(); + for &i in &self.pan.visible { + let h = &self.pan.hills[i]; + if h.is_track_point() || h.flags & FLAG_DUPLIC != 0 { + continue; + } + let hidden = h.flags & FLAG_HIDDEN != 0; + if hidden && !self.show_hidden { + continue; + } + v.push(Marker { + x: h.x, + y: h.y, + name: h.name.as_str().into(), + idx: i, + known: self.known.contains(&i), + hidden, + }); + } + Rc::new(v) + } + + fn track_pts(&self) -> Rc> { + let v = self + .track + .iter() + .filter(|&&i| self.pan.hills[i].flags & FLAG_VISIBLE != 0) + .map(|&i| { + let h = &self.pan.hills[i]; + // C: max(scale * track_width / (dist_m * 10), 1) - distance in + // meters. (Was wrongly divided by km here, inflating widths + // 1000x and breaking stroke triangulation.) Clamped high so + // near-camera points can't produce absurd strokes. + let dist_m = self.pan.get_real_distance(i).max(1.0); + let width = + (self.pan.get_scale() * self.track_width / (dist_m * 10.0)).clamp(1.0, 1000.0); + TrackPt { + x: h.x, + y: h.y, + hidden: h.flags & FLAG_HIDDEN != 0, + width, + } + }) + .collect(); + Rc::new(v) + } + + /// recompute viewport metrics from window size, nav and compass state + fn remeasure(&mut self) { + let nav = if self.core().nav_bar_active() { 340.0 } else { 0.0 }; + let strip = if self.project.is_some() && self.compass_style != CompassStyle::Off { + STRIP_H as f64 + } else { + 0.0 + }; + self.vp_w = (self.win_w - nav).max(100.0); + self.vp_h = (self.win_h - 100.0 - strip).max(100.0); + if self.fit { + self.zoom = self.fit_zoom(); + } + } + + fn refit(&mut self) { + if self.known.is_empty() { + return; + } + if let Err(e) = self.pan.comp_params(&self.known, self.fit_distortion) { + self.status = format!("refit failed: {e}"); + } + self.sync_inputs(); + self.bump(); + } + + fn md_from_state(&self) -> ImageMetaData { + let mut md = ImageMetaData::new(); + if let Some(p) = &self.project { + md.manufacturer = p.md.manufacturer.clone(); + md.model = p.md.model.clone(); + } + md.longitude = self.pan.get_view_long(); + md.latitude = self.pan.get_view_lat(); + md.height = self.pan.get_view_height(); + md.direction = self.pan.get_center_angle(); + md.pitch = self.pan.get_pitch_angle(); + md.roll = self.pan.get_roll_angle(); + md.focal_length_35mm = self.focal_35mm(); + md.projection_type = self.pan.projection_type as i32; + md.k0 = self.pan.parms.k0; + md.k1 = self.pan.parms.k1; + md.x0 = self.pan.parms.x0; + md + } + + fn open_image(&mut self, path: PathBuf) -> Task { + let task = async move { + let img = image::open(&path).map_err(|e| format!("{}: {e}", path.display()))?; + let (w, h) = (img.width(), img.height()); + let rgba = img.to_rgba8(); + let handle = cosmic::widget::image::Handle::from_rgba(w, h, rgba.into_raw()); + let mut md = ImageMetaData::new(); + md.load_image(&path.to_string_lossy()); + Ok(LoadedImage { path, w, h, handle, md }) + }; + cosmic::task::future(async move { Message::ImageLoaded(Box::new(task.await)) }) + } + + fn picker_entries(&self, p: &Picker) -> Vec<(String, usize)> { + let mut idxs: Vec = match p.mode { + PickerMode::Viewpoint => (0..self.pan.hills.len()).collect(), + PickerMode::FindPeak { .. } => self.pan.close.clone(), + }; + let q = p.query.to_uppercase(); + let mut entries: Vec<(String, usize)> = idxs + .drain(..) + .filter(|&i| { + let h = &self.pan.hills[i]; + h.flags & (FLAG_DUPLIC | FLAG_TRACK_POINT) == 0 + && (q.is_empty() || h.name.to_uppercase().contains(&q)) + }) + .map(|i| (format!("{} ({} m)", self.pan.hills[i].name, self.pan.hills[i].height as i64), i)) + .collect(); + entries.sort_by(|a, b| a.0.cmp(&b.0)); + entries.truncate(300); + entries + } + + fn update_vp(&mut self, msg: VpMsg) -> Task { + match msg { + VpMsg::ZoomAt { factor, img } => { + let new_zoom = (self.zoom * factor).clamp(0.02, 16.0); + let z = self.zoom as f64; + let nz = new_zoom as f64; + self.panv = Vector::new( + self.panv.x + (img.x as f64 * (z - nz)) as f32, + self.panv.y + (img.y as f64 * (z - nz)) as f32, + ); + self.zoom = new_zoom; + self.fit = false; + } + VpMsg::PanBy { delta } => { + self.panv = Vector::new(self.panv.x + delta.x, self.panv.y + delta.y); + self.fit = false; + } + VpMsg::PanTo { pan } => { + self.panv = pan; + self.fit = false; + } + VpMsg::MoveMarker { idx, img } => { + if idx < self.pan.hills.len() { + self.pan.hills[idx].x = img.x as f64; + self.pan.hills[idx].y = img.y as f64; + self.bump(); + } + } + VpMsg::DragEnd { had_marker } => { + if had_marker { + self.refit(); + } + } + VpMsg::ToggleKnown(hit) => { + if let Some(i) = hit { + if let Some(pos) = self.known.iter().position(|&k| k == i) { + self.known.remove(pos); + } else { + self.known.push(i); + } + self.refit(); + } + } + VpMsg::ContextAt { img, hit } => { + self.context = Some(ContextState { img_x: img.x as f64, img_y: img.y as f64, hit }); + } + VpMsg::Hover(h) => { + if h != self.hover { + self.hover = h; + self.hover_info = h.map(|i| { + let hm = &self.pan.hills[i]; + format!( + "{} ({} m), distance {:.2} km", + hm.name, + hm.height as i64, + self.pan.get_real_distance(i) / 1000.0 + ) + }); + } + } + } + Task::none() + } +} + +impl Application for Gipfel { + type Executor = cosmic::executor::Default; + type Flags = Flags; + type Message = Message; + const APP_ID: &'static str = "org.gipfel.App"; + + fn core(&self) -> &Core { + &self.core + } + + fn core_mut(&mut self) -> &mut Core { + &mut self.core + } + + fn init(core: Core, flags: Flags) -> (Self, Task) { + let mut pan = Panorama::new(); + let mut status = String::from("no hill database loaded"); + for c in data_candidates(flags.data.as_deref()) { + if c.is_file() { + match pan.load_data(&c.to_string_lossy()) { + Ok(()) => { + status = format!("{} hills loaded", pan.hills.len()); + break; + } + Err(e) => status = format!("cannot read {}: {e}", c.display()), + } + } + } + + let mut app = Gipfel { + core, + pan, + project: None, + known: Vec::new(), + track: Vec::new(), + zoom: 1.0, + fit: true, + panv: Vector::ZERO, + hover: None, + ctrl: false, + version: 0, + win_w: 1280.0, + win_h: 800.0, + vp_w: 0.0, + vp_h: 0.0, + lat_in: String::new(), + lon_in: String::new(), + height_in: String::new(), + k0_in: String::new(), + k1_in: String::new(), + x0_in: String::new(), + profile_in: String::new(), + profiles: gipfel_core::profile::list().into_iter().map(|(n, _)| n).collect(), + profile_overwrite: false, + track_width: 200.0, + show_hidden: false, + label_font: viewport::DEFAULT_LABEL_FONT, + fit_distortion: false, + compass_style: CompassStyle::Aligned, + status, + picker: None, + context: None, + hover_info: None, + about: false, + }; + app.sync_inputs(); + + let task = match flags.image { + Some(p) => app.open_image(p), + None => Task::none(), + }; + (app, task) + } + + fn header_start(&self) -> Vec> { + let mut kb = std::collections::HashMap::new(); + use cosmic::widget::menu::key_bind::{KeyBind, Modifier}; + let kb = &mut kb; + for (ch, action) in [ + ("o", Action::Open), + ("s", Action::Save), + ("d", Action::Dump), + ("t", Action::LoadTrack), + ("v", Action::ChooseViewpoint), + ("f", Action::Fit), + ("q", Action::Quit), + ] { + kb.insert( + KeyBind { modifiers: vec![Modifier::Ctrl], key: Key::Character(ch.into()) }, + action, + ); + } + + let file = menu::Tree::with_children( + cosmic::widget::RcElementWrapper::new(menu::root("File").into()), + menu::items( + kb, + vec![ + menu::Item::Button("Open Image…", None, Action::Open), + menu::Item::Button("Save Image…", None, Action::Save), + menu::Item::Button("Screen Dump…", None, Action::Dump), + menu::Item::Button("Stitch Cylindrical TIFF…", None, Action::Stitch), + menu::Item::Divider, + menu::Item::Button("Load GPX Track…", None, Action::LoadTrack), + menu::Item::Divider, + menu::Item::Button("Quit", None, Action::Quit), + ], + ), + ); + let view = menu::Tree::with_children( + cosmic::widget::RcElementWrapper::new(menu::root("View").into()), + menu::items( + kb, + vec![ + menu::Item::Button("Choose Viewpoint…", None, Action::ChooseViewpoint), + menu::Item::CheckBox( + "Show Hidden", + None, + self.show_hidden, + Action::ShowHidden, + ), + menu::Item::Divider, + menu::Item::Button("Zoom In", None, Action::ZoomIn), + menu::Item::Button("Zoom Out", None, Action::ZoomOut), + menu::Item::Button("Fit", None, Action::Fit), + menu::Item::Divider, + menu::Item::CheckBox( + "Compass strip", + None, + self.compass_style == CompassStyle::Aligned, + Action::CompassAligned, + ), + menu::Item::CheckBox( + "Compass ruler", + None, + self.compass_style == CompassStyle::Ruler, + Action::CompassRuler, + ), + menu::Item::CheckBox( + "Compass off", + None, + self.compass_style == CompassStyle::Off, + Action::CompassOff, + ), + ], + ), + ); + let projection = menu::Tree::with_children( + cosmic::widget::RcElementWrapper::new(menu::root("Projection").into()), + menu::items( + kb, + vec![ + menu::Item::CheckBox( + "Normal", + None, + self.pan.projection_type == ProjectionType::Rectilinear, + Action::ProjNormal, + ), + menu::Item::CheckBox( + "Panoramic", + None, + self.pan.projection_type == ProjectionType::Cylindrical, + Action::ProjPanoramic, + ), + ], + ), + ); + let help = menu::Tree::with_children( + cosmic::widget::RcElementWrapper::new(menu::root("Help").into()), + menu::items(kb, vec![menu::Item::Button("About gipfel", None, Action::About)]), + ); + + vec![ + widget::nav_bar_toggle() + .active(self.core().nav_bar_active()) + .on_toggle(Message::ToggleNav) + .into(), + menu::bar(vec![file, view, projection, help]) + .item_height(menu::ItemHeight::Dynamic(40)) + .item_width(menu::ItemWidth::Uniform(260)) + .spacing(4.0) + .into(), + ] + } + + fn subscription(&self) -> Subscription { + event::listen_with(|event, _status, _id| match event { + Event::Keyboard(keyboard::Event::ModifiersChanged(m)) => { + Some(Message::Event(Event::Keyboard(keyboard::Event::ModifiersChanged(m)))) + } + Event::Keyboard(keyboard::Event::KeyPressed { key, modifiers, .. }) => { + if modifiers.control() { + if let Key::Character(c) = &key { + match c.as_str() { + "o" | "O" => return Some(Message::OpenImage), + "s" | "S" => return Some(Message::SaveImage), + "d" | "D" => return Some(Message::DumpImage), + "t" | "T" => return Some(Message::LoadTrack), + "v" | "V" => return Some(Message::ChooseViewpoint), + "f" | "F" => return Some(Message::Fit), + "q" | "Q" => return Some(Message::Quit), + _ => {} + } + } + } + match key { + Key::Named(cosmic::iced::keyboard::key::Named::Escape) => { + Some(Message::PickerClose) + } + _ => None, + } + } + Event::Window(WinEvent::FileDropped(paths)) => Some(Message::FileDropped(paths)), + Event::Window(WinEvent::Resized(size)) => { + Some(Message::WinResized(size.width as f64, size.height as f64)) + } + _ => None, + }) + } + + fn update(&mut self, message: Message) -> Task { + use cosmic::dialog::file_chooser::FileFilter; + + match message { + Message::Vp(m) => return self.update_vp(m), + Message::Event(Event::Keyboard(keyboard::Event::ModifiersChanged(m))) => { + self.ctrl = m.control(); + } + Message::Event(_) => {} + Message::ToggleNav => { + self.core_mut().nav_bar_toggle(); + self.remeasure(); + } + Message::WinResized(w, h) => { + self.win_w = w; + self.win_h = h; + self.remeasure(); + } + Message::FileDropped(paths) => { + for p in paths { + let ext = p.extension().and_then(|e| e.to_str()).unwrap_or(""); + if matches!(ext.to_ascii_lowercase().as_str(), "jpg" | "jpeg" | "png") { + return self.open_image(p); + } else if ext.eq_ignore_ascii_case("csv") || ext.eq_ignore_ascii_case("dat") { + // treat as track file + let task = async move { + hill::load(&p.to_string_lossy()) + .map(|mut hs| { + for h in hs.iter_mut() { + h.flags |= FLAG_TRACK_POINT; + } + hs + }) + .map_err(|e| e.to_string()) + }; + return cosmic::task::future(async move { + Message::TrackLoaded(Box::new(task.await)) + }); + } + } + } + Message::OpenImage => { + return cosmic::task::future(async move { + let filter = + FileFilter::new("Images").glob("*.jpg").glob("*.jpeg").glob("*.png"); + Message::OpenPath(pick_open("Open Image", filter).await) + }); + } + Message::OpenPath(Some(p)) => return self.open_image(p), + Message::OpenPath(None) => {} + Message::ImageLoaded(res) => match *res { + Ok(li) => { + let title = li.path.file_name().map(|s| s.to_string_lossy().into_owned()).unwrap_or_default(); + self.set_header_title(title); + let (w, h) = (li.w, li.h); + self.project = Some(Project { + path: li.path, + handle: li.handle, + w, + h, + md: li.md.clone(), + }); + let md = li.md; + let have = gipfel_core::metadata::apply_image_metadata(&mut self.pan, &md, w, h); + self.known.clear(); + self.track.clear(); + self.sync_inputs(); + self.fit = true; + self.remeasure(); + self.zoom = self.fit_zoom(); + self.panv = Vector::ZERO; + self.bump(); + // prefill the profile name the auto-match used (C + // get_distortion_profile_name) and refresh the store list + let auto_profile = gipfel_core::metadata::distortion_profile_name(&md); + self.profile_in = auto_profile.clone().unwrap_or_default(); + self.profile_overwrite = false; + self.refresh_profiles(); + self.status = if have { + format!("loaded (gipfel info found)") + } else { + "loaded (no gipfel info: set viewpoint)".into() + }; + if md.k0.is_nan() { + if let Some(name) = + auto_profile.filter(|n| self.profiles.iter().any(|p| p == n)) + { + self.status = format!("{} - distortion profile '{name}'", self.status); + } + } + // the first frame after the image appears can be drawn + // before the texture is fully uploaded; a late nudge + // guarantees the photo shows without user interaction + return Task::perform( + async { tokio::time::sleep(std::time::Duration::from_millis(150)).await }, + |_| cosmic::Action::App(Message::Redraw), + ); + } + Err(e) => self.status = e, + }, + Message::SaveImage => { + let Some(p) = &self.project else { + self.status = "no image loaded".into(); + return Task::none(); + }; + let _ = p; + return cosmic::task::future(async move { + let filter = FileFilter::new("JPEG").glob("*.jpg"); + Message::SavePath(pick_save("Save Image", filter).await) + }); + } + Message::SavePath(out) => { + let (Some(out), Some(p)) = (out, &self.project) else { return Task::none() }; + let md = self.md_from_state(); + let src = p.path.clone(); + return cosmic::task::future(async move { + match md.save_image(&src.to_string_lossy(), &out.to_string_lossy()) { + Ok(()) => Message::ImageSaved(Some(out)), + Err(e) => { + let _ = e; + Message::ImageSaved(None) + } + } + }); + } + Message::ImageSaved(Some(p)) => { + self.status = format!("saved {}", p.display()); + } + Message::ImageSaved(None) => {} + Message::DumpImage => { + let Some(p) = &self.project else { + self.status = "no image loaded".into(); + return Task::none(); + }; + let path = p.path.clone(); + // the dump is at image resolution; screen-pixel overlay sizes + // must be divided by the current zoom to look like the screen + let zoom = self.zoom.max(0.02) as f64; + let markers: Vec = self + .markers() + .iter() + .map(|m| dump::DumpMarker { + x: m.x, + y: m.y, + name: m.name.to_string(), + known: m.known, + hidden: m.hidden, + }) + .collect(); + let track: Vec = self + .track_pts() + .iter() + .map(|t| dump::DumpTrackPt { x: t.x, y: t.y, hidden: t.hidden, width: t.width }) + .collect(); + let font_size = self.label_font; + let font_family = cosmic::config::interface_font().family; + return cosmic::task::future(async move { + let res: Result = async { + let filter = FileFilter::new("Images").glob("*.png").glob("*.jpg"); + let out = + pick_save("Screen Dump", filter).await.ok_or_else(|| "cancelled".to_string())?; + let img = image::open(&path) + .map(|i| i.to_rgba8()) + .map_err(|e| format!("{e}"))?; + dump::render(&img, &markers, &track, &out, font_size, 1.0 / zoom, &font_family)?; + Ok(out) + } + .await; + Message::ImageDumped(Box::new(res)) + }); + } + Message::ImageDumped(res) => { + self.status = match *res { + Ok(p) => format!("dumped {}", p.display()), + Err(e) => e, + }; + } + Message::StitchTiff => { + let Some(p) = &self.project else { + self.status = "no image loaded".into(); + return Task::none(); + }; + // cylindrical output spanning the lens' actual field of view + let edge = self.view_angle_limit(); // half-angle, radians + if !(edge > 0.0 && edge < 1.5) { + self.status = "cannot determine lens field of view".into(); + return Task::none(); + } + let center = self.pan.get_center_angle(); + let (from, to) = (center - edge.to_degrees(), center + edge.to_degrees()); + // output scale = the source's focal scale (pixels per radian + // at image center): the cylindrical image is 1:1 in the middle + // and slightly downsampled towards the edges (where the + // rectilinear source has increasing angular stretch). Using + // the source width instead would upsample ~5-10% and blur. + let radius = self.pan.get_scale(); + let w_out = (2.0 * edge * radius).round_ties_even().max(16.0) as u32; + // cover the source's vertical extent (may be off-center + // with pitch); output rows are symmetric about elevation 0 + let e_top = elev_at_img_y(&self.pan, -(p.h as f64 / 2.0)); + let e_bot = elev_at_img_y(&self.pan, p.h as f64 / 2.0); + let e_max = e_top.abs().max(e_bot.abs()).min(85f64.to_radians()); + let h_out = (2.0 * radius * e_max.tan()).round_ties_even().max(16.0) as u32; + let path = p.path.clone(); + let pan = self.pan.params_snapshot(); + return cosmic::task::future(async move { + let res: Result = async { + let filter = FileFilter::new("TIFF images").glob("*.tif").glob("*.tiff"); + let out = pick_save("Stitch Cylindrical TIFF", filter) + .await + .ok_or_else(|| "cancelled".to_string())?; + let out_for_task = out.clone(); + tokio::task::spawn_blocking(move || -> Result<(), String> { + let mut src = stitch::StitchSource::load(&path.to_string_lossy())?; + // use the live (possibly adjusted) parameters + src.pan = pan; + let img = stitch::resample( + &[src], + stitch::ScanMode::Bicubic, + w_out, + h_out, + from, + to, + ); + stitch::write_tiff(&out_for_task.to_string_lossy(), &img, true) + }) + .await + .map_err(|e| e.to_string())??; + Ok(out) + } + .await; + Message::StitchDone(Box::new(res)) + }); + } + Message::StitchDone(res) => { + self.status = match *res { + Ok(p) => format!("stitched {}", p.display()), + Err(e) => e, + }; + } + Message::LoadTrack => { + return cosmic::task::future(async move { + let res: Result, String> = async { + let filter = FileFilter::new("Track files") + .glob("*.gpx") + .glob("*.csv") + .glob("*.dat"); + let path = pick_open("Load GPX Track", filter) + .await + .ok_or_else(|| "cancelled".to_string())?; + let hs = load_track_file(&path).map_err(|e| e.to_string())?; + Ok(hs) + } + .await; + Message::TrackLoaded(Box::new(res)) + }); + } + Message::TrackLoaded(res) => match *res { + Ok(hills) => { + // remap known hills by name across the index shift caused by + // removing the old track + let names: Vec = self.known.iter().map(|&i| self.pan.hills[i].name.clone()).collect(); + self.pan.remove_hills(FLAG_TRACK_POINT); + let n = hills.len(); + self.pan.add_hills(&hills); + let base = self.pan.hills.len() - n; + self.track = (base..base + n).collect(); + self.known = names + .iter() + .filter_map(|n| self.pan.get_pos(n)) + .collect(); + self.status = format!("track: {n} points"); + self.bump(); + } + Err(e) => self.status = e, + }, + Message::ChooseViewpoint => { + self.picker = Some(Picker { + title: "Choose Viewpoint".into(), + query: String::new(), + mode: PickerMode::Viewpoint, + }); + } + Message::ContextFindPeak(x, y) => { + self.context = None; + self.picker = Some(Picker { + title: "Find Peak".into(), + query: String::new(), + mode: PickerMode::FindPeak { img_x: x, img_y: y }, + }); + } + Message::PickerQuery(q) => { + if let Some(p) = &mut self.picker { + p.query = q; + } + } + Message::PickerPick(i) => { + if let Some(p) = self.picker.take() { + match p.mode { + PickerMode::Viewpoint => { + self.pan.set_viewpoint(Some(i)); + self.sync_inputs(); + self.bump(); + self.status = format!("viewpoint: {}", self.pan.hills[i].name); + } + PickerMode::FindPeak { img_x, img_y } => { + if !self.known.contains(&i) { + self.known.push(i); + } + self.pan.hills[i].flags |= FLAG_VISIBLE; + if !self.pan.visible.contains(&i) { + self.pan.visible.push(i); + } + self.pan.hills[i].x = img_x; + self.pan.hills[i].y = img_y; + self.refit(); + } + } + } + } + Message::PickerClose => { + self.picker = None; + self.about = false; + self.context = None; + } + Message::ContextToggleHidden(i) => { + self.context = None; + self.pan.hills[i].flags ^= FLAG_HIDDEN; + self.bump(); + } + Message::Surface(action) => { + return cosmic::task::message(cosmic::Action::Surface(action)); + } + Message::Redraw => {} + Message::Center(a) => { + self.pan.set_center_angle(a); + self.bump(); + } + Message::CompassStyle(s) => { + self.compass_style = s; + self.remeasure(); + self.bump(); + } + Message::Pitch(a) => { + self.pan.set_pitch_angle(a); + self.bump(); + } + Message::Roll(a) => { + self.pan.set_roll_angle(a); + self.bump(); + } + Message::LabelFont(s) => { + self.label_font = s; + self.bump(); + } + Message::Focal(f) => { + self.set_focal_35mm(f); + self.bump(); + } + Message::Visibility(v) => { + // the slider shows more visibility to the right; the core + // stores a height/distance threshold (higher = fewer hills) + self.pan.set_height_dist_ratio(0.11 - v); + self.bump(); + } + Message::TrackWidth(w) => self.track_width = w, + Message::Lat(s) => self.lat_in = s, + Message::Lon(s) => self.lon_in = s, + Message::Height(s) => self.height_in = s, + Message::ViewApply => { + if let Ok(v) = self.lat_in.parse() { + self.pan.set_view_lat(v); + } + if let Ok(v) = self.lon_in.parse() { + self.pan.set_view_long(v); + } + if let Ok(v) = self.height_in.parse() { + self.pan.set_view_height(v); + } + self.bump(); + } + Message::K0(s) => self.k0_in = s, + Message::K1(s) => self.k1_in = s, + Message::X0(s) => self.x0_in = s, + Message::DistortionApply => { + let k0 = self.k0_in.parse().unwrap_or(0.0); + let k1 = self.k1_in.parse().unwrap_or(0.0); + let x0 = self.x0_in.parse().unwrap_or(0.0); + self.pan.set_distortion_params(k0, k1, x0); + self.bump(); + } + Message::ProfileName(s) => { + self.profile_in = s; + self.profile_overwrite = false; + } + Message::ProfileSave => { + let name = self.profile_in.trim().to_string(); + if name.is_empty() { + self.status = "enter a profile name".into(); + } else { + let p = gipfel_core::profile::Profile { + k0: self.pan.parms.k0, + k1: self.pan.parms.k1, + x0: self.pan.parms.x0, + }; + let exists = self.profiles.iter().any(|n| *n == name); + if exists && !self.profile_overwrite { + self.profile_overwrite = true; + self.status = + format!("profile '{name}' exists - press Save again to overwrite"); + } else { + match gipfel_core::profile::save(&name, p, true) { + Ok(()) => { + self.status = format!("saved profile '{name}'"); + self.profile_overwrite = false; + self.refresh_profiles(); + } + Err(e) => self.status = format!("save failed: {e}"), + } + } + } + } + Message::ProfileLoad => { + let name = self.profile_in.trim().to_string(); + if !name.is_empty() { + self.apply_profile(&name); + } + } + Message::ProfileSelect(i) => { + if let Some(name) = self.profiles.get(i).cloned() { + self.profile_in = name.clone(); + self.apply_profile(&name); + } + } + Message::Projection(t) => { + self.pan.set_projection(t); + self.bump(); + } + Message::ToggleShowHidden => { + self.show_hidden = !self.show_hidden; + self.bump(); + } + Message::FitDistortion(v) => { + self.fit_distortion = v; + // fitting a newly unlocked calibration takes effect at once + if v && self.known.len() >= 4 { + self.refit(); + } + } + Message::Fit => { + self.fit = true; + self.zoom = self.fit_zoom(); + self.panv = Vector::ZERO; + } + Message::ZoomIn => { + return self.update_vp(VpMsg::ZoomAt { factor: 1.25, img: cosmic::iced::Point::ORIGIN }); + } + Message::ZoomOut => { + return self.update_vp(VpMsg::ZoomAt { factor: 0.8, img: cosmic::iced::Point::ORIGIN }); + } + Message::About(b) => self.about = b, + Message::Quit => { + if let Some(id) = self.core.main_window_id() { + return cosmic::iced::window::close(id); + } + } + } + Task::none() + } + + /// The settings panel lives in the framework-managed left nav slot + /// (standard padding, condensed-mode handling, collapse state). + fn nav_bar(&self) -> Option>> { + if !self.core().nav_bar_active() { + return None; + } + let close: Element = widget::button::icon( + widget::icon::from_name("pan-start-symbolic").size(16), + ) + .on_press(Message::ToggleNav) + .into(); + let panel: Element = widget::container( + widget::column::with_children([ + widget::row::with_children([ + widget::text::heading("gipfel").into(), + widget::container(close) + .width(Length::Fill) + .align_x(cosmic::iced::alignment::Horizontal::Right) + .into(), + ]) + .align_y(cosmic::iced::alignment::Vertical::Center) + .spacing(8) + .into(), + widget::scrollable(self.sidebar()).into(), + ]) + .spacing(10), + ) + .padding([10, 12]) + .width(Length::Fixed(340.0)) + .height(Length::Fill) + .class(cosmic::theme::Container::Secondary) + .into(); + Some(panel.map(cosmic::Action::App)) + } + + fn view(&self) -> Element<'_, Message> { + let content: Element = match &self.project { + Some(p) => { + let vp = ViewportWidget:: { + handle: p.handle.clone(), + img_w: p.w as f32, + img_h: p.h as f32, + zoom: self.zoom, + pan: self.panv, + markers: self.markers(), + track: self.track_pts(), + hover: self.hover, + ctrl: self.ctrl, + font_size: self.label_font, + version: self.version, + _msg: PhantomData, + }; + let vp: Element = Element::new(vp); + // standard popover context menu (like cosmic-files): opens in + // its own surface, never shifts the main content + let vp: Element = + cosmic::widget::context_menu(vp, self.context_menu_tree()) + .window_id( + self.core + .main_window_id() + .unwrap_or(cosmic::iced::window::Id::RESERVED), + ) + .on_surface_action(Message::Surface) + .into(); + if self.compass_style == CompassStyle::Off { + vp + } else { + let fov = match self.pan.projection_type { + ProjectionType::Rectilinear => ( + 2.0 * ((p.w as f64 * 0.5) / self.pan.parms.scale).atan() + ) + .to_degrees(), + ProjectionType::Cylindrical => { + (p.w as f64 / self.pan.parms.scale).to_degrees() + } + }; + let strip: Element = Element::new(Compass { + style: self.compass_style, + center_deg: self.pan.get_center_angle(), + zoom: self.zoom, + pan_x: self.panv.x, + parms: self.pan.parms, + ptype: self.pan.projection_type, + fov_deg: fov, + max_view_angle: self.pan.max_view_angle, + on_change: Rc::new(Message::Center), + }); + widget::column::with_children([vp, strip]).into() + } + } + None => widget::container( + widget::column::with_children([ + widget::text::title2("gipfel").into(), + widget::text("Open a JPEG to start (File ▸ Open Image…, Ctrl+O).").into(), + widget::button::suggested("Open Image…").on_press(Message::OpenImage).into(), + ]) + .spacing(12) + .align_x(cosmic::iced::alignment::Horizontal::Center), + ) + .into(), + }; + + let main: Element = content; + + let status_text = self.hover_info.as_deref().unwrap_or(&self.status); + let status: Element = widget::container(widget::text::body(status_text.to_string())) + .padding([4, 10]) + .width(Length::Fill) + .style(|_: &Theme| widget::container::Style { + background: Some(cosmic::iced::Background::Color(Color::from_rgba(0.0, 0.0, 0.0, 0.15))), + ..Default::default() + }) + .into(); + + let body: Element = widget::column::with_children([main, status]).into(); + + if let Some(picker) = &self.picker { + return self.with_dialog(body, self.view_picker(picker)); + } + if self.about { + return self.with_dialog(body, self.view_about()); + } + body + } +} + +// helpers split out to keep view() readable + +/// settings row with the label above a full-width, easy-grip slider +fn slider_item( + label: &str, + value: String, + slider: impl Into>, +) -> Element<'static, Message> { + let header: Element<'static, Message> = widget::row::with_children([ + widget::text::body(label.to_string()).into(), + widget::container(widget::text::body(value)) + .width(Length::Fill) + .align_x(cosmic::iced::alignment::Horizontal::Right) + .into(), + ]) + .into(); + let body: Element<'static, Message> = slider.into(); + widget::column::with_children([header, body]) + .spacing(2) + .into() +} + +/// Async file chooser: returns the chosen path, or None if cancelled. +async fn pick_open( + title: &str, + filter: cosmic::dialog::file_chooser::FileFilter, +) -> Option { + cosmic::dialog::file_chooser::open::Dialog::new() + .title(title.to_string()) + .filter(filter) + .open_file() + .await + .ok()? + .url() + .to_file_path() + .ok() +} + +/// Async save dialog: returns the target path, or None if cancelled. +async fn pick_save( + title: &str, + filter: cosmic::dialog::file_chooser::FileFilter, +) -> Option { + cosmic::dialog::file_chooser::save::Dialog::new() + .title(title.to_string()) + .filter(filter) + .save_file() + .await + .ok()? + .url() + .and_then(|u| u.to_file_path().ok()) +} + +impl Gipfel { + fn with_dialog<'a>(&self, body: Element<'a, Message>, dialog: Element<'a, Message>) -> Element<'a, Message> { + widget::column::with_children([ + widget::container(dialog).width(Length::Fill).into(), + body, + ]) + .into() + } + + fn view_picker<'a>(&self, p: &'a Picker) -> Element<'a, Message> { + let entries = self.picker_entries(p); + let mut list = widget::column::with_capacity(entries.len()); + let mut empty = true; + for (label, idx) in entries { + empty = false; + list = list.push( + widget::button::text(label) + .on_press(Message::PickerPick(idx)) + .width(Length::Fill), + ); + } + if empty { + list = list.push(widget::text("no matches")); + } + + widget::dialog() + .title(p.title.clone()) + .body("Type to search, click to select.") + .control( + widget::text_input("search…", &p.query) + .on_input(Message::PickerQuery), + ) + .control({ + let el: Element = list.into(); + widget::scrollable(el).height(Length::Fixed(320.0)) + }) + .secondary_action(widget::button::text("Cancel").on_press(Message::PickerClose)) + .is_overlay(false) + .into() + } + + fn context_menu_tree(&self) -> Option>> { + let ctx = self.context.as_ref()?; + let mut items: Vec> = Vec::new(); + // placing a marker only makes sense on empty image space, not on an + // existing peak marker + if ctx.hit.is_none() { + items.push(menu::Item::Button( + "Find Peak…".to_string(), + None, + CtxAction::FindPeak(ctx.img_x, ctx.img_y), + )); + } + if let Some(i) = ctx.hit { + let h = &self.pan.hills[i]; + items.push(menu::Item::Button( + format!( + "{} {} ({} m)", + if h.flags & FLAG_HIDDEN != 0 { "Unhide" } else { "Hide" }, + h.name, + h.height as i64 + ), + None, + CtxAction::ToggleHidden(i), + )); + } + if items.is_empty() { + return None; + } + Some(menu::items(&std::collections::HashMap::new(), items)) + } + + fn view_about<'a>(&self) -> Element<'a, Message> { + widget::dialog() + .title("About gipfel") + .body("gipfel 0.5.0 (COSMIC preview)") + .control(widget::text("Annotate mountain panoramas.")) + .control(widget::text("Port of gipfel by Johannes Hofmann - mathematical models unchanged.")) .primary_action(widget::button::suggested("Close").on_press(Message::About(false))) + .is_overlay(false) + .into() + } + + fn sidebar(&self) -> Element<'_, Message> { + let has_img = self.project.is_some(); + let focal = self.focal_35mm().clamp(1.0, 300.0); + let pitch = self.pan.get_pitch_angle().clamp(-40.0, 20.0); + let roll = self.pan.get_roll_angle().clamp(-10.0, 10.0); + let vis = self.pan.height_dist_ratio.clamp(0.01, 0.1); + let vis_ui = 0.11 - vis; // slider: right = more hills visible + let label_font = self.label_font; + let center = self.pan.get_center_angle().rem_euclid(360.0); + + let vp_name = self + .pan + .view_name + .clone() + .unwrap_or_else(|| "manual".to_string()); + + let sec_viewpoint = widget::settings::section() + .title("Viewpoint") + .add(widget::settings::item( + "Position", + widget::row::with_children([ + widget::text(vp_name).into(), + widget::button::text("Choose…").on_press(Message::ChooseViewpoint).into(), + ]) + .spacing(8), + )) + .add(widget::settings::item( + "Latitude", + widget::text_input("°", &self.lat_in) + .on_input(Message::Lat) + .on_submit(|_| Message::ViewApply), + )) + .add(widget::settings::item( + "Longitude", + widget::text_input("°", &self.lon_in) + .on_input(Message::Lon) + .on_submit(|_| Message::ViewApply), + )) + .add(widget::settings::item( + "Height (m)", + widget::text_input("m", &self.height_in) + .on_input(Message::Height) + .on_submit(|_| Message::ViewApply), + )); + + let sec_orientation = widget::settings::section() + .title("Orientation") + .add(widget::settings::item( + "Center (deg)", + widget::text(format!("{center:.1}°")), + )) + .add(slider_item( + "Focal length (35mm)", + format!("{focal:.1} mm"), + widget::slider(1.0..=300.0, focal, Message::Focal) + .step(0.01) + .height(30.0), + )) + .add(slider_item( + "Pitch (deg)", + format!("{pitch:.1}°"), + widget::slider(-40.0..=20.0, pitch, Message::Pitch) + .step(0.01) + .height(30.0), + )) + .add(slider_item( + "Roll (deg)", + format!("{roll:.1}°"), + widget::slider(-10.0..=10.0, roll, Message::Roll) + .step(0.01) + .height(30.0), + )); + + let sec_display = widget::settings::section() + .title("Display") + .add(slider_item( + "Visibility", + format!("{:.0}%", (vis_ui - 0.01) / 0.09 * 100.0), + widget::slider(0.01..=0.1, vis_ui, Message::Visibility) + .step(0.001) + .height(30.0), + )) + .add(slider_item( + "Track width", + format!("{:.0} m at 1 km", self.track_width), + widget::slider(1.0..=500.0, self.track_width, Message::TrackWidth) + .step(1.0) + .height(30.0), + )) + .add(slider_item( + "Label size", + format!("{label_font:.0} px"), + widget::slider(8.0f32..=28.0, label_font, Message::LabelFont) + .step(1.0) + .height(30.0), + )) + .add(widget::settings::item( + "Show hidden", + widget::toggler(self.show_hidden).on_toggle(|_| Message::ToggleShowHidden), + )) + .add(widget::row::with_children([ + widget::button::text("Zoom −").on_press(Message::ZoomOut).into(), + widget::button::text("Fit").on_press(Message::Fit).into(), + widget::button::text("Zoom +").on_press(Message::ZoomIn).into(), + widget::text(format!("{:.0}%", self.zoom * 100.0)).into(), + ]) + .spacing(8)); + + let sec_projection = widget::settings::section().title("Projection").add( + widget::row::with_children([ + widget::radio( + "Normal", + ProjectionType::Rectilinear, + Some(self.pan.projection_type), + Message::Projection, + ) + .into(), + widget::radio( + "Panoramic", + ProjectionType::Cylindrical, + Some(self.pan.projection_type), + Message::Projection, + ) + .into(), + ]) + .spacing(16), + ); + + // distortion fitting is locked unless enough markers are placed + let can_fit = self.known.len() >= 4; + let fit_cb = widget::checkbox(self.fit_distortion); + let fit_cb = if can_fit { + fit_cb.on_toggle(Message::FitDistortion) + } else { + fit_cb.on_toggle_maybe(None as Option Message>) + }; + let mut fit_col = vec![fit_cb.into()]; + if !can_fit { + fit_col.push(widget::text("position 4 or more markers to enable").size(11).into()); + } + + // fixed-width label column: every control starts at the same x and + // therefore shares the same width (settings::item sizes labels to + // their text, which made the profile field narrower than k0/k1/x0) + let mut sec_distortion = widget::settings::section() + .title("Distortion") + .add(widget::settings::item_row(vec![ + widget::text("Fit distortion").width(110).into(), + widget::column::with_children(fit_col).spacing(2).into(), + ])) + .add(widget::settings::item_row(vec![ + widget::text("k0").width(110).into(), + widget::text_input("0", &self.k0_in) + .on_input(Message::K0) + .on_submit(|_| Message::DistortionApply) + .into(), + ])) + .add(widget::settings::item_row(vec![ + widget::text("k1").width(110).into(), + widget::text_input("0", &self.k1_in) + .on_input(Message::K1) + .on_submit(|_| Message::DistortionApply) + .into(), + ])) + .add(widget::settings::item_row(vec![ + widget::text("x0").width(110).into(), + widget::text_input("0", &self.x0_in) + .on_input(Message::X0) + .on_submit(|_| Message::DistortionApply) + .into(), + ])) + .add(widget::settings::item_row(vec![ + widget::text("profile").width(110).into(), + widget::text_input("Camera_Model_50.00_mm", &self.profile_in) + .on_input(Message::ProfileName) + .on_submit(|_| Message::ProfileSave) + .into(), + ])) + .add( + widget::row::with_children([ + widget::button::suggested("Save profile") + .on_press(Message::ProfileSave) + .into(), + widget::button::text("Load profile") + .on_press(Message::ProfileLoad) + .into(), + ]) + .spacing(8), + ); + if !self.profiles.is_empty() { + // stored profiles (shared with the C version); selecting one + // loads its parameters + let win = self + .core + .main_window_id() + .unwrap_or(cosmic::iced::window::Id::RESERVED); + sec_distortion = sec_distortion.add(widget::settings::item_row(vec![ + widget::text("stored").width(110).into(), + widget::dropdown::popup_dropdown( + self.profiles.clone(), + None, + Message::ProfileSelect, + win, + Message::Surface, + |m| m, + ) + .into(), + ])); + } + + let mut col = widget::column::with_children([ + sec_viewpoint.into(), + sec_orientation.into(), + sec_projection.into(), + sec_distortion.into(), + sec_display.into(), + ]) + .spacing(12) + .padding(12); + + if !has_img { + col = col.push(widget::text("Load an image to enable the controls.")); + } + col.into() + } +} + +fn main() -> cosmic::iced::Result { + let mut flags = Flags { data: None, image: None }; + let mut args = std::env::args().skip(1); + while let Some(a) = args.next() { + match a.as_str() { + "--data" => flags.data = args.next().map(PathBuf::from), + s if !s.starts_with('-') && flags.image.is_none() => flags.image = Some(PathBuf::from(s)), + _ => {} + } + } + cosmic::app::run::( + cosmic::app::Settings::default().size(cosmic::iced::Size::new(1280.0, 800.0)), + flags, + ) +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn elev_inverse_matches_projection() { + let mut pan = Panorama::new(); + pan.set_projection(ProjectionType::Rectilinear); + pan.set_center_angle(0.0); + pan.set_pitch_angle(0.0); + pan.set_scale(200.0); + // undistorted rectilinear: y = scale * tan(pitch - elev) + for &y in &[-100.0, -30.0, 0.0, 45.0, 120.0] { + let e = elev_at_img_y(&pan, y); + let (_, y_back) = pan.project_point(0.0, e).unwrap(); + assert!((y_back - y).abs() < 0.01, "y={y} -> elev={e} -> {y_back}"); + } + // with pitch the band is asymmetric but still invertible + pan.set_pitch_angle(10.0); + let p = pan.parms.a_pitch; + let top = elev_at_img_y(&pan, -100.0); + let bot = elev_at_img_y(&pan, 100.0); + assert!(top > p && bot < p, "top={top} p={p} bot={bot}"); + let (_, y_back) = pan.project_point(0.0, top).unwrap(); + assert!((y_back + 100.0).abs() < 0.01); + } +} diff --git a/crates/gipfel-app/src/viewport.rs b/crates/gipfel-app/src/viewport.rs new file mode 100644 index 0000000..3001fc4 --- /dev/null +++ b/crates/gipfel-app/src/viewport.rs @@ -0,0 +1,659 @@ +//! The gipfel viewport: a custom widget that draws the photo, hill markers, +//! decluttered labels, GPS track, and handles all canvas interaction. +//! +//! Marker/label placement mirrors the original GipfelWidget::draw/set_labels +//! (image-space coordinates; the screen transform is applied at draw time). + +use std::rc::Rc; +use std::time::Instant; + +use cosmic::iced::advanced::image::{FilterMethod, Handle, Image as CoreImage}; +use cosmic::iced::advanced::renderer::Style as RendererStyle; +use cosmic::iced::advanced::text::{ + self as adv_text, paragraph::Plain, LineHeight, Renderer as _, +}; +use cosmic::iced::advanced::widget::{tree, Tree, Widget as WidgetTrait}; +use cosmic::iced::advanced::{Clipboard, Layout, Shell}; +use cosmic::iced::mouse::{self, Button, Cursor, Interaction, ScrollDelta}; +use cosmic::iced::widget::canvas::{Frame, LineCap, Path, Stroke}; +use cosmic::iced::advanced::graphics::geometry::Style as GeoStyle; +use cosmic::iced::{ + alignment, Color, Event, Length, Pixels, Point, Rectangle, Radians, Size, Vector, +}; + +use cosmic::Renderer as R; +use cosmic::Theme; + +pub const DEFAULT_LABEL_FONT: f32 = 12.0; +const CROSS_SIZE: f64 = 4.0; + +#[derive(Clone, Debug)] +pub struct Marker { + /// image-space coordinates, origin at image center (like C's m->x/m->y) + pub x: f64, + pub y: f64, + pub name: Rc, + pub idx: usize, + pub known: bool, + pub hidden: bool, +} + +#[derive(Clone, Debug)] +pub struct TrackPt { + pub x: f64, + pub y: f64, + pub hidden: bool, + pub width: f64, +} + +#[derive(Debug, Clone)] +pub enum VpMsg { + /// zoom by `factor`, keeping image point `img` under the cursor + ZoomAt { factor: f32, img: Point }, + PanBy { delta: Vector }, + PanTo { pan: Vector }, + MoveMarker { idx: usize, img: Point }, + DragEnd { had_marker: bool }, + /// hill index of the toggled marker (None = clicked empty space) + ToggleKnown(Option), + ContextAt { img: Point, hit: Option }, + /// hill index of the hovered marker + Hover(Option), +} + +#[derive(Clone)] +struct LabelLayout { + version: u64, + font_size: f32, + zoom: f32, + entries: Vec<(f32, f32)>, +} + +#[derive(Default)] +struct DragState { + drag: Drag, + last_click: Option<(Instant, Point)>, + label_cache: std::cell::RefCell>, + /// Paragraph primitives must outlive the frame: the renderer stores only a + /// Weak reference, so keep strong refs here until the next draw. + paragraphs: std::cell::RefCell::Paragraph>>>, +} + +#[derive(Default, Clone)] +enum Drag { + #[default] + None, + Marker(usize), + Pan { start: Point, pan0: Vector }, +} + +pub struct ViewportWidget + Send + 'static> { + pub handle: Handle, + pub img_w: f32, + pub img_h: f32, + pub zoom: f32, + pub pan: Vector, + /// hill markers in visible order; label layout depends on this order + pub markers: Rc>, + pub track: Rc>, + pub hover: Option, + pub ctrl: bool, + /// label font size in pixels (UI-tunable) + pub font_size: f32, + /// bump whenever marker positions change to invalidate the label cache + pub version: u64, + pub _msg: std::marker::PhantomData M>, +} + +impl + Send + 'static> Clone for ViewportWidget { + fn clone(&self) -> Self { + ViewportWidget { + handle: self.handle.clone(), + img_w: self.img_w, + img_h: self.img_h, + zoom: self.zoom, + pan: self.pan, + markers: self.markers.clone(), + track: self.track.clone(), + hover: self.hover, + ctrl: self.ctrl, + font_size: self.font_size, + version: self.version, + _msg: std::marker::PhantomData, + } + } +} + +impl + Send + 'static> ViewportWidget { + /// line height: paragraph metrics (cosmic LineHeight default = 1.3x) plus + /// a 2px gap, so stacked labels never touch + fn line_h(&self) -> f32 { + self.font_size * 1.3 + 2.0 + } + + fn center(bounds: Rectangle, pan: Vector) -> Point { + Point::new( + bounds.x + bounds.width / 2.0 + pan.x, + bounds.y + bounds.height / 2.0 + pan.y, + ) + } + + fn to_img(&self, center: Point, p: Point) -> Point { + Point::new((p.x - center.x) / self.zoom, (p.y - center.y) / self.zoom) + } + + /// hit-test in image space; tolerance scales with zoom (~6 screen px). + /// Returns the *hill index* of the hit marker. + fn hit(&self, center: Point, p: Point) -> Option { + let img = self.to_img(center, p); + let tol = (6.0 / self.zoom).clamp(2.0, 20.0); + self.markers + .iter() + .find(|m| (m.x as f32 - img.x).abs() < tol && (m.y as f32 - img.y).abs() < tol) + .map(|m| m.idx) + } +} + +impl + Send + 'static> WidgetTrait for ViewportWidget { + fn size(&self) -> Size { + Size::new(Length::Fill, Length::Fill) + } + + fn tag(&self) -> tree::Tag { + tree::Tag::of::() + } + + fn state(&self) -> tree::State { + tree::State::new(DragState::default()) + } + + fn layout( + &mut self, + _tree: &mut Tree, + _renderer: &R, + limits: &cosmic::iced::advanced::layout::Limits, + ) -> cosmic::iced::advanced::layout::Node { + cosmic::iced::advanced::layout::Node::new(limits.max()) + } + + fn draw( + &self, + tree: &Tree, + renderer: &mut R, + _theme: &Theme, + _style: &RendererStyle, + layout: Layout<'_>, + _cursor: Cursor, + _viewport: &Rectangle, + ) { + use cosmic::iced::advanced::image::Renderer as _; + use cosmic::iced::advanced::Renderer as _; + + let bounds = layout.bounds(); + let center = Self::center(bounds, self.pan); + let z = self.zoom; + + // 1. photo + let dest = Rectangle::new( + Point::new(center.x - self.img_w * z / 2.0, center.y - self.img_h * z / 2.0), + Size::new(self.img_w * z, self.img_h * z), + ); + renderer.draw_image( + CoreImage { + handle: self.handle.clone(), + filter_method: FilterMethod::Linear, + rotation: Radians(0.0), + border_radius: 0.0.into(), + opacity: 1.0, + snap: false, + }, + dest, + bounds, + ); + + let layout_entries = self.label_layout(tree, renderer); + + // The wgpu renderer draws geometry before images within one layer; + // record the overlay in a fresh layer so it lands on top of the photo. + renderer.start_layer(bounds); + self.draw_overlay(renderer, tree, bounds, center, &layout_entries); + renderer.end_layer(); + } + + fn update( + &mut self, + tree: &mut Tree, + event: &Event, + layout: Layout<'_>, + cursor: Cursor, + _renderer: &R, + _clipboard: &mut dyn Clipboard, + shell: &mut Shell<'_, M>, + _viewport: &Rectangle, + ) { + let bounds = layout.bounds(); + let center = Self::center(bounds, self.pan); + let state = tree.state.downcast_mut::(); + + match event { + Event::Mouse(mouse::Event::WheelScrolled { delta }) => { + let Some(p) = cursor.position_over(bounds) else { return }; + let (dx, dy) = match delta { + ScrollDelta::Pixels { x, y } => (*x, *y), + ScrollDelta::Lines { x, y } => (*x * 30.0, *y * 30.0), + }; + if self.ctrl { + shell.publish(M::from(VpMsg::ZoomAt { + factor: (1.0 + dy * 0.0015).clamp(0.5, 2.0), + img: self.to_img(center, p), + })); + } else { + shell.publish(M::from(VpMsg::PanBy { delta: Vector::new(-dx, -dy) })); + } + shell.request_redraw(); + } + Event::Mouse(mouse::Event::ButtonPressed(Button::Left)) => { + let Some(p) = cursor.position_over(bounds) else { return }; + let now = Instant::now(); + let dbl = state + .last_click + .is_some_and(|(t, q)| { + t.elapsed().as_millis() < 400 + && (q.x - p.x).abs() < 6.0 + && (q.y - p.y).abs() < 6.0 + }); + state.last_click = Some((now, p)); + + if dbl { + state.drag = Drag::None; + state.last_click = None; + shell.publish(M::from(VpMsg::ToggleKnown(self.hit(center, p)))); + return; + } + + if let Some(idx) = self.hit(center, p) { + if self.markers.iter().any(|m| m.idx == idx && m.known) { + state.drag = Drag::Marker(idx); + return; + } + } + state.drag = Drag::Pan { start: p, pan0: self.pan }; + } + Event::Mouse(mouse::Event::ButtonReleased(Button::Left)) => { + let had_marker = matches!(state.drag, Drag::Marker(_)); + state.drag = Drag::None; + shell.publish(M::from(VpMsg::DragEnd { had_marker })); + } + Event::Mouse(mouse::Event::ButtonPressed(Button::Middle)) => { + if let Some(p) = cursor.position_over(bounds) { + shell.publish(M::from(VpMsg::ToggleKnown(self.hit(center, p)))); + } + } + Event::Mouse(mouse::Event::ButtonPressed(Button::Right)) => { + if let Some(p) = cursor.position_over(bounds) { + shell.publish(M::from(VpMsg::ContextAt { + img: self.to_img(center, p), + hit: self.hit(center, p), + })); + } + } + Event::Mouse(mouse::Event::CursorMoved { .. }) => { + match state.drag.clone() { + Drag::Marker(idx) => { + if let Some(p) = cursor.position() { + shell.publish(M::from(VpMsg::MoveMarker { + idx, + img: self.to_img(center, p), + })); + shell.request_redraw(); + } + } + Drag::Pan { start, pan0 } => { + if let Some(p) = cursor.position() { + shell.publish(M::from(VpMsg::PanTo { + pan: pan0 + Vector::new(p.x - start.x, p.y - start.y), + })); + shell.request_redraw(); + } + } + Drag::None => { + match cursor.position_over(bounds) { + Some(p) => shell.publish(M::from(VpMsg::Hover(self.hit(center, p)))), + // cursor moved off the viewport (e.g. into the + // sidebar): clear any stale hover info + None => shell.publish(M::from(VpMsg::Hover(None))), + } + } + } + } + _ => {} + } + } + + fn mouse_interaction( + &self, + _tree: &Tree, + layout: Layout<'_>, + cursor: Cursor, + _viewport: &Rectangle, + _renderer: &R, + ) -> Interaction { + let center = Self::center(layout.bounds(), self.pan); + match cursor.position_over(layout.bounds()) { + Some(p) if self.hit(center, p).is_some() => Interaction::Pointer, + Some(_) => Interaction::Grab, + None => Interaction::None, + } + } +} + +impl + Send + 'static> ViewportWidget { + fn draw_overlay( + &self, + renderer: &mut R, + tree: &Tree, + bounds: Rectangle, + center: Point, + layout_entries: &[(f32, f32)], + ) { + use cosmic::iced::advanced::graphics::geometry::Renderer as _; + use cosmic::iced::advanced::text::Renderer as _; + + let z = self.zoom; + let to_screen = |p: Point| Point::new(center.x + p.x * z, center.y + p.y * z); + let layout_entries = layout_entries; + let mut frame = Frame::with_bounds(renderer, bounds); + + // 2. GPS track + for w in self.track.windows(2) { + let (a, b) = (&w[0], &w[1]); + let pa = to_screen(Point::new(a.x as f32, a.y as f32)); + let pb = to_screen(Point::new(b.x as f32, b.y as f32)); + let color = if a.hidden { + Color::from_rgb(0.2, 0.4, 1.0) + } else { + Color::from_rgb(1.0, 0.2, 0.2) + }; + let stroke = Stroke { + width: (a.width as f32 * z).max(1.0), + style: GeoStyle::Solid(color), + line_cap: LineCap::Round, + ..Default::default() + }; + let path = Path::new(|bld| { + bld.move_to(pa); + bld.line_to(pb); + }); + frame.stroke(&path, stroke); + } + + // 3. hill markers: crosses + leaders + // C's draw() skips markers outside the image (fabs(x) > w/2 || ...); + // must match label_layout's filter, otherwise off-image markers get + // drawn (in the letterbox area) without ever being laid out. + let in_img = |m: &Marker| { + m.x.abs() * 2.0 <= self.img_w as f64 && m.y.abs() * 2.0 <= self.img_h as f64 + }; + let yellow = Color::from_rgb(1.0, 0.9, 0.2); + for (i, m) in self.markers.iter().enumerate() { + if !in_img(m) { + continue; + } + let p = to_screen(Point::new(m.x as f32, m.y as f32)); + if !bounds.contains(p) { + continue; + } + let (lx, ly) = layout_entries[i]; + let d = CROSS_SIZE as f32; + let label_line = p.y + ly - self.line_h(); + let stroke = Stroke { + width: 1.0, + style: GeoStyle::Solid(yellow), + ..Default::default() + }; + let path = Path::new(|bld| { + bld.move_to(Point::new(p.x - d, p.y)); + bld.line_to(Point::new(p.x + d, p.y)); + bld.move_to(Point::new(p.x, p.y)); + bld.line_to(Point::new(p.x, label_line)); + bld.move_to(Point::new(p.x, label_line)); + bld.line_to(Point::new(p.x + lx, label_line)); + }); + frame.stroke(&path, stroke); + } + + // known-hill flags (pole + pennant + circle) + let flag_color = Color::from_rgb(1.0, 0.25, 0.25); + for m in self.markers.iter().filter(|m| m.known) { + if !in_img(m) { + continue; + } + let p = to_screen(Point::new(m.x as f32, m.y as f32)); + if !bounds.contains(p) { + continue; + } + let stroke = Stroke { + width: 1.5, + style: GeoStyle::Solid(flag_color), + ..Default::default() + }; + let pole = Path::new(|bld| { + bld.move_to(p); + bld.line_to(Point::new(p.x, p.y - 16.0)); + bld.circle(p, 3.0); + }); + let pennant = Path::new(|bld| { + bld.move_to(Point::new(p.x, p.y - 16.0)); + bld.line_to(Point::new(p.x + 10.0, p.y - 13.0)); + bld.line_to(Point::new(p.x, p.y - 10.0)); + bld.close(); + }); + frame.stroke(&pole, stroke.clone()); + frame.fill(&pennant, flag_color); + } + + renderer.draw_geometry(frame.into_geometry()); + + // 4. labels (screen space, constant font size) + let font = renderer.default_font(); + let size = Pixels(self.font_size); + // the renderer keeps only a Weak ref to each paragraph, so hold strong + // refs in the widget state until the next frame replaces them + let state = tree.state.downcast_ref::(); + { + let mut paras = state.paragraphs.borrow_mut(); + paras.clear(); + paras.reserve(self.markers.len()); + for m in self.markers.iter() { + paras.push(self.make_paragraph(font, size, &m.name)); + } + } + let paras = state.paragraphs.borrow(); + + for (i, m) in self.markers.iter().enumerate() { + if !in_img(m) { + continue; + } + let p = to_screen(Point::new(m.x as f32, m.y as f32)); + if !bounds.contains(p) { + continue; + } + let (_, ly) = layout_entries[i]; + let pos = Point::new(p.x + 2.0, p.y + ly - self.line_h()); + let is_hover = self.hover == Some(m.idx); + + let para = ¶s[i]; + let mb = para.min_bounds(); + + let mut bg = Frame::with_bounds(renderer, bounds); + let rect = Rectangle::new( + Point::new(p.x + 1.0, pos.y - 1.0), + Size::new(mb.width + 2.0, mb.height + 1.0), + ); + let path = Path::new(|bld| bld.rectangle(rect.position(), rect.size())); + bg.fill( + &path, + if is_hover { + Color::from_rgba(1.0, 0.9, 0.2, 0.95) + } else { + Color::from_rgba(0.0, 0.0, 0.0, 0.55) + }, + ); + renderer.draw_geometry(bg.into_geometry()); + + let color = if is_hover { + Color::BLACK + } else if m.hidden { + Color::from_rgb(0.4, 0.6, 1.0) + } else { + Color::WHITE + }; + renderer.fill_paragraph( + para.raw(), + Point::new(p.x + 2.0, pos.y), + color, + bounds, + ); + } + } + + fn make_paragraph(&self, font: ::Font, size: Pixels, text: &str) -> Plain<::Paragraph> { + Plain::<::Paragraph>::new(adv_text::Text { + content: text.to_string(), + bounds: Size::new(600.0, 40.0), + size, + line_height: LineHeight::default(), + font, + align_x: alignment::Horizontal::Left.into(), + align_y: alignment::Vertical::Top.into(), + shaping: adv_text::Shaping::Basic, + wrapping: adv_text::Wrapping::None, + ellipsize: adv_text::Ellipsize::None, + }) + } + + /// C set_labels port: place labels above crosses, pushing colliding ones up. + /// Markers outside the image are skipped entirely (C: fabs(x) > w/2 filter), + /// so off-image peaks never push visible labels up. + fn label_layout(&self, tree: &Tree, renderer: &mut R) -> Vec<(f32, f32)> { + let state = tree.state.downcast_ref::(); + if let Some(l) = state.label_cache.borrow().as_ref() { + if l.version == self.version + && l.font_size == self.font_size + && l.zoom == self.zoom + && l.entries.len() == self.markers.len() + { + return l.entries.clone(); + } + } + + let font = renderer.default_font(); + let size = Pixels(self.font_size); + let line_h = self.line_h(); + + let widths: Vec = self + .markers + .iter() + .map(|m| self.make_paragraph(font, size, &m.name).min_bounds().width) + .collect(); + + // C: fabs(m->x) > img->w()/2 || fabs(m->y) > img->h()/2 -> skip + // (marker coords are center-relative, like C's m->x/m->y) + let in_img: Vec = self + .markers + .iter() + .map(|m| m.x.abs() * 2.0 <= self.img_w as f64 && m.y.abs() * 2.0 <= self.img_h as f64) + .collect(); + + // Labels are drawn at constant size in *screen* space, so collisions + // depend on the on-screen distance between markers: run the layout on + // zoom-scaled coordinates (C had no zoom, image space == screen space). + let z = self.zoom; + let sx: Vec = self.markers.iter().map(|m| m.x as f32 * z).collect(); + let sy: Vec = self.markers.iter().map(|m| m.y as f32 * z).collect(); + + let entries = layout_label_y(&sx, &sy, &widths, &in_img, line_h); + *state.label_cache.borrow_mut() = Some(LabelLayout { + version: self.version, + font_size: self.font_size, + zoom: self.zoom, + entries: entries.clone(), + }); + entries + } +} + +/// Shared label collision-avoidance layout, used for both the on-screen +/// viewport and the screen dump so the dump shows exactly the same label +/// placement (plan.md M3: same code path). +/// +/// `xs`/`ys` are marker positions (center-relative) in one uniform pixel +/// space, `widths` the label widths and `line_h` the line height, all in +/// that same space. Markers must be ordered by x (the sweep breaks on the +/// first non-overlapping predecessor). Returns `(label_width + 1, y_offset)` +/// per marker; the label is drawn at `(marker.x + 2, marker.y + y_offset - +/// line_h)`. +pub fn layout_label_y( + xs: &[f32], + ys: &[f32], + widths: &[f32], + in_img: &[bool], + line_h: f32, +) -> Vec<(f32, f32)> { + let n = xs.len(); + let mut label_x = vec![0.0f32; n]; + let mut label_y = vec![0.0f32; n]; + let cross = CROSS_SIZE as f32; + let overlap = |a1: f32, l1: f32, a2: f32, l2: f32| a1 <= a2 + l2 && a1 + l1 >= a2; + + for i in 0..n { + if !in_img[i] { + continue; + } + label_x[i] = widths[i] + 1.0; + label_y[i] = 0.0; + + let mut colliding: Vec = Vec::new(); + for j in (0..i).rev() { + if !in_img[j] { + continue; + } + if overlap(xs[i], label_x[i], xs[j] - cross, label_x[j]) { + colliding.push(j); + } else { + break; + } + } + colliding.sort_by(|&a, &b| { + let ya = ys[a] + label_y[a]; + let yb = ys[b] + label_y[b]; + yb.total_cmp(&ya) + }); + // C does a single overwrite pass here, which can leave overlaps + // when a push moves the label into a band that was already + // checked. Iterate to a fixed point instead (converges quickly: + // each round only moves the label up, bounded by the stack). + for _round in 0..64 { + let mut pushed = false; + for &j in &colliding { + if overlap( + ys[i] + label_y[i] - line_h, + line_h, + ys[j] + label_y[j] - line_h, + line_h, + ) || overlap(ys[i] + label_y[i] - line_h, line_h, ys[j] - cross, 2.0 * cross) { + let ny = (ys[j] + label_y[j] - ys[i] - line_h - 2.0).round(); + if ny < label_y[i] { + label_y[i] = ny; + pushed = true; + } + } + } + if !pushed { + break; + } + } + } + + (0..n).map(|i| (label_x[i], label_y[i])).collect() +} diff --git a/crates/gipfel-cli/Cargo.toml b/crates/gipfel-cli/Cargo.toml new file mode 100644 index 0000000..9f0f96d --- /dev/null +++ b/crates/gipfel-cli/Cargo.toml @@ -0,0 +1,13 @@ +[package] +name = "gipfel-cli" +version = "0.1.0" +edition = "2021" +description = "gipfel command line tools (export-position, export-hills)" +license = "GPL-2.0-or-later" + +[[bin]] +name = "gipfel-cli" +path = "src/main.rs" + +[dependencies] +gipfel-core = { path = "../gipfel-core" } diff --git a/crates/gipfel-cli/src/main.rs b/crates/gipfel-cli/src/main.rs new file mode 100644 index 0000000..f0690a8 --- /dev/null +++ b/crates/gipfel-cli/src/main.rs @@ -0,0 +1,371 @@ +//! gipfel-cli: headless command line tools, mirroring the C binary's +//! `-p/--export-position`, `-e/--export-hills` and stitch modes. + +use gipfel_core::hill::FLAG_EXPORT; +use gipfel_core::metadata::{apply_image_metadata, image_dimensions, ImageMetaData}; +use gipfel_core::panorama::Panorama; +use gipfel_core::stitch::{self, ScanMode}; + +fn usage() { + eprintln!( + "usage: gipfel-cli [options]\n\ + \n\ + commands:\n\ + \x20 export-position print position CSV line\n\ + \x20 export-hills [-d ] [-e ] export hill positions\n\ + \x20 info show parsed metadata\n\ + \x20 stitch [-t out.tif [-4]] [-j out.jpg] [-w N]\n\ + \x20 [-h N] [-r from,to] [-b] img... cylindrical panorama\n\ + \x20 profiles list distortion profiles\n\ + \x20 profile-save save a profile\n" + ); +} + +fn export_position(img: &str) -> i32 { + // C reads metadata directly (no image decode); missing values print as nan + let mut md = ImageMetaData::new(); + md.load_image(img); + + // C: printf(",%s,,%f,%f,%d\n", img_file, latitude, longitude, (int)rint(height)) + println!( + ",{},,{},{},{}", + img, + fmt_f(md.latitude), + fmt_f(md.longitude), + md.height.round_ties_even() as i64 + ); + 0 +} + +fn fmt_f(v: f64) -> String { + if v.is_nan() { + "nan".to_string() + } else { + format!("{:.6}", v) + } +} + +fn export_hills(img: &str, data: &str, export_file: Option<&str>) -> i32 { + let Some((w, h)) = image_dimensions(img) else { + eprintln!("could not read image {}", img); + return 1; + }; + + let mut md = ImageMetaData::new(); + md.load_image(img); + + let mut pan = Panorama::new(); + let have_info = apply_image_metadata(&mut pan, &md, w, h); + + if let Err(e) = pan.load_data(data) { + eprintln!("Could not load datafile {}: {}", data, e); + return 1; + } + pan.set_height_dist_ratio(0.07); // C default "visibility" + + if !have_info { + eprintln!("No gipfel info available for {}.", img); + return 0; + } + + // -e: restrict export to the hills from the given file + if let Some(file) = export_file { + match gipfel_core::hill::load(file) { + Ok(mut extra) => { + for m in extra.iter_mut() { + m.flags |= FLAG_EXPORT; + } + pan.add_hills(&extra); + } + Err(_) => return 1, + } + } + + println!("#\n# name\theight\tx\ty\tdistance\tflags\n#"); + + for row in pan.export_hill_rows(w as i64, h as i64, export_file.is_some()) { + println!("{}", row); + } + + if export_file.is_some() { + for m in pan.hills.iter_mut() { + m.flags &= !FLAG_EXPORT; + } + } + + 0 +} + +/// Port of the C stitch mode: resample the given images into one +/// cylindrical panorama. Flags follow the C binary: -w/-h output size, +/// -r from,to azimuth range in degrees, -b bicubic, -4 16-bit TIFF, +/// -j/-t output file (tiff wins if both given, as in C). +fn stitch(imgs: &[String], opts: &StitchOpts) -> i32 { + if imgs.is_empty() || (opts.tiff.is_none() && opts.jpeg.is_none()) { + usage(); + return 2; + } + let mut sources = Vec::new(); + for img in imgs { + match stitch::StitchSource::load(img) { + Ok(s) => sources.push(s), + Err(e) => { + eprintln!("{e}"); + return 1; + } + } + } + let mode = if opts.bicubic { ScanMode::Bicubic } else { ScanMode::Nearest }; + let out = stitch::resample(&sources, mode, opts.w, opts.h, opts.from, opts.to); + let hit = out.hit.iter().filter(|h| **h).count(); + eprintln!( + "stitched {}x{}: {}/{} pixels covered", + out.w, out.h, hit, out.hit.len() + ); + let res = if let Some(path) = &opts.tiff { + stitch::write_tiff(path, &out, opts.bits16) + } else { + stitch::write_jpeg(opts.jpeg.as_ref().unwrap(), &out, 90) + }; + match res { + Ok(()) => 0, + Err(e) => { + eprintln!("{e}"); + 1 + } + } +} + +struct StitchOpts { + w: u32, + h: u32, + from: f64, + to: f64, + bicubic: bool, + bits16: bool, + tiff: Option, + jpeg: Option, +} + +/// C defaults: stitch_w=2000, stitch_h=500, stitch_from=0.0, stitch_to=380.0 +fn stitch_defaults() -> StitchOpts { + StitchOpts { + w: 2000, + h: 500, + from: 0.0, + to: 380.0, + bicubic: false, + bits16: false, + tiff: None, + jpeg: None, + } +} + +fn stitch_cmd(args: &[String]) -> i32 { + let mut o = stitch_defaults(); + let mut imgs = Vec::new(); + let mut i = 0; + while i < args.len() { + let a = args[i].as_str(); + let need = |i: &mut usize| -> Option { + if *i + 1 < args.len() { + *i += 1; + Some(args[*i].clone()) + } else { + None + } + }; + let ok = match a { + "-w" => match need(&mut i) { + Some(v) => { + o.w = v.parse().unwrap_or(2000); + true + } + None => false, + }, + "-h" => match need(&mut i) { + Some(v) => { + o.h = v.parse().unwrap_or(500); + true + } + None => false, + }, + "-r" => match need(&mut i) { + Some(v) => { + // C: atof(optarg) [, atof(after comma)] — atof reads the + // numeric prefix before the comma + let (first, second) = v.split_once(',').unwrap_or((v.as_str(), "")); + o.from = first.trim().parse().unwrap_or(0.0); + if !second.trim().is_empty() { + o.to = second.trim().parse().unwrap_or(380.0); + } + true + } + None => false, + }, + "-j" => match need(&mut i) { + Some(v) => { + o.jpeg = Some(v); + true + } + None => false, + }, + "-t" => match need(&mut i) { + Some(v) => { + o.tiff = Some(v); + true + } + None => false, + }, + "-b" => { + o.bicubic = true; + true + } + "-4" => { + o.bits16 = true; + true + } + _ => { + if a.starts_with('-') { + false + } else { + imgs.push(a.to_string()); + true + } + } + }; + if !ok { + usage(); + return 2; + } + i += 1; + } + stitch(&imgs, &o) +} + +fn info(img: &str) -> i32 { + let mut md = ImageMetaData::new(); + md.load_image(img); + println!("manufacturer: {:?}", md.manufacturer); + println!("model: {:?}", md.model); + println!("longitude: {}", md.longitude); + println!("latitude: {}", md.latitude); + println!("height: {}", md.height); + println!("direction: {}", md.direction); + println!("pitch: {}", md.pitch); + println!("roll: {}", md.roll); + println!("focal_length: {}", md.focal_length); + println!("focal_length_35mm: {}", md.focal_length_35mm); + println!("projection type: {}", md.projection_type); + println!("k0: {}", md.k0); + println!("k1: {}", md.k1); + println!("x0: {}", md.x0); + 0 +} + +fn profiles() -> i32 { + let all = gipfel_core::profile::list(); + if all.is_empty() { + eprintln!("no profiles in {}", gipfel_core::profile::prefs_path().display()); + } + for (name, p) in all { + println!("{name}\tk0={}\tk1={}\tx0={}", p.k0, p.k1, p.x0); + } + 0 +} + +fn profile_save(name: &str, k0: f64, k1: f64, x0: f64) -> i32 { + let p = gipfel_core::profile::Profile { k0, k1, x0 }; + match gipfel_core::profile::save(name, p, true) { + Ok(()) => { + eprintln!("saved profile {name}"); + 0 + } + Err(e) => { + eprintln!("{e}"); + 1 + } + } +} + +fn profile_save_cmd(args: &[String]) -> i32 { + if args.len() != 4 { + usage(); + return 2; + } + let (Ok(k0), Ok(k1), Ok(x0)) = ( + args[1].parse::(), + args[2].parse::(), + args[3].parse::(), + ) else { + usage(); + return 2; + }; + profile_save(&args[0], k0, k1, x0) +} + +fn main() { + let args: Vec = std::env::args().collect(); + let Some(cmd) = args.get(1) else { + usage(); + std::process::exit(2); + }; + + let code = match cmd.as_str() { + "export-position" => { + if args.len() != 3 { + usage(); + 2 + } else { + export_position(&args[2]) + } + } + "export-hills" => { + if args.len() < 3 { + usage(); + 2 + } else { + let img = &args[2]; + let mut data = "gipfel.dat".to_string(); + let mut export_file: Option = None; + let mut i = 3; + while i < args.len() { + match args[i].as_str() { + "-d" => { + if i + 1 < args.len() { + data = args[i + 1].clone(); + i += 1; + } + } + "-e" => { + if i + 1 < args.len() { + export_file = Some(args[i + 1].clone()); + i += 1; + } + } + _ => {} + } + i += 1; + } + export_hills(img, &data, export_file.as_deref()) + } + } + "info" => { + if args.len() != 3 { + usage(); + 2 + } else { + info(&args[2]) + } + } + "profiles" => profiles(), + "profile-save" => profile_save_cmd(&args[2..]), + "stitch" => stitch_cmd(&args[2..]), + _ => { + usage(); + 2 + } + }; + + std::process::exit(code); +} diff --git a/crates/gipfel-core/Cargo.toml b/crates/gipfel-core/Cargo.toml new file mode 100644 index 0000000..d8ddfd9 --- /dev/null +++ b/crates/gipfel-core/Cargo.toml @@ -0,0 +1,13 @@ +[package] +name = "gipfel-core" +version = "0.1.0" +edition = "2021" +description = "Port of gipfel's computational core (hills, panorama, projections, least-squares fit, metadata)" +license = "GPL-2.0-or-later" + +[dependencies] +kamadak-exif.workspace = true +quick-xml.workspace = true +image.workspace = true +tiff.workspace = true +dirs.workspace = true diff --git a/crates/gipfel-core/examples/genjpg.rs b/crates/gipfel-core/examples/genjpg.rs new file mode 100644 index 0000000..ae12b94 --- /dev/null +++ b/crates/gipfel-core/examples/genjpg.rs @@ -0,0 +1,24 @@ +// creates a test JPEG with a gipfel comment: genjpg +use gipfel_core::metadata::ImageMetaData; + +fn main() { + let out = std::env::args().nth(1).unwrap(); + let (w, h) = (1600u32, 900u32); + let img = image::DynamicImage::new_rgb8(w, h).grayscale(); + img.save(&out).unwrap(); + + let mut md = ImageMetaData::new(); + md.longitude = 11.15; + md.latitude = 47.36; + md.height = 2100.0; + md.direction = 20.0; + md.pitch = 0.0; + md.roll = 1.0; + md.focal_length_35mm = 35.0; + md.projection_type = 0; + md.k0 = 0.0; + md.k1 = 0.0; + md.x0 = 0.0; + md.save_image(&out, &out).unwrap(); + println!("wrote {}", out); +} diff --git a/crates/gipfel-core/src/hill.rs b/crates/gipfel-core/src/hill.rs new file mode 100644 index 0000000..a904948 --- /dev/null +++ b/crates/gipfel-core/src/hill.rs @@ -0,0 +1,224 @@ +//! Hill database: data structures, CSV loading, duplicate detection. +//! +//! Faithful port of gipfel's `Hill.cxx`. + +use std::fs::File; +use std::io::{BufRead, BufReader}; + +pub const FLAG_DUPLIC: u32 = 1 << 0; +pub const FLAG_KNOWN: u32 = 1 << 1; +pub const FLAG_HIDDEN: u32 = 1 << 2; +pub const FLAG_VISIBLE: u32 = 1 << 3; +pub const FLAG_TRACK_POINT: u32 = 1 << 4; +pub const FLAG_EXPORT: u32 = 1 << 5; + +pub const PI: f64 = std::f64::consts::PI; +pub const DEG2RAD: f64 = PI / 180.0; + +#[derive(Clone, Debug)] +pub struct Hill { + pub name: String, + /// latitude, radians + pub phi: f64, + /// longitude, radians + pub lam: f64, + /// height, meters + pub height: f64, + /// azimuth from viewpoint, radians + pub alph: f64, + /// peak angle over horizon, radians + pub a_elev: f64, + /// distance from viewpoint, fraction of earth radius + pub dist: f64, + pub x: f64, + pub y: f64, + pub label_x: f64, + pub label_y: f64, + pub flags: u32, +} + +impl Hill { + pub fn new(name: &str, phi: f64, lam: f64, height: f64) -> Self { + Hill { + name: name.to_string(), + phi, + lam, + height, + alph: 0.0, + a_elev: 0.0, + dist: 0.0, + x: 0.0, + y: 0.0, + label_x: 0.0, + label_y: 0.0, + flags: 0, + } + } + + pub fn is_track_point(&self) -> bool { + self.flags & FLAG_TRACK_POINT != 0 + } +} + +/// `atof()` semantics: parse the longest numeric prefix, 0.0 on failure. +fn atof(s: &str) -> f64 { + let s = s.trim(); + if s.is_empty() { + return 0.0; + } + if let Ok(v) = s.parse::() { + return v; + } + for i in (1..s.len()).rev() { + if let Ok(v) = s[..i].parse::() { + return v; + } + } + 0.0 +} + +/// Load a hill database. Two line formats are accepted: +/// - `SHORT,NAME,DESC,LAT,LON,HEIGHT` (6 fields) +/// - `LAT,LON,HEIGHT` (3 fields, track point) +pub fn load(path: &str) -> std::io::Result> { + let file = File::open(path)?; + let mut hills = Vec::new(); + + // the database is Latin-1; decode to Unicode code points + for line in BufReader::new(file).split(b'\n') { + let raw = line?; + let line: String = raw.iter().map(|b| *b as char).collect(); + // strsep() semantics: split into at most 10 fields, keeping empties + let vals: Vec<&str> = line.split(',').take(10).collect(); + let n = vals.len(); + + if n == 6 { + // standard format including name and description + let phi = atof(vals[3]) * DEG2RAD; + let lam = atof(vals[4]) * DEG2RAD; + let height = atof(vals[5]); + hills.push(Hill::new(vals[1], phi, lam, height)); + } else if n == 3 { + // track point format + let phi = atof(vals[0]) * DEG2RAD; + let lam = atof(vals[1]) * DEG2RAD; + let height = atof(vals[2]); + let mut h = Hill::new("", phi, lam, height); + // NOTE: original does not set TRACK_POINT here + h.flags |= 0; + hills.push(h); + } + } + + Ok(hills) +} + +/// Mark hills that are (near-)exact duplicates of an earlier hill. +/// Faithful port of `Hills::mark_duplicates`. +/// +/// Unlike the C original (which sorts the pointer array in place), this +/// must NOT reorder `hills`: in this port, hill indices are stable handles +/// (track lists, known-hill sets, visible lists). Sorting an index +/// permutation instead keeps the semantics identical while preserving +/// insertion order. +pub fn mark_duplicates(hills: &mut [Hill], dist: f64) { + // descending phi (mirrors the C qsort comparator order) + let mut order: Vec = (0..hills.len()).collect(); + order.sort_by(|&a, &b| { + hills[b] + .phi + .partial_cmp(&hills[a].phi) + .unwrap_or(std::cmp::Ordering::Equal) + }); + + let len = hills.len(); + for pos in 0..len { + let i = order[pos]; + if hills[i].is_track_point() { + continue; + } + + let mut j = pos + 1; + while j < len && (hills[order[j]].phi - hills[i].phi).abs() <= dist { + let jj = order[j]; + if hills[jj].flags & FLAG_DUPLIC == 0 + && (hills[jj].lam - hills[i].lam).abs() <= dist + && (hills[jj].height - hills[i].height).abs() <= 50.0 + { + hills[jj].flags |= FLAG_DUPLIC; + } + j += 1; + } + } +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn atof_semantics() { + assert_eq!(atof("12.5"), 12.5); + assert_eq!(atof(""), 0.0); + assert_eq!(atof("abc"), 0.0); + assert_eq!(atof("12abc"), 12.0); + assert_eq!(atof(" -3.5x"), -3.5); + } + + #[test] + fn load_formats() { + let dir = std::env::temp_dir().join("gipfel-core-test"); + std::fs::create_dir_all(&dir).unwrap(); + let path = dir.join("hills.csv"); + std::fs::write( + &path, + "# comment line\n\ + ZSP,Zugspitze,Germany,47.421,10.9855,2962\n\ + FOO,Foo Peak,,47.5,11.0,1000\n\ + 47.6,11.1,1234\n\ + bad,line\n", + ) + .unwrap(); + + let hills = load(path.to_str().unwrap()).unwrap(); + assert_eq!(hills.len(), 3); + assert_eq!(hills[0].name, "Zugspitze"); + assert!((hills[0].phi - 47.421 * DEG2RAD).abs() < 1e-12); + assert_eq!(hills[0].height, 2962.0); + assert_eq!(hills[2].name, ""); + let _ = std::fs::remove_file(&path); + } + + #[test] + fn duplicates() { + let mut hills = vec![ + Hill::new("A", 47.0, 11.0, 1000.0), + Hill::new("A2", 47.000001, 11.000001, 1010.0), // dup of A + Hill::new("B", 47.5, 11.5, 1000.0), + Hill::new("C", 47.0, 11.0, 1200.0), // >50m diff -> not dup + ]; + mark_duplicates(&mut hills, 0.00001); + assert_eq!(hills.iter().filter(|h| h.flags & FLAG_DUPLIC != 0).count(), 1); + } + + #[test] + fn mark_duplicates_preserves_order() { + // indices are stable handles in the port: mark_duplicates must not + // reorder the vec (regression: an in-place phi sort broke the app's + // track index range when add_hills ran it on the full hill list) + let mut hills = vec![ + Hill::new("south", 46.0, 11.0, 1000.0), + Hill::new("north", 48.0, 11.0, 1000.0), + Hill::new("mid", 47.0, 11.0, 1000.0), + Hill::new("mid-dup", 47.000001, 11.0, 1010.0), + ]; + mark_duplicates(&mut hills, 0.00001); + assert_eq!(hills[0].name, "south"); + assert_eq!(hills[1].name, "north"); + assert_eq!(hills[2].name, "mid"); + assert_eq!(hills[3].name, "mid-dup"); + // descending phi visits mid-dup first, so 'mid' is the one marked + assert_eq!(hills[2].flags & FLAG_DUPLIC, FLAG_DUPLIC); + assert_eq!(hills[3].flags & FLAG_DUPLIC, 0); + } +} diff --git a/crates/gipfel-core/src/lib.rs b/crates/gipfel-core/src/lib.rs new file mode 100644 index 0000000..2e2fac9 --- /dev/null +++ b/crates/gipfel-core/src/lib.rs @@ -0,0 +1,23 @@ +//! gipfel-core: headless port of the gipfel computational core. +//! +//! Modules: +//! - [`hill`]: hill database (CSV load, duplicate marking) +//! - [`panorama`]: viewpoint, geodesy, refraction, occlusion, visibility +//! - [`projection`]: rectilinear/cylindrical projections with symbolic +//! derivatives generated by maxima (see `codegen/`) +//! - [`lsq`]: Levenberg-Marquardt fit of the view parameters (system GSL) +//! - [`metadata`]: `gipfel:` JPEG comment + EXIF read/write + +pub mod hill; +pub mod lsq; +pub mod metadata; +pub mod panorama; +pub mod profile; +pub mod projection; +pub mod stitch; +pub mod track; + +pub use hill::Hill; +pub use metadata::ImageMetaData; +pub use panorama::Panorama; +pub use projection::{Projection, ProjectionType, ViewParams}; diff --git a/crates/gipfel-core/src/lsq.rs b/crates/gipfel-core/src/lsq.rs new file mode 100644 index 0000000..6a9d6c7 --- /dev/null +++ b/crates/gipfel-core/src/lsq.rs @@ -0,0 +1,392 @@ +//! Non-linear least squares fit of the view parameters. +//! +//! Faithful port of `ProjectionLSQ.cxx`. The C++ original uses GSL's +//! `gsl_multifit_fdfsolver_lmsder`; this is a self-contained Levenberg- +//! Marquardt with the same scaled damping (`JᵀJ + λ·diag(JᵀJ)`) and the +//! analytic Jacobian produced by the maxima-generated derivatives. For the +//! well-conditioned gipfel fit (residuals reach zero at the optimum) it +//! converges to the same minimum; the golden tests verify agreement with the +//! C/GSL output. + +use crate::hill::Hill; +use crate::projection::{Projection, ViewParams}; + +const MAX_ITER: usize = 200; + +/// Number of fitted parameters at a given level. +fn np_for_level(level: usize) -> usize { + match level { + 0 => 2, + 1 => 4, + _ => 7, + } +} + +struct FitData<'a> { + proj: &'a dyn Projection, + hills: &'a [Hill], + known: &'a [usize], + level: usize, + /// parameters at solver start; entries above `level` keep these values + parms0: ViewParams, +} + +impl FitData<'_> { + fn params_from_x(&self, x: &[f64]) -> ViewParams { + // C starts from the initial params and only overwrites the entries + // that are actually fitted at this level. + let mut p = self.parms0; + p.a_center = x[0]; + p.a_pitch = x[1]; + if self.level > 0 { + p.a_roll = x[2]; + p.scale = x[3]; + } + if self.level > 1 { + p.k0 = x[4]; + p.k1 = x[5]; + p.x0 = x[6]; + } + p + } + + fn args(&self, p: &ViewParams, i: usize) -> crate::projection::ProjArgs { + let m = &self.hills[self.known[i]]; + crate::projection::ProjArgs::from_params(p, m.alph, m.a_elev) + } + + /// Residual vector: 2 entries (x,y) per known hill, length `2*known.len()`. + fn residuals(&self, x: &[f64], f: &mut [f64]) { + let p = self.params_from_x(x); + for i in 0..self.known.len() { + let m = &self.hills[self.known[i]]; + let a = self.args(&p, i); + f[2 * i] = self.proj.mac_x(a) - m.x; + f[2 * i + 1] = self.proj.mac_y(a) - m.y; + } + } + + /// Jacobian, row-major `m`×`np` (`m = 2*known.len()`). + fn jacobian(&self, x: &[f64], j: &mut [f64]) { + let p = self.params_from_x(x); + let np = np_for_level(self.level); + let row = |r: usize, c: usize, v: f64, j: &mut [f64]| j[r * np + c] = v; + for i in 0..self.known.len() { + let a = self.args(&p, i); + + row(2 * i, 0, self.proj.mac_x_dc_view(a), j); + row(2 * i + 1, 0, self.proj.mac_y_dc_view(a), j); + row(2 * i, 1, self.proj.mac_x_dc_pitch(a), j); + row(2 * i + 1, 1, self.proj.mac_y_dc_pitch(a), j); + + if np > 2 { + row(2 * i, 2, self.proj.mac_x_dc_roll(a), j); + row(2 * i + 1, 2, self.proj.mac_y_dc_roll(a), j); + row(2 * i, 3, self.proj.mac_x_dscale(a), j); + row(2 * i + 1, 3, self.proj.mac_y_dscale(a), j); + } + + if np > 4 { + row(2 * i, 4, self.proj.mac_x_dk0(a), j); + row(2 * i + 1, 4, self.proj.mac_y_dk0(a), j); + row(2 * i, 5, self.proj.mac_x_dk1(a), j); + row(2 * i + 1, 5, self.proj.mac_y_dk1(a), j); + row(2 * i, 6, self.proj.mac_x_dx0(a), j); + row(2 * i + 1, 6, self.proj.mac_y_dx0(a), j); + } + } + } +} + +fn dot(a: &[f64], b: &[f64]) -> f64 { + a.iter().zip(b).map(|(x, y)| x * y).sum() +} + +/// Solve `A·δ = b` in place (A row-major `n`×`n`, overwritten) by Gaussian +/// elimination with partial pivoting. Returns `None` if the matrix is +/// (numerically) singular. +fn gauss_solve(a: &mut [f64], b: &[f64], n: usize) -> Option> { + let mut x = b.to_vec(); + for col in 0..n { + // partial pivot + let mut piv = col; + let mut maxv = a[col * n + col].abs(); + for r in (col + 1)..n { + let v = a[r * n + col].abs(); + if v > maxv { + maxv = v; + piv = r; + } + } + if maxv < 1e-300 { + return None; + } + if piv != col { + for c in 0..n { + a.swap(col * n + c, piv * n + c); + } + x.swap(col, piv); + } + let diag = a[col * n + col]; + for r in (col + 1)..n { + let factor = a[r * n + col] / diag; + a[r * n + col] = 0.0; + for c in (col + 1)..n { + a[r * n + c] -= factor * a[col * n + c]; + } + x[r] -= factor * x[col]; + } + } + for r in (0..n).rev() { + let mut s = x[r]; + for c in (r + 1)..n { + s -= a[r * n + c] * x[c]; + } + x[r] = s / a[r * n + r]; + } + Some(x) +} + +/// One Levenberg-Marquardt fit at a single parameter level. On return, `x` +/// holds the fitted parameters (length `np`). +fn lm(data: &FitData, x: &mut [f64], np: usize) { + let m = 2 * data.known.len(); + let mut f = vec![0.0f64; m]; + let mut fnew = vec![0.0f64; m]; + let mut jac = vec![0.0f64; m * np]; + let mut a = vec![0.0f64; np * np]; + let mut ad = vec![0.0f64; np * np]; + let mut g = vec![0.0f64; np]; + let mut neg_g = vec![0.0f64; np]; + let mut xnew = vec![0.0f64; np]; + + data.residuals(x, &mut f); + let mut cost = 0.5 * dot(&f, &f); + let mut lambda: f64 = 1e-3; + + for _ in 0..MAX_ITER { + if !cost.is_finite() { + break; + } + data.jacobian(x, &mut jac); + + // A = JᵀJ (symmetric), g = Jᵀf + for i in 0..np { + for k in i..np { + let mut s = 0.0; + for r in 0..m { + s += jac[r * np + i] * jac[r * np + k]; + } + a[i * np + k] = s; + a[k * np + i] = s; + } + let mut s = 0.0; + for r in 0..m { + s += jac[r * np + i] * f[r]; + } + g[i] = s; + } + + // Solve the damped normal equations (JᵀJ + λ·diag(JᵀJ))·δ = -Jᵀf, + // inflating λ until the system is non-singular. + let mut delta: Option> = None; + for _ in 0..50 { + ad.copy_from_slice(&a); + for i in 0..np { + let d = a[i * np + i]; + ad[i * np + i] = if d > 1e-300 { d * (1.0 + lambda) } else { lambda.max(1e-300) }; + } + for i in 0..np { + neg_g[i] = -g[i]; + } + let mut ad_try = ad.clone(); + if let Some(d) = gauss_solve(&mut ad_try, &neg_g, np) { + delta = Some(d); + break; + } + lambda *= 10.0; + } + let delta = match delta { + Some(d) => d, + None => break, + }; + + for i in 0..np { + xnew[i] = x[i] + delta[i]; + } + data.residuals(&xnew, &mut fnew); + let cost_new = 0.5 * dot(&fnew, &fnew); + + if cost_new < cost { + let improvement = cost - cost_new; + x.copy_from_slice(&xnew); + std::mem::swap(&mut f, &mut fnew); + cost = cost_new; + lambda = (lambda * 0.1).max(1e-12); + // Converged when the cost can no longer improve meaningfully. + if improvement <= 1e-18 * cost { + break; + } + } else { + lambda *= 10.0; + if lambda > 1e16 { + break; + } + } + } +} + +/// Run the Levenberg-Marquardt fit for one parameter level and write the +/// result back into `parms`. Mirrors `ProjectionLSQ::lsq`. +fn lsq(proj: &dyn Projection, hills: &[Hill], known: &[usize], parms: &mut ViewParams, level: usize) { + let np = np_for_level(level); + let data = FitData { + proj, + hills, + known, + level, + parms0: *parms, + }; + + let mut x = vec![ + parms.a_center, + parms.a_pitch, + parms.a_roll, + parms.scale, + parms.k0, + parms.k1, + parms.x0, + ]; + x.truncate(np); + + lm(&data, &mut x, np); + + parms.a_center = x[0]; + parms.a_pitch = x[1]; + if level > 0 { + parms.a_roll = x[2]; + parms.scale = x[3]; + } + if level > 1 { + parms.k0 = x[4]; + parms.k1 = x[5]; + parms.x0 = x[6]; + } +} + +/// Fit the view parameters to the known hills. +/// Mirrors `ProjectionLSQ::comp_params`. `known` indexes into `hills`. +/// With `fit_distortion == false` the lens distortion parameters +/// (k0/k1/x0) are held at their current values and only the viewpoint +/// parameters are fitted, even for 4+ known hills. +pub fn comp_params( + proj: &dyn Projection, + hills: &[Hill], + known: &[usize], + parms: &mut ViewParams, + fit_distortion: bool, +) -> Result<(), String> { + let n = known.len(); + if n < 1 { + return Err("Please position at least 1 hill".to_string()); + } + + let m1 = &hills[known[0]]; + let mut new_parms = *parms; + new_parms.a_center = m1.alph; + new_parms.a_pitch = 0.0; + + if n > 1 { + let m2 = &hills[known[1]]; + let scale_tmp = (m2.x - m1.x).abs() / (m2.alph - m1.alph).abs(); + if scale_tmp.is_nan() || scale_tmp < 50.0 { + return Err(format!("Could not determine initial scale value ({scale_tmp})")); + } + new_parms.scale = scale_tmp; + new_parms.a_roll = 0.0; + } + + if n > 3 && fit_distortion { + new_parms.k0 = 0.0; + new_parms.k1 = 0.0; + new_parms.x0 = 0.0; + eprintln!("Performing calibration"); + } + + let level = if n == 1 { 0 } else { 1 }; + lsq(proj, hills, known, &mut new_parms, level); + if n >= 4 && fit_distortion { + lsq(proj, hills, known, &mut new_parms, 2); + } + + if new_parms.scale.is_nan() || new_parms.scale < 50.0 { + return Err("Could not determine reasonable view parameters".to_string()); + } + + *parms = new_parms; + Ok(()) +} + +#[cfg(test)] +mod tests { + use super::*; + use crate::hill::Hill; + use crate::projection::{Projection, ProjectionType, ViewParams, projection_from_type}; + + /// Five observed hills projected with known parameters (incl. distortion). + fn fixture() -> (Box, Vec, Vec, ViewParams) { + let proj = projection_from_type(ProjectionType::Rectilinear); + let truth = ViewParams { + a_center: 0.3, + scale: 2000.0, + a_pitch: 0.05, + a_roll: 0.0, + k0: -0.0012, + k1: 2e-5, + x0: 0.0, + }; + let mut hills = Vec::new(); + let mut known = Vec::new(); + for i in 0..5 { + let alph = 0.15 + i as f64 * 0.07; + let elev = 0.02 + i as f64 * 0.015; + let (x, y) = proj.get_coordinates(alph, elev, &truth); + let mut h = Hill::new(&format!("h{i}"), 0.0, 0.0, 1000.0); + h.alph = alph; + h.a_elev = elev; + h.x = x; + h.y = y; + hills.push(h); + known.push(i); + } + (proj, hills, known, truth) + } + + #[test] + fn locked_distortion_stays_untouched() { + let (proj, hills, known, truth) = fixture(); + let mut p = truth; + p.a_center -= 0.01; + p.scale *= 1.02; + comp_params(&*proj, &hills, &known, &mut p, false).unwrap(); + // 5 known hills, but the distortion params must not move + assert_eq!(p.k0, truth.k0); + assert_eq!(p.k1, truth.k1); + assert_eq!(p.x0, truth.x0); + // the viewpoint params did fit + assert!((p.a_center - truth.a_center).abs() < 1e-6); + } + + #[test] + fn unlocked_distortion_fit_recovers() { + let (proj, hills, known, truth) = fixture(); + let mut p = truth; + p.k0 = 0.0; + p.k1 = 0.0; + p.a_center -= 0.01; + p.scale *= 1.02; + comp_params(&*proj, &hills, &known, &mut p, true).unwrap(); + assert!((p.k0 - truth.k0).abs() < 1e-5, "k0={}", p.k0); + assert!((p.k1 - truth.k1).abs() < 1e-7, "k1={}", p.k1); + } +} diff --git a/crates/gipfel-core/src/metadata.rs b/crates/gipfel-core/src/metadata.rs new file mode 100644 index 0000000..e2a888f --- /dev/null +++ b/crates/gipfel-core/src/metadata.rs @@ -0,0 +1,697 @@ +//! Image metadata: the `gipfel:` JPEG comment (read/write, byte-compatible +//! with the C version) plus EXIF fallback (GPS position, focal length). +//! +//! Faithful port of gipfel's `ImageMetaData.cxx`. + +use std::fs::File; +use std::io::{BufReader, Write}; + +use exif::{In, Tag, Value}; + +pub const NAN: f64 = f64::NAN; + +#[derive(Clone, Debug)] +pub struct ImageMetaData { + pub manufacturer: Option, + pub model: Option, + pub longitude: f64, + pub latitude: f64, + pub height: f64, + pub direction: f64, + /// camera pitch, degrees (written as the "nick" key in the gipfel: + /// comment for backwards compatibility with the C original) + pub pitch: f64, + /// camera roll, degrees (written as the "tilt" key for compatibility) + pub roll: f64, + pub k0: f64, + pub k1: f64, + pub x0: f64, + pub focal_length: f64, + pub focal_length_35mm: f64, + pub projection_type: i32, +} + +impl Default for ImageMetaData { + fn default() -> Self { + ImageMetaData { + manufacturer: None, + model: None, + longitude: NAN, + latitude: NAN, + height: NAN, + direction: NAN, + pitch: NAN, + roll: NAN, + k0: NAN, + k1: NAN, + x0: NAN, + focal_length: NAN, + focal_length_35mm: NAN, + projection_type: 0, + } + } +} + +impl ImageMetaData { + pub fn new() -> Self { + Self::default() + } + + /// Read JPEG comment first, then fill remaining gaps from EXIF. + pub fn load_image(&mut self, path: &str) { + *self = ImageMetaData::new(); + if let Some(com) = read_jpeg_comment(path) { + self.parse_gipfel_comment(&com); + } + self.load_image_exif(path); + } + + /// Write the gipfel comment into a copy of `in_img`, saved as `out_img`. + /// Writes atomically via a temp file, like the C original. + pub fn save_image(&self, in_img: &str, out_img: &str) -> std::io::Result<()> { + let comment = self.format_comment(); + let tmp = format!("{}.tmp", out_img); + write_jpeg_comment(in_img, &tmp, comment.as_bytes())?; + std::fs::rename(&tmp, out_img)?; + Ok(()) + } + + pub fn format_comment(&self) -> String { + format!( + "gipfel: longitude {}, latitude {}, height {}, direction {}, nick {}, tilt {}, focal_length_35mm {}, projection type {}, k0 {}, k1 {}, x0 {}", + fmt_f64(self.longitude), + fmt_f64(self.latitude), + fmt_f64(self.height), + fmt_f64(self.direction), + fmt_f64(self.pitch), + fmt_f64(self.roll), + fmt_f64(self.focal_length_35mm), + self.projection_type, + fmt_f64(self.k0), + fmt_f64(self.k1), + fmt_f64(self.x0), + ) + } + + /// sscanf(gipfel_format_scan, ...) equivalent; returns true when at + /// least 8 fields were matched (the C requirement). + pub fn parse_gipfel_comment(&mut self, com: &[u8]) -> bool { + let s = String::from_utf8_lossy(com); + let mut sc = Scanner { s: s.as_ref() }; + + if !sc.lit("gipfel: longitude ") { + return false; + } + let Some(lo) = sc.f64() else { return false }; + if !sc.lit(", latitude ") { return false; } + let la = sc.f64(); + if la.is_none() { return false; } + if !sc.lit(", height ") { return false; } + let he = sc.f64(); + if he.is_none() { return false; } + if !sc.lit(", direction ") { return false; } + let dir = sc.f64(); + if dir.is_none() { return false; } + if !sc.lit(", nick ") { return false; } + let pit = sc.f64(); + if pit.is_none() { return false; } + if !sc.lit(", tilt ") { return false; } + let rol = sc.f64(); + if rol.is_none() { return false; } + if !sc.lit(", focal_length_35mm ") { return false; } + let fr = sc.f64(); + if fr.is_none() { return false; } + if !sc.lit(", projection type ") { return false; } + let pt = sc.i32(); + + // 8 matches required: through projection type + let Some(pt) = pt else { return false }; + + self.longitude = lo; + self.latitude = la.unwrap(); + self.height = he.unwrap(); + self.direction = dir.unwrap(); + self.pitch = pit.unwrap(); + self.roll = rol.unwrap(); + self.focal_length_35mm = fr.unwrap(); + self.projection_type = pt; + + // n >= 10: k0 and k1 were matched; x0 defaults to 0.0 when only 10 + if let Some(k0) = sc.f64_after(", k0 ") { + if let Some(k1) = sc.f64_after(", k1 ") { + self.k0 = k0; + self.k1 = k1; + self.x0 = sc.f64_after(", x0 ").unwrap_or(0.0); + } + } + + true + } + + fn load_image_exif(&mut self, path: &str) { + let Ok(file) = File::open(path) else { return }; + let mut reader = BufReader::new(file); + let Ok(exif) = exif::Reader::new().read_from_container(&mut reader) else { + return; + }; + + if self.manufacturer.is_none() { + if let Some(f) = exif.get_field(Tag::Make, In::PRIMARY) { + let s = field_string(f); + if !s.is_empty() { + self.manufacturer = Some(s); + } + } + } + if self.model.is_none() { + if let Some(f) = exif.get_field(Tag::Model, In::PRIMARY) { + let s = field_string(f); + if !s.is_empty() { + self.model = Some(s); + } + } + } + + if self.focal_length.is_nan() { + exif_set_value(self, &exif, Tag::FocalLength, In::PRIMARY, |md, v| md.focal_length = v); + } + if self.focal_length_35mm.is_nan() { + exif_set_value( + self, + &exif, + Tag::FocalLengthIn35mmFilm, + In::PRIMARY, + |md, v| md.focal_length_35mm = v, + ); + } + if self.longitude.is_nan() { + exif_set_coordinate(self, &exif, Tag::GPSLongitude, In::PRIMARY, |md, v| md.longitude = v); + if !self.longitude.is_nan() && exif_ref_is(&exif, Tag::GPSLongitudeRef, "W") { + self.longitude = -self.longitude; + } + } + if self.latitude.is_nan() { + exif_set_coordinate(self, &exif, Tag::GPSLatitude, In::PRIMARY, |md, v| md.latitude = v); + if !self.latitude.is_nan() && exif_ref_is(&exif, Tag::GPSLatitudeRef, "S") { + self.latitude = -self.latitude; + } + } + if self.height.is_nan() { + exif_set_value(self, &exif, Tag::GPSAltitude, In::PRIMARY, |md, v| md.height = v); + // GPSAltitudeRef: 0 = above sea level, 1 = below + if !self.height.is_nan() { + if let Some(f) = exif.get_field(Tag::GPSAltitudeRef, In::PRIMARY) { + if let Value::Byte(b) = &f.value { + if b.first() == Some(&1) { + self.height = -self.height; + } + } + } + } + } + } +} + +/// True if an EXIF ASCII reference tag (GPSLatitudeRef etc.) equals `val`. +fn exif_ref_is(exif: &exif::Exif, tag: Tag, val: &str) -> bool { + exif.get_field(tag, In::PRIMARY) + .and_then(|f| match &f.value { + Value::Ascii(a) => a.first().map(|v| String::from_utf8_lossy(v).trim() == val), + _ => None, + }) + .unwrap_or(false) +} + +fn fmt_f64(v: f64) -> String { + if v.is_nan() { + // C printf("%f", nan) prints "nan"; keep it parseable by both + "nan".to_string() + } else { + format!("{:.6}", v) + } +} + +fn field_string(f: &exif::Field) -> String { + let join = |v: &[String]| v.join(" "); + match &f.value { + Value::Ascii(a) => a + .first() + .map(|v| String::from_utf8_lossy(v).trim().to_string()) + .unwrap_or_default(), + Value::Rational(r) => join(r.iter().map(|v| v.to_string()).collect::>().as_slice()), + Value::Short(v) => join(v.iter().map(|v| v.to_string()).collect::>().as_slice()), + Value::Long(v) => join(v.iter().map(|v| v.to_string()).collect::>().as_slice()), + Value::SShort(v) => join(v.iter().map(|v| v.to_string()).collect::>().as_slice()), + Value::SLong(v) => join(v.iter().map(|v| v.to_string()).collect::>().as_slice()), + _ => String::new(), + } +} + +/// exifSetValue: assign when the tag exists and its value is >= 0. +fn exif_set_value(md: &mut ImageMetaData, exif: &exif::Exif, tag: Tag, inp: In, set: F) +where + F: FnOnce(&mut ImageMetaData, f64), +{ + if let Some(f) = exif.get_field(tag, inp) { + if let Some(v) = field_tofloat(f, 0) { + if v >= 0.0 { + set(md, v); + } + } + } +} + +/// exifSetCoordinate: degrees + minutes/60 + seconds/3600, each component +/// only used when >= 0 (faithful to the C original, incl. ignoring N/S refs). +fn exif_set_coordinate(md: &mut ImageMetaData, exif: &exif::Exif, tag: Tag, inp: In, set: F) +where + F: FnOnce(&mut ImageMetaData, f64), +{ + let Some(f) = exif.get_field(tag, inp) else { + return; + }; + + let mut dest: Option = None; + for i in 0..3 { + if let Some(v) = field_tofloat(f, i) { + if v >= 0.0 { + dest = Some(match i { + 0 => v, + 1 => dest.unwrap_or(0.0) + v / 60.0, + _ => dest.unwrap_or(0.0) + v / 3600.0, + }); + } + } + } + if let Some(v) = dest { + set(md, v); + } +} + +fn field_tofloat(f: &exif::Field, idx: usize) -> Option { + match &f.value { + Value::Rational(r) => r.get(idx).map(|v| v.num as f64 / v.denom as f64), + Value::Short(v) => v.get(idx).map(|v| *v as f64), + Value::Long(v) => v.get(idx).map(|v| *v as f64), + Value::SShort(v) => v.get(idx).map(|v| *v as f64), + Value::SLong(v) => v.get(idx).map(|v| *v as f64), + _ => None, + } +} + +// --------------------------------------------------------------------------- +// strtod-like scanner for the sscanf format +// --------------------------------------------------------------------------- + +struct Scanner<'a> { + s: &'a str, +} + +impl Scanner<'_> { + fn skip_ws(&mut self) { + self.s = self.s.trim_start(); + } + + fn lit(&mut self, lit: &str) -> bool { + if self.s.starts_with(lit) { + self.s = &self.s[lit.len()..]; + true + } else { + false + } + } + + /// %lf: skip whitespace, parse the longest numeric prefix. + fn f64(&mut self) -> Option { + self.skip_ws(); + for i in (1..=self.s.len()).rev() { + if let Ok(v) = self.s[..i].parse::() { + self.s = &self.s[i..]; + return Some(v); + } + } + None + } + + fn f64_after(&mut self, lit: &str) -> Option { + if self.lit(lit) { + self.f64() + } else { + None + } + } + + /// %d + fn i32(&mut self) -> Option { + self.skip_ws(); + for i in (1..=self.s.len()).rev() { + if let Ok(v) = self.s[..i].parse::() { + self.s = &self.s[i..]; + return Some(v); + } + } + None + } +} + +// --------------------------------------------------------------------------- +// JPEG COM segment reading / writing +// --------------------------------------------------------------------------- + +fn marker_len_ok(marker: u8) -> bool { + // standalone markers without length field + marker == 0x01 || (0xD0..=0xD9).contains(&marker) +} + +/// Return the payload of the first COM segment before SOS, if any. +pub fn read_jpeg_comment(path: &str) -> Option> { + let bytes = std::fs::read(path).ok()?; + if bytes.len() < 2 || bytes[0] != 0xFF || bytes[1] != 0xD8 { + return None; + } + + let mut i = 2; + while i + 1 < bytes.len() { + if bytes[i] != 0xFF { + return None; + } + let marker = bytes[i + 1]; + if marker == 0xFF { + i += 1; + continue; + } + if marker == 0xDA || marker == 0xD9 { + return None; + } + if marker_len_ok(marker) { + i += 2; + continue; + } + if i + 4 > bytes.len() { + return None; + } + let len = (u16::from_be_bytes([bytes[i + 2], bytes[i + 3]]) as usize).checked_add(2)?; + if i + len > bytes.len() { + return None; + } + if marker == 0xFE { + return Some(bytes[i + 4..i + 2 + len - 2].to_vec()); + } + i += len; + } + None +} + +/// Copy `in_img` to `out_img`, dropping all COM segments and inserting a +/// fresh COM segment right after SOI. +pub fn write_jpeg_comment(in_img: &str, out_img: &str, comment: &[u8]) -> std::io::Result<()> { + let bytes = std::fs::read(in_img)?; + if bytes.len() < 2 || bytes[0] != 0xFF || bytes[1] != 0xD8 { + return Err(std::io::Error::new(std::io::ErrorKind::InvalidInput, "not a JPEG file")); + } + if comment.len() + 2 > 0xFFFF { + return Err(std::io::Error::new(std::io::ErrorKind::InvalidInput, "comment too long")); + } + + let mut out: Vec = Vec::with_capacity(bytes.len() + comment.len() + 4); + out.extend_from_slice(&bytes[..2]); // SOI + + // Copy the SOI-adjacent header region (APPn/COM/DPI/etc.), dropping COM + // segments, and insert our fresh COM after the last APPn segment: + // that is where exiv2 (and thus the C gipfel) expects it. + let mut i = 2; + let mut insert_at = out.len(); + let mut copy_rest_from = bytes.len(); + while i + 1 < bytes.len() { + if bytes[i] != 0xFF { + break; + } + let marker = bytes[i + 1]; + if marker == 0xFF { + i += 1; + continue; + } + if marker == 0xDA || marker == 0xD9 { + copy_rest_from = i; + break; + } + if marker_len_ok(marker) { + out.push(0xFF); + out.push(marker); + i += 2; + continue; + } + if i + 4 > bytes.len() { + break; + } + let len = match u16::from_be_bytes([bytes[i + 2], bytes[i + 3]]) { + l => (l as usize).checked_add(2).unwrap_or(usize::MAX), + }; + if len < 2 || i + len > bytes.len() { + break; + } + if marker != 0xFE { + out.extend_from_slice(&bytes[i..i + len]); + if (0xE0..=0xEF).contains(&marker) { + insert_at = out.len(); + } + } + i += len; + } + copy_rest_from = copy_rest_from.min(i.max(2)); + + let mut com: Vec = Vec::with_capacity(comment.len() + 4); + com.push(0xFF); + com.push(0xFE); + com.extend_from_slice(&((comment.len() + 2) as u16).to_be_bytes()); + com.extend_from_slice(comment); + out.splice(insert_at..insert_at, com); + + out.extend_from_slice(&bytes[copy_rest_from..]); + + let mut f = File::create(out_img)?; + f.write_all(&out)?; + f.flush() +} + +/// Read raw image dimensions (width, height) of a JPEG. +pub fn image_dimensions(path: &str) -> Option<(u32, u32)> { + image::image_dimensions(path).ok() +} + +/// Apply image metadata to a panorama like `GipfelWidget::load_image` does. +/// Returns false when the gipfel info is incomplete (C: `have_gipfel_info`). +pub fn apply_image_metadata( + pan: &mut crate::panorama::Panorama, + md: &ImageMetaData, + img_w: u32, + img_h: u32, +) -> bool { + use crate::projection::ProjectionType; + + let mut have_info = true; + + pan.set_view_long(md.longitude); + pan.set_view_lat(md.latitude); + pan.set_view_height(md.height); + + pan.set_projection(match md.projection_type { + 1 => ProjectionType::Cylindrical, + _ => ProjectionType::Rectilinear, + }); + + if md.direction.is_nan() { + pan.set_center_angle(0.0); + have_info = false; + } else { + pan.set_center_angle(md.direction); + } + + if md.pitch.is_nan() { + pan.set_pitch_angle(0.0); + have_info = false; + } else { + pan.set_pitch_angle(md.pitch); + } + + if md.roll.is_nan() { + pan.set_roll_angle(0.0); + have_info = false; + } else { + pan.set_roll_angle(md.roll); + } + + // set_focal_length_35mm: scale = fl * max(w, h) / 35 + let w = img_w.max(img_h); + if md.focal_length_35mm.is_nan() || md.focal_length_35mm == 0.0 { + pan.set_scale(35.0 * w as f64 / 35.0); + have_info = false; + } else { + pan.set_scale(md.focal_length_35mm * w as f64 / 35.0); + } + + // distortion parameters, C ordering: + // 1. gipfel data in JPEG comment + // 2. matching distortion profile (camera + focal length) + // 3. 0.0, 0.0, 0.0 + let (mut k0, mut k1, mut x0) = (md.k0, md.k1, md.x0); + if k0.is_nan() { + if let Some(p) = md + .manufacturer + .as_deref() + .zip(md.model.as_deref()) + .and_then(|(mfr, model)| crate::profile::profile_name(mfr, model, md.focal_length)) + .and_then(|name| crate::profile::load(&name).ok()) + { + k0 = p.k0; + k1 = p.k1; + x0 = p.x0; + } + } + if k0.is_nan() { + k0 = 0.0; + } + if k1.is_nan() { + k1 = 0.0; + } + if x0.is_nan() { + x0 = 0.0; + } + pan.set_distortion_params(k0, k1, x0); + + have_info +} + +/// The profile name auto-matched from this image's camera and focal length +/// (C `get_distortion_profile_name`). Suggested whenever the camera is +/// known — [`apply_image_metadata`] only *applies* it when the JPEG +/// comment carries no distortion params. +pub fn distortion_profile_name(md: &ImageMetaData) -> Option { + md.manufacturer + .as_deref() + .zip(md.model.as_deref()) + .and_then(|(mfr, model)| crate::profile::profile_name(mfr, model, md.focal_length)) +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn jpeg_comment_write_read() { + let dir = std::env::temp_dir().join("gipfel-core-mdtest"); + std::fs::create_dir_all(&dir).unwrap(); + let img = dir.join("t.jpg"); + let out = dir.join("o.jpg"); + + image::DynamicImage::new_rgb8(64, 48) + .grayscale() + .save(&img) + .unwrap(); + + let comment = b"gipfel: longitude 1.000000, latitude 2.000000, height 3.000000, direction 4.000000, nick 5.000000, tilt 6.000000, focal_length_35mm 7.000000, projection type 0, k0 0.000000, k1 0.000000, x0 0.000000"; + write_jpeg_comment(img.to_str().unwrap(), out.to_str().unwrap(), comment).unwrap(); + + // no COM before, one after + assert!(read_jpeg_comment(img.to_str().unwrap()).is_none()); + assert_eq!(read_jpeg_comment(out.to_str().unwrap()).unwrap(), comment); + + // rewriting replaces the existing comment + let out2 = dir.join("o2.jpg"); + write_jpeg_comment(out.to_str().unwrap(), out2.to_str().unwrap(), b"short").unwrap(); + assert_eq!(read_jpeg_comment(out2.to_str().unwrap()).unwrap(), b"short"); + + // image still decodes + assert!(image::image_dimensions(out.to_str().unwrap()).is_ok()); + + let _ = std::fs::remove_dir_all(&dir); + } + + #[test] + fn distortion_precedence_and_legacy_import() { + use crate::panorama::Panorama; + use crate::profile; + + let dir = std::env::temp_dir().join(format!("gipfel-prec-{}", std::process::id())); + let _ = std::fs::remove_dir_all(&dir); + let data = dir.join("data"); + let home = dir.join("home"); + std::fs::create_dir_all(&data).unwrap(); + + let old = (std::env::var("HOME").ok(), std::env::var("XDG_DATA_HOME").ok()); + std::env::set_var("HOME", &home); + std::env::set_var("XDG_DATA_HOME", &data); + + // legacy FLTK prefs file (C version's location and format) + let legacy = home + .join(".fltk") + .join("Johannes.HofmannATgmx.de") + .join("gipfel") + .join("DistortionProfiles.prefs"); + std::fs::create_dir_all(legacy.parent().unwrap()).unwrap(); + std::fs::write( + &legacy, + "; FLTK preferences file format 1.0\n; vendor: Johannes.HofmannATgmx.de\n\ + ; application: gipfel/DistortionProfiles\n\n[.]\n\n\n\ + [./Canon_EOS_5D_50.00_mm]\n\nk0:-0.01\nk1:0.001\nx0:3\n", + ) + .unwrap(); + + // first access imports into the XDG data dir + assert!(!profile::prefs_path().exists()); + let all = profile::list(); + assert!(all.iter().any(|(n, p)| n == "Canon_EOS_5D_50.00_mm" && p.k0 == -0.01)); + assert!(profile::prefs_path().starts_with(&data)); + assert!(profile::prefs_path().exists()); + + let mut md = ImageMetaData::new(); + md.manufacturer = Some("Canon".into()); + md.model = Some("EOS_5D".into()); + md.focal_length = 50.0; + let mut pan = Panorama::new(); + + // 1. comment wins when present + md.k0 = -0.9; + md.k1 = 0.0; + md.x0 = 0.0; + apply_image_metadata(&mut pan, &md, 1000, 800); + assert!((pan.parms.k0 + 0.9).abs() < 1e-12, "comment should win"); + + // 2. profile used when comment absent + md.k0 = f64::NAN; + md.k1 = f64::NAN; + md.x0 = f64::NAN; + apply_image_metadata(&mut pan, &md, 1000, 800); + assert!( + (pan.parms.k0 + 0.01).abs() < 1e-12, + "profile should apply: k0={}", + pan.parms.k0 + ); + // name is suggested regardless of the comment (C save dialog) + assert_eq!( + distortion_profile_name(&md).as_deref(), + Some("Canon_EOS_5D_50.00_mm") + ); + md.k0 = 0.0; + assert_eq!( + distortion_profile_name(&md).as_deref(), + Some("Canon_EOS_5D_50.00_mm") + ); + + // 3. zero when neither + md.k0 = f64::NAN; + md.manufacturer = None; + apply_image_metadata(&mut pan, &md, 1000, 800); + assert_eq!(pan.parms.k0, 0.0); + + match old.0 { + Some(h) => std::env::set_var("HOME", h), + None => std::env::remove_var("HOME"), + } + match old.1 { + Some(h) => std::env::set_var("XDG_DATA_HOME", h), + None => std::env::remove_var("XDG_DATA_HOME"), + } + let _ = std::fs::remove_dir_all(&dir); + } +} 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, + /// indices into `hills`, sorted by alpha descending (rebuilt by update_angles) + pub alpha_order: Vec, + /// indices into `hills`, sorted by alpha descending + pub close: Vec, + /// indices into `hills` + pub visible: Vec, + + pub view_name: Option, + /// 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, +} + +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) { + 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 { + 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 = (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 + /// + 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 { + 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)); + } +} diff --git a/crates/gipfel-core/src/profile.rs b/crates/gipfel-core/src/profile.rs new file mode 100644 index 0000000..e8bda68 --- /dev/null +++ b/crates/gipfel-core/src/profile.rs @@ -0,0 +1,362 @@ +//! Lens distortion profiles. +//! +//! Profiles live in the standard per-user data directory +//! (`$XDG_DATA_HOME/gipfel/distortion-profiles.prefs`, i.e. +//! `~/.local/share/gipfel/…`). The file format is the FLTK preferences +//! format used by the C version (`[./name]` groups, `key:value` entries, +//! `%g` values), so the legacy file +//! `~/.fltk/Johannes.HofmannATgmx.de/gipfel/DistortionProfiles.prefs` +//! imports automatically into the new location on first use. +//! +//! File layout (as written by FLTK): +//! ```text +//! ; FLTK preferences file format 1.0 +//! ; vendor: Johannes.HofmannATgmx.de +//! ; application: gipfel/DistortionProfiles +//! +//! [.] +//! +//! [./Canon_EOS_5D_50.00_mm] +//! +//! k0:-0.00123 +//! k1:2.5e-05 +//! x0:0.5 +//! ``` + +use std::path::{Path, PathBuf}; + +#[derive(Clone, Copy, Debug, PartialEq)] +pub struct Profile { + pub k0: f64, + pub k1: f64, + pub x0: f64, +} + +#[derive(Debug)] +pub enum Error { + NotFound, + Exists, + BadName, + Io(std::io::Error), +} + +impl std::fmt::Display for Error { + fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { + match self { + Error::NotFound => write!(f, "profile not found"), + Error::Exists => write!(f, "profile already exists"), + Error::BadName => write!(f, "invalid profile name"), + Error::Io(e) => write!(f, "{e}"), + } + } +} + +/// The profile store in the standard data directory +/// (`$XDG_DATA_HOME/gipfel/distortion-profiles.prefs`). +pub fn prefs_path() -> PathBuf { + dirs::data_dir() + .unwrap_or_else(|| PathBuf::from(std::env::var("HOME").unwrap_or_default())) + .join("gipfel") + .join("distortion-profiles.prefs") +} + +/// Legacy FLTK preferences file used by the C version. +fn fltk_prefs_path() -> PathBuf { + let home = std::env::var("HOME").unwrap_or_default(); + Path::new(&home) + .join(".fltk") + .join("Johannes.HofmannATgmx.de") + .join("gipfel") + .join("DistortionProfiles.prefs") +} + +/// Path to the profile store; imports legacy FLTK profiles on first use. +fn ensure_store() -> PathBuf { + let path = prefs_path(); + if !path.exists() { + let legacy = fltk_prefs_path(); + if legacy.exists() { + if let Some(dir) = path.parent() { + let _ = std::fs::create_dir_all(dir); + } + // same format, so a copy suffices + let _ = std::fs::copy(&legacy, &path); + } + } + path +} + +/// C `get_distortion_profile_name`: `"%s_%s_%.2f_mm"`. `None` when the +/// camera is unknown or the focal length is unusable (C would build a name +/// containing "nan" that can never match a stored profile). +pub fn profile_name(manufacturer: &str, model: &str, focal_mm: f64) -> Option { + if manufacturer.is_empty() || model.is_empty() || !focal_mm.is_finite() || focal_mm <= 0.0 { + return None; + } + let name = format!("{manufacturer}_{model}_{focal_mm:.2}_mm"); + sanitize_name(&name).ok().map(|_| name) +} + +/// Reject names that would break the prefs file structure (FLTK group paths +/// use '/' as separator; control chars/newlines corrupt the file). +pub fn sanitize_name(name: &str) -> Result<(), Error> { + if name.trim().is_empty() || name.len() > 255 { + return Err(Error::BadName); + } + if name + .chars() + .any(|c| c == '/' || c == '\\' || c.is_control()) + { + return Err(Error::BadName); + } + Ok(()) +} + +/// Parse a prefs file into ordered (name, profile) pairs. Tolerant: +/// comments (`;`), the root group `[.]`, blank lines and unknown keys are +/// ignored; values are parsed with Rust's `f64` (superset of C `%g`). +fn parse(text: &str) -> Vec<(String, Profile)> { + let mut out: Vec<(String, Profile)> = Vec::new(); + let mut cur: Option = None; + for line in text.lines() { + let line = line.trim_end(); + if line.is_empty() || line.starts_with(';') { + continue; + } + if let Some(group) = line.strip_prefix('[').and_then(|s| s.strip_suffix(']')) { + cur = match group.strip_prefix("./") { + Some(name) => { + out.push(( + name.to_string(), + Profile { + k0: f64::NAN, + k1: f64::NAN, + x0: f64::NAN, + }, + )); + Some(out.len() - 1) + } + None => None, // "[.]" root group + }; + continue; + } + let Some(idx) = cur else { continue }; + if let Some((key, val)) = line.split_once(':') { + let Ok(v) = val.trim().parse::() else { + continue; + }; + let p = &mut out[idx].1; + match key.trim() { + "k0" => p.k0 = v, + "k1" => p.k1 = v, + "x0" => p.x0 = v, + _ => {} + } + } + } + out.retain(|(_, p)| p.k0.is_finite() && p.k1.is_finite() && p.x0.is_finite()); + out +} + +/// Serialize in FLTK's layout so the C version can read our file. +fn serialize(profiles: &[(String, Profile)]) -> String { + let mut s = String::from( + "; FLTK preferences file format 1.0\n\ + ; vendor: Johannes.HofmannATgmx.de\n\ + ; application: gipfel/DistortionProfiles\n\ + \n\ + [.]\n\ + \n", + ); + for (name, p) in profiles { + s.push_str(&format!( + "\n[./{name}]\n\nk0:{}\nk1:{}\nx0:{}\n", + fmt_g(p.k0), + fmt_g(p.k1), + fmt_g(p.x0) + )); + } + s +} + +/// Closest Rust equivalent of C `%g`: scientific notation for magnitudes +/// outside [1e-4, 1e15), exponent zero-padded to two digits ("2.5e-05"). +/// Unlike C we keep the shortest round-trip mantissa (more precise; the +/// FLTK reader uses strtod and reads either fine). +fn fmt_g(v: f64) -> String { + if v != 0.0 && (v.abs() < 1e-4 || v.abs() >= 1e15) { + if let Some((mant, exp)) = format!("{v:e}").split_once('e') { + let neg = exp.starts_with('-'); + let digits = exp.trim_start_matches(['-', '+']); + let digits = if digits.len() < 2 { format!("0{digits}") } else { digits.to_string() }; + let sign = if neg { '-' } else { '+' }; + return format!("{mant}e{sign}{digits}"); + } + } + if v == v.trunc() && v.abs() < 1e16 { + format!("{}", v as i64) + } else { + format!("{v}") + } +} + +fn read_all(path: &Path) -> Result, Error> { + match std::fs::read_to_string(path) { + Ok(text) => Ok(parse(&text)), + Err(e) if e.kind() == std::io::ErrorKind::NotFound => Ok(Vec::new()), + Err(e) => Err(Error::Io(e)), + } +} + +fn write_all(path: &Path, profiles: &[(String, Profile)]) -> Result<(), Error> { + if let Some(dir) = path.parent() { + std::fs::create_dir_all(dir).map_err(Error::Io)?; + } + // atomic: temp file + rename + let tmp = path.with_extension("prefs.tmp"); + std::fs::write(&tmp, serialize(profiles)).map_err(Error::Io)?; + std::fs::rename(&tmp, path).map_err(Error::Io)?; + Ok(()) +} + +/// Load a profile by exact name from the default store. +pub fn load(name: &str) -> Result { + load_from(&ensure_store(), name) +} + +pub fn load_from(path: &Path, name: &str) -> Result { + read_all(path)? + .into_iter() + .find(|(n, _)| n == name) + .map(|(_, p)| p) + .ok_or(Error::NotFound) +} + +/// Save a profile. With `force == false` an existing profile of that name is +/// not overwritten (C `save_distortion_params(name, force)`). +pub fn save(name: &str, profile: Profile, force: bool) -> Result<(), Error> { + save_to(&ensure_store(), name, profile, force) +} + +pub fn save_to(path: &Path, name: &str, profile: Profile, force: bool) -> Result<(), Error> { + sanitize_name(name)?; + let mut all = read_all(path)?; + if let Some(pos) = all.iter().position(|(n, _)| n == name) { + if !force { + return Err(Error::Exists); + } + all[pos].1 = profile; + } else { + all.push((name.to_string(), profile)); + } + write_all(path, &all) +} + +/// All profiles in the default store, in file order. +pub fn list() -> Vec<(String, Profile)> { + read_all(&ensure_store()).unwrap_or_default() +} + +pub fn list_from(path: &Path) -> Vec<(String, Profile)> { + read_all(path).unwrap_or_default() +} + +#[cfg(test)] +mod tests { + use super::*; + + /// Exact bytes as written by the C version (verified against FLTK 1.4). + const FLTK_FIXTURE: &str = "; FLTK preferences file format 1.0\n\ + ; vendor: Johannes.HofmannATgmx.de\n\ + ; application: gipfel/DistortionProfiles\n\ + \n\ + [.]\n\ + \n\ + \n\ + [./Canon_EOS_5D_50.00_mm]\n\ + \n\ + k0:-0.00123\n\ + k1:2.5e-05\n\ + x0:0.5\n\ + \n\ + [./Nikon D850_35.00_mm]\n\ + \n\ + k0:0.01\n\ + k1:0\n\ + x0:0\n"; + + #[test] + fn parses_fltk_layout() { + let all = parse(FLTK_FIXTURE); + assert_eq!(all.len(), 2); + assert_eq!(all[0].0, "Canon_EOS_5D_50.00_mm"); + assert_eq!(all[0].1.k0, -0.00123); + assert_eq!(all[0].1.k1, 2.5e-5); + assert_eq!(all[0].1.x0, 0.5); + assert_eq!(all[1].0, "Nikon D850_35.00_mm"); + assert_eq!(all[1].1.k0, 0.01); + assert_eq!(all[1].1.k1, 0.0); + } + + #[test] + fn serialize_replicates_fltk_layout() { + let all = parse(FLTK_FIXTURE); + assert_eq!(serialize(&all), FLTK_FIXTURE); + } + + #[test] + fn fmt_g_matches_c_percent_g() { + assert_eq!(fmt_g(2.5e-5), "2.5e-05"); + assert_eq!(fmt_g(-0.00123), "-0.00123"); + assert_eq!(fmt_g(0.0001), "0.0001"); + assert_eq!(fmt_g(1e15), "1e+15"); + assert_eq!(fmt_g(0.0), "0"); + assert_eq!(fmt_g(0.01), "0.01"); + } + + #[test] + fn profile_name_matches_c() { + assert_eq!( + profile_name("Canon", "EOS_5D", 50.0).as_deref(), + Some("Canon_EOS_5D_50.00_mm") + ); + assert_eq!( + profile_name("Nikon", "D850", 34.666).as_deref(), + Some("Nikon_D850_34.67_mm") + ); + assert!(profile_name("", "D850", 35.0).is_none()); + assert!(profile_name("Canon", "D850", f64::NAN).is_none()); + assert!(profile_name("Can/on", "D850", 35.0).is_none()); + } + + #[test] + fn roundtrip_and_overwrite() { + let dir = std::env::temp_dir().join(format!("gipfel-prof-{}", std::process::id())); + let path = dir.join("DistortionProfiles.prefs"); + let _ = std::fs::remove_dir_all(&dir); + + let p = Profile { k0: -0.001, k1: 2e-5, x0: 0.0 }; + save_to(&path, "Cam_A_24.00_mm", p, false).unwrap(); + // exists -> error without force + assert!(matches!( + save_to(&path, "Cam_A_24.00_mm", p, false), + Err(Error::Exists) + )); + save_to(&path, "Cam_B_50.00_mm", Profile { k0: 0.1, k1: 0.0, x0: 1.0 }, false).unwrap(); + // overwrite first, keep second + let p2 = Profile { k0: -0.002, k1: 2e-5, x0: 0.0 }; + save_to(&path, "Cam_A_24.00_mm", p2, true).unwrap(); + + let all = list_from(&path); + assert_eq!(all.len(), 2); + assert_eq!(all[0].1, p2); + assert_eq!(load_from(&path, "Cam_B_50.00_mm").unwrap().k0, 0.1); + assert!(matches!(load_from(&path, "nope"), Err(Error::NotFound))); + + // bad names rejected + assert!(save_to(&path, "a/b", p, true).is_err()); + assert!(save_to(&path, "with\nnewline", p, true).is_err()); + + std::fs::remove_dir_all(&dir).unwrap(); + } +} diff --git a/crates/gipfel-core/src/projection/cylindrical.rs b/crates/gipfel-core/src/projection/cylindrical.rs new file mode 100644 index 0000000..9adf79c --- /dev/null +++ b/crates/gipfel-core/src/projection/cylindrical.rs @@ -0,0 +1,86 @@ +//! Cylindrical projection. + +mod funcs { + #![allow(non_snake_case, unused_parens, unused_variables, dead_code, clippy::all)] + + // math shims: the maxima-generated code uses C-style calls + #[inline(always)] + fn sin(x: f64) -> f64 { + x.sin() + } + #[inline(always)] + fn cos(x: f64) -> f64 { + x.cos() + } + #[inline(always)] + fn tan(x: f64) -> f64 { + x.tan() + } + #[inline(always)] + fn sqrt(x: f64) -> f64 { + x.sqrt() + } + #[inline(always)] + fn pow(x: f64, y: f64) -> f64 { + x.powf(y) + } + #[inline(always)] + fn exp(x: f64) -> f64 { + x.exp() + } + #[inline(always)] + fn fabs(x: f64) -> f64 { + x.abs() + } + #[inline(always)] + fn atan2(y: f64, x: f64) -> f64 { + y.atan2(x) + } + #[inline(always)] + fn sec(x: f64) -> f64 { + 1.0 / x.cos() + } + + include!("generated/cylindrical_funcs.rs"); +} + +use super::{ProjArgs, Projection}; + +pub struct ProjectionCylindrical; + +impl Projection for ProjectionCylindrical { + fn get_view_angle(&self) -> f64 { + 6.2831853 /* 360 deg */ + } + + fn mac_x(&self, a: ProjArgs) -> f64 { + funcs::mac_x(a.c_view, a.c_pitch, a.c_roll, a.scale, a.k0, a.k1, a.x0, a.m_view, a.m_elev) + } + fn mac_y(&self, a: ProjArgs) -> f64 { + funcs::mac_y(a.c_view, a.c_pitch, a.c_roll, a.scale, a.k0, a.k1, a.x0, a.m_view, a.m_elev) + } + fn mac_x_dc_view(&self, a: ProjArgs) -> f64 { + funcs::mac_x_dc_view(a.c_view, a.c_pitch, a.c_roll, a.scale, a.k0, a.k1, a.x0, a.m_view, a.m_elev) + } + fn mac_x_dc_pitch(&self, a: ProjArgs) -> f64 { + funcs::mac_x_dc_pitch(a.c_view, a.c_pitch, a.c_roll, a.scale, a.k0, a.k1, a.x0, a.m_view, a.m_elev) + } + fn mac_x_dc_roll(&self, a: ProjArgs) -> f64 { + funcs::mac_x_dc_roll(a.c_view, a.c_pitch, a.c_roll, a.scale, a.k0, a.k1, a.x0, a.m_view, a.m_elev) + } + fn mac_x_dscale(&self, a: ProjArgs) -> f64 { + funcs::mac_x_dscale(a.c_view, a.c_pitch, a.c_roll, a.scale, a.k0, a.k1, a.x0, a.m_view, a.m_elev) + } + fn mac_y_dc_view(&self, a: ProjArgs) -> f64 { + funcs::mac_y_dc_view(a.c_view, a.c_pitch, a.c_roll, a.scale, a.k0, a.k1, a.x0, a.m_view, a.m_elev) + } + fn mac_y_dc_pitch(&self, a: ProjArgs) -> f64 { + funcs::mac_y_dc_pitch(a.c_view, a.c_pitch, a.c_roll, a.scale, a.k0, a.k1, a.x0, a.m_view, a.m_elev) + } + fn mac_y_dc_roll(&self, a: ProjArgs) -> f64 { + funcs::mac_y_dc_roll(a.c_view, a.c_pitch, a.c_roll, a.scale, a.k0, a.k1, a.x0, a.m_view, a.m_elev) + } + fn mac_y_dscale(&self, a: ProjArgs) -> f64 { + funcs::mac_y_dscale(a.c_view, a.c_pitch, a.c_roll, a.scale, a.k0, a.k1, a.x0, a.m_view, a.m_elev) + } +} diff --git a/crates/gipfel-core/src/projection/generated/cylindrical_funcs.rs b/crates/gipfel-core/src/projection/generated/cylindrical_funcs.rs new file mode 100644 index 0000000..220a7d8 --- /dev/null +++ b/crates/gipfel-core/src/projection/generated/cylindrical_funcs.rs @@ -0,0 +1,10 @@ +pub fn mac_x ( c_view: f64, c_pitch: f64, c_roll: f64, scale: f64, k0: f64, k1: f64, x0: f64, m_view: f64, m_elev: f64 ) -> f64 { let mut _1 : f64 = 0.0; let mut _2 : f64 = 0.0; _1 = tan(c_pitch) ; _2 = tan(m_elev) ; return (cos(c_roll)*(m_view-c_view)-(sin(c_roll)*(_2-_1))/(_1*_2+1.0))*scale ; } +pub fn mac_y ( c_view: f64, c_pitch: f64, c_roll: f64, scale: f64, k0: f64, k1: f64, x0: f64, m_view: f64, m_elev: f64 ) -> f64 { let mut _1 : f64 = 0.0; let mut _2 : f64 = 0.0; _1 = tan(c_pitch) ; _2 = tan(m_elev) ; return ((-sin(c_roll)*(m_view-c_view))-(cos(c_roll)*(_2-_1))/(_1*_2+1.0))*scale ; } +pub fn mac_x_dc_view ( c_view: f64, c_pitch: f64, c_roll: f64, scale: f64, k0: f64, k1: f64, x0: f64, m_view: f64, m_elev: f64 ) -> f64 { return -cos(c_roll)*scale ; } +pub fn mac_x_dc_pitch ( c_view: f64, c_pitch: f64, c_roll: f64, scale: f64, k0: f64, k1: f64, x0: f64, m_view: f64, m_elev: f64 ) -> f64 { let mut _1 : f64 = 0.0; let mut _2 : f64 = 0.0; let mut _3 : f64 = 0.0; let mut _4 : f64 = 0.0; let mut _5 : f64 = 0.0; _1 = pow(sec(c_pitch),2.0) ; _2 = sin(c_roll) ; _3 = tan(m_elev) ; _4 = tan(c_pitch) ; _5 = _4*_3+1.0 ; return ((_1*_2*_3*(_3-_4))/pow(_5,2.0)+(_1*_2)/_5)*scale ; } +pub fn mac_x_dc_roll ( c_view: f64, c_pitch: f64, c_roll: f64, scale: f64, k0: f64, k1: f64, x0: f64, m_view: f64, m_elev: f64 ) -> f64 { let mut _1 : f64 = 0.0; let mut _2 : f64 = 0.0; _1 = tan(c_pitch) ; _2 = tan(m_elev) ; return ((-sin(c_roll)*(m_view-c_view))-(cos(c_roll)*(_2-_1))/(_1*_2+1.0))*scale ; } +pub fn mac_x_dscale ( c_view: f64, c_pitch: f64, c_roll: f64, scale: f64, k0: f64, k1: f64, x0: f64, m_view: f64, m_elev: f64 ) -> f64 { let mut _1 : f64 = 0.0; let mut _2 : f64 = 0.0; _1 = tan(c_pitch) ; _2 = tan(m_elev) ; return cos(c_roll)*(m_view-c_view)-(sin(c_roll)*(_2-_1))/(_1*_2+1.0) ; } +pub fn mac_y_dc_view ( c_view: f64, c_pitch: f64, c_roll: f64, scale: f64, k0: f64, k1: f64, x0: f64, m_view: f64, m_elev: f64 ) -> f64 { return sin(c_roll)*scale ; } +pub fn mac_y_dc_pitch ( c_view: f64, c_pitch: f64, c_roll: f64, scale: f64, k0: f64, k1: f64, x0: f64, m_view: f64, m_elev: f64 ) -> f64 { let mut _1 : f64 = 0.0; let mut _2 : f64 = 0.0; let mut _3 : f64 = 0.0; let mut _4 : f64 = 0.0; let mut _5 : f64 = 0.0; _1 = pow(sec(c_pitch),2.0) ; _2 = cos(c_roll) ; _3 = tan(m_elev) ; _4 = tan(c_pitch) ; _5 = _4*_3+1.0 ; return ((_1*_2*_3*(_3-_4))/pow(_5,2.0)+(_1*_2)/_5)*scale ; } +pub fn mac_y_dc_roll ( c_view: f64, c_pitch: f64, c_roll: f64, scale: f64, k0: f64, k1: f64, x0: f64, m_view: f64, m_elev: f64 ) -> f64 { let mut _1 : f64 = 0.0; let mut _2 : f64 = 0.0; _1 = tan(c_pitch) ; _2 = tan(m_elev) ; return ((sin(c_roll)*(_2-_1))/(_1*_2+1.0)-cos(c_roll)*(m_view-c_view))*scale ; } +pub fn mac_y_dscale ( c_view: f64, c_pitch: f64, c_roll: f64, scale: f64, k0: f64, k1: f64, x0: f64, m_view: f64, m_elev: f64 ) -> f64 { let mut _1 : f64 = 0.0; let mut _2 : f64 = 0.0; _1 = tan(c_pitch) ; _2 = tan(m_elev) ; return (-sin(c_roll)*(m_view-c_view))-(cos(c_roll)*(_2-_1))/(_1*_2+1.0) ; } diff --git a/crates/gipfel-core/src/projection/generated/rectilinear_funcs.rs b/crates/gipfel-core/src/projection/generated/rectilinear_funcs.rs new file mode 100644 index 0000000..7bd450d --- /dev/null +++ b/crates/gipfel-core/src/projection/generated/rectilinear_funcs.rs @@ -0,0 +1,16 @@ +pub fn mac_x ( c_view: f64, c_pitch: f64, c_roll: f64, scale: f64, k0: f64, k1: f64, x0: f64, m_view: f64, m_elev: f64 ) -> f64 { let mut _1 : f64 = 0.0; let mut _2 : f64 = 0.0; let mut _3 : f64 = 0.0; let mut _4 : f64 = 0.0; let mut _5 : f64 = 0.0; let mut _6 : f64 = 0.0; let mut _7 : f64 = 0.0; let mut _8 : f64 = 0.0; let mut _9 : f64 = 0.0; let mut _10 : f64 = 0.0; let mut _11 : f64 = 0.0; let mut _12 : f64 = 0.0; let mut _13 : f64 = 0.0; let mut _14 : f64 = 0.0; let mut _15 : f64 = 0.0; let mut _16 : f64 = 0.0; _1 = sin(c_roll) ; _2 = sin(c_pitch) ; _3 = sin(m_elev) ; _4 = cos(c_pitch) ; _5 = cos(m_elev) ; _6 = cos(c_view) ; _7 = cos(m_view) ; _8 = sin(c_view) ; _9 = sin(m_view) ; _10 = _8*_9+_6*_7 ; _11 = 1.0/(_4*_5*_10+_2*_3) ; _12 = _2*_5*_10-_4*_3 ; _13 = cos(c_roll) ; _14 = _6*_9-_8*_7 ; _15 = x0+_13*_5*_14*_11+_1*_11*_12 ; _16 = pow(_15,2.0)+pow(_13*_11*_12-_1*_5*_14*_11,2.0) ; return scale*_15*(k0*sqrt(_16)+k1*_16+1.0) ; } +pub fn mac_y ( c_view: f64, c_pitch: f64, c_roll: f64, scale: f64, k0: f64, k1: f64, x0: f64, m_view: f64, m_elev: f64 ) -> f64 { let mut _1 : f64 = 0.0; let mut _2 : f64 = 0.0; let mut _3 : f64 = 0.0; let mut _4 : f64 = 0.0; let mut _5 : f64 = 0.0; let mut _6 : f64 = 0.0; let mut _7 : f64 = 0.0; let mut _8 : f64 = 0.0; let mut _9 : f64 = 0.0; let mut _10 : f64 = 0.0; let mut _11 : f64 = 0.0; let mut _12 : f64 = 0.0; let mut _13 : f64 = 0.0; let mut _14 : f64 = 0.0; let mut _15 : f64 = 0.0; let mut _16 : f64 = 0.0; _1 = cos(c_roll) ; _2 = sin(c_pitch) ; _3 = sin(m_elev) ; _4 = cos(c_pitch) ; _5 = cos(m_elev) ; _6 = cos(c_view) ; _7 = cos(m_view) ; _8 = sin(c_view) ; _9 = sin(m_view) ; _10 = _8*_9+_6*_7 ; _11 = 1.0/(_4*_5*_10+_2*_3) ; _12 = _2*_5*_10-_4*_3 ; _13 = sin(c_roll) ; _14 = _6*_9-_8*_7 ; _15 = _1*_11*_12-_13*_5*_14*_11 ; _16 = pow(x0+_1*_5*_14*_11+_13*_11*_12,2.0)+pow(_15,2.0) ; return _15*scale*(k0*sqrt(_16)+k1*_16+1.0) ; } +pub fn mac_x_dc_view ( c_view: f64, c_pitch: f64, c_roll: f64, scale: f64, k0: f64, k1: f64, x0: f64, m_view: f64, m_elev: f64 ) -> f64 { let mut _1 : f64 = 0.0; let mut _2 : f64 = 0.0; let mut _3 : f64 = 0.0; let mut _4 : f64 = 0.0; let mut _5 : f64 = 0.0; let mut _6 : f64 = 0.0; let mut _7 : f64 = 0.0; let mut _8 : f64 = 0.0; let mut _9 : f64 = 0.0; let mut _10 : f64 = 0.0; let mut _11 : f64 = 0.0; let mut _12 : f64 = 0.0; let mut _13 : f64 = 0.0; let mut _14 : f64 = 0.0; let mut _15 : f64 = 0.0; let mut _16 : f64 = 0.0; let mut _17 : f64 = 0.0; let mut _18 : f64 = 0.0; let mut _19 : f64 = 0.0; let mut _20 : f64 = 0.0; let mut _21 : f64 = 0.0; let mut _22 : f64 = 0.0; let mut _23 : f64 = 0.0; let mut _24 : f64 = 0.0; let mut _25 : f64 = 0.0; _1 = sin(c_roll) ; _2 = sin(c_pitch) ; _3 = sin(m_elev) ; _4 = cos(c_pitch) ; _5 = cos(m_elev) ; _6 = cos(c_view) ; _7 = cos(m_view) ; _8 = sin(c_view) ; _9 = sin(m_view) ; _10 = _8*_9+_6*_7 ; _11 = _4*_5*_10+_2*_3 ; _12 = 1.0/_11 ; _13 = _2*_5*_10-_4*_3 ; _14 = cos(c_roll) ; _15 = _6*_9-_8*_7 ; _16 = x0+_14*_5*_15*_12+_1*_12*_13 ; _17 = _14*_12*_13-_1*_5*_15*_12 ; _18 = 1.0/pow(_11,2.0) ; _19 = pow(_5,2.0) ; _20 = pow(_15,2.0) ; _21 = (-_8*_9)-_6*_7 ; _22 = _14*_5*_21*_12+_2*_1*_5*_15*_12-_4*_14*_19*_20*_18-_4*_1*_5*_15*_18*_13 ; _23 = 2.0*_22*_16+2.0*_17*((-_1*_5*_21*_12)+_2*_14*_5*_15*_12+_4*_1*_19*_20*_18-_4*_14*_5*_15*_18*_13) ; _24 = pow(_17,2.0)+pow(_16,2.0) ; _25 = sqrt(_24) ; return _22*scale*(k0*_25+k1*_24+1.0)+scale*_16*((k0*_23)/(2.0*_25)+k1*_23) ; } +pub fn mac_x_dc_pitch ( c_view: f64, c_pitch: f64, c_roll: f64, scale: f64, k0: f64, k1: f64, x0: f64, m_view: f64, m_elev: f64 ) -> f64 { let mut _1 : f64 = 0.0; let mut _2 : f64 = 0.0; let mut _3 : f64 = 0.0; let mut _4 : f64 = 0.0; let mut _5 : f64 = 0.0; let mut _6 : f64 = 0.0; let mut _7 : f64 = 0.0; let mut _8 : f64 = 0.0; let mut _9 : f64 = 0.0; let mut _10 : f64 = 0.0; let mut _11 : f64 = 0.0; let mut _12 : f64 = 0.0; let mut _13 : f64 = 0.0; let mut _14 : f64 = 0.0; let mut _15 : f64 = 0.0; let mut _16 : f64 = 0.0; let mut _17 : f64 = 0.0; let mut _18 : f64 = 0.0; let mut _19 : f64 = 0.0; let mut _20 : f64 = 0.0; let mut _21 : f64 = 0.0; let mut _22 : f64 = 0.0; let mut _23 : f64 = 0.0; _1 = sin(c_roll) ; _2 = sin(c_pitch) ; _3 = sin(m_elev) ; _4 = cos(c_pitch) ; _5 = cos(m_elev) ; _6 = cos(c_view) ; _7 = cos(m_view) ; _8 = sin(c_view) ; _9 = sin(m_view) ; _10 = _8*_9+_6*_7 ; _11 = _4*_5*_10+_2*_3 ; _12 = 1.0/_11 ; _13 = _2*_5*_10-_4*_3 ; _14 = cos(c_roll) ; _15 = _6*_9-_8*_7 ; _16 = x0+_14*_5*_15*_12+_1*_12*_13 ; _17 = 1.0/pow(_11,2.0) ; _18 = _4*_3-_2*_5*_10 ; _19 = _14*_12*_13-_1*_5*_15*_12 ; _20 = (-_1*_17*_18*_13)-_14*_5*_15*_17*_18+_1 ; _21 = 2.0*_20*_16+2.0*((-_14*_17*_18*_13)+_1*_5*_15*_17*_18+_14)*_19 ; _22 = pow(_19,2.0)+pow(_16,2.0) ; _23 = sqrt(_22) ; return _20*scale*(k0*_23+k1*_22+1.0)+scale*_16*((k0*_21)/(2.0*_23)+k1*_21) ; } +pub fn mac_x_dc_roll ( c_view: f64, c_pitch: f64, c_roll: f64, scale: f64, k0: f64, k1: f64, x0: f64, m_view: f64, m_elev: f64 ) -> f64 { let mut _1 : f64 = 0.0; let mut _2 : f64 = 0.0; let mut _3 : f64 = 0.0; let mut _4 : f64 = 0.0; let mut _5 : f64 = 0.0; let mut _6 : f64 = 0.0; let mut _7 : f64 = 0.0; let mut _8 : f64 = 0.0; let mut _9 : f64 = 0.0; let mut _10 : f64 = 0.0; let mut _11 : f64 = 0.0; let mut _12 : f64 = 0.0; let mut _13 : f64 = 0.0; let mut _14 : f64 = 0.0; let mut _15 : f64 = 0.0; let mut _16 : f64 = 0.0; let mut _17 : f64 = 0.0; let mut _18 : f64 = 0.0; let mut _19 : f64 = 0.0; _1 = sin(c_roll) ; _2 = sin(c_pitch) ; _3 = sin(m_elev) ; _4 = cos(c_pitch) ; _5 = cos(m_elev) ; _6 = cos(c_view) ; _7 = cos(m_view) ; _8 = sin(c_view) ; _9 = sin(m_view) ; _10 = _8*_9+_6*_7 ; _11 = 1.0/(_4*_5*_10+_2*_3) ; _12 = _2*_5*_10-_4*_3 ; _13 = cos(c_roll) ; _14 = _6*_9-_8*_7 ; _15 = x0+_13*_5*_14*_11+_1*_11*_12 ; _16 = _13*_11*_12-_1*_5*_14*_11 ; _17 = 2.0*_16*_15+2.0*((-_13*_5*_14*_11)-_1*_11*_12)*_16 ; _18 = pow(_16,2.0)+pow(_15,2.0) ; _19 = sqrt(_18) ; return _16*scale*(k0*_19+k1*_18+1.0)+scale*_15*((k0*_17)/(2.0*_19)+k1*_17) ; } +pub fn mac_x_dscale ( c_view: f64, c_pitch: f64, c_roll: f64, scale: f64, k0: f64, k1: f64, x0: f64, m_view: f64, m_elev: f64 ) -> f64 { let mut _1 : f64 = 0.0; let mut _2 : f64 = 0.0; let mut _3 : f64 = 0.0; let mut _4 : f64 = 0.0; let mut _5 : f64 = 0.0; let mut _6 : f64 = 0.0; let mut _7 : f64 = 0.0; let mut _8 : f64 = 0.0; let mut _9 : f64 = 0.0; let mut _10 : f64 = 0.0; let mut _11 : f64 = 0.0; let mut _12 : f64 = 0.0; let mut _13 : f64 = 0.0; let mut _14 : f64 = 0.0; let mut _15 : f64 = 0.0; let mut _16 : f64 = 0.0; _1 = sin(c_roll) ; _2 = sin(c_pitch) ; _3 = sin(m_elev) ; _4 = cos(c_pitch) ; _5 = cos(m_elev) ; _6 = cos(c_view) ; _7 = cos(m_view) ; _8 = sin(c_view) ; _9 = sin(m_view) ; _10 = _8*_9+_6*_7 ; _11 = 1.0/(_4*_5*_10+_2*_3) ; _12 = _2*_5*_10-_4*_3 ; _13 = cos(c_roll) ; _14 = _6*_9-_8*_7 ; _15 = x0+_13*_5*_14*_11+_1*_11*_12 ; _16 = pow(_15,2.0)+pow(_13*_11*_12-_1*_5*_14*_11,2.0) ; return _15*(k0*sqrt(_16)+k1*_16+1.0) ; } +pub fn mac_x_dk0 ( c_view: f64, c_pitch: f64, c_roll: f64, scale: f64, k0: f64, k1: f64, x0: f64, m_view: f64, m_elev: f64 ) -> f64 { let mut _1 : f64 = 0.0; let mut _2 : f64 = 0.0; let mut _3 : f64 = 0.0; let mut _4 : f64 = 0.0; let mut _5 : f64 = 0.0; let mut _6 : f64 = 0.0; let mut _7 : f64 = 0.0; let mut _8 : f64 = 0.0; let mut _9 : f64 = 0.0; let mut _10 : f64 = 0.0; let mut _11 : f64 = 0.0; let mut _12 : f64 = 0.0; let mut _13 : f64 = 0.0; let mut _14 : f64 = 0.0; let mut _15 : f64 = 0.0; _1 = sin(c_roll) ; _2 = sin(c_pitch) ; _3 = sin(m_elev) ; _4 = cos(c_pitch) ; _5 = cos(m_elev) ; _6 = cos(c_view) ; _7 = cos(m_view) ; _8 = sin(c_view) ; _9 = sin(m_view) ; _10 = _8*_9+_6*_7 ; _11 = 1.0/(_4*_5*_10+_2*_3) ; _12 = _2*_5*_10-_4*_3 ; _13 = cos(c_roll) ; _14 = _6*_9-_8*_7 ; _15 = x0+_13*_5*_14*_11+_1*_11*_12 ; return _15*sqrt(pow(_15,2.0)+pow(_13*_11*_12-_1*_5*_14*_11,2.0))*scale ; } +pub fn mac_x_dk1 ( c_view: f64, c_pitch: f64, c_roll: f64, scale: f64, k0: f64, k1: f64, x0: f64, m_view: f64, m_elev: f64 ) -> f64 { let mut _1 : f64 = 0.0; let mut _2 : f64 = 0.0; let mut _3 : f64 = 0.0; let mut _4 : f64 = 0.0; let mut _5 : f64 = 0.0; let mut _6 : f64 = 0.0; let mut _7 : f64 = 0.0; let mut _8 : f64 = 0.0; let mut _9 : f64 = 0.0; let mut _10 : f64 = 0.0; let mut _11 : f64 = 0.0; let mut _12 : f64 = 0.0; let mut _13 : f64 = 0.0; let mut _14 : f64 = 0.0; let mut _15 : f64 = 0.0; _1 = sin(c_roll) ; _2 = sin(c_pitch) ; _3 = sin(m_elev) ; _4 = cos(c_pitch) ; _5 = cos(m_elev) ; _6 = cos(c_view) ; _7 = cos(m_view) ; _8 = sin(c_view) ; _9 = sin(m_view) ; _10 = _8*_9+_6*_7 ; _11 = 1.0/(_4*_5*_10+_2*_3) ; _12 = _2*_5*_10-_4*_3 ; _13 = cos(c_roll) ; _14 = _6*_9-_8*_7 ; _15 = x0+_13*_5*_14*_11+_1*_11*_12 ; return _15*(pow(_15,2.0)+pow(_13*_11*_12-_1*_5*_14*_11,2.0))*scale ; } +pub fn mac_x_dx0 ( c_view: f64, c_pitch: f64, c_roll: f64, scale: f64, k0: f64, k1: f64, x0: f64, m_view: f64, m_elev: f64 ) -> f64 { let mut _1 : f64 = 0.0; let mut _2 : f64 = 0.0; let mut _3 : f64 = 0.0; let mut _4 : f64 = 0.0; let mut _5 : f64 = 0.0; let mut _6 : f64 = 0.0; let mut _7 : f64 = 0.0; let mut _8 : f64 = 0.0; let mut _9 : f64 = 0.0; let mut _10 : f64 = 0.0; let mut _11 : f64 = 0.0; let mut _12 : f64 = 0.0; let mut _13 : f64 = 0.0; let mut _14 : f64 = 0.0; let mut _15 : f64 = 0.0; let mut _16 : f64 = 0.0; let mut _17 : f64 = 0.0; _1 = sin(c_roll) ; _2 = sin(c_pitch) ; _3 = sin(m_elev) ; _4 = cos(c_pitch) ; _5 = cos(m_elev) ; _6 = cos(c_view) ; _7 = cos(m_view) ; _8 = sin(c_view) ; _9 = sin(m_view) ; _10 = _8*_9+_6*_7 ; _11 = 1.0/(_4*_5*_10+_2*_3) ; _12 = _2*_5*_10-_4*_3 ; _13 = cos(c_roll) ; _14 = _6*_9-_8*_7 ; _15 = x0+_13*_5*_14*_11+_1*_11*_12 ; _16 = pow(_15,2.0)+pow(_13*_11*_12-_1*_5*_14*_11,2.0) ; _17 = sqrt(_16) ; return scale*(k0*_17+k1*_16+1.0)+scale*_15*((k0*_15)/_17+2.0*k1*_15) ; } +pub fn mac_y_dc_view ( c_view: f64, c_pitch: f64, c_roll: f64, scale: f64, k0: f64, k1: f64, x0: f64, m_view: f64, m_elev: f64 ) -> f64 { let mut _1 : f64 = 0.0; let mut _2 : f64 = 0.0; let mut _3 : f64 = 0.0; let mut _4 : f64 = 0.0; let mut _5 : f64 = 0.0; let mut _6 : f64 = 0.0; let mut _7 : f64 = 0.0; let mut _8 : f64 = 0.0; let mut _9 : f64 = 0.0; let mut _10 : f64 = 0.0; let mut _11 : f64 = 0.0; let mut _12 : f64 = 0.0; let mut _13 : f64 = 0.0; let mut _14 : f64 = 0.0; let mut _15 : f64 = 0.0; let mut _16 : f64 = 0.0; let mut _17 : f64 = 0.0; let mut _18 : f64 = 0.0; let mut _19 : f64 = 0.0; let mut _20 : f64 = 0.0; let mut _21 : f64 = 0.0; let mut _22 : f64 = 0.0; let mut _23 : f64 = 0.0; let mut _24 : f64 = 0.0; let mut _25 : f64 = 0.0; _1 = cos(c_roll) ; _2 = sin(c_pitch) ; _3 = sin(m_elev) ; _4 = cos(c_pitch) ; _5 = cos(m_elev) ; _6 = cos(c_view) ; _7 = cos(m_view) ; _8 = sin(c_view) ; _9 = sin(m_view) ; _10 = _8*_9+_6*_7 ; _11 = _4*_5*_10+_2*_3 ; _12 = 1.0/_11 ; _13 = _2*_5*_10-_4*_3 ; _14 = sin(c_roll) ; _15 = _6*_9-_8*_7 ; _16 = _1*_12*_13-_14*_5*_15*_12 ; _17 = 1.0/pow(_11,2.0) ; _18 = pow(_5,2.0) ; _19 = pow(_15,2.0) ; _20 = (-_8*_9)-_6*_7 ; _21 = (-_14*_5*_20*_12)+_2*_1*_5*_15*_12+_4*_14*_18*_19*_17-_4*_1*_5*_15*_17*_13 ; _22 = x0+_1*_5*_15*_12+_14*_12*_13 ; _23 = 2.0*(_1*_5*_20*_12+_2*_14*_5*_15*_12-_4*_1*_18*_19*_17-_4*_14*_5*_15*_17*_13)*_22+2.0*_16*_21 ; _24 = pow(_22,2.0)+pow(_16,2.0) ; _25 = sqrt(_24) ; return _21*scale*(k0*_25+k1*_24+1.0)+_16*scale*((k0*_23)/(2.0*_25)+k1*_23) ; } +pub fn mac_y_dc_pitch ( c_view: f64, c_pitch: f64, c_roll: f64, scale: f64, k0: f64, k1: f64, x0: f64, m_view: f64, m_elev: f64 ) -> f64 { let mut _1 : f64 = 0.0; let mut _2 : f64 = 0.0; let mut _3 : f64 = 0.0; let mut _4 : f64 = 0.0; let mut _5 : f64 = 0.0; let mut _6 : f64 = 0.0; let mut _7 : f64 = 0.0; let mut _8 : f64 = 0.0; let mut _9 : f64 = 0.0; let mut _10 : f64 = 0.0; let mut _11 : f64 = 0.0; let mut _12 : f64 = 0.0; let mut _13 : f64 = 0.0; let mut _14 : f64 = 0.0; let mut _15 : f64 = 0.0; let mut _16 : f64 = 0.0; let mut _17 : f64 = 0.0; let mut _18 : f64 = 0.0; let mut _19 : f64 = 0.0; let mut _20 : f64 = 0.0; let mut _21 : f64 = 0.0; let mut _22 : f64 = 0.0; let mut _23 : f64 = 0.0; _1 = cos(c_roll) ; _2 = sin(c_pitch) ; _3 = sin(m_elev) ; _4 = cos(c_pitch) ; _5 = cos(m_elev) ; _6 = cos(c_view) ; _7 = cos(m_view) ; _8 = sin(c_view) ; _9 = sin(m_view) ; _10 = _8*_9+_6*_7 ; _11 = _4*_5*_10+_2*_3 ; _12 = 1.0/_11 ; _13 = _2*_5*_10-_4*_3 ; _14 = sin(c_roll) ; _15 = _6*_9-_8*_7 ; _16 = _1*_12*_13-_14*_5*_15*_12 ; _17 = 1.0/pow(_11,2.0) ; _18 = _4*_3-_2*_5*_10 ; _19 = (-_1*_17*_18*_13)+_14*_5*_15*_17*_18+_1 ; _20 = x0+_1*_5*_15*_12+_14*_12*_13 ; _21 = 2.0*((-_14*_17*_18*_13)-_1*_5*_15*_17*_18+_14)*_20+2.0*_19*_16 ; _22 = pow(_20,2.0)+pow(_16,2.0) ; _23 = sqrt(_22) ; return _19*scale*(k0*_23+k1*_22+1.0)+_16*scale*((k0*_21)/(2.0*_23)+k1*_21) ; } +pub fn mac_y_dc_roll ( c_view: f64, c_pitch: f64, c_roll: f64, scale: f64, k0: f64, k1: f64, x0: f64, m_view: f64, m_elev: f64 ) -> f64 { let mut _1 : f64 = 0.0; let mut _2 : f64 = 0.0; let mut _3 : f64 = 0.0; let mut _4 : f64 = 0.0; let mut _5 : f64 = 0.0; let mut _6 : f64 = 0.0; let mut _7 : f64 = 0.0; let mut _8 : f64 = 0.0; let mut _9 : f64 = 0.0; let mut _10 : f64 = 0.0; let mut _11 : f64 = 0.0; let mut _12 : f64 = 0.0; let mut _13 : f64 = 0.0; let mut _14 : f64 = 0.0; let mut _15 : f64 = 0.0; let mut _16 : f64 = 0.0; let mut _17 : f64 = 0.0; let mut _18 : f64 = 0.0; let mut _19 : f64 = 0.0; let mut _20 : f64 = 0.0; _1 = cos(c_roll) ; _2 = sin(c_pitch) ; _3 = sin(m_elev) ; _4 = cos(c_pitch) ; _5 = cos(m_elev) ; _6 = cos(c_view) ; _7 = cos(m_view) ; _8 = sin(c_view) ; _9 = sin(m_view) ; _10 = _8*_9+_6*_7 ; _11 = 1.0/(_4*_5*_10+_2*_3) ; _12 = _2*_5*_10-_4*_3 ; _13 = sin(c_roll) ; _14 = _6*_9-_8*_7 ; _15 = _1*_11*_12-_13*_5*_14*_11 ; _16 = (-_1*_5*_14*_11)-_13*_11*_12 ; _17 = x0+_1*_5*_14*_11+_13*_11*_12 ; _18 = 2.0*_15*_17+2.0*_16*_15 ; _19 = pow(_17,2.0)+pow(_15,2.0) ; _20 = sqrt(_19) ; return _16*scale*(k0*_20+k1*_19+1.0)+_15*scale*((k0*_18)/(2.0*_20)+k1*_18) ; } +pub fn mac_y_dscale ( c_view: f64, c_pitch: f64, c_roll: f64, scale: f64, k0: f64, k1: f64, x0: f64, m_view: f64, m_elev: f64 ) -> f64 { let mut _1 : f64 = 0.0; let mut _2 : f64 = 0.0; let mut _3 : f64 = 0.0; let mut _4 : f64 = 0.0; let mut _5 : f64 = 0.0; let mut _6 : f64 = 0.0; let mut _7 : f64 = 0.0; let mut _8 : f64 = 0.0; let mut _9 : f64 = 0.0; let mut _10 : f64 = 0.0; let mut _11 : f64 = 0.0; let mut _12 : f64 = 0.0; let mut _13 : f64 = 0.0; let mut _14 : f64 = 0.0; let mut _15 : f64 = 0.0; let mut _16 : f64 = 0.0; _1 = cos(c_roll) ; _2 = sin(c_pitch) ; _3 = sin(m_elev) ; _4 = cos(c_pitch) ; _5 = cos(m_elev) ; _6 = cos(c_view) ; _7 = cos(m_view) ; _8 = sin(c_view) ; _9 = sin(m_view) ; _10 = _8*_9+_6*_7 ; _11 = 1.0/(_4*_5*_10+_2*_3) ; _12 = _2*_5*_10-_4*_3 ; _13 = sin(c_roll) ; _14 = _6*_9-_8*_7 ; _15 = _1*_11*_12-_13*_5*_14*_11 ; _16 = pow(x0+_1*_5*_14*_11+_13*_11*_12,2.0)+pow(_15,2.0) ; return _15*(k0*sqrt(_16)+k1*_16+1.0) ; } +pub fn mac_y_dk0 ( c_view: f64, c_pitch: f64, c_roll: f64, scale: f64, k0: f64, k1: f64, x0: f64, m_view: f64, m_elev: f64 ) -> f64 { let mut _1 : f64 = 0.0; let mut _2 : f64 = 0.0; let mut _3 : f64 = 0.0; let mut _4 : f64 = 0.0; let mut _5 : f64 = 0.0; let mut _6 : f64 = 0.0; let mut _7 : f64 = 0.0; let mut _8 : f64 = 0.0; let mut _9 : f64 = 0.0; let mut _10 : f64 = 0.0; let mut _11 : f64 = 0.0; let mut _12 : f64 = 0.0; let mut _13 : f64 = 0.0; let mut _14 : f64 = 0.0; let mut _15 : f64 = 0.0; _1 = cos(c_roll) ; _2 = sin(c_pitch) ; _3 = sin(m_elev) ; _4 = cos(c_pitch) ; _5 = cos(m_elev) ; _6 = cos(c_view) ; _7 = cos(m_view) ; _8 = sin(c_view) ; _9 = sin(m_view) ; _10 = _8*_9+_6*_7 ; _11 = 1.0/(_4*_5*_10+_2*_3) ; _12 = _2*_5*_10-_4*_3 ; _13 = sin(c_roll) ; _14 = _6*_9-_8*_7 ; _15 = _1*_11*_12-_13*_5*_14*_11 ; return _15*scale*sqrt(pow(x0+_1*_5*_14*_11+_13*_11*_12,2.0)+pow(_15,2.0)) ; } +pub fn mac_y_dk1 ( c_view: f64, c_pitch: f64, c_roll: f64, scale: f64, k0: f64, k1: f64, x0: f64, m_view: f64, m_elev: f64 ) -> f64 { let mut _1 : f64 = 0.0; let mut _2 : f64 = 0.0; let mut _3 : f64 = 0.0; let mut _4 : f64 = 0.0; let mut _5 : f64 = 0.0; let mut _6 : f64 = 0.0; let mut _7 : f64 = 0.0; let mut _8 : f64 = 0.0; let mut _9 : f64 = 0.0; let mut _10 : f64 = 0.0; let mut _11 : f64 = 0.0; let mut _12 : f64 = 0.0; let mut _13 : f64 = 0.0; let mut _14 : f64 = 0.0; let mut _15 : f64 = 0.0; _1 = cos(c_roll) ; _2 = sin(c_pitch) ; _3 = sin(m_elev) ; _4 = cos(c_pitch) ; _5 = cos(m_elev) ; _6 = cos(c_view) ; _7 = cos(m_view) ; _8 = sin(c_view) ; _9 = sin(m_view) ; _10 = _8*_9+_6*_7 ; _11 = 1.0/(_4*_5*_10+_2*_3) ; _12 = _2*_5*_10-_4*_3 ; _13 = sin(c_roll) ; _14 = _6*_9-_8*_7 ; _15 = _1*_11*_12-_13*_5*_14*_11 ; return _15*scale*(pow(x0+_1*_5*_14*_11+_13*_11*_12,2.0)+pow(_15,2.0)) ; } +pub fn mac_y_dx0 ( c_view: f64, c_pitch: f64, c_roll: f64, scale: f64, k0: f64, k1: f64, x0: f64, m_view: f64, m_elev: f64 ) -> f64 { let mut _1 : f64 = 0.0; let mut _2 : f64 = 0.0; let mut _3 : f64 = 0.0; let mut _4 : f64 = 0.0; let mut _5 : f64 = 0.0; let mut _6 : f64 = 0.0; let mut _7 : f64 = 0.0; let mut _8 : f64 = 0.0; let mut _9 : f64 = 0.0; let mut _10 : f64 = 0.0; let mut _11 : f64 = 0.0; let mut _12 : f64 = 0.0; let mut _13 : f64 = 0.0; let mut _14 : f64 = 0.0; let mut _15 : f64 = 0.0; let mut _16 : f64 = 0.0; _1 = cos(c_roll) ; _2 = sin(c_pitch) ; _3 = sin(m_elev) ; _4 = cos(c_pitch) ; _5 = cos(m_elev) ; _6 = cos(c_view) ; _7 = cos(m_view) ; _8 = sin(c_view) ; _9 = sin(m_view) ; _10 = _8*_9+_6*_7 ; _11 = 1.0/(_4*_5*_10+_2*_3) ; _12 = _2*_5*_10-_4*_3 ; _13 = sin(c_roll) ; _14 = _6*_9-_8*_7 ; _15 = _1*_11*_12-_13*_5*_14*_11 ; _16 = x0+_1*_5*_14*_11+_13*_11*_12 ; return _15*((_16*k0)/sqrt(pow(_16,2.0)+pow(_15,2.0))+2.0*k1*_16)*scale ; } diff --git a/crates/gipfel-core/src/projection/mod.rs b/crates/gipfel-core/src/projection/mod.rs new file mode 100644 index 0000000..65b74e5 --- /dev/null +++ b/crates/gipfel-core/src/projection/mod.rs @@ -0,0 +1,174 @@ +//! Projections: image coordinates as a function of azimuth/elevation angle, +//! with symbolically pre-generated derivatives (see `codegen/`). +//! +//! Faithful port of gipfel's `ProjectionLSQ.{H,cxx}`, `ProjectionRectilinear.*` +//! and `ProjectionCylindrical.*`. + +pub mod cylindrical; +pub mod rectilinear; + +pub const NAN: f64 = f64::NAN; + +/// The 7 fitted view parameters (radians where angular). +#[derive(Clone, Copy, Debug, Default, PartialEq)] +pub struct ViewParams { + /// azimuth of image center, radians + pub a_center: f64, + /// focal scale (pixels) + pub scale: f64, + /// pitch angle of image center (camera nose up), radians + pub a_pitch: f64, + /// camera roll, radians + pub a_roll: f64, + /// radial distortion coefficient 0 + pub k0: f64, + /// radial distortion coefficient 1 + pub k1: f64, + /// distortion center x, pixels + pub x0: f64, +} + +impl ViewParams { + pub fn new() -> Self { + ViewParams { + a_center: 0.0, + scale: 3500.0, + a_pitch: 0.0, + a_roll: 0.0, + k0: 0.0, + k1: 0.0, + x0: 0.0, + } + } +} + +/// The 9 arguments every generated macro function takes. +#[derive(Clone, Copy, Debug)] +pub struct ProjArgs { + pub c_view: f64, + pub c_pitch: f64, + pub c_roll: f64, + pub scale: f64, + pub k0: f64, + pub k1: f64, + pub x0: f64, + pub m_view: f64, + pub m_elev: f64, +} + +impl ProjArgs { + pub fn from_params(p: &ViewParams, m_view: f64, m_elev: f64) -> Self { + ProjArgs { + c_view: p.a_center, + c_pitch: p.a_pitch, + c_roll: p.a_roll, + scale: p.scale, + k0: p.k0, + k1: p.k1, + x0: p.x0, + m_view, + m_elev, + } + } +} + +#[derive(Clone, Copy, Debug, PartialEq, Eq)] +pub enum ProjectionType { + Rectilinear = 0, + Cylindrical = 1, +} + +/// A projection plus its symbolic partial derivatives w.r.t. the view +/// parameters. Methods default to `NAN`, matching the C++ base class. +pub trait Projection: Send + Sync { + /// half-width of the visible azimuth range, radians + fn get_view_angle(&self) -> f64; + + fn mac_x(&self, a: ProjArgs) -> f64 { + let _ = a; + NAN + } + fn mac_y(&self, a: ProjArgs) -> f64 { + let _ = a; + NAN + } + fn mac_x_dc_view(&self, a: ProjArgs) -> f64 { + let _ = a; + NAN + } + fn mac_x_dc_pitch(&self, a: ProjArgs) -> f64 { + let _ = a; + NAN + } + fn mac_x_dc_roll(&self, a: ProjArgs) -> f64 { + let _ = a; + NAN + } + fn mac_x_dscale(&self, a: ProjArgs) -> f64 { + let _ = a; + NAN + } + fn mac_x_dk0(&self, a: ProjArgs) -> f64 { + let _ = a; + NAN + } + fn mac_x_dk1(&self, a: ProjArgs) -> f64 { + let _ = a; + NAN + } + fn mac_x_dx0(&self, a: ProjArgs) -> f64 { + let _ = a; + NAN + } + fn mac_y_dc_view(&self, a: ProjArgs) -> f64 { + let _ = a; + NAN + } + fn mac_y_dc_pitch(&self, a: ProjArgs) -> f64 { + let _ = a; + NAN + } + fn mac_y_dc_roll(&self, a: ProjArgs) -> f64 { + let _ = a; + NAN + } + fn mac_y_dscale(&self, a: ProjArgs) -> f64 { + let _ = a; + NAN + } + fn mac_y_dk0(&self, a: ProjArgs) -> f64 { + let _ = a; + NAN + } + fn mac_y_dk1(&self, a: ProjArgs) -> f64 { + let _ = a; + NAN + } + fn mac_y_dx0(&self, a: ProjArgs) -> f64 { + let _ = a; + NAN + } + + /// Project azimuth/elevation to image coordinates (relative to image center), + /// normalizing azimuth around the view center like the C++ original. + fn get_coordinates(&self, alph: f64, elev: f64, parms: &ViewParams) -> (f64, f64) { + let mut alph = alph; + let pi = crate::hill::PI; + + if alph - parms.a_center > pi { + alph -= 2.0 * pi; + } else if alph - parms.a_center < -pi { + alph += 2.0 * pi; + } + + let args = ProjArgs::from_params(parms, alph, elev); + (self.mac_x(args), self.mac_y(args)) + } +} + +pub fn projection_from_type(t: ProjectionType) -> Box { + match t { + ProjectionType::Rectilinear => Box::new(rectilinear::ProjectionRectilinear), + ProjectionType::Cylindrical => Box::new(cylindrical::ProjectionCylindrical), + } +} diff --git a/crates/gipfel-core/src/projection/rectilinear.rs b/crates/gipfel-core/src/projection/rectilinear.rs new file mode 100644 index 0000000..a7037d6 --- /dev/null +++ b/crates/gipfel-core/src/projection/rectilinear.rs @@ -0,0 +1,104 @@ +//! Rectilinear (gnomonic) projection with lens distortion. + +mod funcs { + #![allow(non_snake_case, unused_parens, unused_variables, dead_code, clippy::all)] + + // math shims: the maxima-generated code uses C-style calls + #[inline(always)] + fn sin(x: f64) -> f64 { + x.sin() + } + #[inline(always)] + fn cos(x: f64) -> f64 { + x.cos() + } + #[inline(always)] + fn tan(x: f64) -> f64 { + x.tan() + } + #[inline(always)] + fn sqrt(x: f64) -> f64 { + x.sqrt() + } + #[inline(always)] + fn pow(x: f64, y: f64) -> f64 { + x.powf(y) + } + #[inline(always)] + fn exp(x: f64) -> f64 { + x.exp() + } + #[inline(always)] + fn fabs(x: f64) -> f64 { + x.abs() + } + #[inline(always)] + fn atan2(y: f64, x: f64) -> f64 { + y.atan2(x) + } + #[inline(always)] + fn sec(x: f64) -> f64 { + 1.0 / x.cos() + } + + include!("generated/rectilinear_funcs.rs"); +} + +use super::{ProjArgs, Projection}; + +pub struct ProjectionRectilinear; + +impl Projection for ProjectionRectilinear { + fn get_view_angle(&self) -> f64 { + 1.0471976 /* 60 deg */ + } + + fn mac_x(&self, a: ProjArgs) -> f64 { + funcs::mac_x(a.c_view, a.c_pitch, a.c_roll, a.scale, a.k0, a.k1, a.x0, a.m_view, a.m_elev) + } + fn mac_y(&self, a: ProjArgs) -> f64 { + funcs::mac_y(a.c_view, a.c_pitch, a.c_roll, a.scale, a.k0, a.k1, a.x0, a.m_view, a.m_elev) + } + fn mac_x_dc_view(&self, a: ProjArgs) -> f64 { + funcs::mac_x_dc_view(a.c_view, a.c_pitch, a.c_roll, a.scale, a.k0, a.k1, a.x0, a.m_view, a.m_elev) + } + fn mac_x_dc_pitch(&self, a: ProjArgs) -> f64 { + funcs::mac_x_dc_pitch(a.c_view, a.c_pitch, a.c_roll, a.scale, a.k0, a.k1, a.x0, a.m_view, a.m_elev) + } + fn mac_x_dc_roll(&self, a: ProjArgs) -> f64 { + funcs::mac_x_dc_roll(a.c_view, a.c_pitch, a.c_roll, a.scale, a.k0, a.k1, a.x0, a.m_view, a.m_elev) + } + fn mac_x_dscale(&self, a: ProjArgs) -> f64 { + funcs::mac_x_dscale(a.c_view, a.c_pitch, a.c_roll, a.scale, a.k0, a.k1, a.x0, a.m_view, a.m_elev) + } + fn mac_x_dk0(&self, a: ProjArgs) -> f64 { + funcs::mac_x_dk0(a.c_view, a.c_pitch, a.c_roll, a.scale, a.k0, a.k1, a.x0, a.m_view, a.m_elev) + } + fn mac_x_dk1(&self, a: ProjArgs) -> f64 { + funcs::mac_x_dk1(a.c_view, a.c_pitch, a.c_roll, a.scale, a.k0, a.k1, a.x0, a.m_view, a.m_elev) + } + fn mac_x_dx0(&self, a: ProjArgs) -> f64 { + funcs::mac_x_dx0(a.c_view, a.c_pitch, a.c_roll, a.scale, a.k0, a.k1, a.x0, a.m_view, a.m_elev) + } + fn mac_y_dc_view(&self, a: ProjArgs) -> f64 { + funcs::mac_y_dc_view(a.c_view, a.c_pitch, a.c_roll, a.scale, a.k0, a.k1, a.x0, a.m_view, a.m_elev) + } + fn mac_y_dc_pitch(&self, a: ProjArgs) -> f64 { + funcs::mac_y_dc_pitch(a.c_view, a.c_pitch, a.c_roll, a.scale, a.k0, a.k1, a.x0, a.m_view, a.m_elev) + } + fn mac_y_dc_roll(&self, a: ProjArgs) -> f64 { + funcs::mac_y_dc_roll(a.c_view, a.c_pitch, a.c_roll, a.scale, a.k0, a.k1, a.x0, a.m_view, a.m_elev) + } + fn mac_y_dscale(&self, a: ProjArgs) -> f64 { + funcs::mac_y_dscale(a.c_view, a.c_pitch, a.c_roll, a.scale, a.k0, a.k1, a.x0, a.m_view, a.m_elev) + } + fn mac_y_dk0(&self, a: ProjArgs) -> f64 { + funcs::mac_y_dk0(a.c_view, a.c_pitch, a.c_roll, a.scale, a.k0, a.k1, a.x0, a.m_view, a.m_elev) + } + fn mac_y_dk1(&self, a: ProjArgs) -> f64 { + funcs::mac_y_dk1(a.c_view, a.c_pitch, a.c_roll, a.scale, a.k0, a.k1, a.x0, a.m_view, a.m_elev) + } + fn mac_y_dx0(&self, a: ProjArgs) -> f64 { + funcs::mac_y_dx0(a.c_view, a.c_pitch, a.c_roll, a.scale, a.k0, a.k1, a.x0, a.m_view, a.m_elev) + } +} diff --git a/crates/gipfel-core/src/stitch.rs b/crates/gipfel-core/src/stitch.rs new file mode 100644 index 0000000..a311961 --- /dev/null +++ b/crates/gipfel-core/src/stitch.rs @@ -0,0 +1,323 @@ +//! Cylindrical panorama resampling — port of the C++ `Stitch.cxx` / +//! `ScanImage.cxx` stitch mode. +//! +//! The output is a cylindrical (equirectangular-in-azimuth, tan-in-elevation) +//! image spanning `view_start..view_end` degrees of azimuth. Each output +//! pixel is mapped back through every source image's fitted projection; the +//! first source that covers the direction contributes its (resampled) color. +//! +//! All sampling arithmetic reproduces the C implementation bit-for-bit: +//! colors are handled on the 0..=65025 scale (`value * 255`), nearest uses +//! `rint` (round-half-to-even), bicubic uses the same Keys cubic and only +//! rounds after the final vertical interpolation. + +use crate::metadata::{ImageMetaData, apply_image_metadata}; +use crate::panorama::Panorama; +use image::RgbImage; + +#[derive(Clone, Copy, PartialEq, Eq, Debug)] +pub enum ScanMode { + Nearest, + Bicubic, +} + +/// C MAX_VALUE: clamp bound on the 0..=65025 (value*255) scale. +pub const MAX_VALUE: i64 = 65025; + +fn pixel(img: &RgbImage, x: i64, y: i64) -> Option<[i64; 3]> { + if x < 0 || y < 0 || x >= img.width() as i64 || y >= img.height() as i64 { + return None; + } + let p = img.get_pixel(x as u32, y as u32); + Some([p[0] as i64 * 255, p[1] as i64 * 255, p[2] as i64 * 255]) +} + +/// C ScanImage::get_pixel_nearest: NaN-guard, then round to nearest pixel. +fn get_pixel_nearest(img: &RgbImage, x: f64, y: f64) -> Option<[i64; 3]> { + if x.is_nan() || y.is_nan() { + return None; + } + pixel(img, x.round_ties_even() as i64, y.round_ties_even() as i64) +} + +/// C ScanImage interp_cubic (Keys cubic convolution kernel, evaluated in +/// Bernstein-like power basis). +#[inline] +fn interp_cubic(x: f64, x2: f64, x3: f64, v: &[f64; 4]) -> f64 { + let a0 = v[3] - v[2] - v[0] + v[1]; + let a1 = v[0] - v[1] - a0; + let a2 = v[2] - v[0]; + let a3 = v[1]; + a0 * x3 + a1 * x2 + a2 * x + a3 +} + +/// C ScanImage::get_pixel_bicubic. Any of the 16 taps outside the source +/// makes the pixel unset (C returns 1), matching the "leave black" behavior. +fn get_pixel_bicubic(img: &RgbImage, x: f64, y: f64) -> Option<[i64; 3]> { + if x.is_nan() || y.is_nan() { + return None; + } + let fl_x = x.floor(); + let fl_y = y.floor(); + let dx = x - fl_x; + let dx2 = dx * dx; + let dx3 = dx2 * dx; + let dy = y - fl_y; + let dy2 = dy * dy; + let dy3 = dy2 * dy; + + let mut c1 = [[0f64; 4]; 3]; + for iy in 0..4usize { + let mut c = [[0f64; 4]; 3]; + for ix in 0..4usize { + let p = pixel(img, fl_x as i64 + ix as i64 - 1, fl_y as i64 + iy as i64 - 1)?; + for l in 0..3 { + c[l][ix] = p[l] as f64; + } + } + for l in 0..3 { + c1[l][iy] = interp_cubic(dx, dx2, dx3, &c[l]); + } + } + Some([ + interp_cubic(dy, dy2, dy3, &c1[0]).round_ties_even() as i64, + interp_cubic(dy, dy2, dy3, &c1[1]).round_ties_even() as i64, + interp_cubic(dy, dy2, dy3, &c1[2]).round_ties_even() as i64, + ]) +} + +/// C ScanImage::get_pixel(mode, x, y). Values on the 0..=65025 scale; +/// bicubic may overshoot (clamped by the caller, as in C). +pub fn get_pixel(img: &RgbImage, mode: ScanMode, x: f64, y: f64) -> Option<[i64; 3]> { + match mode { + ScanMode::Nearest => get_pixel_nearest(img, x, y), + ScanMode::Bicubic => get_pixel_bicubic(img, x, y), + } +} + +/// One input photo: its fitted projection plus its decoded pixels. +pub struct StitchSource { + pub pan: Panorama, + pub img: RgbImage, +} + +impl StitchSource { + /// Load a JPEG and apply its gipfel metadata (comment or EXIF), like + /// `GipfelWidget::load_image` does for the stitch mode. + pub fn load(path: &str) -> Result { + let img = image::open(path) + .map(|i| i.to_rgb8()) + .map_err(|e| format!("{}: {e}", path))?; + let mut md = ImageMetaData::new(); + md.load_image(path); + let mut pan = Panorama::new(); + apply_image_metadata(&mut pan, &md, img.width(), img.height()); + Ok(Self { pan, img }) + } + + /// C GipfelWidget::get_pixel: forward-project (azimuth, elevation) + /// through this image's projection and sample the source. + pub fn get_pixel(&self, mode: ScanMode, a_alph: f64, a_elev: f64) -> Option<[i64; 3]> { + let (px, py) = self.pan.project_point(a_alph, a_elev)?; + get_pixel( + &self.img, + mode, + px + self.img.width() as f64 / 2.0, + py + self.img.height() as f64 / 2.0, + ) + } +} + +/// The resampled cylindrical panorama. Pixel values are RGB triplets on the +/// C 0..=65025 scale (clamped); `hit` marks pixels some source contributed +/// to (unset pixels are black in the C output). +pub struct StitchedImage { + pub w: u32, + pub h: u32, + pub pixels: Vec<[u16; 3]>, + pub hit: Vec, +} + +/// Port of `Stitch::resample`: output column x maps linearly to azimuth, +/// output row y maps to elevation via `atan((y_off - y) / radius)`. +pub fn resample( + sources: &[StitchSource], + mode: ScanMode, + w: u32, + h: u32, + view_start_deg: f64, + view_end_deg: f64, +) -> StitchedImage { + let view_start = view_start_deg.to_radians(); + let view_end = view_end_deg.to_radians(); + let mut pixels = vec![[0u16; 3]; (w as usize) * (h as usize)]; + let mut hit = vec![false; (w as usize) * (h as usize)]; + if !(view_end > view_start) || w == 0 || h == 0 { + return StitchedImage { + w, + h, + pixels, + hit, + }; + } + let step_view = (view_end - view_start) / w as f64; + let y_off = (h / 2) as f64; + let radius = w as f64 / (view_end - view_start); + + for y in 0..h { + let a_elev = ((y_off - y as f64) / radius).atan(); + for x in 0..w { + let a_view = view_start + x as f64 * step_view; + for src in sources { + let Some(rgb) = src.get_pixel(mode, a_view, a_elev) else { + continue; + }; + let clamped = [ + rgb[0].clamp(0, MAX_VALUE) as u16, + rgb[1].clamp(0, MAX_VALUE) as u16, + rgb[2].clamp(0, MAX_VALUE) as u16, + ]; + let i = (y as usize) * (w as usize) + (x as usize); + pixels[i] = clamped; + hit[i] = true; + break; // C: first contributing image wins + } + } + } + + StitchedImage { + w, + h, + pixels, + hit, + } +} + +/// Write an RGBA TIFF like the C `TIFFOutputImage`: 8 bit (value/255, alpha +/// 255) or 16 bit (raw 0..=65025, alpha 65025), deflate-compressed. +pub fn write_tiff(path: &str, img: &StitchedImage, bits16: bool) -> Result<(), String> { + use tiff::encoder::colortype::{RGBA16, RGBA8}; + use tiff::encoder::{Compression, TiffEncoder}; + + let n = (img.w as usize) * (img.h as usize); + let file = std::fs::File::create(path).map_err(|e| format!("{}: {e}", path))?; + let mut enc = TiffEncoder::new(file) + .map_err(|e| format!("{}: {e}", path))? + .with_compression(Compression::Deflate(Default::default())); + if bits16 { + let mut buf: Vec = Vec::with_capacity(n * 4); + for p in &img.pixels { + buf.extend_from_slice(&[p[0], p[1], p[2], MAX_VALUE as u16]); + } + enc.write_image::(img.w, img.h, &buf) + .map_err(|e| format!("{}: {e}", path))?; + } else { + let mut buf: Vec = Vec::with_capacity(n * 4); + for p in &img.pixels { + buf.extend_from_slice(&[(p[0] / 255) as u8, (p[1] / 255) as u8, (p[2] / 255) as u8, 255]); + } + enc.write_image::(img.w, img.h, &buf) + .map_err(|e| format!("{}: {e}", path))?; + } + Ok(()) +} + +/// Write a JPEG (quality 90, like the C `JPEGOutputImage`). +pub fn write_jpeg(path: &str, img: &StitchedImage, quality: u8) -> Result<(), String> { + use image::codecs::jpeg::JpegEncoder; + use image::ImageEncoder; + + let rgb: Vec = img + .pixels + .iter() + .flat_map(|p| [(p[0] / 255) as u8, (p[1] / 255) as u8, (p[2] / 255) as u8]) + .collect(); + let enc = JpegEncoder::new_with_quality( + std::io::BufWriter::new(std::fs::File::create(path).map_err(|e| format!("{}: {e}", path))?), + quality, + ); + enc.write_image(&rgb, img.w, img.h, image::ExtendedColorType::Rgb8) + .map_err(|e| format!("{}: {e}", path))?; + Ok(()) +} + +#[cfg(test)] +mod tests { + use super::*; + + fn solid(v: u8, w: u32, h: u32) -> RgbImage { + RgbImage::from_pixel(w, h, image::Rgb([v, v, v])) + } + + #[test] + fn nearest_rounds_half_to_even() { + let img = solid(100, 4, 4); + // 1.5 -> 2 (ties to even), 2.5 -> 2 + assert_eq!(get_pixel_nearest(&img, 1.5, 0.0), Some([100 * 255; 3])); + assert_eq!(get_pixel_nearest(&img, 2.5, 0.0), Some([100 * 255; 3])); + assert_eq!(get_pixel_nearest(&img, -0.4, 0.0), Some([100 * 255; 3])); + assert_eq!(get_pixel_nearest(&img, -0.6, 0.0), None); + assert_eq!(get_pixel_nearest(&img, f64::NAN, 0.0), None); + } + + #[test] + fn bicubic_constant_and_bounds() { + let img = solid(100, 4, 4); + // constant image: cubic reproduces the constant exactly + assert_eq!(get_pixel_bicubic(&img, 1.3, 1.7), Some([100 * 255; 3])); + // any tap outside the 4x4 neighborhood -> unset + assert_eq!(get_pixel_bicubic(&img, 0.0, 2.0), None); // fl_x-1 = -1 + assert_eq!(get_pixel_bicubic(&img, 2.0, 3.5), None); // fl_y+2 = 5 + assert_eq!(get_pixel_bicubic(&img, 1.0, 1.0), Some([100 * 255; 3])); + } + + #[test] + fn bicubic_matches_keys_cubic_on_ramp() { + // Keys cubic (a=-1/2): on a linear ramp the kernel yields + // a + b*(2 dx^3 - 3 dx^2 + 2 dx), not the linear value. + let mut img = RgbImage::new(8, 8); + for x in 0..8 { + for y in 0..8 { + img.put_pixel(x, y, image::Rgb([(x as u8) * 30, 0, 0])); + } + } + let got = get_pixel_bicubic(&img, 3.25, 4.0).unwrap(); + let a = 3.0 * 30.0 * 255.0; + let b = 30.0 * 255.0; + let dx = 0.25f64; + let want = (a + b * (2.0 * dx * dx * dx - 3.0 * dx * dx + 2.0 * dx)).round_ties_even(); + assert_eq!(got[0], want as i64); + } + + #[test] + fn bicubic_exact_at_integer_positions() { + // at integer sample coordinates (dx = dy = 0) the cubic kernel + // reproduces the center tap exactly -> a 1:1 resample is lossless + let mut img = RgbImage::new(8, 8); + for x in 0..8 { + for y in 0..8 { + img.put_pixel(x, y, image::Rgb([(x * 7 + y * 3) as u8, (y * 5) as u8, 9])); + } + } + for y in 1..6 { + for x in 1..6 { + let p = img.get_pixel(x, y).0; + assert_eq!( + get_pixel_bicubic(&img, x as f64, y as f64).unwrap(), + [p[0] as i64 * 255, p[1] as i64 * 255, p[2] as i64 * 255] + ); + } + } + } + + #[test] + fn resample_rejects_empty_range() { + let img = StitchSource { + pan: Panorama::new(), + img: solid(50, 8, 8), + }; + let out = resample(&[img], ScanMode::Nearest, 4, 4, 10.0, 10.0); + assert_eq!(out.hit.len(), 16); + assert!(out.hit.iter().all(|h| !h)); + } +} diff --git a/crates/gipfel-core/src/track.rs b/crates/gipfel-core/src/track.rs new file mode 100644 index 0000000..5450aa8 --- /dev/null +++ b/crates/gipfel-core/src/track.rs @@ -0,0 +1,248 @@ +//! GPS track loading (GPX). +//! +//! Minimal streaming parser built directly on quick-xml: we only need the +//! `trkpt` positions, and real-world files (Garmin etc.) are full of vendor +//! `` blocks that stricter parsers choke on. Waypoints and +//! routes are ignored; only track points are read. +//! +//! Security posture for untrusted files: +//! - quick-xml is used with `default-features = false` (no serde, no encoding +//! detection); it never fetches external DTDs/entities and expands only the +//! five predefined XML entities, so XXE and entity-expansion ("billion +//! laughs") attacks are structurally impossible; unknown entity references +//! make the attribute fail to parse and the point is skipped +//! - input size is capped (MAX_GPX_BYTES) before anything is parsed +//! - the parser only reads lat/lon attributes and text; all other +//! elements (including vendor extensions) are ignored without being +//! interpreted + +use crate::hill::{Hill, FLAG_TRACK_POINT}; +use quick_xml::events::Event; + +/// Reject absurdly large track files (a 64 MiB GPX is already ~500k points). +pub const MAX_GPX_BYTES: u64 = 64 * 1024 * 1024; + +/// Load track points from a GPX file (all tracks/segments, in file order). +/// +/// Consecutive duplicate positions are dropped; the result is decimated with +/// a uniform stride to at most `max_points` (GPX logs easily reach 10k +/// points, which would bloat the panorama's hill list). +pub fn load_gpx(path: &str, max_points: usize) -> Result, String> { + let size = std::fs::metadata(path).map_err(|e| format!("{path}: {e}"))?; + if size.len() > MAX_GPX_BYTES { + return Err(format!("{path}: file too large ({} bytes)", size.len())); + } + let data = std::fs::read_to_string(path).map_err(|e| format!("{path}: {e}"))?; + let mut reader = quick_xml::Reader::from_str(&data); + reader.config_mut().trim_text(true); + + let mut pts: Vec<(f64, f64, f64)> = Vec::new(); // (lat, lon, ele) degrees + let mut depth: usize = 0; + let mut trkpt_depth: Option = None; + let mut in_ele = false; + let (mut lat, mut lon, mut ele) = (f64::NAN, f64::NAN, 0.0); + + fn push(pts: &mut Vec<(f64, f64, f64)>, lat: f64, lon: f64, ele: f64) { + if lat.is_nan() || lon.is_nan() { + return; + } + if let Some(&(lata, lona, _)) = pts.last() { + if lata == lat && lona == lon { + return; // consecutive duplicate fix + } + } + pts.push((lat, lon, ele)); + } + + loop { + match reader.read_event().map_err(|e| format!("{path}: {e}"))? { + Event::Eof => break, + Event::Empty(e) => { + if e.local_name().as_ref() == b"trkpt" { + let (la, lo) = read_latlon(&e); + lat = la; + lon = lo; + ele = 0.0; + push(&mut pts, lat, lon, ele); + } + } + Event::Start(e) => { + depth += 1; + match e.local_name().as_ref() { + b"trkpt" => { + trkpt_depth = Some(depth); + ele = 0.0; + let (la, lo) = read_latlon(&e); + lat = la; + lon = lo; + } + // accept only as a direct child of , so + // elevation-like tags inside vendor extensions are ignored + b"ele" => { + in_ele = trkpt_depth == Some(depth - 1); + if in_ele { + ele = 0.0; + } + } + _ => {} + } + } + Event::Text(t) => { + if in_ele { + if let Ok(s) = t.decode() { + if let Ok(v) = s.trim().parse() { + ele = v; + } + } + } + } + Event::End(e) => { + match e.local_name().as_ref() { + b"ele" => in_ele = false, + b"trkpt" => { + push(&mut pts, lat, lon, ele); + trkpt_depth = None; + } + _ => {} + } + depth = depth.saturating_sub(1); + } + _ => {} + } + } + + if pts.is_empty() { + return Err(format!("{path}: no track points")); + } + + let stride = pts.len().div_ceil(max_points.max(1)); + Ok(pts + .iter() + .step_by(stride) + .map(|&(lat, lon, ele)| { + let mut h = Hill::new("track", lat.to_radians(), lon.to_radians(), ele); + h.flags |= FLAG_TRACK_POINT; + h + }) + .collect()) +} + +fn read_latlon(e: &quick_xml::events::BytesStart) -> (f64, f64) { + let (mut lat, mut lon) = (f64::NAN, f64::NAN); + for a in e.attributes().flatten() { + if let Ok(v) = a.normalized_value(quick_xml::XmlVersion::Implicit1_0) { + match a.key.local_name().as_ref() { + b"lat" => lat = v.trim().parse().unwrap_or(f64::NAN), + b"lon" => lon = v.trim().parse().unwrap_or(f64::NAN), + _ => {} + } + } + } + (lat, lon) +} + +#[cfg(test)] +mod tests { + use super::*; + + fn tmp(name: &str, content: &str) -> std::path::PathBuf { + let dir = std::env::temp_dir().join(format!("gipfel-gpx-test-{name}")); + std::fs::create_dir_all(&dir).unwrap(); + let p = dir.join("t.gpx"); + std::fs::write(&p, content).unwrap(); + p + } + + #[test] + fn gpx_track_roundtrip() { + let gpx = r#" + + tour + 1000 + 1100 + 1100 + 1200 + +"#; + let p = tmp("rt", gpx); + let hills = load_gpx(p.to_str().unwrap(), 100).unwrap(); + assert_eq!(hills.len(), 3); // duplicate dropped + assert!(hills.iter().all(|h| h.is_track_point())); + assert!((hills[1].phi - 47.31f64.to_radians()).abs() < 1e-12); + assert!((hills[2].height - 1200.0).abs() < 1e-9); + } + + #[test] + fn gpx_vendor_extensions_ignored() { + // Garmin-style files are full of vendor extension blocks; they must + // not confuse the parser (and must not contribute elevation). + let gpx = r#" + + + T9999 + + + 1000 + 120 + + 1100 + + +"#; + let p = tmp("ext", gpx); + let hills = load_gpx(p.to_str().unwrap(), 100).unwrap(); + assert_eq!(hills.len(), 2); + assert!((hills[0].height - 1000.0).abs() < 1e-9); // not 9999 + } + + #[test] + fn gpx_decimation() { + let mut pts = String::new(); + for i in 0..1000 { + pts.push_str(&format!( + "1", + 47.0 + i as f64 * 1e-5, + 11.0 + i as f64 * 1e-5 + )); + } + let gpx = format!( + r#"{pts}"# + ); + let p = tmp("dec", &gpx); + let hills = load_gpx(p.to_str().unwrap(), 100).unwrap(); + assert!(hills.len() <= 100); + } + + #[test] + fn gpx_bom_is_tolerated() { + // Garmin Desktop App writes a UTF-8 BOM; it must not break parsing. + let gpx = format!( + "\u{feff}" + ); + let p = tmp("bom", &gpx); + let hills = load_gpx(p.to_str().unwrap(), 100).unwrap(); + assert_eq!(hills.len(), 1); + } + + #[test] + fn gpx_entity_injection_is_inert() { + // XXE / entity-expansion attempts: quick-xml only expands the five + // predefined entities; unknown references make the attribute fail to + // parse, so the point is simply skipped (never expanded, never fetched). + let gpx = r#" + ]> + + + 1000 + 1100 + 1200 + +"#; + let p = tmp("ent", gpx); + let hills = load_gpx(p.to_str().unwrap(), 100).unwrap(); + // only the clean point survives; the entity-bearing ones are skipped + assert_eq!(hills.len(), 1); + assert!((hills[0].phi - 47.32f64.to_radians()).abs() < 1e-12); + } +} diff --git a/crates/gipfel-core/tests/golden.rs b/crates/gipfel-core/tests/golden.rs new file mode 100644 index 0000000..fb12e40 --- /dev/null +++ b/crates/gipfel-core/tests/golden.rs @@ -0,0 +1,511 @@ +//! Golden-value tests against the original C++ implementation. +//! +//! Goldens are generated by `tests/gen_golden.sh` (needs the built C++ tree). +//! Run with `cargo test --test golden` (unit tests: `cargo test --lib`). + +use gipfel_core::hill; +use gipfel_core::metadata::ImageMetaData; +use gipfel_core::panorama::Panorama; +use gipfel_core::projection::{ + cylindrical::ProjectionCylindrical, rectilinear::ProjectionRectilinear, ProjArgs, Projection, ProjectionType, +}; + +fn golden(name: &str) -> Vec { + let path = format!( + "{}/../../tests/golden/{}", + env!("CARGO_MANIFEST_DIR"), + name + ); + std::fs::read_to_string(&path) + .unwrap_or_else(|_| panic!("missing golden file {} (run tests/gen_golden.sh)", path)) + .lines() + .map(|l| l.to_string()) + .collect() +} + +fn data_path() -> Option { + if let Ok(p) = std::env::var("GIPFEL_DATA") { + if std::path::Path::new(&p).exists() { + return Some(p); + } + } + let home = std::env::var("HOME").unwrap_or_default(); + let p = format!("{}/gipfel/gipfel.dat", home); + std::path::Path::new(&p).exists().then_some(p) +} + +/// relative-or-absolute closeness (libm differences between glibc and Rust) +fn close(a: f64, b: f64) -> bool { + if a.is_nan() && b.is_nan() { + return true; + } + if a.is_infinite() || b.is_infinite() { + return a == b; + } + (a - b).abs() <= 1e-9f64.max(a.abs() * 1e-9) +} + +#[test] +fn generated_functions_match_c() { + let golden = golden("funcs.txt"); + let rect = ProjectionRectilinear; + let cyl = ProjectionCylindrical; + + assert_eq!(golden.len(), 200); + + for (k, line) in golden.iter().enumerate() { + let kk = k as f64; + let c_view = -3.0 + 0.03 * kk; + let c_pitch = 0.10 * (0.11 * kk).sin(); + let c_roll = 0.05 * (0.07 * kk).cos(); + let scale = 2000.0 + 37.0 * (kk % 19.0); + let k0 = 0.001 * (0.31 * kk).sin(); + let k1 = 0.0001 * (0.17 * kk).cos(); + let x0 = 50.0 * (0.13 * kk).sin(); + let m_view = c_view + 0.5 * (0.23 * kk).sin(); + let m_elev = 0.05 * (0.19 * kk).cos(); + + let a = ProjArgs { c_view, c_pitch, c_roll, scale, k0, k1, x0, m_view, m_elev }; + + let mine: Vec = vec![ + rect.mac_x(a), + rect.mac_y(a), + rect.mac_x_dc_view(a), + rect.mac_x_dc_pitch(a), + rect.mac_x_dc_roll(a), + rect.mac_x_dscale(a), + rect.mac_x_dk0(a), + rect.mac_x_dk1(a), + rect.mac_x_dx0(a), + rect.mac_y_dc_view(a), + rect.mac_y_dc_pitch(a), + rect.mac_y_dc_roll(a), + rect.mac_y_dscale(a), + rect.mac_y_dk0(a), + rect.mac_y_dk1(a), + rect.mac_y_dx0(a), + cyl.mac_x(a), + cyl.mac_y(a), + cyl.mac_x_dc_view(a), + cyl.mac_x_dc_pitch(a), + cyl.mac_x_dc_roll(a), + cyl.mac_x_dscale(a), + cyl.mac_y_dc_view(a), + cyl.mac_y_dc_pitch(a), + cyl.mac_y_dc_roll(a), + cyl.mac_y_dscale(a), + ]; + + let theirs: Vec = line.split_whitespace().map(|t| t.parse().unwrap()).collect(); + assert_eq!(mine.len(), theirs.len(), "line {}", k); + + for (i, (m, t)) in mine.iter().zip(theirs.iter()).enumerate() { + assert!(close(*m, *t), "k={} fn#{}: rust {} != c {}", k, i, m, t); + } + } +} + +/// derivatives of the generated code vs. central finite differences +#[test] +fn derivatives_match_finite_differences() { + let rect = ProjectionRectilinear; + let base = ProjArgs { + c_view: 0.3, + c_pitch: 0.05, + c_roll: 0.02, + scale: 3000.0, + k0: 0.001, + k1: 0.0002, + x0: 12.0, + m_view: 0.5, + m_elev: 0.04, + }; + + let h = 1e-6f64; + + // (arg index, analytic derivative fn) + let dx: [(usize, &dyn Fn(ProjArgs) -> f64); 7] = [ + (0, &|a| rect.mac_x_dc_view(a)), + (1, &|a| rect.mac_x_dc_pitch(a)), + (2, &|a| rect.mac_x_dc_roll(a)), + (3, &|a| rect.mac_x_dscale(a)), + (4, &|a| rect.mac_x_dk0(a)), + (5, &|a| rect.mac_x_dk1(a)), + (6, &|a| rect.mac_x_dx0(a)), + ]; + let dy: [(usize, &dyn Fn(ProjArgs) -> f64); 7] = [ + (0, &|a| rect.mac_y_dc_view(a)), + (1, &|a| rect.mac_y_dc_pitch(a)), + (2, &|a| rect.mac_y_dc_roll(a)), + (3, &|a| rect.mac_y_dscale(a)), + (4, &|a| rect.mac_y_dk0(a)), + (5, &|a| rect.mac_y_dk1(a)), + (6, &|a| rect.mac_y_dx0(a)), + ]; + + let set = |a: ProjArgs, idx: usize, v: f64| -> ProjArgs { + let mut b = a; + match idx { + 0 => b.c_view = v, + 1 => b.c_pitch = v, + 2 => b.c_roll = v, + 3 => b.scale = v, + 4 => b.k0 = v, + 5 => b.k1 = v, + _ => b.x0 = v, + } + b + }; + + for (idx, f) in dx { + let plus = rect.mac_x(set(base, idx, base_val(&base, idx) + h)); + let minus = rect.mac_x(set(base, idx, base_val(&base, idx) - h)); + let fd = (plus - minus) / (2.0 * h); + let an = f(base); + assert!( + (an - fd).abs() <= 1e-4f64.max(an.abs() * 1e-5), + "dx deriv arg {}: analytic {} finite-diff {}", + idx, + an, + fd + ); + } + for (idx, f) in dy { + let plus = rect.mac_y(set(base, idx, base_val(&base, idx) + h)); + let minus = rect.mac_y(set(base, idx, base_val(&base, idx) - h)); + let fd = (plus - minus) / (2.0 * h); + let an = f(base); + assert!( + (an - fd).abs() <= 1e-4f64.max(an.abs() * 1e-5), + "dy deriv arg {}: analytic {} finite-diff {}", + idx, + an, + fd + ); + } +} + +fn base_val(a: &ProjArgs, idx: usize) -> f64 { + match idx { + 0 => a.c_view, + 1 => a.c_pitch, + 2 => a.c_roll, + 3 => a.scale, + 4 => a.k0, + 5 => a.k1, + _ => a.x0, + } +} + +/// full fit (geodesy + projection + GSL LM) against the C++ fitdump +fn fit_golden(projtype: ProjectionType, file: &str) { + let Some(data) = data_path() else { + eprintln!("gipfel.dat not found, skipping"); + return; + }; + let golden = golden(file); + + let mut pan = Panorama::new(); + pan.load_data(&data).unwrap(); + pan.set_projection(projtype); + + pan.set_view_lat(47.36); + pan.set_view_long(11.15); + pan.set_view_height(2100.0); + + // ground truth + pan.set_center_angle(20.0); + pan.set_pitch_angle(0.0); + pan.set_roll_angle(1.0); + pan.set_scale(3000.0); + + // snapshot first 6 visible hills as observations + let mut obs: Vec = pan.visible[..6].iter().map(|&i| pan.hills[i].clone()).collect(); + let n = obs.len(); + let known: Vec = (0..n).collect(); + // the C harness copies a_elev explicitly (C copy ctor omits it) + for (k, &i) in pan.visible[..6].iter().enumerate() { + obs[k].a_elev = pan.hills[i].a_elev; + } + + // perturb initial guess + pan.set_center_angle(pan.get_center_angle() + 0.5); + pan.set_pitch_angle(0.3); + pan.set_roll_angle(0.5); + pan.set_scale(2800.0); + + let ret = pan.comp_params_with(&obs, &known, true); + + let g_ret: i32 = golden[0]["ret ".len()..].parse().unwrap(); + assert_eq!(ret.is_ok(), g_ret == 0); + + let g: Vec = golden[1..8] + .iter() + .map(|l| l.split_whitespace().nth(1).unwrap().parse().unwrap()) + .collect(); + + let mine = [ + pan.parms.a_center, + pan.parms.a_pitch, + pan.parms.a_roll, + pan.parms.scale, + pan.parms.k0, + pan.parms.k1, + pan.parms.x0, + ]; + for (m, t) in mine.iter().zip(g.iter()) { + // solver results: allow small drift, exact zeros must stay ~zero + assert!( + (m - t).abs() <= 1e-6f64.max(m.abs() * 1e-6), + "fit param {} vs {}", + m, + t + ); + } + + // final coordinates + let proj = gipfel_core::projection::projection_from_type(projtype); + for line in &golden[8..] { + let parts: Vec<&str> = line.split_whitespace().collect(); + if parts.len() != 4 { + continue; + } + let i: usize = parts[1].parse().unwrap(); + let tx: f64 = parts[2].parse().unwrap(); + let ty: f64 = parts[3].parse().unwrap(); + let (x, y) = proj.get_coordinates(obs[i].alph, obs[i].a_elev, &pan.parms); + assert!(close(x, tx), "hill {} x: {} vs {}", i, x, tx); + assert!(close(y, ty), "hill {} y: {} vs {}", i, y, ty); + } +} + +#[test] +fn fit_rectilinear_matches_c() { + fit_golden(ProjectionType::Rectilinear, "fit_rect.txt"); +} + +#[test] +fn fit_cylindrical_matches_c() { + fit_golden(ProjectionType::Cylindrical, "fit_cyl.txt"); +} + +#[test] +fn real_database_loads() { + let Some(data) = data_path() else { + eprintln!("gipfel.dat not found, skipping"); + return; + }; + let hills = hill::load(&data).unwrap(); + assert!(hills.len() > 21000, "loaded {}", hills.len()); + assert!(hills.iter().any(|h| h.name == "ZUGSPITZE")); +} + +#[test] +fn viewpoint_geometry() { + let Some(data) = data_path() else { + eprintln!("gipfel.dat not found, skipping"); + return; + }; + let mut pan = Panorama::new(); + pan.load_data(&data).unwrap(); + let idx = pan.get_pos("ZUGSPITZE").expect("ZUGSPITZE in db"); + pan.set_viewpoint(Some(idx)); + + // NOTE: set_viewpoint re-sorts hills by azimuth; look the name up again + let idx = pan.get_pos("ZUGSPITZE").unwrap(); + let m = &pan.hills[idx]; + assert!((m.height - 2963.0).abs() < 1.0); + // distance to self is zero + assert!(pan.distance(m.phi, m.lam).abs() < 1e-9); + // a far hill: Mont Blanc (~45.83 N, 6.86 E) + let mb = pan.get_pos("MONT BLANC").expect("MONT BLANC in db"); + let d = pan.distance(pan.hills[mb].phi, pan.hills[mb].lam); + let km = d * 6371000.785 / 1000.0; + assert!(km > 340.0 && km < 380.0, "Mont Blanc distance {} km", km); +} + +#[test] +fn export_hills_matches_c() { + let Some(data) = data_path() else { + eprintln!("gipfel.dat not found, skipping"); + return; + }; + let jpg = format!( + "{}/../../tests/golden/test.jpg", + env!("CARGO_MANIFEST_DIR") + ); + let golden = golden("export.txt"); + + let (w, h) = gipfel_core::metadata::image_dimensions(&jpg).expect("test.jpg readable"); + let mut md = ImageMetaData::new(); + md.load_image(&jpg); + + let mut pan = Panorama::new(); + assert!(gipfel_core::metadata::apply_image_metadata(&mut pan, &md, w, h)); + pan.load_data(&data).unwrap(); + pan.set_height_dist_ratio(0.07); + + let mut out = vec![ + "#".to_string(), + "# name\theight\tx\ty\tdistance\tflags".to_string(), + "#".to_string(), + ]; + out.extend(pan.export_hill_rows(w as i64, h as i64, false)); + + assert_eq!(out.len(), golden.len()); + for (a, b) in out.iter().zip(golden.iter()) { + assert_eq!(a, b); + } +} + +#[test] +fn metadata_comment_roundtrip_and_c_format() { + let mut md = ImageMetaData::new(); + md.longitude = 11.1533; + md.latitude = 47.3601; + md.height = 2100.5; + md.direction = 20.0; + md.pitch = 0.3; + md.roll = 1.0; + md.focal_length_35mm = 35.0; + md.projection_type = 1; + md.k0 = 0.001; + md.k1 = -0.0002; + md.x0 = 12.0; + + let s = md.format_comment(); + // exact C printf("%f") compatibility + assert_eq!( + s, + "gipfel: longitude 11.153300, latitude 47.360100, height 2100.500000, direction 20.000000, nick 0.300000, tilt 1.000000, focal_length_35mm 35.000000, projection type 1, k0 0.001000, k1 -0.000200, x0 12.000000" + ); + + let mut md2 = ImageMetaData::new(); + assert!(md2.parse_gipfel_comment(s.as_bytes())); + assert!((md2.longitude - 11.1533).abs() < 1e-9); + assert!((md2.latitude - 47.3601).abs() < 1e-9); + assert!((md2.height - 2100.5).abs() < 1e-9); + assert!((md2.direction - 20.0).abs() < 1e-9); + assert_eq!(md2.projection_type, 1); + assert!((md2.k0 - 0.001).abs() < 1e-9); + assert!((md2.k1 + 0.0002).abs() < 1e-9); + assert!((md2.x0 - 12.0).abs() < 1e-9); + + // short form (8 fields, pre-distortion files): k0/k1/x0 stay NaN + let short = b"gipfel: longitude 1.000000, latitude 2.000000, height 3.000000, direction 4.000000, nick 5.000000, tilt 6.000000, focal_length_35mm 7.000000, projection type 0"; + let mut md3 = ImageMetaData::new(); + assert!(md3.parse_gipfel_comment(short)); + assert!(md3.k0.is_nan()); + assert_eq!(md3.projection_type, 0); + + // trailing garbage after x0 must not break parsing (sscanf semantics) + let mut md4 = ImageMetaData::new(); + let mut s4 = s.clone(); + s4.push_str(" some other comment junk"); + assert!(md4.parse_gipfel_comment(s4.as_bytes())); + assert!((md4.x0 - 12.0).abs() < 1e-9); +} + +/// The least-squares fit must run on a small stack: the LM solver keeps its +/// working buffers (residuals, Jacobian, normal equations) on the heap, so a +/// 128 KiB thread stack is plenty. Guards against regressions that move those +/// buffers back onto the stack. +#[test] +fn fit_survives_small_stack() { + let Some(data) = data_path() else { + eprintln!("gipfel.dat not found, skipping"); + return; + }; + std::thread::Builder::new() + .stack_size(128 * 1024) + .spawn(move || { + let mut pan = Panorama::new(); + pan.load_data(&data).unwrap(); + pan.set_view_lat(47.36); + pan.set_view_long(11.15); + pan.set_view_height(2100.0); + pan.set_center_angle(20.0); + pan.set_pitch_angle(0.0); + pan.set_roll_angle(1.0); + pan.set_scale(3000.0); + let known: Vec = pan.visible[..4].to_vec(); + // perturb so the solver actually iterates + pan.set_center_angle(20.5); + pan.set_scale(2800.0); + pan.comp_params(&known, true).unwrap(); + }) + .unwrap() + .join() + .unwrap(); +} + +/// Stitch resample against the C Stitch::resample + ScanImage pipeline. +fn stitch_golden( + name: &str, + proj: ProjectionType, + center: f64, + pitch: f64, + roll: f64, + scale: f64, + k0: f64, + k1: f64, + bicubic: bool, +) { + use gipfel_core::stitch::{self, ScanMode}; + + let ppm = format!( + "{}/../../tests/golden/stitch.ppm", + env!("CARGO_MANIFEST_DIR") + ); + let img = image::open(&ppm) + .unwrap_or_else(|e| panic!("missing {}: {e}", ppm)) + .to_rgb8(); + + let mut pan = Panorama::new(); + pan.set_projection(proj); + pan.set_center_angle(center); + pan.set_pitch_angle(pitch); + pan.set_roll_angle(roll); + pan.set_scale(scale); + pan.set_distortion_params(k0, k1, 0.0); + + let src = stitch::StitchSource { pan, img }; + let mode = if bicubic { ScanMode::Bicubic } else { ScanMode::Nearest }; + let out = stitch::resample(&[src], mode, 64, 48, 10.0, 30.0); + + let mut checked = 0; + for line in golden(name) { + let f: Vec<&str> = line.split_whitespace().collect(); + let (x, y, hit): (usize, usize, i32) = (f[0].parse().unwrap(), f[1].parse().unwrap(), f[2].parse().unwrap()); + let want = [ + f[3].parse::().unwrap(), + f[4].parse::().unwrap(), + f[5].parse::().unwrap(), + ]; + let i = y * 64 + x; + assert_eq!(out.hit[i], hit == 1, "{name} pixel ({x},{y}) hit flag"); + assert_eq!(out.pixels[i], want, "{name} pixel ({x},{y})"); + checked += 1; + } + assert!(checked > 1000); +} + +#[test] +fn stitch_rect_nearest_matches_c() { + stitch_golden("stitch_rect_near.txt", ProjectionType::Rectilinear, 20.0, 2.0, 0.0, 200.0, 0.0, 0.0, false); +} + +#[test] +fn stitch_rect_bicubic_matches_c() { + stitch_golden("stitch_rect_bic.txt", ProjectionType::Rectilinear, 20.0, 2.0, 0.0, 200.0, 0.0, 0.0, true); +} + +#[test] +fn stitch_cyl_bicubic_matches_c() { + stitch_golden("stitch_cyl_bic.txt", ProjectionType::Cylindrical, 20.0, 2.0, 0.0, 200.0, 0.0, 0.0, true); +} + +#[test] +fn stitch_distorted_bicubic_matches_c() { + stitch_golden("stitch_dist_bic.txt", ProjectionType::Rectilinear, 20.0, 2.0, 3.0, 200.0, -0.3, 0.1, true); +} -- cgit v1.2.3