oxiproj-transformations 0.1.2

//! SIMD-batched Helmert and affine transform application.
//!
//! Provides SIMD-accelerated (feature `simd`, using the `wide` crate) and scalar
//! fallback implementations of:
//! - 3-parameter (translation-only) Helmert batch transform ([`apply_helmert_3_batch`])
//! - 7-parameter Helmert batch transform, **position-vector** convention
//!   ([`apply_helmert_7_batch_position_vector`])
//! - 2-D affine batch transform ([`apply_affine_2d_batch`])
//!
//! All functions use **Structure-of-Arrays** (SoA) layout — separate `&[f64]` slices
//! for each coordinate axis — which maps naturally onto 4-wide SIMD registers.
//! The scalar fallback (compiled when the `simd` feature is absent) is
//! otherwise identical in behaviour, enabling use in `no_std` or constrained targets.
//!
//! ## Input constraints
//! All input coordinate slices must have the same length `n`. The output slices
//! must each have length `≥ n`.
//!
//! ## Rotation convention — read this before choosing a 7-parameter function
//!
//! This module's 7-parameter batch transform,
//! [`apply_helmert_7_batch_position_vector`], uses the **position-vector**
//! convention with a **dimensionless** scale (pass `scale_ppm * 1e-6`):
//! ```text
//! R = (1 + scale) × [[  1,  −rz,  ry ],
//!                    [ rz,    1,  −rx ],
//!                    [−ry,   rx,    1 ]]
//! ```
//! The sibling module [`crate::helmert_batch`] provides a *different* batch
//! Helmert entry point, [`crate::helmert_batch::helmert7_batch_scalar`], which
//! uses the **coordinate-frame** convention (the transpose of the matrix
//! above — every off-diagonal rotation term has the *opposite* sign) with a
//! **parts-per-million** scale (`scale_ppm`, not pre-divided by 1e-6).
//! Passing the same `(rx, ry, rz)` values to the "other" function silently
//! flips the sign of every rotation term — a classic geodesy footgun.
//! **Always call the convention-suffixed name** for the one you mean; do not
//! use the deprecated bare `apply_helmert_7_batch` alias in new code.

// ── SIMD implementation ─────────────────────────────────────────────────────
#[cfg(feature = "simd")]
mod inner {
    #![allow(clippy::too_many_arguments)]
    use wide::f64x4;

    /// Apply a 3-parameter (translation-only) Helmert transform to N coordinate triples.
    ///
    /// `xs`, `ys`, `zs` — input slices (all the same length `n`).
    /// `xs_out`, `ys_out`, `zs_out` — output slices (must be at least length `n`).
    ///
    /// Processes 4 points per iteration using `f64x4`; a scalar tail handles
    /// any remaining points when `n` is not a multiple of 4.
    pub fn apply_helmert_3_batch(
        tx: f64,
        ty: f64,
        tz: f64,
        xs: &[f64],
        ys: &[f64],
        zs: &[f64],
        xs_out: &mut [f64],
        ys_out: &mut [f64],
        zs_out: &mut [f64],
    ) {
        let n = xs.len();
        let tx4 = f64x4::splat(tx);
        let ty4 = f64x4::splat(ty);
        let tz4 = f64x4::splat(tz);
        let chunks = n / 4;
        for i in 0..chunks {
            let b = i * 4;
            let x4 = f64x4::new([xs[b], xs[b + 1], xs[b + 2], xs[b + 3]]);
            let y4 = f64x4::new([ys[b], ys[b + 1], ys[b + 2], ys[b + 3]]);
            let z4 = f64x4::new([zs[b], zs[b + 1], zs[b + 2], zs[b + 3]]);
            let rx: [f64; 4] = (x4 + tx4).into();
            let ry: [f64; 4] = (y4 + ty4).into();
            let rz: [f64; 4] = (z4 + tz4).into();
            xs_out[b..b + 4].copy_from_slice(&rx);
            ys_out[b..b + 4].copy_from_slice(&ry);
            zs_out[b..b + 4].copy_from_slice(&rz);
        }
        // Scalar tail
        for i in (chunks * 4)..n {
            xs_out[i] = xs[i] + tx;
            ys_out[i] = ys[i] + ty;
            zs_out[i] = zs[i] + tz;
        }
    }

    /// Apply a 7-parameter Helmert transform (position-vector convention) to N triples.
    ///
    /// `rot_x`, `rot_y`, `rot_z` — rotations in **radians** (linearised).
    /// `scale` — dimensionless scale factor (e.g. pass `scale_ppm * 1e-6`).
    ///
    /// The linearised rotation matrix (position-vector convention):
    /// ```text
    /// R = (1 + scale) × [[  1,  −rz,  ry ],
    ///                    [ rz,    1,  −rx ],
    ///                    [−ry,   rx,    1 ]]
    /// ```
    /// Transformed point: `X′ = R·X + T`
    pub fn apply_helmert_7_batch(
        tx: f64,
        ty: f64,
        tz: f64,
        rot_x: f64,
        rot_y: f64,
        rot_z: f64,
        scale: f64,
        xs: &[f64],
        ys: &[f64],
        zs: &[f64],
        xs_out: &mut [f64],
        ys_out: &mut [f64],
        zs_out: &mut [f64],
    ) {
        let s1 = 1.0 + scale;
        // Pre-compute the 9 matrix elements as SIMD splatted scalars.
        let m00 = f64x4::splat(s1);
        let m01 = f64x4::splat(-rot_z * s1);
        let m02 = f64x4::splat(rot_y * s1);
        let m10 = f64x4::splat(rot_z * s1);
        let m11 = f64x4::splat(s1);
        let m12 = f64x4::splat(-rot_x * s1);
        let m20 = f64x4::splat(-rot_y * s1);
        let m21 = f64x4::splat(rot_x * s1);
        let m22 = f64x4::splat(s1);
        let tx4 = f64x4::splat(tx);
        let ty4 = f64x4::splat(ty);
        let tz4 = f64x4::splat(tz);
        let n = xs.len();
        let chunks = n / 4;
        for i in 0..chunks {
            let b = i * 4;
            let x4 = f64x4::new([xs[b], xs[b + 1], xs[b + 2], xs[b + 3]]);
            let y4 = f64x4::new([ys[b], ys[b + 1], ys[b + 2], ys[b + 3]]);
            let z4 = f64x4::new([zs[b], zs[b + 1], zs[b + 2], zs[b + 3]]);
            let out_x: [f64; 4] = (m00 * x4 + m01 * y4 + m02 * z4 + tx4).into();
            let out_y: [f64; 4] = (m10 * x4 + m11 * y4 + m12 * z4 + ty4).into();
            let out_z: [f64; 4] = (m20 * x4 + m21 * y4 + m22 * z4 + tz4).into();
            xs_out[b..b + 4].copy_from_slice(&out_x);
            ys_out[b..b + 4].copy_from_slice(&out_y);
            zs_out[b..b + 4].copy_from_slice(&out_z);
        }
        // Scalar tail
        for i in (chunks * 4)..n {
            let x = xs[i];
            let y = ys[i];
            let z = zs[i];
            xs_out[i] = s1 * x - rot_z * s1 * y + rot_y * s1 * z + tx;
            ys_out[i] = rot_z * s1 * x + s1 * y - rot_x * s1 * z + ty;
            zs_out[i] = -rot_y * s1 * x + rot_x * s1 * y + s1 * z + tz;
        }
    }

    /// Apply a 2-D affine transform to N coordinate pairs.
    ///
    /// Transform: `x′ = a00·x + a01·y + b0`,  `y′ = a10·x + a11·y + b1`
    ///
    /// `xs`, `ys` — input slices (same length `n`).
    /// `xs_out`, `ys_out` — output slices (must be at least length `n`).
    pub fn apply_affine_2d_batch(
        a00: f64,
        a01: f64,
        b0: f64,
        a10: f64,
        a11: f64,
        b1: f64,
        xs: &[f64],
        ys: &[f64],
        xs_out: &mut [f64],
        ys_out: &mut [f64],
    ) {
        let a00v = f64x4::splat(a00);
        let a01v = f64x4::splat(a01);
        let b0v = f64x4::splat(b0);
        let a10v = f64x4::splat(a10);
        let a11v = f64x4::splat(a11);
        let b1v = f64x4::splat(b1);
        let n = xs.len();
        let chunks = n / 4;
        for i in 0..chunks {
            let b = i * 4;
            let x4 = f64x4::new([xs[b], xs[b + 1], xs[b + 2], xs[b + 3]]);
            let y4 = f64x4::new([ys[b], ys[b + 1], ys[b + 2], ys[b + 3]]);
            let out_x: [f64; 4] = (a00v * x4 + a01v * y4 + b0v).into();
            let out_y: [f64; 4] = (a10v * x4 + a11v * y4 + b1v).into();
            xs_out[b..b + 4].copy_from_slice(&out_x);
            ys_out[b..b + 4].copy_from_slice(&out_y);
        }
        // Scalar tail
        for i in (chunks * 4)..n {
            xs_out[i] = a00 * xs[i] + a01 * ys[i] + b0;
            ys_out[i] = a10 * xs[i] + a11 * ys[i] + b1;
        }
    }
}

// ── Scalar fallback ─────────────────────────────────────────────────────────
#[cfg(not(feature = "simd"))]
mod inner {
    #![allow(clippy::too_many_arguments)]
    /// Scalar fallback for 3-parameter Helmert batch transform.
    pub fn apply_helmert_3_batch(
        tx: f64,
        ty: f64,
        tz: f64,
        xs: &[f64],
        ys: &[f64],
        zs: &[f64],
        xs_out: &mut [f64],
        ys_out: &mut [f64],
        zs_out: &mut [f64],
    ) {
        for i in 0..xs.len() {
            xs_out[i] = xs[i] + tx;
            ys_out[i] = ys[i] + ty;
            zs_out[i] = zs[i] + tz;
        }
    }

    /// Scalar fallback for 7-parameter Helmert batch transform.
    pub fn apply_helmert_7_batch(
        tx: f64,
        ty: f64,
        tz: f64,
        rot_x: f64,
        rot_y: f64,
        rot_z: f64,
        scale: f64,
        xs: &[f64],
        ys: &[f64],
        zs: &[f64],
        xs_out: &mut [f64],
        ys_out: &mut [f64],
        zs_out: &mut [f64],
    ) {
        let s1 = 1.0 + scale;
        for i in 0..xs.len() {
            let x = xs[i];
            let y = ys[i];
            let z = zs[i];
            xs_out[i] = s1 * x - rot_z * s1 * y + rot_y * s1 * z + tx;
            ys_out[i] = rot_z * s1 * x + s1 * y - rot_x * s1 * z + ty;
            zs_out[i] = -rot_y * s1 * x + rot_x * s1 * y + s1 * z + tz;
        }
    }

    /// Scalar fallback for 2-D affine batch transform.
    pub fn apply_affine_2d_batch(
        a00: f64,
        a01: f64,
        b0: f64,
        a10: f64,
        a11: f64,
        b1: f64,
        xs: &[f64],
        ys: &[f64],
        xs_out: &mut [f64],
        ys_out: &mut [f64],
    ) {
        for i in 0..xs.len() {
            xs_out[i] = a00 * xs[i] + a01 * ys[i] + b0;
            ys_out[i] = a10 * xs[i] + a11 * ys[i] + b1;
        }
    }
}

pub use inner::{apply_affine_2d_batch, apply_helmert_3_batch};

/// Apply a 7-parameter Helmert transform, **position-vector** convention, to
/// N coordinate triples.
///
/// `rot_x`, `rot_y`, `rot_z` — rotations in **radians** (linearised).
/// `scale` — dimensionless scale factor (e.g. pass `scale_ppm * 1e-6`).
///
/// The linearised rotation matrix (position-vector convention):
/// ```text
/// R = (1 + scale) × [[  1,  −rz,  ry ],
///                    [ rz,    1,  −rx ],
///                    [−ry,   rx,    1 ]]
/// ```
/// Transformed point: `X′ = R·X + T`
///
/// This is the explicitly-named, recommended entry point for this
/// convention; see the module-level "Rotation convention" section for how it
/// differs from [`crate::helmert_batch::helmert7_batch_scalar`]
/// (coordinate-frame convention, ppm scale).
#[allow(clippy::too_many_arguments)]
pub fn apply_helmert_7_batch_position_vector(
    tx: f64,
    ty: f64,
    tz: f64,
    rot_x: f64,
    rot_y: f64,
    rot_z: f64,
    scale: f64,
    xs: &[f64],
    ys: &[f64],
    zs: &[f64],
    xs_out: &mut [f64],
    ys_out: &mut [f64],
    zs_out: &mut [f64],
) {
    inner::apply_helmert_7_batch(
        tx, ty, tz, rot_x, rot_y, rot_z, scale, xs, ys, zs, xs_out, ys_out, zs_out,
    )
}

/// Deprecated alias for [`apply_helmert_7_batch_position_vector`].
///
/// This name is ambiguous: [`crate::helmert_batch::helmert7_batch_scalar`] is
/// a near-identically-purposed public batch Helmert API that uses the
/// *opposite* (coordinate-frame) rotation convention with a ppm scale.
/// Passing the same rotation values to the two functions silently flips the
/// sign of every rotation term. Call
/// [`apply_helmert_7_batch_position_vector`] (identical behaviour, explicit
/// name) instead.
#[deprecated(
    since = "0.1.2",
    note = "ambiguous name shared with a differently-conventioned sibling API (helmert_batch::helmert7_batch_scalar, coordinate-frame + ppm scale); call `apply_helmert_7_batch_position_vector` instead for the identical position-vector + dimensionless-scale behaviour, spelled out explicitly"
)]
#[allow(clippy::too_many_arguments)]
pub fn apply_helmert_7_batch(
    tx: f64,
    ty: f64,
    tz: f64,
    rot_x: f64,
    rot_y: f64,
    rot_z: f64,
    scale: f64,
    xs: &[f64],
    ys: &[f64],
    zs: &[f64],
    xs_out: &mut [f64],
    ys_out: &mut [f64],
    zs_out: &mut [f64],
) {
    apply_helmert_7_batch_position_vector(
        tx, ty, tz, rot_x, rot_y, rot_z, scale, xs, ys, zs, xs_out, ys_out, zs_out,
    )
}

#[cfg(test)]
mod tests {
    #![allow(clippy::too_many_arguments)]
    use super::{
        apply_affine_2d_batch, apply_helmert_3_batch, apply_helmert_7_batch_position_vector,
    };

    fn assert_close(a: f64, b: f64, tol: f64, label: &str) {
        assert!(
            (a - b).abs() < tol,
            "{label}: got {a:.15e}, expected {b:.15e}, diff {:.3e}",
            (a - b).abs()
        );
    }

    // ── Scalar reference implementations ────────────────────────────────────
    fn h3_ref(
        tx: f64,
        ty: f64,
        tz: f64,
        xs: &[f64],
        ys: &[f64],
        zs: &[f64],
    ) -> (Vec<f64>, Vec<f64>, Vec<f64>) {
        let n = xs.len();
        let mut xo = vec![0.0; n];
        let mut yo = vec![0.0; n];
        let mut zo = vec![0.0; n];
        for i in 0..n {
            xo[i] = xs[i] + tx;
            yo[i] = ys[i] + ty;
            zo[i] = zs[i] + tz;
        }
        (xo, yo, zo)
    }

    fn h7_ref(
        tx: f64,
        ty: f64,
        tz: f64,
        rx: f64,
        ry: f64,
        rz: f64,
        scale: f64,
        xs: &[f64],
        ys: &[f64],
        zs: &[f64],
    ) -> (Vec<f64>, Vec<f64>, Vec<f64>) {
        let n = xs.len();
        let s1 = 1.0 + scale;
        let mut xo = vec![0.0; n];
        let mut yo = vec![0.0; n];
        let mut zo = vec![0.0; n];
        for i in 0..n {
            let x = xs[i];
            let y = ys[i];
            let z = zs[i];
            xo[i] = s1 * x - rz * s1 * y + ry * s1 * z + tx;
            yo[i] = rz * s1 * x + s1 * y - rx * s1 * z + ty;
            zo[i] = -ry * s1 * x + rx * s1 * y + s1 * z + tz;
        }
        (xo, yo, zo)
    }

    fn aff2d_ref(
        a00: f64,
        a01: f64,
        b0: f64,
        a10: f64,
        a11: f64,
        b1: f64,
        xs: &[f64],
        ys: &[f64],
    ) -> (Vec<f64>, Vec<f64>) {
        let n = xs.len();
        let mut xo = vec![0.0; n];
        let mut yo = vec![0.0; n];
        for i in 0..n {
            xo[i] = a00 * xs[i] + a01 * ys[i] + b0;
            yo[i] = a10 * xs[i] + a11 * ys[i] + b1;
        }
        (xo, yo)
    }

    // ── Helmert 3-parameter tests ────────────────────────────────────────────

    #[test]
    fn helmert3_batch_eight_points() {
        let xs = [1.0_f64, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0];
        let ys = [10.0_f64, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0];
        let zs = [100.0_f64, 200.0, 300.0, 400.0, 500.0, 600.0, 700.0, 800.0];
        let (tx, ty, tz) = (5.0_f64, 10.0, 20.0);
        let n = xs.len();
        let mut xs_out = vec![0.0; n];
        let mut ys_out = vec![0.0; n];
        let mut zs_out = vec![0.0; n];
        apply_helmert_3_batch(
            tx,
            ty,
            tz,
            &xs,
            &ys,
            &zs,
            &mut xs_out,
            &mut ys_out,
            &mut zs_out,
        );
        let (xr, yr, zr) = h3_ref(tx, ty, tz, &xs, &ys, &zs);
        for i in 0..n {
            assert_close(xs_out[i], xr[i], 1e-12, &format!("h3x[{i}]"));
            assert_close(ys_out[i], yr[i], 1e-12, &format!("h3y[{i}]"));
            assert_close(zs_out[i], zr[i], 1e-12, &format!("h3z[{i}]"));
        }
    }

    #[test]
    fn helmert3_batch_empty() {
        let mut xs_out: Vec<f64> = Vec::new();
        let mut ys_out: Vec<f64> = Vec::new();
        let mut zs_out: Vec<f64> = Vec::new();
        apply_helmert_3_batch(
            1.0,
            2.0,
            3.0,
            &[],
            &[],
            &[],
            &mut xs_out,
            &mut ys_out,
            &mut zs_out,
        );
        assert!(xs_out.is_empty());
        assert!(ys_out.is_empty());
        assert!(zs_out.is_empty());
    }

    #[test]
    fn helmert3_batch_len_three() {
        let xs = [100.0_f64, 200.0, 300.0];
        let ys = [400.0_f64, 500.0, 600.0];
        let zs = [700.0_f64, 800.0, 900.0];
        let (tx, ty, tz) = (1.5_f64, 2.5, 3.5);
        let n = 3;
        let mut xs_out = vec![0.0; n];
        let mut ys_out = vec![0.0; n];
        let mut zs_out = vec![0.0; n];
        apply_helmert_3_batch(
            tx,
            ty,
            tz,
            &xs,
            &ys,
            &zs,
            &mut xs_out,
            &mut ys_out,
            &mut zs_out,
        );
        let (xr, yr, zr) = h3_ref(tx, ty, tz, &xs, &ys, &zs);
        for i in 0..n {
            assert_close(xs_out[i], xr[i], 1e-12, &format!("h3_3x[{i}]"));
            assert_close(ys_out[i], yr[i], 1e-12, &format!("h3_3y[{i}]"));
            assert_close(zs_out[i], zr[i], 1e-12, &format!("h3_3z[{i}]"));
        }
    }

    #[test]
    fn helmert3_batch_len_seven() {
        let n = 7;
        let xs: Vec<f64> = (1..=n).map(|i| i as f64 * 100.0).collect();
        let ys: Vec<f64> = (1..=n).map(|i| i as f64 * 200.0).collect();
        let zs: Vec<f64> = (1..=n).map(|i| i as f64 * 300.0).collect();
        let (tx, ty, tz) = (7.7_f64, 8.8, 9.9);
        let mut xs_out = vec![0.0; n];
        let mut ys_out = vec![0.0; n];
        let mut zs_out = vec![0.0; n];
        apply_helmert_3_batch(
            tx,
            ty,
            tz,
            &xs,
            &ys,
            &zs,
            &mut xs_out,
            &mut ys_out,
            &mut zs_out,
        );
        let (xr, yr, zr) = h3_ref(tx, ty, tz, &xs, &ys, &zs);
        for i in 0..n {
            assert_close(xs_out[i], xr[i], 1e-12, &format!("h3_7x[{i}]"));
            assert_close(ys_out[i], yr[i], 1e-12, &format!("h3_7y[{i}]"));
            assert_close(zs_out[i], zr[i], 1e-12, &format!("h3_7z[{i}]"));
        }
    }

    // ── Helmert 7-parameter tests ────────────────────────────────────────────

    // Approximate WGS84 → ETRS89 parameters (EUREF Permanent Network, epoch 1989.0)
    const TX: f64 = 0.565; // m
    const TY: f64 = 0.052; // m
    const TZ: f64 = -0.500; // m
    const RX: f64 = 0.891e-6; // rad (linearised)
    const RY: f64 = 5.390e-6; // rad
    const RZ: f64 = -8.712e-6; // rad
    const SC: f64 = -0.0088e-6; // dimensionless (0.0088 ppb)

    #[test]
    fn helmert7_batch_eight_points_vs_scalar() {
        let n = 8;
        let xs: Vec<f64> = (0..n).map(|i| 3_000_000.0 + i as f64 * 100_000.0).collect();
        let ys: Vec<f64> = (0..n).map(|i| 1_000_000.0 + i as f64 * 50_000.0).collect();
        let zs: Vec<f64> = (0..n).map(|i| 5_000_000.0 + i as f64 * 80_000.0).collect();
        let mut xs_out = vec![0.0; n];
        let mut ys_out = vec![0.0; n];
        let mut zs_out = vec![0.0; n];
        apply_helmert_7_batch_position_vector(
            TX,
            TY,
            TZ,
            RX,
            RY,
            RZ,
            SC,
            &xs,
            &ys,
            &zs,
            &mut xs_out,
            &mut ys_out,
            &mut zs_out,
        );
        let (xr, yr, zr) = h7_ref(TX, TY, TZ, RX, RY, RZ, SC, &xs, &ys, &zs);
        for i in 0..n {
            assert_close(xs_out[i], xr[i], 1e-9, &format!("h7x[{i}]"));
            assert_close(ys_out[i], yr[i], 1e-9, &format!("h7y[{i}]"));
            assert_close(zs_out[i], zr[i], 1e-9, &format!("h7z[{i}]"));
        }
    }

    #[test]
    fn helmert7_batch_empty() {
        let mut xs_out: Vec<f64> = Vec::new();
        let mut ys_out: Vec<f64> = Vec::new();
        let mut zs_out: Vec<f64> = Vec::new();
        apply_helmert_7_batch_position_vector(
            1.0,
            2.0,
            3.0,
            0.0,
            0.0,
            0.0,
            0.0,
            &[],
            &[],
            &[],
            &mut xs_out,
            &mut ys_out,
            &mut zs_out,
        );
        assert!(xs_out.is_empty());
        assert!(ys_out.is_empty());
        assert!(zs_out.is_empty());
    }

    #[test]
    fn helmert7_batch_len_three() {
        let n = 3;
        let xs = vec![4_000_000.0_f64, 3_500_000.0, 3_000_000.0];
        let ys = vec![500_000.0_f64, 800_000.0, 1_200_000.0];
        let zs = vec![4_800_000.0_f64, 4_900_000.0, 5_000_000.0];
        let (tx, ty, tz, rx, ry, rz, sc) = (1.0, -2.0, 3.0, 1e-7, 2e-7, -1e-7, 1.5e-9);
        let mut xs_out = vec![0.0; n];
        let mut ys_out = vec![0.0; n];
        let mut zs_out = vec![0.0; n];
        apply_helmert_7_batch_position_vector(
            tx,
            ty,
            tz,
            rx,
            ry,
            rz,
            sc,
            &xs,
            &ys,
            &zs,
            &mut xs_out,
            &mut ys_out,
            &mut zs_out,
        );
        let (xr, yr, zr) = h7_ref(tx, ty, tz, rx, ry, rz, sc, &xs, &ys, &zs);
        for i in 0..n {
            assert_close(xs_out[i], xr[i], 1e-9, &format!("h7_3x[{i}]"));
            assert_close(ys_out[i], yr[i], 1e-9, &format!("h7_3y[{i}]"));
            assert_close(zs_out[i], zr[i], 1e-9, &format!("h7_3z[{i}]"));
        }
    }

    #[test]
    fn helmert7_batch_len_seven() {
        let n = 7;
        let xs: Vec<f64> = (0..n).map(|i| 3_000_000.0 + i as f64 * 100_000.0).collect();
        let ys: Vec<f64> = (0..n).map(|i| 1_000_000.0 + i as f64 * 50_000.0).collect();
        let zs: Vec<f64> = (0..n).map(|i| 5_000_000.0 + i as f64 * 80_000.0).collect();
        let (tx, ty, tz, rx, ry, rz, sc) = (10.0, -5.0, 3.0, 1e-7, 2e-7, -3e-7, 2e-9);
        let mut xs_out = vec![0.0; n];
        let mut ys_out = vec![0.0; n];
        let mut zs_out = vec![0.0; n];
        apply_helmert_7_batch_position_vector(
            tx,
            ty,
            tz,
            rx,
            ry,
            rz,
            sc,
            &xs,
            &ys,
            &zs,
            &mut xs_out,
            &mut ys_out,
            &mut zs_out,
        );
        let (xr, yr, zr) = h7_ref(tx, ty, tz, rx, ry, rz, sc, &xs, &ys, &zs);
        for i in 0..n {
            assert_close(xs_out[i], xr[i], 1e-9, &format!("h7_7x[{i}]"));
            assert_close(ys_out[i], yr[i], 1e-9, &format!("h7_7y[{i}]"));
            assert_close(zs_out[i], zr[i], 1e-9, &format!("h7_7z[{i}]"));
        }
    }

    /// Regression test for the "two near-identically-named public batch
    /// Helmert APIs use opposite rotation conventions" finding: the
    /// deprecated bare `apply_helmert_7_batch` alias must still exist and
    /// must be byte-for-byte identical in behaviour to
    /// [`apply_helmert_7_batch_position_vector`] (it just forwards to it),
    /// so existing callers are not broken by the rename/deprecation.
    #[test]
    #[allow(deprecated)]
    fn deprecated_bare_alias_matches_explicit_position_vector_fn() {
        let n = 5;
        let xs: Vec<f64> = (0..n).map(|i| 3_000_000.0 + i as f64 * 100_000.0).collect();
        let ys: Vec<f64> = (0..n).map(|i| 1_000_000.0 + i as f64 * 50_000.0).collect();
        let zs: Vec<f64> = (0..n).map(|i| 5_000_000.0 + i as f64 * 80_000.0).collect();

        let mut xs_alias = vec![0.0; n];
        let mut ys_alias = vec![0.0; n];
        let mut zs_alias = vec![0.0; n];
        super::apply_helmert_7_batch(
            TX,
            TY,
            TZ,
            RX,
            RY,
            RZ,
            SC,
            &xs,
            &ys,
            &zs,
            &mut xs_alias,
            &mut ys_alias,
            &mut zs_alias,
        );

        let mut xs_explicit = vec![0.0; n];
        let mut ys_explicit = vec![0.0; n];
        let mut zs_explicit = vec![0.0; n];
        apply_helmert_7_batch_position_vector(
            TX,
            TY,
            TZ,
            RX,
            RY,
            RZ,
            SC,
            &xs,
            &ys,
            &zs,
            &mut xs_explicit,
            &mut ys_explicit,
            &mut zs_explicit,
        );

        for i in 0..n {
            // Both call paths delegate to the same underlying implementation
            // with the same inputs, so the results must be bit-identical;
            // a tiny tolerance is used only to avoid relying on exact `f64`
            // equality semantics in the test itself.
            assert_close(xs_alias[i], xs_explicit[i], 1e-15, &format!("alias_x[{i}]"));
            assert_close(ys_alias[i], ys_explicit[i], 1e-15, &format!("alias_y[{i}]"));
            assert_close(zs_alias[i], zs_explicit[i], 1e-15, &format!("alias_z[{i}]"));
        }
    }

    /// Regression test proving the two sibling batch Helmert APIs
    /// (`batch_helmert::apply_helmert_7_batch_position_vector`, this module,
    /// and `helmert_batch::helmert7_batch_scalar`) do in fact use opposite
    /// rotation-sign conventions for the same input rotation values, i.e.
    /// that the module docs' warning is accurate and would catch a future
    /// accidental convention unification/regression.
    #[test]
    fn position_vector_and_coordinate_frame_conventions_disagree_in_sign() {
        use crate::helmert_batch::{helmert7_batch_scalar, Helmert7Params};

        // A pure Z-axis rotation applied to a point offset only along Y: the
        // two conventions' rotation matrices differ only in the sign of
        // their off-diagonal terms, so this isolates a clean sign flip on
        // the X output (rather than a value that could coincidentally
        // cancel out for both conventions, as symmetric rx/ry/rz/xyz inputs
        // can).
        let (tx, ty, tz) = (0.0, 0.0, 0.0);
        let (rx, ry, rz) = (0.0, 0.0, 1e-3);
        let scale = 0.0;
        let xs = [0.0_f64];
        let ys = [1_000_000.0_f64];
        let zs = [0.0_f64];

        let mut pv_x = [0.0_f64];
        let mut pv_y = [0.0_f64];
        let mut pv_z = [0.0_f64];
        apply_helmert_7_batch_position_vector(
            tx, ty, tz, rx, ry, rz, scale, &xs, &ys, &zs, &mut pv_x, &mut pv_y, &mut pv_z,
        );

        let cf_params = Helmert7Params {
            tx,
            ty,
            tz,
            rx,
            ry,
            rz,
            scale_ppm: scale * 1e6,
        };
        let cf = helmert7_batch_scalar(&[[xs[0], ys[0], zs[0]]], &cf_params);

        // Position-vector: x' = x - rz*y + ry*z = -1e-3 * 1_000_000 = -1000.
        // Coordinate-frame: x' = x + rz*y - ry*z = +1e-3 * 1_000_000 = +1000.
        assert_close(pv_x[0], -1000.0, 1e-9, "pv_x");
        assert_close(cf[0][0], 1000.0, 1e-9, "cf_x");
        // The two conventions must disagree (opposite-signed off-diagonal
        // terms), proving they are genuinely different and not accidentally
        // aliased to the same behaviour.
        assert_close(pv_x[0], -cf[0][0], 1e-9, "pv_x == -cf_x");
    }

    // ── Affine 2-D tests ─────────────────────────────────────────────────────

    #[test]
    fn affine2d_batch_identity() {
        let n = 8;
        let xs: Vec<f64> = (0..n).map(|i| i as f64 * 100.0).collect();
        let ys: Vec<f64> = (0..n).map(|i| i as f64 * 200.0).collect();
        let mut xs_out = vec![0.0; n];
        let mut ys_out = vec![0.0; n];
        apply_affine_2d_batch(
            1.0,
            0.0,
            0.0,
            0.0,
            1.0,
            0.0,
            &xs,
            &ys,
            &mut xs_out,
            &mut ys_out,
        );
        for i in 0..n {
            assert_close(xs_out[i], xs[i], 1e-12, &format!("aff2d_id_x[{i}]"));
            assert_close(ys_out[i], ys[i], 1e-12, &format!("aff2d_id_y[{i}]"));
        }
    }

    #[test]
    fn affine2d_batch_translation() {
        let n = 8;
        let xs: Vec<f64> = (0..n).map(|i| i as f64 * 1000.0).collect();
        let ys: Vec<f64> = (0..n).map(|i| i as f64 * 500.0).collect();
        let (b0, b1) = (300.0_f64, 700.0);
        let mut xs_out = vec![0.0; n];
        let mut ys_out = vec![0.0; n];
        apply_affine_2d_batch(
            1.0,
            0.0,
            b0,
            0.0,
            1.0,
            b1,
            &xs,
            &ys,
            &mut xs_out,
            &mut ys_out,
        );
        let (xr, yr) = aff2d_ref(1.0, 0.0, b0, 0.0, 1.0, b1, &xs, &ys);
        for i in 0..n {
            assert_close(xs_out[i], xr[i], 1e-12, &format!("aff2d_tx_x[{i}]"));
            assert_close(ys_out[i], yr[i], 1e-12, &format!("aff2d_tx_y[{i}]"));
        }
    }

    #[test]
    fn affine2d_batch_scale_rotate_vs_scalar() {
        let n = 8;
        let theta = core::f64::consts::FRAC_PI_4; // 45 degrees
        let (a00, a01) = (theta.cos(), -theta.sin());
        let (a10, a11) = (theta.sin(), theta.cos());
        let xs: Vec<f64> = (0..n).map(|i| i as f64 * 250.0).collect();
        let ys: Vec<f64> = (0..n).map(|i| i as f64 * 150.0).collect();
        let mut xs_out = vec![0.0; n];
        let mut ys_out = vec![0.0; n];
        apply_affine_2d_batch(
            a00,
            a01,
            0.0,
            a10,
            a11,
            0.0,
            &xs,
            &ys,
            &mut xs_out,
            &mut ys_out,
        );
        let (xr, yr) = aff2d_ref(a00, a01, 0.0, a10, a11, 0.0, &xs, &ys);
        for i in 0..n {
            assert_close(xs_out[i], xr[i], 1e-12, &format!("aff2d_rot_x[{i}]"));
            assert_close(ys_out[i], yr[i], 1e-12, &format!("aff2d_rot_y[{i}]"));
        }
    }

    #[test]
    fn affine2d_batch_empty() {
        let mut xs_out: Vec<f64> = Vec::new();
        let mut ys_out: Vec<f64> = Vec::new();
        apply_affine_2d_batch(
            1.0,
            0.0,
            0.0,
            0.0,
            1.0,
            0.0,
            &[],
            &[],
            &mut xs_out,
            &mut ys_out,
        );
        assert!(xs_out.is_empty());
        assert!(ys_out.is_empty());
    }

    #[test]
    fn affine2d_batch_len_three() {
        let n = 3;
        let xs = vec![1.0_f64, 2.0, 3.0];
        let ys = vec![4.0_f64, 5.0, 6.0];
        let (a00, a01, b0, a10, a11, b1) = (2.0_f64, 1.0, 5.0, -1.0, 3.0, -2.0);
        let mut xs_out = vec![0.0; n];
        let mut ys_out = vec![0.0; n];
        apply_affine_2d_batch(
            a00,
            a01,
            b0,
            a10,
            a11,
            b1,
            &xs,
            &ys,
            &mut xs_out,
            &mut ys_out,
        );
        let (xr, yr) = aff2d_ref(a00, a01, b0, a10, a11, b1, &xs, &ys);
        for i in 0..n {
            assert_close(xs_out[i], xr[i], 1e-12, &format!("aff2d_3x[{i}]"));
            assert_close(ys_out[i], yr[i], 1e-12, &format!("aff2d_3y[{i}]"));
        }
    }

    #[test]
    fn affine2d_batch_len_seven() {
        let n = 7;
        let xs: Vec<f64> = (0..n).map(|i| i as f64 * 300.0).collect();
        let ys: Vec<f64> = (0..n).map(|i| i as f64 * 400.0).collect();
        let (a00, a01, b0, a10, a11, b1) = (1.5_f64, 0.5, 10.0, -0.5, 1.5, -10.0);
        let mut xs_out = vec![0.0; n];
        let mut ys_out = vec![0.0; n];
        apply_affine_2d_batch(
            a00,
            a01,
            b0,
            a10,
            a11,
            b1,
            &xs,
            &ys,
            &mut xs_out,
            &mut ys_out,
        );
        let (xr, yr) = aff2d_ref(a00, a01, b0, a10, a11, b1, &xs, &ys);
        for i in 0..n {
            assert_close(xs_out[i], xr[i], 1e-12, &format!("aff2d_7x[{i}]"));
            assert_close(ys_out[i], yr[i], 1e-12, &format!("aff2d_7y[{i}]"));
        }
    }
}