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.

1 files changed, 323 insertions, 0 deletions

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<Self, String> {
+        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<bool>,
+}
+
+/// 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<u16> = 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::<RGBA16>(img.w, img.h, &buf)
+            .map_err(|e| format!("{}: {e}", path))?;
+    } else {
+        let mut buf: Vec<u8> = 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::<RGBA8>(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<u8> = 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));
+    }
+}