oxiproj-engine 0.1.0

//! The `Pj` operation wrapper: prepare/finalize around an inner operation.
//!
//! Ported from PROJ 9.8.0 `src/fwd.cpp` (`fwd_prepare`/`fwd_finalize`) and
//! `src/inv.cpp` (`inv_prepare`/`inv_finalize`).

use oxiproj_core::{Coord, Ellipsoid, IoUnits, Lp, Operation, ProjError, ProjResult, Xy};

/// Latitude tolerance for the prepare-stage domain check (radians).
const PJ_EPS_LAT: f64 = 1e-12;

/// A prepared coordinate operation with input/output unit handling.
///
/// Wraps an inner [`Operation`] and applies the unit conversions, prime-meridian
/// and central-meridian shifts, and axis bookkeeping that PROJ performs around
/// the raw projection math.
#[derive(Debug)]
pub struct Pj {
    /// The inner operation (projection or conversion) driven via 4D.
    pub operation: Box<dyn Operation>,
    /// The reference ellipsoid.
    pub ellipsoid: Ellipsoid,
    /// Central meridian (radians).
    pub lam0: f64,
    /// Latitude of origin (radians).
    pub phi0: f64,
    /// False easting (meters, pre-scale).
    pub x0: f64,
    /// False northing (meters, pre-scale).
    pub y0: f64,
    /// False up (meters, pre-scale).
    pub z0: f64,
    /// Scale factor at the central meridian.
    pub k0: f64,
    /// Horizontal input-unit-to-meter factor.
    pub to_meter: f64,
    /// Horizontal meter-to-output-unit factor.
    pub fr_meter: f64,
    /// Vertical input-unit-to-meter factor.
    pub vto_meter: f64,
    /// Vertical meter-to-output-unit factor.
    pub vfr_meter: f64,
    /// Prime-meridian offset from Greenwich (radians).
    pub from_greenwich: f64,
    /// Allow longitudes outside the normal range (`+over`).
    pub over: bool,
    /// Geocentric latitude flag (`+geoc`).
    pub geoc: bool,
    /// Whether this is a pass-through geographic (lat/long) operation.
    pub is_latlong: bool,
    /// Input units (forward direction).
    pub left: IoUnits,
    /// Output units (forward direction).
    pub right: IoUnits,
    /// Whether the operation runs inverted in its containing pipeline.
    pub inverted: bool,
    /// Skip prepare/finalize entirely (used by the top-level pipeline driver).
    pub bypass_prepare_finalize: bool,
}

impl Pj {
    /// Forward transform with prepare/finalize.
    ///
    /// Ported from `src/fwd.cpp` (`fwd_prepare` -> operation -> `fwd_finalize`).
    fn forward_impl(&self, c: Coord) -> ProjResult<Coord> {
        // --- fwd_prepare ---
        let v = c.v();
        if !v[0].is_finite() || !v[1].is_finite() || !v[2].is_finite() {
            return Ok(Coord::error());
        }

        let prepared = if self.left == IoUnits::Radians {
            let mut lam = v[0];
            let mut phi = v[1];
            let z = v[2];
            let t = v[3];

            let tt = phi.abs() - oxiproj_core::M_HALFPI;
            if tt > PJ_EPS_LAT {
                return Err(ProjError::InvalidCoord);
            }
            if !(-10.0..=10.0).contains(&lam) {
                return Err(ProjError::InvalidCoord);
            }
            // Clamp latitude to ±π/2 (PROJ src/fwd.cpp).
            phi = phi.clamp(-oxiproj_core::M_HALFPI, oxiproj_core::M_HALFPI);

            if self.geoc && self.ellipsoid.es != 0.0 {
                phi = (self.ellipsoid.one_es * phi.tan()).atan();
            }
            if !self.over {
                lam = oxiproj_core::adjlon(lam);
            }
            // TODO(v0.2): datum shift in prepare (hgridshift / helmert / cart).
            // TODO(v0.2): vgridshift.
            lam = (lam - self.from_greenwich) - self.lam0;
            if !self.over {
                lam = oxiproj_core::adjlon(lam);
            }
            Coord::new(lam, phi, z, t)
        } else {
            // TODO(v0.2): cartesian helmert in prepare.
            c
        };

        // --- operation ---
        let mid = self.operation.forward_4d(prepared)?;

        // --- fwd_finalize ---
        let mv = mid.v();
        let mut x = mv[0];
        let mut y = mv[1];
        let mut z = mv[2];
        let t = mv[3];
        match self.right {
            IoUnits::Classic | IoUnits::Projected => {
                if self.right == IoUnits::Classic {
                    x *= self.ellipsoid.a;
                    y *= self.ellipsoid.a;
                }
                x = self.fr_meter * (x + self.x0);
                y = self.fr_meter * (y + self.y0);
                z = self.vfr_meter * (z + self.z0);
            }
            IoUnits::Radians => {
                z = self.vfr_meter * (z + self.z0);
                if !self.over {
                    x = oxiproj_core::adjlon(x);
                }
            }
            IoUnits::Cartesian => {
                x *= self.fr_meter;
                y *= self.fr_meter;
                z *= self.fr_meter;
                // TODO(v0.2): is_geocent cartesian handling.
            }
            IoUnits::Whatever | IoUnits::Degrees => {}
        }
        // TODO(v0.2): axisswap.
        Ok(Coord::new(x, y, z, t))
    }

    /// Inverse transform with prepare/finalize.
    ///
    /// Ported from `src/inv.cpp` (`inv_prepare` -> operation -> `inv_finalize`).
    fn inverse_impl(&self, c: Coord) -> ProjResult<Coord> {
        // --- inv_prepare ---
        let v = c.v();
        if !v[0].is_finite() || !v[1].is_finite() || !v[2].is_finite() {
            return Err(ProjError::OutsideProjectionDomain);
        }
        // TODO(v0.2): axisswap.
        let t = v[3];
        let prepared = match self.right {
            IoUnits::Whatever | IoUnits::Degrees => c,
            IoUnits::Cartesian => {
                let x = self.to_meter * v[0];
                let y = self.to_meter * v[1];
                let z = self.to_meter * v[2];
                // TODO(v0.2): is_geocent cartesian handling.
                Coord::new(x, y, z, t)
            }
            IoUnits::Projected | IoUnits::Classic => {
                let mut x = self.to_meter * v[0] - self.x0;
                let mut y = self.to_meter * v[1] - self.y0;
                let z = self.vto_meter * v[2] - self.z0;
                if self.right == IoUnits::Classic {
                    x *= self.ellipsoid.ra;
                    y *= self.ellipsoid.ra;
                }
                Coord::new(x, y, z, t)
            }
            IoUnits::Radians => {
                let z = self.vto_meter * v[2] - self.z0;
                Coord::new(v[0], v[1], z, t)
            }
        };

        // --- operation ---
        let mid = self.operation.inverse_4d(prepared)?;

        // --- inv_finalize ---
        if self.left == IoUnits::Radians {
            let mv = mid.v();
            let mut lam = mv[0];
            let phi = mv[1];
            let z = mv[2];
            let t = mv[3];
            lam = lam + self.from_greenwich + self.lam0;
            if !self.over {
                lam = oxiproj_core::adjlon(lam);
            }
            let phi = if self.geoc && self.ellipsoid.es != 0.0 {
                let limit = oxiproj_core::M_HALFPI - 1e-9;
                if phi.abs() < limit {
                    (self.ellipsoid.rone_es * phi.tan()).atan()
                } else {
                    phi
                }
            } else {
                phi
            };
            // TODO(v0.2): vgridshift / hgridshift / helmert.
            Ok(Coord::new(lam, phi, z, t))
        } else {
            Ok(mid)
        }
    }

    /// Run the operation in the forward direction, honoring `inverted`.
    pub fn forward(&self, c: Coord) -> ProjResult<Coord> {
        if self.bypass_prepare_finalize {
            if self.inverted {
                return self.operation.inverse_4d(c);
            } else {
                return self.operation.forward_4d(c);
            }
        }
        if self.inverted {
            self.inverse_impl(c)
        } else {
            self.forward_impl(c)
        }
    }

    /// Run the operation in the inverse direction, honoring `inverted`.
    pub fn inverse(&self, c: Coord) -> ProjResult<Coord> {
        if self.bypass_prepare_finalize {
            if self.inverted {
                return self.operation.forward_4d(c);
            } else {
                return self.operation.inverse_4d(c);
            }
        }
        if self.inverted {
            self.forward_impl(c)
        } else {
            self.inverse_impl(c)
        }
    }

    /// Compute map-scale factors (Tissot indicatrix) at a geographic coordinate.
    ///
    /// `coord` must be in geographic input units — radians if `self.left == IoUnits::Radians`,
    /// degrees if `self.left == IoUnits::Degrees`. The coordinate layout is
    /// `Coord::new(longitude, latitude, z, t)`.
    ///
    /// Returns the [`oxiproj_core::Factors`] struct describing Tissot indicatrix parameters
    /// at the given point.
    ///
    /// Ported from PROJ 9.8.0 `src/factors.cpp` (`proj_factors` → `pj_factors` → `pj_deriv`).
    pub fn factors(&self, coord: Coord) -> ProjResult<oxiproj_core::Factors> {
        let v = coord.v();
        // Convert lon/lat to radians
        let (phi, lam_abs) = if self.left == IoUnits::Radians {
            (v[1], v[0])
        } else {
            // Assume degrees
            (v[1].to_radians(), v[0].to_radians())
        };

        // Longitude relative to central meridian (as done in fwd_prepare)
        let lam = lam_abs - self.from_greenwich - self.lam0;

        let es = self.ellipsoid.es;
        let one_es = self.ellipsoid.one_es;

        // The closure calls the raw inner operation's forward in normalized space.
        let op = &*self.operation;

        oxiproj_core::Factors::compute(phi, lam, es, one_es, |lam_f, phi_f| {
            let c = Coord::new(lam_f, phi_f, 0.0, 0.0);
            let r = op.forward_4d(c)?;
            let rv = r.v();
            Ok((rv[0], rv[1]))
        })
    }
}

impl Operation for Pj {
    fn forward_4d(&self, c: Coord) -> ProjResult<Coord> {
        self.forward(c)
    }
    fn inverse_4d(&self, c: Coord) -> ProjResult<Coord> {
        self.inverse(c)
    }
    fn forward_2d(&self, lp: Lp) -> ProjResult<Xy> {
        let c = self.forward(Coord::new(lp.lam, lp.phi, 0.0, 0.0))?;
        Ok(c.xy())
    }
    fn inverse_2d(&self, xy: Xy) -> ProjResult<Lp> {
        let c = self.inverse(Coord::new(xy.x, xy.y, 0.0, 0.0))?;
        Ok(c.lp())
    }
}

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

    #[derive(Debug)]
    struct Id;
    impl Operation for Id {
        fn forward_2d(&self, lp: Lp) -> ProjResult<Xy> {
            Ok(Xy::new(lp.lam, lp.phi))
        }
        fn inverse_2d(&self, xy: Xy) -> ProjResult<Lp> {
            Ok(Lp::new(xy.x, xy.y))
        }
    }

    fn unit_latlong() -> Pj {
        Pj {
            operation: Box::new(Id),
            ellipsoid: Ellipsoid::named("WGS84").unwrap(),
            lam0: 0.0,
            phi0: 0.0,
            x0: 0.0,
            y0: 0.0,
            z0: 0.0,
            k0: 1.0,
            to_meter: 1.0,
            fr_meter: 1.0,
            vto_meter: 1.0,
            vfr_meter: 1.0,
            from_greenwich: 0.0,
            over: false,
            geoc: false,
            is_latlong: true,
            left: IoUnits::Radians,
            right: IoUnits::Radians,
            inverted: false,
            bypass_prepare_finalize: false,
        }
    }

    #[test]
    fn latlong_round_trip() {
        let pj = unit_latlong();
        let lam = 0.2;
        let phi = 0.5;
        let fwd = pj.forward(Coord::new(lam, phi, 0.0, 0.0)).unwrap();
        let inv = pj.inverse(fwd).unwrap();
        let r = inv.v();
        assert!((r[0] - lam).abs() < 1e-12);
        assert!((r[1] - phi).abs() < 1e-12);
    }

    #[test]
    fn non_finite_forward_yields_error_coord() {
        let pj = unit_latlong();
        let out = pj.forward(Coord::new(f64::NAN, 0.0, 0.0, 0.0)).unwrap();
        assert!(out.is_error());
    }

    #[test]
    fn out_of_domain_latitude_rejected() {
        let pj = unit_latlong();
        let res = pj.forward(Coord::new(0.0, 2.0, 0.0, 0.0));
        assert!(res.is_err());
    }
}