//! Numeric parity against tomopy (M1 golden-data gate).
//!
//! Loads checked-in golden data generated by `tools/gen_tomopy_golden.py`
//! (tomopy 1.15.3): a fixed Shepp-Logan phantom, its parallel-beam sinogram,
//! and tomopy's own `fbp`/`gridrec` reconstructions of that sinogram. tomoxide
//! reconstructs the SAME sinogram and is compared to those references. Both
//! sides therefore invert identical measured data, isolating the reconstruction
//! algorithm from the projector. The golden `.npy` is committed, so the test
//! runs offline — only regenerating the golden needs a tomopy install.
//!
//! Orientation: tomopy's reconstruction grid is a vertical (row-axis) flip of
//! tomoxide's — the two differ only by the grid↔detector handedness, a single
//! fixed `iy → n-1-iy` re-indexing (`ROW_FLIP` below), the same for every
//! algorithm.
//!
//! FBP reference: tomopy's own `fbp` is NOT a usable numeric reference on this
//! data — it correlates only r≈0.47 with its own phantom (in every orientation)
//! and r≈0.48 with its own gridrec, with a large DC offset. tomopy documents
//! `gridrec` as its analytic method; `fbp` is a basic path. So gridrec is gated
//! against tomopy's gridrec (true cross-implementation parity), and FBP is gated
//! against the phantom — i.e. tomoxide must recover the phantom from tomopy's
//! foreign sinogram, which is a stronger check than the self-round-trip because
//! it crosses projector conventions.

use ndarray::{Array1, Array2, Array3, Axis};
use ndarray_npy::read_npy;
use tomoxide::{
    recon, Algorithm, Angles, CpuBackend, FilterName, Geometry, Layout, ReconParams, Tomo,
};

const FIXTURES: &str = concat!(env!("CARGO_MANIFEST_DIR"), "/tests/fixtures");

/// tomopy's grid is a vertical flip of tomoxide's (grid↔detector handedness).
fn row_flip(a: &Array2<f32>) -> Array2<f32> {
    let n = a.dim().0;
    Array2::from_shape_fn((n, n), |(iy, ix)| a[[n - 1 - iy, ix]])
}

/// Pearson correlation over a centered disk, amplitude-scale invariant.
fn pearson_disk(a: &Array2<f32>, b: &Array2<f32>, n: usize, radius_frac: f32) -> f32 {
    let c = (n as f32 - 1.0) / 2.0;
    let r2 = (radius_frac * n as f32 / 2.0).powi(2);
    let (mut xs, mut ys) = (Vec::new(), Vec::new());
    for iy in 0..n {
        for ix in 0..n {
            let (dy, dx) = (iy as f32 - c, ix as f32 - c);
            if dx * dx + dy * dy <= r2 {
                xs.push(a[[iy, ix]]);
                ys.push(b[[iy, ix]]);
            }
        }
    }
    let nn = xs.len() as f32;
    let mx = xs.iter().sum::<f32>() / nn;
    let my = ys.iter().sum::<f32>() / nn;
    let (mut sxy, mut sxx, mut syy) = (0.0f32, 0.0f32, 0.0f32);
    for (&x, &y) in xs.iter().zip(ys.iter()) {
        let (dx, dy) = (x - mx, y - my);
        sxy += dx * dy;
        sxx += dx * dx;
        syy += dy * dy;
    }
    sxy / (sxx.sqrt() * syy.sqrt())
}

/// Reconstruct the golden sinogram with tomoxide, aligned to tomopy's grid.
fn tomoxide_recon(algorithm: Algorithm) -> (Array2<f32>, usize) {
    let sino: Array3<f32> = read_npy(format!("{FIXTURES}/sino.npy")).unwrap(); // (nang,1,ndet)
    let angles: Array1<f32> = read_npy(format!("{FIXTURES}/angles.npy")).unwrap();
    let n = sino.dim().2;

    let tomo = Tomo::new(sino, Layout::Projection);
    let geom = Geometry::parallel(Angles(angles.to_vec()), n, 1, 1.0);
    let params = ReconParams {
        num_gridx: Some(n),
        // Pin ramp: this test checks algorithm recovery against the phantom,
        // calibrated for the sharp filter (the default is parzen).
        filter_name: FilterName::Ramp,
        ..Default::default()
    };
    let vol = recon::recon(&tomo, &geom, algorithm, &params, &CpuBackend::new()).unwrap();
    let slice = vol.array.index_axis(Axis(0), 0).to_owned();
    (row_flip(&slice), n)
}

#[test]
fn gridrec_matches_tomopy_gridrec() {
    let golden: Array2<f32> = read_npy(format!("{FIXTURES}/tomopy_gridrec.npy")).unwrap();
    let (rec, n) = tomoxide_recon(Algorithm::Gridrec);
    let r = pearson_disk(&rec, &golden, n, 0.85);
    eprintln!("gridrec vs tomopy gridrec: Pearson r = {r:.4}");
    assert!(
        r >= 0.95,
        "tomoxide gridrec diverges from tomopy gridrec: r = {r:.4}"
    );
}

#[test]
fn fbp_recovers_phantom_from_tomopy_sinogram() {
    let phantom: Array2<f32> = read_npy(format!("{FIXTURES}/phantom.npy")).unwrap();
    let (rec, n) = tomoxide_recon(Algorithm::Fbp);
    let r = pearson_disk(&rec, &phantom, n, 0.85);
    eprintln!("fbp (tomopy sinogram) vs phantom: Pearson r = {r:.4}");
    assert!(
        r >= 0.85,
        "tomoxide FBP fails to recover the phantom from tomopy's sinogram: r = {r:.4}"
    );
}
