path: root/crates/gipfel-core/tests/golden.rs

//! 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<String> {
    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<String> {
    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<f64> = 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<f64> = 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<hill::Hill> = pan.visible[..6].iter().map(|&i| pan.hills[i].clone()).collect();
    let n = obs.len();
    let known: Vec<usize> = (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<f64> = 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<usize> = 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::<u16>().unwrap(),
            f[4].parse::<u16>().unwrap(),
            f[5].parse::<u16>().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);
}