diff options
Diffstat (limited to 'crates/gipfel-core/tests')
| -rw-r--r-- | crates/gipfel-core/tests/golden.rs | 511 |
1 files changed, 511 insertions, 0 deletions
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<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); +} |