oxiproj-transformations 0.1.2

//! Horner polynomial transformation — ported from PROJ 9.8.0 `src/transformations/horner.cpp`.

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

fn real_coeff_count(deg: usize) -> usize {
    (deg + 1) * (deg + 2) / 2
}

fn complex_coeff_count(deg: usize) -> usize {
    2 * deg + 2
}

fn parse_coefs(s: &str, expected: usize) -> ProjResult<Vec<f64>> {
    let vals: Result<Vec<f64>, _> = s.split(',').map(|t| t.trim().parse::<f64>()).collect();
    let vals = vals.map_err(|_| ProjError::IllegalArgValue)?;
    if vals.len() != expected {
        return Err(ProjError::IllegalArgValue);
    }
    Ok(vals)
}

fn parse_origin(s: &str) -> ProjResult<[f64; 2]> {
    let parts: Vec<&str> = s.split(',').collect();
    if parts.len() != 2 {
        return Err(ProjError::IllegalArgValue);
    }
    let u = parts[0]
        .trim()
        .parse::<f64>()
        .map_err(|_| ProjError::IllegalArgValue)?;
    let v = parts[1]
        .trim()
        .parse::<f64>()
        .map_err(|_| ProjError::IllegalArgValue)?;
    Ok([u, v])
}

/// Evaluate the 2D real Horner polynomial for both E and N simultaneously.
/// Returns (E, N). Ported exactly from PROJ `double_real_horner_eval`.
fn double_real_horner_eval(
    order: usize,
    cx: &[f64], // fwd_u or inv_u
    cy: &[f64], // fwd_v or inv_v
    e: f64,
    n: f64,
    order_offset: usize,
) -> (f64, f64) {
    let sz = (order + 1) * (order + 2) / 2;
    let mut idx_x = sz;
    let mut idx_y = sz;

    idx_y -= 1;
    let mut big_n = cy[idx_y];
    idx_x -= 1;
    let mut big_e = cx[idx_x];

    let mut r = order;
    while r > order_offset {
        idx_y -= 1;
        let mut u = cy[idx_y];
        idx_x -= 1;
        let mut v = cx[idx_x];

        let mut c = order;
        loop {
            if c < r {
                break;
            }
            idx_y -= 1;
            u = n * u + cy[idx_y];
            idx_x -= 1;
            v = e * v + cx[idx_x];
            if c == 0 {
                break;
            }
            c -= 1;
        }

        big_n = e * big_n + u;
        big_e = n * big_e + v;

        r -= 1;
    }

    (big_e, big_n)
}

/// Evaluate a single 1D real Horner polynomial.
/// Ported exactly from PROJ `single_real_horner_eval`.
fn single_real_horner_eval(order: usize, cx: &[f64], x: f64, order_offset: usize) -> f64 {
    let sz = order + 1;
    let mut idx = sz;

    idx -= 1;
    let mut u = cx[idx];

    let mut r = order;
    while r > order_offset {
        idx -= 1;
        u = x * u + cx[idx];
        r -= 1;
    }

    u
}

/// Evaluate the complex Horner polynomial.
/// Returns (E, N). Ported exactly from PROJ `complex_horner_eval`.
fn complex_horner_eval(order: usize, c: &[f64], e: f64, n: f64, order_offset: usize) -> (f64, f64) {
    let sz = 2 * order + 2;
    let cbeg_idx = order_offset * 2;
    let mut idx = sz;

    idx -= 1;
    let mut big_e = c[idx];
    idx -= 1;
    let mut big_n = c[idx];

    while idx > cbeg_idx {
        idx -= 1;
        let w = n * big_e + e * big_n + c[idx];
        idx -= 1;
        big_n = n * big_n - e * big_e + c[idx];
        big_e = w;
    }

    (big_e, big_n)
}

#[derive(Debug)]
struct Horner {
    degree: usize,
    range: f64,
    inverse_tolerance: f64,
    fwd_origin: [f64; 2],
    inv_origin: [f64; 2],
    fwd_u: Vec<f64>,
    fwd_v: Vec<f64>,
    inv_u: Vec<f64>,
    inv_v: Vec<f64>,
    fwd_c: Vec<f64>,
    inv_c: Vec<f64>,
    uneg: bool,
    vneg: bool,
    has_inv: bool,
    complex_mode: bool,
}

impl Operation for Horner {
    fn forward_4d(&self, coord: Coord) -> ProjResult<Coord> {
        let v = coord.v();
        let mut e = v[0] - self.fwd_origin[0];
        let mut n = v[1] - self.fwd_origin[1];
        let z = v[2];
        let t = v[3];

        if self.complex_mode {
            if self.uneg {
                e = -e;
            }
            if self.vneg {
                n = -n;
            }
        }

        if e.abs() > self.range || n.abs() > self.range {
            return Err(ProjError::OutsideProjectionDomain);
        }

        let (out_e, out_n) = if self.complex_mode {
            complex_horner_eval(self.degree, &self.fwd_c, e, n, 0)
        } else {
            double_real_horner_eval(self.degree, &self.fwd_u, &self.fwd_v, e, n, 0)
        };

        Ok(Coord::new(out_e, out_n, z, t))
    }

    fn inverse_4d(&self, coord: Coord) -> ProjResult<Coord> {
        let v = coord.v();
        let z = v[2];
        let t = v[3];

        if self.has_inv {
            let mut e = v[0] - self.inv_origin[0];
            let mut n = v[1] - self.inv_origin[1];

            if self.complex_mode {
                if self.uneg {
                    e = -e;
                }
                if self.vneg {
                    n = -n;
                }
            }

            if e.abs() > self.range || n.abs() > self.range {
                return Err(ProjError::OutsideProjectionDomain);
            }

            let (out_e, out_n) = if self.complex_mode {
                complex_horner_eval(self.degree, &self.inv_c, e, n, 0)
            } else {
                double_real_horner_eval(self.degree, &self.inv_u, &self.inv_v, e, n, 0)
            };

            Ok(Coord::new(out_e, out_n, z, t))
        } else if self.complex_mode {
            Err(ProjError::NoInverseOp)
        } else {
            // Newton's iterative inverse for real mode
            let de = v[0] - self.fwd_u[0];
            let dn = v[1] - self.fwd_v[0];
            let mut x0 = 0.0_f64;
            let mut y0 = 0.0_f64;
            let mut converged = false;

            for _ in 0..32 {
                let (mb, mc) =
                    double_real_horner_eval(self.degree, &self.fwd_u, &self.fwd_v, x0, y0, 1);
                let ma = single_real_horner_eval(self.degree, &self.fwd_u, x0, 1);
                let md = single_real_horner_eval(self.degree, &self.fwd_v, y0, 1);

                let det = ma * md - mb * mc;
                if det.abs() < f64::EPSILON {
                    return Err(ProjError::NoConvergence);
                }
                let idet = 1.0 / det;
                let x = idet * (md * de - mb * dn);
                let y = idet * (ma * dn - mc * de);

                if (x - x0).abs() < self.inverse_tolerance
                    && (y - y0).abs() < self.inverse_tolerance
                {
                    converged = true;
                    x0 = x;
                    y0 = y;
                    break;
                }

                x0 = x;
                y0 = y;
            }

            if !converged {
                return Err(ProjError::NoConvergence);
            }

            Ok(Coord::new(
                x0 + self.fwd_origin[0],
                y0 + self.fwd_origin[1],
                z,
                t,
            ))
        }
    }

    fn has_inverse(&self) -> bool {
        self.has_inv || !self.complex_mode
    }
}

/// Construct the `horner` transformation.
///
/// Ported from PROJ `horner.cpp`. Supports real 2D polynomial and complex
/// polynomial modes, with optional explicit inverse coefficients.
pub fn new(p: &crate::TransParams) -> ProjResult<crate::TransBuild> {
    let params = p.params;

    let degree = match params.get_int("deg") {
        Some(d) if d > 0 => d as usize,
        Some(_) => return Err(ProjError::IllegalArgValue),
        None => return Err(ProjError::MissingArg),
    };

    let range = params.get_f64("range").unwrap_or(500_000.0);
    let inverse_tolerance = params.get_f64("inv_tolerance").unwrap_or(0.001);

    let fwd_origin = match params.get_str("fwd_origin") {
        Some(s) => parse_origin(s)?,
        None => [0.0, 0.0],
    };

    let uneg = params.get_bool("uneg");
    let vneg = params.get_bool("vneg");

    let complex_mode = params.exists("fwd_c");

    let (fwd_u, fwd_v, fwd_c, inv_u, inv_v, inv_c, has_inv) = if complex_mode {
        let n_complex = complex_coeff_count(degree);

        let fwd_c = match params.get_str("fwd_c") {
            Some(s) => parse_coefs(s, n_complex)?,
            None => return Err(ProjError::MissingArg),
        };

        let (inv_c, has_inv) = match params.get_str("inv_c") {
            Some(s) => (parse_coefs(s, n_complex)?, true),
            None => (Vec::new(), false),
        };

        (
            Vec::new(),
            Vec::new(),
            fwd_c,
            Vec::new(),
            Vec::new(),
            inv_c,
            has_inv,
        )
    } else {
        let n_real = real_coeff_count(degree);

        let fwd_u = match params.get_str("fwd_u") {
            Some(s) => parse_coefs(s, n_real)?,
            None => return Err(ProjError::MissingArg),
        };
        let fwd_v = match params.get_str("fwd_v") {
            Some(s) => parse_coefs(s, n_real)?,
            None => return Err(ProjError::MissingArg),
        };

        let has_inv = params.exists("inv_u") && params.exists("inv_v");
        let inv_u = if has_inv {
            match params.get_str("inv_u") {
                Some(s) => parse_coefs(s, n_real)?,
                None => Vec::new(),
            }
        } else {
            Vec::new()
        };
        let inv_v = if has_inv {
            match params.get_str("inv_v") {
                Some(s) => parse_coefs(s, n_real)?,
                None => Vec::new(),
            }
        } else {
            Vec::new()
        };

        (fwd_u, fwd_v, Vec::new(), inv_u, inv_v, Vec::new(), has_inv)
    };

    let inv_origin = match params.get_str("inv_origin") {
        Some(s) => parse_origin(s)?,
        None => {
            if has_inv {
                [0.0, 0.0]
            } else {
                fwd_origin
            }
        }
    };

    let op = Horner {
        degree,
        range,
        inverse_tolerance,
        fwd_origin,
        inv_origin,
        fwd_u,
        fwd_v,
        inv_u,
        inv_v,
        fwd_c,
        inv_c,
        uneg,
        vneg,
        has_inv,
        complex_mode,
    };

    Ok(crate::TransBuild::new(
        Box::new(op),
        IoUnits::Whatever,
        IoUnits::Whatever,
    ))
}

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

    struct MockParams {
        strings: HashMap<&'static str, String>,
        ints: HashMap<&'static str, i64>,
        floats: HashMap<&'static str, f64>,
        bools: HashMap<&'static str, bool>,
    }

    impl MockParams {
        fn new() -> Self {
            Self {
                strings: HashMap::new(),
                ints: HashMap::new(),
                floats: HashMap::new(),
                bools: HashMap::new(),
            }
        }
        fn with_str(mut self, k: &'static str, v: &str) -> Self {
            self.strings.insert(k, v.to_string());
            self
        }
        fn with_int(mut self, k: &'static str, v: i64) -> Self {
            self.ints.insert(k, v);
            self
        }
        fn with_float(mut self, k: &'static str, v: f64) -> Self {
            self.floats.insert(k, v);
            self
        }
        #[allow(dead_code)]
        fn with_flag(mut self, k: &'static str) -> Self {
            self.bools.insert(k, true);
            self
        }
    }

    impl crate::TransParamLookup for MockParams {
        fn get_dms(&self, key: &str) -> Option<f64> {
            self.floats.get(key).copied()
        }
        fn get_f64(&self, key: &str) -> Option<f64> {
            self.floats.get(key).copied()
        }
        fn get_int(&self, key: &str) -> Option<i64> {
            self.ints.get(key).copied()
        }
        fn get_str(&self, key: &str) -> Option<&str> {
            self.strings.get(key).map(|s| s.as_str())
        }
        fn get_bool(&self, key: &str) -> bool {
            *self.bools.get(key).unwrap_or(&false)
        }
        fn exists(&self, key: &str) -> bool {
            self.strings.contains_key(key)
                || self.ints.contains_key(key)
                || self.floats.contains_key(key)
                || self.bools.contains_key(key)
        }
    }

    fn make_params<'a>(
        mp: &'a MockParams,
        ell: &'a oxiproj_core::Ellipsoid,
    ) -> crate::TransParams<'a> {
        crate::TransParams {
            ellipsoid: ell,
            params: mp,
            registry: None,
        }
    }

    #[test]
    fn missing_deg_is_error() {
        let mp = MockParams::new()
            .with_str("fwd_u", "0,1,0")
            .with_str("fwd_v", "0,1,0");
        let ell = oxiproj_core::Ellipsoid::named("WGS84").unwrap();
        let pp = make_params(&mp, &ell);
        assert!(new(&pp).is_err());
    }

    #[test]
    fn degree1_identity_real() {
        // degree=1: 3 coefficients per polynomial
        // For the 2D Horner scheme, identity E=e, N=n uses:
        // fwd_u = [0, 1, 0] => E = e*1 + n*0 + 0 = e
        // fwd_v = [0, 1, 0] => N = n*1 + e*0 + 0 = n
        let mp = MockParams::new()
            .with_int("deg", 1)
            .with_str("fwd_u", "0,1,0")
            .with_str("fwd_v", "0,1,0");
        let ell = oxiproj_core::Ellipsoid::named("WGS84").unwrap();
        let pp = make_params(&mp, &ell);
        let b = new(&pp).expect("build failed");
        let out = b
            .operation
            .forward_4d(Coord::new(100.0, 200.0, 0.0, 0.0))
            .expect("fwd failed");
        let v = out.v();
        assert!((v[0] - 100.0).abs() < 1e-9, "x = {}", v[0]);
        assert!((v[1] - 200.0).abs() < 1e-9, "y = {}", v[1]);
    }

    #[test]
    fn coefficient_count_validation() {
        // degree=2 needs 6 coefficients, but we provide only 3
        let mp = MockParams::new()
            .with_int("deg", 2)
            .with_str("fwd_u", "0,1,0")
            .with_str("fwd_v", "0,1,0");
        let ell = oxiproj_core::Ellipsoid::named("WGS84").unwrap();
        let pp = make_params(&mp, &ell);
        assert!(new(&pp).is_err());
    }

    #[test]
    fn explicit_inverse_round_trip() {
        // degree=1 with explicit inv_u and inv_v (identity)
        // forward then inverse should return original coord
        let mp = MockParams::new()
            .with_int("deg", 1)
            .with_str("fwd_u", "0,1,0")
            .with_str("fwd_v", "0,1,0")
            .with_str("inv_u", "0,1,0")
            .with_str("inv_v", "0,1,0");
        let ell = oxiproj_core::Ellipsoid::named("WGS84").unwrap();
        let pp = make_params(&mp, &ell);
        let b = new(&pp).expect("build failed");
        let fwd = b
            .operation
            .forward_4d(Coord::new(50.0, 100.0, 5.0, 0.0))
            .expect("fwd failed");
        let back = b.operation.inverse_4d(fwd).expect("inv failed");
        let v = back.v();
        assert!((v[0] - 50.0).abs() < 1e-9, "x = {}", v[0]);
        assert!((v[1] - 100.0).abs() < 1e-9, "y = {}", v[1]);
    }

    #[test]
    fn iterative_inverse_round_trip() {
        // Real mode without explicit inverse: iterative Newton's method
        let mp = MockParams::new()
            .with_int("deg", 1)
            .with_str("fwd_u", "0,1,0")
            .with_str("fwd_v", "0,1,0");
        let ell = oxiproj_core::Ellipsoid::named("WGS84").unwrap();
        let pp = make_params(&mp, &ell);
        let b = new(&pp).expect("build failed");
        let fwd = b
            .operation
            .forward_4d(Coord::new(75.0, 150.0, 3.0, 0.0))
            .expect("fwd failed");
        let back = b.operation.inverse_4d(fwd).expect("inv failed");
        let v = back.v();
        assert!((v[0] - 75.0).abs() < 1e-6, "x = {}", v[0]);
        assert!((v[1] - 150.0).abs() < 1e-6, "y = {}", v[1]);
    }

    #[test]
    fn range_check_error() {
        let mp = MockParams::new()
            .with_int("deg", 1)
            .with_str("fwd_u", "0,1,0")
            .with_str("fwd_v", "0,1,0")
            .with_float("range", 100.0);
        let ell = oxiproj_core::Ellipsoid::named("WGS84").unwrap();
        let pp = make_params(&mp, &ell);
        let b = new(&pp).expect("build failed");
        // 200.0 > range=100.0 => error
        let err = b
            .operation
            .forward_4d(Coord::new(200.0, 50.0, 0.0, 0.0))
            .err();
        assert_eq!(err, Some(ProjError::OutsideProjectionDomain));
    }
}