oxiproj-engine 0.1.1

//! Distortion-field rasters: grids of Tissot indicatrix values over a geographic extent.

use crate::pj::Pj;
use oxiproj_core::{Coord, IoUnits, ProjError, TissotParams};

/// Whether Tissot values came from exact AD or numeric finite-difference.
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub enum TissotMethod {
    Exact,
    Numeric,
    Failed,
}

/// One cell in the distortion raster.
#[derive(Debug, Clone)]
pub struct TissotCell {
    pub lon: f64,
    pub lat: f64,
    pub semi_major: f64,
    pub semi_minor: f64,
    pub angular_distortion: f64,
    pub areal_scale: f64,
    pub meridian_scale: f64,
    pub parallel_scale: f64,
    pub meridian_convergence: f64,
    pub is_conformal: bool,
    pub is_equal_area: bool,
    pub method: TissotMethod,
}

/// A distortion-field raster: Tissot parameters on a regular lon/lat grid.
pub struct DistortionRaster {
    pub lon_min: f64,
    pub lon_max: f64,
    pub lat_min: f64,
    pub lat_max: f64,
    pub ncols: usize,
    pub nrows: usize,
    /// Flattened row-major grid of TissotCell (len = nrows * ncols).
    pub cells: Vec<TissotCell>,
}

impl DistortionRaster {
    /// Compute a distortion raster for the given projection over the given extent.
    ///
    /// All lon/lat values are in degrees. The raster has:
    /// - `ncols = ((lon_max - lon_min) / step_deg).ceil() as usize + 1`
    /// - `nrows = ((lat_max - lat_min) / step_deg).ceil() as usize + 1`
    pub fn compute(
        pj: &Pj,
        lon_min: f64,
        lat_min: f64,
        lon_max: f64,
        lat_max: f64,
        step_deg: f64,
    ) -> Self {
        let ncols = ((lon_max - lon_min) / step_deg).ceil() as usize + 1;
        let nrows = ((lat_max - lat_min) / step_deg).ceil() as usize + 1;
        let use_radians = pj.left == IoUnits::Radians;

        let mut cells = Vec::with_capacity(nrows * ncols);

        for r in 0..nrows {
            let lat = lat_min + r as f64 * step_deg;
            for c in 0..ncols {
                let lon = lon_min + c as f64 * step_deg;

                let coord = if use_radians {
                    Coord::new(lon.to_radians(), lat.to_radians(), 0.0, 0.0)
                } else {
                    Coord::new(lon, lat, 0.0, 0.0)
                };

                let cell = compute_cell(pj, lon, lat, coord);
                cells.push(cell);
            }
        }

        DistortionRaster {
            lon_min,
            lon_max,
            lat_min,
            lat_max,
            ncols,
            nrows,
            cells,
        }
    }

    /// Get a cell by row and column index.
    pub fn get(&self, row: usize, col: usize) -> Option<&TissotCell> {
        if row < self.nrows && col < self.ncols {
            Some(&self.cells[row * self.ncols + col])
        } else {
            None
        }
    }
}

fn failed_cell(lon: f64, lat: f64) -> TissotCell {
    TissotCell {
        lon,
        lat,
        semi_major: f64::NAN,
        semi_minor: f64::NAN,
        angular_distortion: f64::NAN,
        areal_scale: f64::NAN,
        meridian_scale: f64::NAN,
        parallel_scale: f64::NAN,
        meridian_convergence: f64::NAN,
        is_conformal: false,
        is_equal_area: false,
        method: TissotMethod::Failed,
    }
}

fn tissot_to_cell(lon: f64, lat: f64, t: TissotParams, method: TissotMethod) -> TissotCell {
    TissotCell {
        lon,
        lat,
        semi_major: t.semi_major,
        semi_minor: t.semi_minor,
        angular_distortion: t.angular_distortion,
        areal_scale: t.areal_scale,
        meridian_scale: t.meridian_scale,
        parallel_scale: t.parallel_scale,
        meridian_convergence: t.meridian_convergence,
        is_conformal: t.is_conformal,
        is_equal_area: t.is_equal_area,
        method,
    }
}

fn compute_cell(pj: &Pj, lon: f64, lat: f64, coord: Coord) -> TissotCell {
    match pj.factors_exact(coord) {
        Ok(fe) => {
            let t = TissotParams::from_factors(&fe);
            tissot_to_cell(lon, lat, t, TissotMethod::Exact)
        }
        Err(ProjError::UnsupportedOperation) => {
            // Fall back to numeric finite-difference
            match pj.factors(coord) {
                Ok(f) => {
                    // Factors struct fields: h (meridian scale), k (parallel scale),
                    // s (areal scale), conv (meridian convergence), theta_prime
                    let t = TissotParams::from_exact(f.h, f.k, f.s, f.theta_prime, f.conv);
                    tissot_to_cell(lon, lat, t, TissotMethod::Numeric)
                }
                Err(_) => failed_cell(lon, lat),
            }
        }
        Err(_) => failed_cell(lon, lat),
    }
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::create::create;

    #[test]
    fn mercator_3x3_equator() {
        let pj = create("+proj=merc +datum=WGS84").unwrap();
        let raster = DistortionRaster::compute(&pj, -5.0, -5.0, 5.0, 5.0, 5.0);
        assert_eq!(raster.ncols, 3);
        assert_eq!(raster.nrows, 3);

        // All cells must produce a valid (non-Failed) result.
        // Exact is returned when an AD-capable projection is registered;
        // Numeric is returned via finite-difference fallback.
        // Both are correct; the raster must never be Failed.
        for cell in &raster.cells {
            assert_ne!(
                cell.method,
                TissotMethod::Failed,
                "expected non-Failed at lon={} lat={}",
                cell.lon,
                cell.lat
            );
        }

        // Center cell (row=1, col=1) at (0,0): conformal, angular_distortion ≈ 0
        let center = raster.get(1, 1).expect("center cell");
        assert!(
            center.angular_distortion.abs() < 1e-6,
            "angular_distortion at origin should be ~0, got {}",
            center.angular_distortion
        );
        assert!(
            (center.semi_major - center.semi_minor).abs() < 1e-4,
            "semi_major ≈ semi_minor for conformal at origin, got {} vs {}",
            center.semi_major,
            center.semi_minor
        );
    }

    #[test]
    fn europe_5x5_raster() {
        let pj = create("+proj=merc +datum=WGS84").unwrap();
        let raster = DistortionRaster::compute(&pj, -10.0, 35.0, 30.0, 75.0, 10.0);
        assert_eq!(raster.ncols, 5);
        assert_eq!(raster.nrows, 5);

        for cell in &raster.cells {
            assert_ne!(
                cell.method,
                TissotMethod::Failed,
                "cell at lon={} lat={} should not fail",
                cell.lon,
                cell.lat
            );
            // Mercator is conformal so angular_distortion ≈ 0, but general check: >= 0
            assert!(
                cell.angular_distortion >= 0.0,
                "angular_distortion should be >= 0 at lon={} lat={}",
                cell.lon,
                cell.lat
            );
        }
    }

    #[test]
    fn get_out_of_bounds_returns_none() {
        let pj = create("+proj=merc +datum=WGS84").unwrap();
        let raster = DistortionRaster::compute(&pj, -5.0, -5.0, 5.0, 5.0, 5.0);
        assert!(raster.get(10, 10).is_none());
        assert!(raster.get(0, 0).is_some());
    }
}