rwarp 0.1.1

//! Coordinate transformers for the warp pipeline.
//!
//! Maps GDAL's transformer hierarchy (gdaltransformer.cpp, commit a3b7b01d3e):
//!
//!   Transformer trait          ← GDALTransformerFunc callback signature
//!   GenImgProjTransformer      ← GDALGenImgProjTransform
//!
//! Default gdalwarp architecture (from gdalwarp_lib.cpp line 3260):
//!
//!   ApproxTransformer {                    ← GDALApproxTransform (max_error=0.125)
//!     wraps: GenImgProjTransformer {       ← GDALGenImgProjTransform
//!       step 1: dst pixel → dst geo         (apply dst geotransform)
//!       step 2: dst geo → src geo           (PROJ coordinate transform)
//!       step 3: src geo → src pixel         (apply inv src geotransform)
//!     }
//!   }
//!
//! The ApproxTransformer wraps the GenImgProjTransform and not just
//! the reprojection step. The error threshold of 0.125 is in the output
//! space of GenImgProjTransform = source pixel coordinates (since the
//! warp kernel always calls with bDstToSrc=TRUE).
//!
//! Evidence from gdalwarp_lib.cpp:
//!   line 2841: pfnTransformer = GDALGenImgProjTransform
//!   line 3260: hTransformArg = GDALCreateApproxTransformer(
//!                  GDALGenImgProjTransform, hTransformArg, dfErrorThreshold)
//!   line 3263: pfnTransformer = GDALApproxTransform
//!   line 1598: dfErrorThreshold = 0.125 (default)
//!
//! Note: There is ALSO a "REPROJECTION_APPROX_ERROR_IN_DST_SRS_UNIT" option
//! that wraps just the reprojection inside GenImgProjTransformer, but this
//! is NEVER used by default gdalwarp. It's a separate code path for other
//! callers.
//!
//! ## Cell-centre convention
//!
//! GDAL's warp kernel constructs destination pixel coordinates as
//! `iDstX + 0.5` (cell centres, not top-left corners). The high-level
//! [`transform_scanline`] function handles this automatically. The raw
//! [`Transformer::transform`] method receives whatever coordinates the
//! caller provides — integer pixel indices give top-left corner coordinates,
//! not centres. This is a silent half-pixel shift that is easy to miss.
//!
//! ## Axis order
//!
//! [`proj::Proj::new_known_crs`] normalises axis order to easting/northing
//! (lon/lat) for EPSG codes. Raw PROJ strings or WKT2 definitions with
//! explicit lat/lon axis order may not be normalized. Prefer EPSG codes.
//!
//! ## Two PROJ objects
//!
//! Because the `proj` crate's `Proj::convert()` always transforms in the
//! forward direction of the pipeline it was constructed with, we create
//! separate `Proj` instances for `dst→src` and `src→dst`. This doubles
//! PROJ initialisation cost but is negligible for tile-serving use cases
//! where the transformer is created once and reused across tiles.

use vaster::inv_geotransform;

// ---------------------------------------------------------------------------
// Transformer trait
// ---------------------------------------------------------------------------

/// Trait for coordinate transformers.
///
/// Matches the semantics of GDAL's GDALTransformerFunc:
/// - Takes arrays of x, y coordinates (modified in place)
/// - Returns per-point success flags
/// - `dst_to_src=true`: input is dest pixel coords, output is src pixel coords
/// - `dst_to_src=false`: input is src pixel coords, output is dest pixel coords
///
/// For warp, the kernel always calls with dst_to_src=true.
///
/// The trait is intentionally minimal — it maps coordinate arrays without
/// knowing what data values are. This makes it the right seam for future
/// transformer plugins: curvilinear grids, displacement fields, GCP
/// polynomials, and others all implement this trait and compose unchanged
/// with `ApproxTransformer`, `compute_source_window`, and `warp_resample`.
///
/// For vector-quantity warps (velocity fields, wind, ocean currents) the
/// optional [`Transformer::jacobian`] method provides the local frame
/// rotation needed to correctly reproject `(u, v)` components. Scalar
/// resampling via `warp_resample` does not call `jacobian`.
pub trait Transformer {
    fn transform(
        &self,
        dst_to_src: bool,
        x: &mut [f64],
        y: &mut [f64],
    ) -> Vec<bool>;

    /// Local Jacobian of the transform at a single point.
    ///
    /// Returns the 2×2 matrix:
    /// ```text
    /// [[∂x_out/∂x_in,  ∂x_out/∂y_in],
    ///  [∂y_out/∂x_in,  ∂y_out/∂y_in]]
    /// ```
    /// evaluated at `(x, y)` in the input coordinate space.
    ///
    /// Used by vector-quantity-preserving warp kernels to rotate `(u, v)`
    /// fields from the source index-grid frame into the destination CRS
    /// frame. This is the generalized form of the north-arrow correction:
    /// "north" (or any reference direction) rotates relative to the display
    /// grid as you move across a projection, and `(u, v)` components stored
    /// in index-grid coordinates need the same correction when reprojected.
    ///
    /// Scalar resampling via `warp_resample` does not call this method.
    ///
    /// The default implementation uses central finite differences (step 1e-5).
    /// Implementors with analytic Jacobians (e.g. from `proj_factors`) should
    /// override for accuracy and performance.
    fn jacobian(&self, dst_to_src: bool, x: f64, y: f64) -> Option<[[f64; 2]; 2]> {
        let h = 1e-5_f64;

        let mut xp = [x + h, x - h];
        let mut yp = [y,     y    ];
        let okx = self.transform(dst_to_src, &mut xp, &mut yp);
        if !okx[0] || !okx[1] { return None; }
        let dxdx = (xp[0] - xp[1]) / (2.0 * h);
        let dydx = (yp[0] - yp[1]) / (2.0 * h);

        let mut xq = [x,     x    ];
        let mut yq = [y + h, y - h];
        let oky = self.transform(dst_to_src, &mut xq, &mut yq);
        if !oky[0] || !oky[1] { return None; }
        let dxdy = (xq[0] - xq[1]) / (2.0 * h);
        let dydy = (yq[0] - yq[1]) / (2.0 * h);

        Some([[dxdx, dxdy], [dydx, dydy]])
    }
}

// ---------------------------------------------------------------------------
// GenImgProjTransformer
// ---------------------------------------------------------------------------

/// The three-step pixel ↔ geo ↔ CRS pipeline.
///
/// Maps GDAL's GDALGenImgProjTransform (gdaltransformer.cpp line ~3100).
///
/// Composes:
///   1. Apply geotransform: pixel → geo (always exact)
///   2. PROJ: CRS → CRS (always exact, per-point)
///   3. Apply inverse geotransform: geo → pixel (always exact)
///
/// The ApproxTransformer wraps this entire struct from outside (in lib.rs),
/// matching gdalwarp's default architecture.
pub struct GenImgProjTransformer {
    // Source grid
    src_gt: [f64; 6],
    src_inv_gt: [f64; 6],

    // Destination grid
    dst_gt: [f64; 6],
    dst_inv_gt: [f64; 6],

    // PROJ coordinate transforms
    proj_dst_to_src: proj::Proj,
    proj_src_to_dst: proj::Proj,
}

impl GenImgProjTransformer {
    pub fn new(
        src_crs: &str,
        src_gt: [f64; 6],
        dst_crs: &str,
        dst_gt: [f64; 6],
    ) -> Result<Self, String> {
        let src_inv_gt = inv_geotransform(&src_gt)
            .ok_or_else(|| "Source geotransform not invertible".to_string())?;
        let dst_inv_gt = inv_geotransform(&dst_gt)
            .ok_or_else(|| "Dest geotransform not invertible".to_string())?;

        let proj_dst_to_src = proj::Proj::new_known_crs(dst_crs, src_crs, None)
            .map_err(|e| format!("PROJ dst→src failed: {}", e))?;
        let proj_src_to_dst = proj::Proj::new_known_crs(src_crs, dst_crs, None)
            .map_err(|e| format!("PROJ src→dst failed: {}", e))?;

        Ok(Self {
            src_gt,
            src_inv_gt,
            dst_gt,
            dst_inv_gt,
            proj_dst_to_src,
            proj_src_to_dst,
        })
    }
}

impl Transformer for GenImgProjTransformer {
    /// Transform coordinates through the three-step pipeline.
    ///
    /// When `dst_to_src = true` (the warp kernel path):
    ///   1. Apply dst geotransform: pixel → geo
    ///   2. PROJ: dst CRS → src CRS
    ///   3. Apply inverse src geotransform: geo → pixel
    fn transform(
        &self,
        dst_to_src: bool,
        x: &mut [f64],
        y: &mut [f64],
    ) -> Vec<bool> {
        let n = x.len();

        if dst_to_src {
            // Step 1: dst pixel → dst geo (line 3113-3142)
            for i in 0..n {
                let pixel = x[i];
                let line = y[i];
                x[i] = self.dst_gt[0] + pixel * self.dst_gt[1] + line * self.dst_gt[2];
                y[i] = self.dst_gt[3] + pixel * self.dst_gt[4] + line * self.dst_gt[5];
            }

            // Step 2: PROJ dst geo → src geo (line 3148-3152)
            let mut success = vec![true; n];
            for i in 0..n {
                match self.proj_dst_to_src.convert((x[i], y[i])) {
                    Ok((rx, ry)) => {
                        x[i] = rx;
                        y[i] = ry;
                    }
                    Err(_) => {
                        success[i] = false;
                    }
                }
            }

            // Step 3: src geo → src pixel (line 3158-3187)
            for i in 0..n {
                if !success[i] {
                    continue;
                }
                let geo_x = x[i];
                let geo_y = y[i];
                x[i] = self.src_inv_gt[0]
                    + geo_x * self.src_inv_gt[1]
                    + geo_y * self.src_inv_gt[2];
                y[i] = self.src_inv_gt[3]
                    + geo_x * self.src_inv_gt[4]
                    + geo_y * self.src_inv_gt[5];
            }

            success
        } else {
            // Reverse: src pixel → src geo → dst geo → dst pixel

            for i in 0..n {
                let pixel = x[i];
                let line = y[i];
                x[i] = self.src_gt[0] + pixel * self.src_gt[1] + line * self.src_gt[2];
                y[i] = self.src_gt[3] + pixel * self.src_gt[4] + line * self.src_gt[5];
            }

            let mut success = vec![true; n];
            for i in 0..n {
                match self.proj_src_to_dst.convert((x[i], y[i])) {
                    Ok((rx, ry)) => {
                        x[i] = rx;
                        y[i] = ry;
                    }
                    Err(_) => {
                        success[i] = false;
                    }
                }
            }

            for i in 0..n {
                if !success[i] {
                    continue;
                }
                let geo_x = x[i];
                let geo_y = y[i];
                x[i] = self.dst_inv_gt[0]
                    + geo_x * self.dst_inv_gt[1]
                    + geo_y * self.dst_inv_gt[2];
                y[i] = self.dst_inv_gt[3]
                    + geo_x * self.dst_inv_gt[4]
                    + geo_y * self.dst_inv_gt[5];
            }

            success
        }
    }
}

// ---------------------------------------------------------------------------
// Convenience: transform scanlines
// ---------------------------------------------------------------------------


/// Build destination pixel coordinates for a scanline and transform to source.
///
/// Matches GDAL's pattern at the top of GWKNearestThread (line 5560):
/// ```text
///   for (iDstX = 0; iDstX < nDstXSize; iDstX++) {
///       padfX[iDstX] = iDstX + 0.5 + nDstXOff;
///       padfY[iDstX] = iDstY + 0.5 + nDstYOff;
///   }
///   pfnTransformer(..., TRUE, nDstXSize, padfX, padfY, padfZ, pabSuccess);
/// ```
pub fn transform_scanline(
    transformer: &impl Transformer,
    dst_y: usize,
    dst_ncol: usize,
    dst_x_off: usize,
    dst_y_off: usize,
) -> (Vec<f64>, Vec<f64>, Vec<bool>) {
    let mut x: Vec<f64> = (0..dst_ncol)
        .map(|col| (col + dst_x_off) as f64 + 0.5)
        .collect();
    let mut y = vec![(dst_y + dst_y_off) as f64 + 0.5; dst_ncol];

    let success = transformer.transform(true, &mut x, &mut y);

    (x, y, success)
}

/// Compute source pixel coordinates for every destination pixel.
pub fn transform_grid(
    transformer: &impl Transformer,
    dst_dim: &[usize; 2],
) -> (Vec<f64>, Vec<f64>) {
    let n = dst_dim[0] * dst_dim[1];
    let mut all_x = Vec::with_capacity(n);
    let mut all_y = Vec::with_capacity(n);

    for row in 0..dst_dim[1] {
        let (sx, sy, ok) = transform_scanline(transformer, row, dst_dim[0], 0, 0);
        for col in 0..dst_dim[0] {
            if ok[col] {
                all_x.push(sx[col]);
                all_y.push(sy[col]);
            } else {
                all_x.push(f64::NAN);
                all_y.push(f64::NAN);
            }
        }
    }

    (all_x, all_y)
}