oxiproj-engine 0.1.1

//! Top-level construction and transformation entry points.
//!
//! Ported from PROJ 9.8.0 `src/4D_api.cpp` (`proj_create`, `proj_trans`,
//! `proj_trans_array`) and `src/pipeline.cpp` (pipeline assembly).

use crate::context::Context;
use crate::params::{parse, ParamList};
use crate::pipeline::Pipeline;
use crate::pj::Pj;
use oxiproj_core::{Coord, Direction, Ellipsoid, IoUnits, ProjError, ProjResult};

/// The keys that, if present, mean a parameter set carries its own ellipsoid.
const ELLIPSOID_KEYS: [&str; 9] = ["R", "a", "ellps", "datum", "rf", "f", "es", "e", "b"];

fn has_ellipsoid_def(pl: &ParamList) -> bool {
    ELLIPSOID_KEYS.iter().any(|k| pl.exists(k))
}

// Pipeline boundary i/o units mirror PROJ 9.8.0 src/internal.cpp
// `pj_left`/`pj_right`, consumed by src/pipeline.cpp
// (`P->left = pj_left(front); P->right = pj_right(back)`). A step that runs
// inverted swaps its raw left/right when reported in forward sense, and the
// `Classic` sentinel is normalized to `Projected`.

/// Effective forward-input units of a single step, mirroring PROJ
/// `pj_left` (src/internal.cpp): an inverted step swaps its raw left/right,
/// and `Classic` is normalized to `Projected`.
fn effective_left(s: &Pj) -> IoUnits {
    let u = if s.inverted { s.right } else { s.left };
    match u {
        IoUnits::Classic => IoUnits::Projected,
        other => other,
    }
}

/// Effective forward-output units of a single step, mirroring PROJ
/// `pj_right` (src/internal.cpp): an inverted step swaps its raw left/right,
/// and `Classic` is normalized to `Projected`.
fn effective_right(s: &Pj) -> IoUnits {
    let u = if s.inverted { s.left } else { s.right };
    match u {
        IoUnits::Classic => IoUnits::Projected,
        other => other,
    }
}

/// Construct a [`Pj`] from a proj-string.
///
/// Handles both single operations and `+proj=pipeline` definitions. Pipelines
/// resolve a global ellipsoid (GRS80 by default) inherited by steps lacking an
/// explicit ellipsoid, and are driven with prepare/finalize bypassed so each
/// step does its own unit handling.
pub fn create(proj_string: &str) -> ProjResult<Pj> {
    let ctx = Context::new();
    let params = parse(proj_string);

    let is_pipeline = params
        .entries
        .iter()
        .any(|(k, v)| k == "proj" && v.as_deref() == Some("pipeline"));

    if is_pipeline {
        create_pipeline(&params, &ctx)
    } else {
        create_single(&params, &ctx)
    }
}

/// Construct a [`Pj`] from a proj-string using a supplied [`Context`].
///
/// Unlike [`create`], this allows the caller to register grid data in the
/// context before construction so that `+nadgrids=` and `+hgridshift` steps
/// can resolve their grids.
pub fn create_with_ctx(proj_string: &str, ctx: &Context) -> ProjResult<Pj> {
    let params = parse(proj_string);

    let is_pipeline = params
        .entries
        .iter()
        .any(|(k, v)| k == "proj" && v.as_deref() == Some("pipeline"));

    if is_pipeline {
        create_pipeline(&params, ctx)
    } else {
        create_single(&params, ctx)
    }
}

fn extract_towgs84_from_params(params: &ParamList) -> Option<String> {
    // 1. Check for explicit +towgs84=... param
    if let Some(s) = params.get_str("towgs84") {
        if !s.is_empty() {
            return Some(s.to_string());
        }
    }
    // 2. Check +datum=<name> and look up in datum table
    if let Some(datum_name) = params.get_str("datum") {
        if !datum_name.is_empty() {
            if let Some(d) = oxiproj_core::find_datum(datum_name) {
                if let Some(rest) = d.defn.strip_prefix("towgs84=") {
                    return Some(rest.to_string());
                }
            }
        }
    }
    None
}

fn build_datum_shift_pipeline(
    params: &ParamList,
    towgs84_str: &str,
    ctx: &Context,
) -> ProjResult<Pj> {
    // Determine source ellipsoid name
    let src_ellps = params
        .get_str("ellps")
        .filter(|s| !s.is_empty())
        .map(|s| s.to_string())
        .or_else(|| {
            params
                .get_str("datum")
                .filter(|s| !s.is_empty())
                .and_then(|d| oxiproj_core::find_datum(d))
                .map(|d| d.ellipse_id.to_string())
        })
        .unwrap_or_else(|| "WGS84".to_string());

    // Parse towgs84 values to detect all-zero (identity)
    let vals: Vec<&str> = towgs84_str.split(',').collect();

    // All-zero means identity: skip datum shift injection
    let all_zero = vals
        .iter()
        .all(|v| matches!(v.trim(), "0" | "0.0" | "-0" | "-0.0"));
    if all_zero {
        return create_single_core(params, ctx);
    }

    // Build helmert params string
    let helmert_params = if vals.len() >= 7 {
        format!(
            "+proj=helmert +x={} +y={} +z={} +rx={} +ry={} +rz={} +s={} +convention=position_vector",
            vals[0].trim(),
            vals[1].trim(),
            vals[2].trim(),
            vals[3].trim(),
            vals[4].trim(),
            vals[5].trim(),
            vals[6].trim()
        )
    } else {
        format!(
            "+proj=helmert +x={} +y={} +z={}",
            vals.first().map(|v| v.trim()).unwrap_or("0"),
            vals.get(1).map(|v| v.trim()).unwrap_or("0"),
            vals.get(2).map(|v| v.trim()).unwrap_or("0")
        )
    };

    // Build the original projection step params (exclude towgs84, datum)
    let other_params: String = params
        .entries
        .iter()
        .filter(|(k, _)| k != "towgs84" && k != "datum")
        .map(|(k, v)| match v {
            Some(val) => format!("+{}={}", k, val),
            None => format!("+{}", k),
        })
        .collect::<Vec<_>>()
        .join(" ");

    // Synthetic pipeline:
    // step 1: cart (geodetic -> ECEF) using src ellipsoid
    // step 2: helmert (datum shift)
    // step 3: cart +inv (ECEF -> geodetic) using WGS84
    // step 4: the original projection
    let pipeline_str = format!(
        "+proj=pipeline \
         +step +proj=cart +ellps={src_ellps} \
         +step {helmert_params} \
         +step +proj=cart +inv +ellps=WGS84 \
         +step {other_params}"
    );

    let pipeline_params = parse(&pipeline_str);
    create_pipeline(&pipeline_params, ctx)
}

fn extract_nadgrids_from_params(params: &ParamList) -> Option<String> {
    if let Some(s) = params.get_str("nadgrids") {
        if !s.is_empty() {
            return Some(s.to_string());
        }
    }
    None
}

fn build_nadgrids_pipeline(params: &ParamList, grid_name: &str, ctx: &Context) -> ProjResult<Pj> {
    // Build the original projection step params (exclude nadgrids, datum)
    let other_params: String = params
        .entries
        .iter()
        .filter(|(k, _)| k != "nadgrids" && k != "datum")
        .map(|(k, v)| match v {
            Some(val) => format!("+{}={}", k, val),
            None => format!("+{}", k),
        })
        .collect::<Vec<_>>()
        .join(" ");

    // Synthetic pipeline:
    // step 1: hgridshift (geographic shift from source datum to WGS84, in radians)
    // step 2: the original projection
    let pipeline_str = format!(
        "+proj=pipeline \
         +step +proj=hgridshift +grids={grid_name} \
         +step {other_params}"
    );

    let pipeline_params = parse(&pipeline_str);
    create_pipeline(&pipeline_params, ctx)
}

fn create_single(params: &ParamList, ctx: &Context) -> ProjResult<Pj> {
    // Inject datum-shift pipeline when towgs84 params are present
    if let Some(towgs84_str) = extract_towgs84_from_params(params) {
        return build_datum_shift_pipeline(params, &towgs84_str, ctx);
    }
    // Inject hgridshift pipeline when nadgrids is present
    if let Some(grid_name) = extract_nadgrids_from_params(params) {
        return build_nadgrids_pipeline(params, &grid_name, ctx);
    }
    create_single_core(params, ctx)
}

fn create_single_core(params: &ParamList, ctx: &Context) -> ProjResult<Pj> {
    let name = params
        .entries
        .iter()
        .find(|(k, _)| k == "proj")
        .and_then(|(_, v)| v.as_deref())
        .ok_or(ProjError::MissingArg)?;

    let ellipsoid = crate::setup::setup_ellipsoid(params)?;
    let mut pj = crate::registry::build_single_op(name, params, ellipsoid, ctx)?;

    if params
        .entries
        .iter()
        .any(|(k, v)| (k == "inv" || k == "inverted") && v.is_none())
    {
        pj.inverted = true;
    }
    Ok(pj)
}

fn create_pipeline(params: &ParamList, ctx: &Context) -> ProjResult<Pj> {
    // Split entries into a global section (before the first `step`) and one
    // vec per subsequent step. The `step` markers are not retained.
    let mut global: Vec<(String, Option<String>)> = Vec::new();
    let mut steps_entries: Vec<Vec<(String, Option<String>)>> = Vec::new();
    let mut seen_step = false;
    for (k, v) in &params.entries {
        if k == "step" && v.is_none() {
            seen_step = true;
            steps_entries.push(Vec::new());
            continue;
        }
        if !seen_step {
            // Drop the pipeline marker from globals; it is never looked up.
            if k == "proj" && v.as_deref() == Some("pipeline") {
                continue;
            }
            global.push((k.clone(), v.clone()));
        } else if let Some(last) = steps_entries.last_mut() {
            last.push((k.clone(), v.clone()));
        }
    }

    let global_pl = ParamList { entries: global };

    // Global ellipsoid: explicit definition, else GRS80.
    //
    // GRS80 default per src/pipeline.cpp set_ellipsoid (a = 6378137,
    // f = 1/298.257222101).
    let global_ellipsoid = if has_ellipsoid_def(&global_pl) {
        crate::setup::setup_ellipsoid(&global_pl)?
    } else {
        Ellipsoid::from_a_rf(6378137.0, 298.257222101)?
    };

    let mut steps: Vec<Pj> = Vec::new();
    for step_vec in steps_entries {
        let step_pl = ParamList { entries: step_vec };
        let name = step_pl.get_str("proj").ok_or(ProjError::MissingArg)?;
        // get_str returns "" for a bare `proj`; that is still missing.
        if name.is_empty() {
            return Err(ProjError::MissingArg);
        }
        let inverted = step_pl
            .entries
            .iter()
            .any(|(k, v)| k == "inv" && v.is_none());
        let omit_fwd = step_pl
            .entries
            .iter()
            .any(|(k, v)| k == "omit_fwd" && v.is_none());
        let omit_inv = step_pl
            .entries
            .iter()
            .any(|(k, v)| k == "omit_inv" && v.is_none());
        let step_ellipsoid = if has_ellipsoid_def(&step_pl) {
            crate::setup::setup_ellipsoid(&step_pl)?
        } else {
            global_ellipsoid
        };
        let mut step_pj = crate::registry::build_single_op(name, &step_pl, step_ellipsoid, ctx)?;
        step_pj.inverted = inverted;
        step_pj.omit_fwd = omit_fwd;
        step_pj.omit_inv = omit_inv;
        steps.push(step_pj);
    }

    if steps.is_empty() {
        return Err(ProjError::MissingArg);
    }

    // Apply pipeline optimization: eliminate redundant/cancelling steps.
    // Run after the empty-pipeline guard (a truly empty pipeline string is an
    // error regardless) but before computing boundary i/o units (those must
    // reflect the post-optimization boundary steps).
    let steps = crate::pipeline_opt::optimize_pipeline(steps);

    if steps.is_empty() {
        // All steps were optimized away (e.g. every step was a complete no-op).
        return Err(ProjError::MissingArg);
    }

    let top_inverted = global_pl
        .entries
        .iter()
        .any(|(k, v)| k == "inv" && v.is_none());

    // Pipeline i/o units mirror PROJ src/pipeline.cpp
    // (`P->left = pj_left(front); P->right = pj_right(back)`): take the
    // EFFECTIVE units of the boundary steps so an inverted boundary step
    // (e.g. `+proj=cart +inv`, whose effective output is geographic radians)
    // propagates correctly to the pipeline's forward-sense units.
    let left = steps
        .first()
        .map(effective_left)
        .unwrap_or(IoUnits::Whatever);
    let right = steps
        .last()
        .map(effective_right)
        .unwrap_or(IoUnits::Whatever);

    Ok(Pj {
        operation: Box::new(Pipeline { steps }),
        ellipsoid: global_ellipsoid,
        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: false,
        left,
        right,
        inverted: top_inverted,
        bypass_prepare_finalize: true,
        omit_fwd: false,
        omit_inv: false,
        ad_proj: None,
        op_name: String::new(),
    })
}

/// Transform a single coordinate in the given direction.
///
/// The `inverted` flag is honored inside [`Pj::forward`]/[`Pj::inverse`], so the
/// direction is not swapped here.
pub fn trans(pj: &Pj, dir: Direction, c: Coord) -> ProjResult<Coord> {
    match dir {
        Direction::Fwd => pj.forward(c),
        Direction::Inv => pj.inverse(c),
        Direction::Ident => Ok(c),
    }
}

/// Transform a slice of coordinates in place.
///
/// A failed coordinate becomes [`Coord::error`]; the loop does not abort early.
pub fn trans_array(pj: &Pj, dir: Direction, coords: &mut [Coord]) -> ProjResult<()> {
    for c in coords.iter_mut() {
        match trans(pj, dir, *c) {
            Ok(r) => *c = r,
            Err(_) => *c = Coord::error(),
        }
    }
    Ok(())
}

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

    #[test]
    fn merc_forward_and_round_trip() {
        let pj = create("+proj=merc +ellps=WGS84").unwrap();
        let input = Coord::new(12.0 * DEG_TO_RAD, 55.0 * DEG_TO_RAD, 0.0, 0.0);
        let fwd = trans(&pj, Direction::Fwd, input).unwrap();
        let f = fwd.v();
        assert!((f[0] - 1335833.8895192828).abs() < 1e-6, "x got {}", f[0]);
        assert!(
            (f[1] - 7_326_837.715_045_549).abs() < 1e-6,
            "y got {}",
            f[1]
        );
        let inv = trans(&pj, Direction::Inv, fwd).unwrap();
        let i = inv.v();
        assert!((i[0] - 12.0 * DEG_TO_RAD).abs() < 1e-9, "lam got {}", i[0]);
        assert!((i[1] - 55.0 * DEG_TO_RAD).abs() < 1e-9, "phi got {}", i[1]);
    }

    #[test]
    fn utm_known_values() {
        let pj = create("+proj=utm +zone=32 +ellps=WGS84").unwrap();
        let at_origin = trans(
            &pj,
            Direction::Fwd,
            Coord::new(9.0 * DEG_TO_RAD, 0.0, 0.0, 0.0),
        )
        .unwrap();
        let o = at_origin.v();
        assert!((o[0] - 500000.0).abs() < 1e-6, "x got {}", o[0]);
        assert!((o[1] - 0.0).abs() < 1e-6, "y got {}", o[1]);

        let p = trans(
            &pj,
            Direction::Fwd,
            Coord::new(12.0 * DEG_TO_RAD, 55.0 * DEG_TO_RAD, 0.0, 0.0),
        )
        .unwrap();
        let pv = p.v();
        assert!(
            (pv[0] - 691_875.632_137_542).abs() < 1e-6,
            "x got {}",
            pv[0]
        );
        assert!(
            (pv[1] - 6_098_907.825_129_169).abs() < 1e-6,
            "y got {}",
            pv[1]
        );
    }

    #[test]
    fn etmerc_central_meridian() {
        let pj = create("+proj=etmerc +lon_0=9 +ellps=WGS84").unwrap();
        let out = trans(
            &pj,
            Direction::Fwd,
            Coord::new(9.0 * DEG_TO_RAD, 50.0 * DEG_TO_RAD, 0.0, 0.0),
        )
        .unwrap();
        let o = out.v();
        assert!(o[0].abs() < 1e-6, "x got {}", o[0]);
        assert!(
            (o[1] - 5_540_847.041_684_148).abs() < 1e-6,
            "y got {}",
            o[1]
        );
    }

    #[test]
    fn pipeline_utm_round_trip() {
        let pj = create(
            "+proj=pipeline +step +proj=utm +zone=32 +ellps=WGS84 +step +proj=utm +zone=32 +ellps=WGS84 +inv",
        )
        .unwrap();
        let out = trans(
            &pj,
            Direction::Fwd,
            Coord::new(12.0 * DEG_TO_RAD, 55.0 * DEG_TO_RAD, 0.0, 0.0),
        )
        .unwrap();
        let o = out.v();
        assert!((o[0] - 12.0 * DEG_TO_RAD).abs() < 1e-9, "lam got {}", o[0]);
        assert!((o[1] - 55.0 * DEG_TO_RAD).abs() < 1e-9, "phi got {}", o[1]);
    }

    #[test]
    fn trans_array_maps_in_place() {
        let pj = create("+proj=merc +ellps=WGS84").unwrap();
        let mut coords = [
            Coord::new(12.0 * DEG_TO_RAD, 55.0 * DEG_TO_RAD, 0.0, 0.0),
            Coord::new(0.0, 0.0, 0.0, 0.0),
        ];
        trans_array(&pj, Direction::Fwd, &mut coords).unwrap();
        assert!((coords[0].v()[0] - 1335833.8895192828).abs() < 1e-6);
        assert!(coords[1].v()[0].abs() < 1e-6);
    }

    /// Build a minimal NTv2 file in memory (uniform zero shift) for test use.
    ///
    /// Covers lon [-1°, 0°] and lat [0°, 1°] — a 2×2 grid with zero shift so
    /// hgridshift is an identity over that region.
    fn make_zero_shift_ntv2() -> Vec<u8> {
        let mut buf = Vec::new();

        // File overview header (11 records × 16 bytes = 176 bytes)
        // Record 0: NUM_OREC
        buf.extend_from_slice(b"NUM_OREC");
        buf.extend_from_slice(&11i32.to_le_bytes());
        buf.extend_from_slice(&[0u8; 4]);
        // Record 1: NUM_SREC
        buf.extend_from_slice(b"NUM_SREC");
        buf.extend_from_slice(&11i32.to_le_bytes());
        buf.extend_from_slice(&[0u8; 4]);
        // Record 2: NUM_FILE (1 subfile)
        buf.extend_from_slice(b"NUM_FILE");
        buf.extend_from_slice(&1u32.to_le_bytes());
        buf.extend_from_slice(&[0u8; 4]);
        // Record 3: GS_TYPE = "SECONDS "
        buf.extend_from_slice(b"GS_TYPE ");
        buf.extend_from_slice(b"SECONDS ");
        // Records 4-10: unused filler (7 × 16 = 112 bytes)
        buf.extend_from_slice(&[0u8; 112]);

        // Sub-grid header (11 records × 16 bytes = 176 bytes)
        // Record 0: SUB_NAME
        buf.extend_from_slice(b"SUB_NAME");
        buf.extend_from_slice(b"TESTGRID");
        // Record 1: PARENT
        buf.extend_from_slice(b"PARENT  ");
        buf.extend_from_slice(b"NONE    ");
        // Record 2: CREATED
        buf.extend_from_slice(b"CREATED ");
        buf.extend_from_slice(b"20240101");
        // Record 3: UPDATED
        buf.extend_from_slice(b"UPDATED ");
        buf.extend_from_slice(b"20240101");
        // Record 4: S_LAT = 0 arcsec (equator)
        buf.extend_from_slice(b"S_LAT   ");
        buf.extend_from_slice(&0.0f64.to_le_bytes());
        // Record 5: N_LAT = 3600 arcsec = 1 degree
        buf.extend_from_slice(b"N_LAT   ");
        buf.extend_from_slice(&3600.0f64.to_le_bytes());
        // Record 6: E_LONG = 0 arcsec (east boundary, positive-westward = 0°W = 0°E)
        buf.extend_from_slice(b"E_LONG  ");
        buf.extend_from_slice(&0.0f64.to_le_bytes());
        // Record 7: W_LONG = 3600 arcsec (west boundary, positive-westward = 1°W = -1°E)
        buf.extend_from_slice(b"W_LONG  ");
        buf.extend_from_slice(&3600.0f64.to_le_bytes());
        // Record 8: LAT_INC = 3600 arcsec (1 degree → 2 rows)
        buf.extend_from_slice(b"LAT_INC ");
        buf.extend_from_slice(&3600.0f64.to_le_bytes());
        // Record 9: LONG_INC = 3600 arcsec (1 degree → 2 cols)
        buf.extend_from_slice(b"LONG_INC");
        buf.extend_from_slice(&3600.0f64.to_le_bytes());
        // Record 10: GS_COUNT = 4 (2×2 grid)
        buf.extend_from_slice(b"GS_COUNT");
        buf.extend_from_slice(&4i32.to_le_bytes());
        buf.extend_from_slice(&[0u8; 4]);

        // Data: 4 cells × 4 f32 = 64 bytes (zero shift everywhere)
        for _ in 0..4 {
            buf.extend_from_slice(&0.0f32.to_le_bytes()); // lat shift
            buf.extend_from_slice(&0.0f32.to_le_bytes()); // lon shift
            buf.extend_from_slice(&0.0f32.to_le_bytes()); // lat accuracy
            buf.extend_from_slice(&0.0f32.to_le_bytes()); // lon accuracy
        }

        buf
    }

    #[test]
    fn nadgrids_pipeline_constructed() {
        // Verify that +nadgrids= on a single proj-string builds a pipeline
        // containing an hgridshift step, and that the grid is used for the shift.
        // We use a zero-shift grid so the round-trip is identity.
        let mut ctx = crate::context::Context::new();
        ctx.register_grid("test.gsb", make_zero_shift_ntv2());

        // merc with nadgrids — the point (lon=-0.5°, lat=0.5°) is inside the grid
        let pj = create_with_ctx("+proj=merc +ellps=WGS84 +nadgrids=test.gsb", &ctx).unwrap();

        // With zero shift, result should match plain merc
        let plain = create("+proj=merc +ellps=WGS84").unwrap();
        let coord = Coord::new(-0.5 * DEG_TO_RAD, 0.5 * DEG_TO_RAD, 0.0, 0.0);
        let shifted = trans(&pj, Direction::Fwd, coord).unwrap();
        let expected = trans(&plain, Direction::Fwd, coord).unwrap();
        let sv = shifted.v();
        let ev = expected.v();
        assert!(
            (sv[0] - ev[0]).abs() < 1e-3,
            "x: got {}, expected {}",
            sv[0],
            ev[0]
        );
        assert!(
            (sv[1] - ev[1]).abs() < 1e-3,
            "y: got {}, expected {}",
            sv[1],
            ev[1]
        );
    }

    #[test]
    fn pipeline_omit_fwd_skips_middle_step() {
        // Pipeline: noop -> axisswap(+omit_fwd) -> noop
        // Forward: middle step omitted => x,y unchanged
        // Inverse: middle step applied => x,y swapped
        let pj = create(
            "+proj=pipeline \
             +step +proj=noop \
             +step +proj=axisswap +order=2,1 +omit_fwd \
             +step +proj=noop",
        )
        .unwrap();

        let input = Coord::new(1.0, 2.0, 3.0, 0.0);

        // Forward: middle step omitted, so output equals input (noop+noop)
        let fwd = trans(&pj, Direction::Fwd, input).unwrap();
        let fv = fwd.v();
        assert!(
            (fv[0] - 1.0).abs() < 1e-12,
            "forward x should be 1.0, got {}",
            fv[0]
        );
        assert!(
            (fv[1] - 2.0).abs() < 1e-12,
            "forward y should be 2.0, got {}",
            fv[1]
        );

        // Inverse: middle step included (axisswap), so x and y are swapped
        let inv = trans(&pj, Direction::Inv, input).unwrap();
        let iv = inv.v();
        assert!(
            (iv[0] - 2.0).abs() < 1e-12,
            "inverse x should be 2.0 (swapped), got {}",
            iv[0]
        );
        assert!(
            (iv[1] - 1.0).abs() < 1e-12,
            "inverse y should be 1.0 (swapped), got {}",
            iv[1]
        );
    }

    #[test]
    fn optimizer_removes_double_omit_step() {
        // A step with both +omit_fwd and +omit_inv is never executed — the
        // pipeline optimizer removes it. The remaining noop steps let coordinates
        // pass through unchanged.
        let pj = create(
            "+proj=pipeline \
             +step +proj=noop \
             +step +proj=noop +omit_fwd +omit_inv \
             +step +proj=noop",
        )
        .expect("pipeline with double-omit middle step");
        let input = Coord::new(1.0, 2.0, 3.0, 0.0);
        let fwd = trans(&pj, Direction::Fwd, input).expect("forward through optimized pipeline");
        let fv = fwd.v();
        assert!(
            (fv[0] - 1.0).abs() < 1e-12,
            "x should pass through unchanged, got {}",
            fv[0]
        );
        assert!(
            (fv[1] - 2.0).abs() < 1e-12,
            "y should pass through unchanged, got {}",
            fv[1]
        );
    }

    #[test]
    fn optimizer_cancels_axisswap_pair() {
        // Two consecutive axisswap(2,1) steps cancel each other (axisswap is
        // self-inverse). After optimization the pipeline reduces to noop only,
        // and the input coordinate passes through unchanged.
        let pj = create(
            "+proj=pipeline \
             +step +proj=noop \
             +step +proj=axisswap +order=2,1 \
             +step +proj=axisswap +order=2,1",
        )
        .expect("pipeline with noop + cancelling axisswap pair");
        let input = Coord::new(3.0, 7.0, 0.0, 0.0);
        let fwd = trans(&pj, Direction::Fwd, input).expect("forward through optimized pipeline");
        let fv = fwd.v();
        assert!(
            (fv[0] - 3.0).abs() < 1e-12,
            "x should be unchanged after axisswap cancellation, got {}",
            fv[0]
        );
        assert!(
            (fv[1] - 7.0).abs() < 1e-12,
            "y should be unchanged after axisswap cancellation, got {}",
            fv[1]
        );
    }

    #[test]
    fn optimizer_all_cancelled_returns_error() {
        // Every step has both +omit_fwd and +omit_inv — all steps are removed by
        // the optimizer, leaving zero steps. This must return an error.
        let result = create(
            "+proj=pipeline \
             +step +proj=noop +omit_fwd +omit_inv \
             +step +proj=noop +omit_fwd +omit_inv",
        );
        assert!(
            result.is_err(),
            "all-cancelled pipeline must return an error"
        );
    }
}