noobase 0.0.5

//! Surface-brightness-conserving exact image reprojection.
//!
//! The caller supplies an input image on its own pixel grid plus a corner
//! field that gives, for every node of the *output* pixel grid, the
//! corresponding location in the input image's pixel coordinates. The
//! reprojector then computes each output pixel by polygon-clipping the
//! input pixels that fall under the quadrilateral defined by the four
//! output corners and taking the area-weighted mean of the contributing
//! input values.
//!
//! Design notes:
//!
//! - **No WCS.** The corner field is plain pixel coordinates. The caller
//!   is expected to derive it from whichever WCS toolkit they prefer
//!   (astropy.wcs, gwcs, etc.) and hand the result in as an ndarray.
//!   Spherical curvature is therefore *not* corrected — the residual
//!   error is the curvature of the input projection over a single
//!   output pixel, which is negligible for any reasonable astronomical
//!   image.
//!
//! - **Pixel coordinate convention.** Integer `(x, y) = (column, row)`
//!   refers to the *center* of pixel `(row=y, column=x)` (matching
//!   astropy). Pixel `(i, j)` therefore occupies the continuous square
//!   `[j - 0.5, j + 0.5] x [i - 0.5, i + 0.5]`. Internally the
//!   reprojector applies a `+0.5` shift to every corner so input pixels
//!   align with integer-aligned unit cells `[j, j + 1] x [i, i + 1]`;
//!   this shift is invisible to the caller.
//!
//! - **Corner indexing.** Output pixel `(i_out, j_out)` is bounded by
//!   `pixel_corners[i_out, j_out]`, `[i_out, j_out + 1]`,
//!   `[i_out + 1, j_out + 1]`, `[i_out + 1, j_out]`, in that winding
//!   order. Adjacent output pixels share corners by construction, which
//!   makes the polygon partition watertight.
//!
//! - **NaN semantics.** Any NaN in the four corners of an output pixel
//!   produces `image = NaN, footprint = 0, weight = 0`. NaN values
//!   inside the contributing input pixels are excluded from both the
//!   numerator and the *valid* denominator (`weight`), but still count
//!   toward the purely geometric `footprint`.
//!
//! - **Surface brightness conservation.** The default reduction is the
//!   area-weighted mean over the valid (non-NaN) input area, which
//!   conserves surface brightness. For total flux, multiply the output
//!   image by `weight` and by the output pixel area in input-pixel
//!   units.
//!
//! - **Footprint vs weight.** The output exposes two coverage arrays.
//!   `footprint` is the pure geometric overlap fraction between the
//!   output pixel and the input image bounds (independent of data
//!   validity). `weight` is the fraction of the output pixel that
//!   contributed *valid* (non-NaN) input. Invariant: `weight <=
//!   footprint`; with no NaN inputs they are equal. The ratio
//!   `weight / footprint` (where `footprint > 0`) is the fraction of
//!   the covered region that was free of NaN.
//!
//! - **Optional error propagation.** When a 1-sigma error image is
//!   supplied alongside `image`, the output carries a propagated
//!   1-sigma error array. The propagation is first-order under the
//!   assumption that input pixel errors are independent: for an output
//!   pixel that averages input pixels with overlap areas `A_k` and
//!   errors `sigma_k` over the valid input area `W = sum(A_k)`, the
//!   output error is `sqrt(sum(A_k^2 * sigma_k^2)) / W`. Supplying the
//!   error image does **not** change `image`, `footprint`, or `weight`
//!   bit-for-bit — propagation is purely additive. Only the
//!   per-output-pixel marginal 1-sigma is produced; adjacent output
//!   pixels that share input pixels have correlated errors, and that
//!   off-diagonal covariance is not represented (the same caveat as
//!   `Spectrum::rebin` on upsampling). An input pixel whose value is
//!   valid but whose error is non-finite or negative poisons that
//!   output pixel's error to `NaN`: the pixel still enters the image
//!   mean, so silently dropping its unknown variance would bias the
//!   output error low. A value-NaN input contributes neither value nor
//!   error, so its error entry is never consulted.

use ndarray::{Array2, ArrayView2, ArrayView3, ArrayViewMut1, Axis};
use rayon::prelude::*;
use thiserror::Error;

use crate::float::Float;
use crate::image::polygon::{Point, clip_quad_against_unit_cell, signed_area};

/// Output of [`reproject_exact`].
#[derive(Debug, Clone, PartialEq)]
pub struct ReprojectExactOutput {
    /// Surface-brightness-conserving reprojection of `image`.
    /// Shape: `(H_out, W_out)`. Pixels that received no valid input
    /// contribution are `NaN`.
    pub image: Array2<f64>,
    /// Pure geometric overlap fraction in `[0, 1]`: the ratio of the
    /// in-bounds input area covered by the output pixel to the area of
    /// the output pixel projected into input-pixel space. Independent
    /// of data validity (NaN input pixels still count toward this).
    /// Shape: `(H_out, W_out)`.
    pub footprint: Array2<f64>,
    /// Valid overlap fraction in `[0, 1]`: same numerator as
    /// `footprint` but restricted to non-NaN input pixels. This is the
    /// effective weight that fed the surface-brightness average.
    /// Always `<= footprint`. Shape: `(H_out, W_out)`.
    pub weight: Array2<f64>,
    /// Propagated 1-sigma error of `image`, in the same units as
    /// `image`. `Some` with shape `(H_out, W_out)` when `error` was
    /// supplied to [`reproject_exact`]; `None` otherwise.
    ///
    /// First-order propagation under independent input-pixel errors:
    /// `sqrt(sum(A_k^2 * sigma_k^2)) / sum(A_k)` over the valid (non-NaN)
    /// input pixels that overlap the output pixel (`A_k` = overlap area,
    /// `sigma_k` = input 1-sigma). An output pixel is `NaN` when it has
    /// no valid contribution (matching `image`) or when a contributing
    /// input pixel has a valid value but a non-finite or negative error
    /// (such a pixel stays in the image mean, so its unknown variance
    /// cannot be silently dropped). Only the marginal per-pixel 1-sigma
    /// is represented; correlations between output pixels sharing input
    /// pixels are not.
    pub error: Option<Array2<f64>>,
}

/// Errors returned by [`reproject_exact`] for ill-shaped inputs.
///
/// Algorithmic edge cases (out-of-footprint output pixels, NaN corners,
/// NaN input pixels, degenerate quadrilaterals) are *not* errors — they
/// fall out of the NaN / zero-weight conventions documented on
/// [`reproject_exact`].
#[derive(Debug, Error, PartialEq)]
pub enum ReprojectError {
    #[error(
        "pixel_corners must have shape (H_out + 1, W_out + 1, 2) with last dim == 2 and both grid dims >= 2; got ({rows}, {cols}, {last})"
    )]
    PixelCornersShape {
        rows: usize,
        cols: usize,
        last: usize,
    },
    #[error(
        "error must have the same shape as image ({image_rows}, {image_cols}); got ({error_rows}, {error_cols})"
    )]
    ErrorShapeMismatch {
        image_rows: usize,
        image_cols: usize,
        error_rows: usize,
        error_cols: usize,
    },
}

/// Surface-brightness-conserving exact reprojection of a 2-D image.
///
/// See the module-level documentation for the coordinate convention,
/// corner indexing rule, NaN semantics, and conservation property.
///
/// # Parameters
///
/// - `image`: input image of shape `(H_in, W_in)`. NaN values are
///   treated as mask (excluded from both numerator and denominator).
///   May be `f32` or `f64`; internally upcast to `f64` for the kernel.
/// - `pixel_corners`: shape `(H_out + 1, W_out + 1, 2)`. The last
///   dimension stores `[x_in, y_in]`, i.e. the location of the
///   corresponding output-pixel corner in the input image's continuous
///   pixel coordinates (integer = pixel center). NaN corners propagate
///   to `(NaN, 0)` in the output for any output pixel that touches
///   them.
/// - `error`: optional 1-sigma error image, same shape and dtype as
///   `image`. When `Some`, the output carries a propagated `error`
///   array; when `None`, `error` is `None` and there is no extra work.
///   Supplying it does not change `image`, `footprint`, or `weight`.
///
/// # Returns
///
/// [`ReprojectExactOutput`] containing `image`, `footprint`, and
/// `weight` (each of shape `(H_out, W_out)`, dtype `f64`), plus `error`
/// which is `Some` of the same shape exactly when `error` was
/// supplied.
///
/// # Errors
///
/// Returns [`ReprojectError::PixelCornersShape`] if `pixel_corners`
/// does not have last dim `2`, or if either of its front dimensions is
/// below `2` (which would produce zero output pixels in that
/// direction). Returns [`ReprojectError::ErrorShapeMismatch`] if
/// `error` is supplied with a shape different from `image`.
pub fn reproject_exact<T: Float>(
    image: ArrayView2<T>,
    pixel_corners: ArrayView3<f64>,
    error: Option<ArrayView2<T>>,
) -> Result<ReprojectExactOutput, ReprojectError> {
    let corner_shape = pixel_corners.shape();
    let rows = corner_shape[0];
    let cols = corner_shape[1];
    let last = corner_shape[2];
    if last != 2 || rows < 2 || cols < 2 {
        return Err(ReprojectError::PixelCornersShape { rows, cols, last });
    }

    // The error image, when present, must align pixel-for-pixel with the
    // input image so the kernel can read both at the same index.
    if let Some(error_view) = error.as_ref() {
        let image_shape = image.shape();
        let error_shape = error_view.shape();
        if image_shape != error_shape {
            return Err(ReprojectError::ErrorShapeMismatch {
                image_rows: image_shape[0],
                image_cols: image_shape[1],
                error_rows: error_shape[0],
                error_cols: error_shape[1],
            });
        }
    }

    // Upcast input image to f64 once, up front. The kernel is f64-only
    // because the polygon math is f64; carrying T into the kernel would
    // either force a generic kernel (more code, no benefit) or repeat
    // the cast inside the hot loop. The error image (if any) is upcast
    // the same way.
    let image_f64: Array2<f64> = image.mapv(|value| value.to_f64().unwrap_or(f64::NAN));
    let error_f64: Option<Array2<f64>> =
        error.map(|view| view.mapv(|value| value.to_f64().unwrap_or(f64::NAN)));

    let height_out = rows - 1;
    let width_out = cols - 1;
    let mut image_out = Array2::<f64>::from_elem((height_out, width_out), f64::NAN);
    let mut footprint = Array2::<f64>::zeros((height_out, width_out));
    let mut weight = Array2::<f64>::zeros((height_out, width_out));

    // Rows are independent: each output row only reads `image_f64`
    // (and `error_f64`) plus the slice of `pixel_corners` between
    // rows `i_out` and `i_out + 1`, and only writes its own row of the
    // outputs. Drive them in parallel via rayon. The error path is
    // branched out so the no-error call allocates nothing extra and
    // does zero per-pixel error work.
    let error_out = match error_f64.as_ref() {
        Some(error_f64_view) => {
            let mut error_out = Array2::<f64>::from_elem((height_out, width_out), f64::NAN);
            image_out
                .axis_iter_mut(Axis(0))
                .into_par_iter()
                .zip(footprint.axis_iter_mut(Axis(0)).into_par_iter())
                .zip(weight.axis_iter_mut(Axis(0)).into_par_iter())
                .zip(error_out.axis_iter_mut(Axis(0)).into_par_iter())
                .enumerate()
                .for_each(
                    |(row_index, (((row_image, row_footprint), row_weight), row_error))| {
                        process_output_row(
                            row_index,
                            row_image,
                            row_footprint,
                            row_weight,
                            Some(row_error),
                            image_f64.view(),
                            Some(error_f64_view.view()),
                            pixel_corners,
                        );
                    },
                );
            Some(error_out)
        }
        None => {
            image_out
                .axis_iter_mut(Axis(0))
                .into_par_iter()
                .zip(footprint.axis_iter_mut(Axis(0)).into_par_iter())
                .zip(weight.axis_iter_mut(Axis(0)).into_par_iter())
                .enumerate()
                .for_each(|(row_index, ((row_image, row_footprint), row_weight))| {
                    process_output_row(
                        row_index,
                        row_image,
                        row_footprint,
                        row_weight,
                        None,
                        image_f64.view(),
                        None,
                        pixel_corners,
                    );
                });
            None
        }
    };

    Ok(ReprojectExactOutput {
        image: image_out,
        footprint,
        weight,
        error: error_out,
    })
}

/// Sequential driver retained for tests. Used by the parallel-vs-
/// sequential equivalence test to confirm the rayon driver does not
/// change the result bitwise.
#[cfg(test)]
fn reproject_exact_sequential<T: Float>(
    image: ArrayView2<T>,
    pixel_corners: ArrayView3<f64>,
    error: Option<ArrayView2<T>>,
) -> Result<ReprojectExactOutput, ReprojectError> {
    let corner_shape = pixel_corners.shape();
    let rows = corner_shape[0];
    let cols = corner_shape[1];
    let last = corner_shape[2];
    if last != 2 || rows < 2 || cols < 2 {
        return Err(ReprojectError::PixelCornersShape { rows, cols, last });
    }
    if let Some(error_view) = error.as_ref() {
        let image_shape = image.shape();
        let error_shape = error_view.shape();
        if image_shape != error_shape {
            return Err(ReprojectError::ErrorShapeMismatch {
                image_rows: image_shape[0],
                image_cols: image_shape[1],
                error_rows: error_shape[0],
                error_cols: error_shape[1],
            });
        }
    }
    let image_f64: Array2<f64> = image.mapv(|value| value.to_f64().unwrap_or(f64::NAN));
    let error_f64: Option<Array2<f64>> =
        error.map(|view| view.mapv(|value| value.to_f64().unwrap_or(f64::NAN)));
    let height_out = rows - 1;
    let width_out = cols - 1;
    let mut image_out = Array2::<f64>::from_elem((height_out, width_out), f64::NAN);
    let mut footprint = Array2::<f64>::zeros((height_out, width_out));
    let mut weight = Array2::<f64>::zeros((height_out, width_out));
    let mut error_out = error_f64
        .as_ref()
        .map(|_| Array2::<f64>::from_elem((height_out, width_out), f64::NAN));
    let error_f64_view = error_f64.as_ref().map(|array| array.view());
    for row_index in 0..height_out {
        let row_error = error_out
            .as_mut()
            .map(|array| array.index_axis_mut(Axis(0), row_index));
        process_output_row(
            row_index,
            image_out.index_axis_mut(Axis(0), row_index),
            footprint.index_axis_mut(Axis(0), row_index),
            weight.index_axis_mut(Axis(0), row_index),
            row_error,
            image_f64.view(),
            error_f64_view,
            pixel_corners,
        );
    }
    Ok(ReprojectExactOutput {
        image: image_out,
        footprint,
        weight,
        error: error_out,
    })
}

/// Compute one output row in place. Factored out of the driver so a
/// parallel driver can reuse the exact same per-row work unit.
///
/// `row_error` and `error_f64` are `Some` together (error requested)
/// or `None` together (no error). When present, the per-output-pixel
/// marginal 1-sigma is written into `row_error` using first-order
/// independent propagation; see [`ReprojectExactOutput::error`].
#[allow(clippy::too_many_arguments)]
pub(crate) fn process_output_row(
    row_index: usize,
    mut row_image: ArrayViewMut1<f64>,
    mut row_footprint: ArrayViewMut1<f64>,
    mut row_weight: ArrayViewMut1<f64>,
    mut row_error: Option<ArrayViewMut1<f64>>,
    image_f64: ArrayView2<f64>,
    error_f64: Option<ArrayView2<f64>>,
    pixel_corners: ArrayView3<f64>,
) {
    let height_in = image_f64.shape()[0];
    let width_in = image_f64.shape()[1];
    let width_out = row_image.len();

    for column_index in 0..width_out {
        // Read the four corners that bound output pixel
        // (row_index, column_index). Winding (CCW or CW) is irrelevant
        // because we take the absolute area; the only requirement is
        // that consecutive indices walk around the quad.
        let corner_top_left = corner_at(pixel_corners, row_index, column_index);
        let corner_top_right = corner_at(pixel_corners, row_index, column_index + 1);
        let corner_bottom_right = corner_at(pixel_corners, row_index + 1, column_index + 1);
        let corner_bottom_left = corner_at(pixel_corners, row_index + 1, column_index);

        let any_nan_corner = corner_top_left[0].is_nan()
            || corner_top_left[1].is_nan()
            || corner_top_right[0].is_nan()
            || corner_top_right[1].is_nan()
            || corner_bottom_right[0].is_nan()
            || corner_bottom_right[1].is_nan()
            || corner_bottom_left[0].is_nan()
            || corner_bottom_left[1].is_nan();
        if any_nan_corner {
            row_image[column_index] = f64::NAN;
            row_footprint[column_index] = 0.0;
            row_weight[column_index] = 0.0;
            continue;
        }

        // Apply the half-pixel shift so input pixel (i, j) occupies the
        // unit cell [j, j+1] x [i, i+1]. Keeping this internal lets
        // callers reason in astropy's convention (integer = center).
        let quad: [Point; 4] = [
            [corner_top_left[0] + 0.5, corner_top_left[1] + 0.5],
            [corner_top_right[0] + 0.5, corner_top_right[1] + 0.5],
            [corner_bottom_right[0] + 0.5, corner_bottom_right[1] + 0.5],
            [corner_bottom_left[0] + 0.5, corner_bottom_left[1] + 0.5],
        ];

        let output_pixel_area = signed_area(&quad).abs();

        // Bounding box of the (shifted) quad in continuous coords.
        let mut x_min_f = quad[0][0];
        let mut x_max_f = quad[0][0];
        let mut y_min_f = quad[0][1];
        let mut y_max_f = quad[0][1];
        for vertex in quad.iter().skip(1) {
            if vertex[0] < x_min_f {
                x_min_f = vertex[0];
            }
            if vertex[0] > x_max_f {
                x_max_f = vertex[0];
            }
            if vertex[1] < y_min_f {
                y_min_f = vertex[1];
            }
            if vertex[1] > y_max_f {
                y_max_f = vertex[1];
            }
        }

        // Clamp to image bounds. Right-open: pixel (W_in - 1) covers
        // [W_in - 1, W_in] after the shift, so the maximum integer cell
        // index we enumerate is W_in - 1 (inclusive).
        if x_max_f <= 0.0
            || y_max_f <= 0.0
            || x_min_f >= width_in as f64
            || y_min_f >= height_in as f64
        {
            row_image[column_index] = f64::NAN;
            row_footprint[column_index] = 0.0;
            row_weight[column_index] = 0.0;
            continue;
        }

        let column_lo = x_min_f.floor().max(0.0) as i64;
        let column_hi = (x_max_f.ceil() as i64).min(width_in as i64);
        let row_lo = y_min_f.floor().max(0.0) as i64;
        let row_hi = (y_max_f.ceil() as i64).min(height_in as i64);

        let mut numerator = 0.0_f64;
        let mut denominator_valid = 0.0_f64;
        let mut denominator_geom = 0.0_f64;
        // Error propagation accumulators, used only when error_f64 is
        // Some. `variance_numerator` collects sum(A_k^2 * sigma_k^2);
        // `error_poisoned` records a valid pixel with an unusable error.
        let mut variance_numerator = 0.0_f64;
        let mut error_poisoned = false;
        for cell_row in row_lo..row_hi {
            for cell_column in column_lo..column_hi {
                // Always compute the geometric overlap first: the
                // footprint accumulator must count all in-bounds area
                // regardless of input data validity.
                let clipped =
                    clip_quad_against_unit_cell(&quad, (cell_column as i32, cell_row as i32));
                if clipped.is_empty() {
                    continue;
                }
                let overlap_area = signed_area(&clipped).abs();
                if overlap_area == 0.0 {
                    continue;
                }
                denominator_geom += overlap_area;
                let pixel_value = image_f64[(cell_row as usize, cell_column as usize)];
                if pixel_value.is_nan() {
                    // NaN-as-mask: drop this input pixel from the
                    // numerator and the valid denominator (weight),
                    // but it has already been counted toward footprint.
                    continue;
                }
                numerator += overlap_area * pixel_value;
                denominator_valid += overlap_area;
                if let Some(error) = error_f64.as_ref() {
                    let sigma = error[(cell_row as usize, cell_column as usize)];
                    if sigma.is_finite() && sigma >= 0.0 {
                        variance_numerator += overlap_area * overlap_area * sigma * sigma;
                    } else {
                        // The value is valid (it has entered the mean) but
                        // its error is unusable. Dropping the term would
                        // bias the output 1-sigma low, so poison this
                        // output pixel's error instead.
                        error_poisoned = true;
                    }
                }
            }
        }

        if denominator_valid > 0.0 {
            row_image[column_index] = numerator / denominator_valid;
        } else {
            row_image[column_index] = f64::NAN;
        }
        if output_pixel_area > 0.0 {
            row_footprint[column_index] = (denominator_geom / output_pixel_area).clamp(0.0, 1.0);
            row_weight[column_index] = (denominator_valid / output_pixel_area).clamp(0.0, 1.0);
        } else {
            row_footprint[column_index] = 0.0;
            row_weight[column_index] = 0.0;
        }
        if let Some(row_error) = row_error.as_mut() {
            // NaN when there is no valid contribution (matching image) or
            // when a contributing pixel had an unusable error; otherwise
            // the propagated marginal 1-sigma. NaN-corner and
            // out-of-bounds output pixels took an early `continue` above
            // and keep the array's initial NaN.
            row_error[column_index] = if denominator_valid > 0.0 && !error_poisoned {
                variance_numerator.sqrt() / denominator_valid
            } else {
                f64::NAN
            };
        }
    }
}

/// Read the `(x, y)` corner at grid node `(row, column)` out of the
/// corner field.
#[inline]
fn corner_at(pixel_corners: ArrayView3<f64>, row: usize, column: usize) -> [f64; 2] {
    [
        pixel_corners[(row, column, 0)],
        pixel_corners[(row, column, 1)],
    ]
}

#[cfg(test)]
mod tests {
    use super::*;
    use ndarray::{Array2, Array3};

    const TOL: f64 = 1e-12;

    fn approx_eq(a: f64, b: f64, tol: f64) -> bool {
        (a - b).abs() <= tol * a.abs().max(b.abs()).max(1.0)
    }

    /// Build an identity corner field: output grid coincides with input
    /// pixel boundaries. Output pixel (i, j) is bounded by input-pixel
    /// corners at integer + 0.5 spacing (centers convention), but the
    /// caller passes them as the boundary nodes, i.e. integer-aligned
    /// half-pixel-offset coordinates.
    fn identity_corners(height_out: usize, width_out: usize) -> Array3<f64> {
        // Output pixel (i, j) should sample input pixel (i, j). In
        // continuous coords, input pixel (i, j) has center (j, i) and
        // edges at j - 0.5 and i - 0.5. So corner (i_node, j_node) is
        // at (j_node - 0.5, i_node - 0.5).
        let mut corners = Array3::<f64>::zeros((height_out + 1, width_out + 1, 2));
        for i_node in 0..=height_out {
            for j_node in 0..=width_out {
                corners[(i_node, j_node, 0)] = j_node as f64 - 0.5;
                corners[(i_node, j_node, 1)] = i_node as f64 - 0.5;
            }
        }
        corners
    }

    #[test]
    fn identity_reprojection_preserves_image_footprint_and_weight_f64() {
        let image: Array2<f64> =
            ndarray::array![[1.0, 2.0, 3.0], [4.0, 5.0, 6.0], [7.0, 8.0, 9.0],];
        let corners = identity_corners(3, 3);
        let output = reproject_exact(image.view(), corners.view(), None).unwrap();
        for i in 0..3 {
            for j in 0..3 {
                assert!(approx_eq(output.image[(i, j)], image[(i, j)], TOL));
                assert!(approx_eq(output.footprint[(i, j)], 1.0, TOL));
                assert!(approx_eq(output.weight[(i, j)], 1.0, TOL));
            }
        }
    }

    #[test]
    fn identity_reprojection_preserves_image_footprint_and_weight_f32() {
        let image: Array2<f32> =
            ndarray::array![[1.0, 2.0, 3.0], [4.0, 5.0, 6.0], [7.0, 8.0, 9.0],];
        let corners = identity_corners(3, 3);
        let output = reproject_exact(image.view(), corners.view(), None).unwrap();
        for i in 0..3 {
            for j in 0..3 {
                assert!(approx_eq(output.image[(i, j)], image[(i, j)] as f64, TOL));
                assert!(approx_eq(output.footprint[(i, j)], 1.0, TOL));
                assert!(approx_eq(output.weight[(i, j)], 1.0, TOL));
            }
        }
    }

    #[test]
    fn constant_image_yields_constant_output_where_weight_positive() {
        let constant_value = 4.2_f64;
        let image: Array2<f64> = Array2::from_elem((5, 5), constant_value);
        // Translate output grid by (0.3, -0.2) in input-pixel units;
        // shape stays the same.
        let mut corners = identity_corners(5, 5);
        for i_node in 0..corners.shape()[0] {
            for j_node in 0..corners.shape()[1] {
                corners[(i_node, j_node, 0)] += 0.3;
                corners[(i_node, j_node, 1)] -= 0.2;
            }
        }
        let output = reproject_exact(image.view(), corners.view(), None).unwrap();
        for i in 0..5 {
            for j in 0..5 {
                // No NaN inputs anywhere, so footprint and weight must
                // be identical wherever there is any coverage.
                assert!(
                    approx_eq(output.footprint[(i, j)], output.weight[(i, j)], TOL),
                    "i={i} j={j} footprint {} != weight {}",
                    output.footprint[(i, j)],
                    output.weight[(i, j)]
                );
                if output.weight[(i, j)] > 0.0 {
                    assert!(
                        approx_eq(output.image[(i, j)], constant_value, TOL),
                        "i={i} j={j} got {} weight {}",
                        output.image[(i, j)],
                        output.weight[(i, j)]
                    );
                }
            }
        }
    }

    #[test]
    fn half_pixel_shift_linear_average_of_two_columns() {
        // Image: 1x4. Shift output by +0.5 in x. Each output pixel
        // should average two adjacent input pixels (1:1 mix).
        let image: Array2<f64> = ndarray::array![[1.0, 2.0, 3.0, 4.0]];
        // Identity corners for 1x4 output.
        let mut corners = identity_corners(1, 4);
        for i_node in 0..corners.shape()[0] {
            for j_node in 0..corners.shape()[1] {
                corners[(i_node, j_node, 0)] += 0.5;
            }
        }
        let output = reproject_exact(image.view(), corners.view(), None).unwrap();
        // Output pixel j has x-extent (j - 0.5 + 0.5, j + 0.5 + 0.5)
        //                              = (j, j + 1).
        // After internal +0.5 shift -> (j + 0.5, j + 1.5).
        // Input pixel j occupies cell [j, j + 1] (after shift).
        // So output pixel 0 overlaps input pixels 0 and 1, each by 0.5.
        // => image_out[0, 0] = (image[0]*0.5 + image[1]*0.5) / 1.0
        //   = 1.5; weight = 1.0.
        assert!(approx_eq(output.image[(0, 0)], 1.5, TOL));
        assert!(approx_eq(output.image[(0, 1)], 2.5, TOL));
        assert!(approx_eq(output.image[(0, 2)], 3.5, TOL));
        // Output pixel 3 overlaps input pixel 3 (half) only; pixel 4
        // does not exist, so footprint = weight = 0.5 and image_out
        // = 4.0. No NaN inputs anywhere, so footprint and weight match.
        assert!(approx_eq(output.image[(0, 3)], 4.0, TOL));
        assert!(approx_eq(output.footprint[(0, 0)], 1.0, TOL));
        assert!(approx_eq(output.footprint[(0, 3)], 0.5, TOL));
        assert!(approx_eq(output.weight[(0, 0)], 1.0, TOL));
        assert!(approx_eq(output.weight[(0, 3)], 0.5, TOL));
    }

    #[test]
    fn rotation_by_90_degrees_preserves_surface_brightness() {
        // Rotate output grid by 90 degrees about the (1, 1) input
        // center. Use a constant image -> output should be the same
        // constant everywhere it is fully covered.
        let constant_value = 7.0_f64;
        let image: Array2<f64> = Array2::from_elem((3, 3), constant_value);
        // For output pixel node (i_node, j_node), instead of mapping
        // 1:1 we rotate -90 degrees: (x_in, y_in) =
        //   (y_node_centered + center_x, -x_node_centered + center_y).
        // With center = (1, 1).
        let mut corners = Array3::<f64>::zeros((4, 4, 2));
        for i_node in 0..4 {
            for j_node in 0..4 {
                let x_orig = j_node as f64 - 0.5;
                let y_orig = i_node as f64 - 0.5;
                let x_centered = x_orig - 1.0;
                let y_centered = y_orig - 1.0;
                let x_rotated = y_centered + 1.0;
                let y_rotated = -x_centered + 1.0;
                corners[(i_node, j_node, 0)] = x_rotated;
                corners[(i_node, j_node, 1)] = y_rotated;
            }
        }
        let output = reproject_exact(image.view(), corners.view(), None).unwrap();
        for i in 0..3 {
            for j in 0..3 {
                assert!(approx_eq(output.image[(i, j)], constant_value, TOL));
                assert!(approx_eq(output.footprint[(i, j)], 1.0, TOL));
                assert!(approx_eq(output.weight[(i, j)], 1.0, TOL));
            }
        }
    }

    #[test]
    fn output_pixels_outside_footprint_are_nan_zero() {
        // 2x2 input, output corner field places the whole output
        // grid far away from the input image.
        let image: Array2<f64> = ndarray::array![[1.0, 2.0], [3.0, 4.0]];
        let mut corners = Array3::<f64>::zeros((3, 3, 2));
        for i_node in 0..3 {
            for j_node in 0..3 {
                corners[(i_node, j_node, 0)] = 100.0 + j_node as f64;
                corners[(i_node, j_node, 1)] = 100.0 + i_node as f64;
            }
        }
        let output = reproject_exact(image.view(), corners.view(), None).unwrap();
        for i in 0..2 {
            for j in 0..2 {
                assert!(output.image[(i, j)].is_nan());
                assert_eq!(output.footprint[(i, j)], 0.0);
                assert_eq!(output.weight[(i, j)], 0.0);
            }
        }
    }

    #[test]
    fn nan_input_separates_footprint_from_weight() {
        // 1x3 image with the middle pixel NaN. Output pixel that spans
        // the middle two pixels in equal proportion should reduce to
        // the non-NaN pixel's value, with weight = 0.5 (only half of
        // the area contributed valid data). Footprint is 1.0 because
        // the output pixel is fully covered by the input grid
        // geometrically — NaN does not reduce footprint.
        let image: Array2<f64> = ndarray::array![[1.0, f64::NAN, 3.0]];
        // Identity 1x3 corners, then shift +0.5 in x: output pixel j
        // spans input pixels j and j+1, 50/50.
        let mut corners = identity_corners(1, 3);
        for i_node in 0..corners.shape()[0] {
            for j_node in 0..corners.shape()[1] {
                corners[(i_node, j_node, 0)] += 0.5;
            }
        }
        let output = reproject_exact(image.view(), corners.view(), None).unwrap();
        // Output pixel 0 spans input 0 (1.0) + input 1 (NaN), 50/50.
        // image = 1.0; footprint = 1.0 (fully geometrically covered);
        // weight = 0.5 (only the non-NaN half contributed).
        assert!(approx_eq(output.image[(0, 0)], 1.0, TOL));
        assert!(approx_eq(output.footprint[(0, 0)], 1.0, TOL));
        assert!(approx_eq(output.weight[(0, 0)], 0.5, TOL));
        // Output pixel 1 spans input 1 (NaN) + input 2 (3.0). Same
        // pattern.
        assert!(approx_eq(output.image[(0, 1)], 3.0, TOL));
        assert!(approx_eq(output.footprint[(0, 1)], 1.0, TOL));
        assert!(approx_eq(output.weight[(0, 1)], 0.5, TOL));
        // Output pixel 2 spans input 2 (3.0) + (off-grid pixel 3, NaN
        // because out of bounds). image = 3.0; footprint = 0.5 (only
        // half the output pixel is inside the input image);
        // weight = 0.5 (the half that is inside is non-NaN).
        assert!(approx_eq(output.image[(0, 2)], 3.0, TOL));
        assert!(approx_eq(output.footprint[(0, 2)], 0.5, TOL));
        assert!(approx_eq(output.weight[(0, 2)], 0.5, TOL));
    }

    #[test]
    fn nan_in_corners_marks_output_nan_zero() {
        let image: Array2<f64> = ndarray::array![[1.0, 2.0], [3.0, 4.0]];
        let mut corners = identity_corners(2, 2);
        corners[(1, 1, 0)] = f64::NAN;
        let output = reproject_exact(image.view(), corners.view(), None).unwrap();
        // The NaN corner is shared by output pixels (0, 0), (0, 1),
        // (1, 0), (1, 1) — all 4. Each must be NaN with both footprint
        // and weight zeroed.
        for i in 0..2 {
            for j in 0..2 {
                assert!(output.image[(i, j)].is_nan());
                assert_eq!(output.footprint[(i, j)], 0.0);
                assert_eq!(output.weight[(i, j)], 0.0);
            }
        }
    }

    #[test]
    fn shape_error_last_dim_not_two() {
        let image: Array2<f64> = Array2::zeros((2, 2));
        let corners = Array3::<f64>::zeros((3, 3, 3));
        let err = reproject_exact(image.view(), corners.view(), None).unwrap_err();
        assert_eq!(
            err,
            ReprojectError::PixelCornersShape {
                rows: 3,
                cols: 3,
                last: 3
            }
        );
    }

    #[test]
    fn parallel_matches_sequential_bitwise_on_moderate_input() {
        // Moderate size so rayon actually splits the work, and a
        // sub-pixel-shifted corner field so every output pixel does
        // real area-weighted averaging.
        let height = 32usize;
        let width = 40usize;
        let mut image: Array2<f64> = Array2::zeros((height, width));
        for i in 0..height {
            for j in 0..width {
                image[(i, j)] = ((i * 13 + j * 7) % 97) as f64 / 11.0 + (i as f64).sin();
            }
        }
        let mut corners = Array3::<f64>::zeros((height + 1, width + 1, 2));
        for i_node in 0..=height {
            for j_node in 0..=width {
                corners[(i_node, j_node, 0)] = j_node as f64 - 0.5 + 0.37;
                corners[(i_node, j_node, 1)] = i_node as f64 - 0.5 - 0.21;
            }
        }
        let parallel = reproject_exact(image.view(), corners.view(), None).unwrap();
        let sequential = reproject_exact_sequential(image.view(), corners.view(), None).unwrap();
        assert_eq!(parallel.image.shape(), sequential.image.shape());
        for i in 0..height {
            for j in 0..width {
                let p_image = parallel.image[(i, j)];
                let s_image = sequential.image[(i, j)];
                if p_image.is_nan() {
                    assert!(s_image.is_nan());
                } else {
                    // Bitwise equality: each output pixel's per-input
                    // accumulation order is fixed (row-major over the
                    // bbox) and lives entirely inside one row task, so
                    // floats accumulate identically regardless of how
                    // rayon schedules rows.
                    assert_eq!(p_image.to_bits(), s_image.to_bits());
                }
                assert_eq!(
                    parallel.footprint[(i, j)].to_bits(),
                    sequential.footprint[(i, j)].to_bits()
                );
                assert_eq!(
                    parallel.weight[(i, j)].to_bits(),
                    sequential.weight[(i, j)].to_bits()
                );
            }
        }
    }

    #[test]
    fn shape_error_front_dims_too_small() {
        let image: Array2<f64> = Array2::zeros((2, 2));
        let corners = Array3::<f64>::zeros((1, 3, 2));
        let err = reproject_exact(image.view(), corners.view(), None).unwrap_err();
        assert_eq!(
            err,
            ReprojectError::PixelCornersShape {
                rows: 1,
                cols: 3,
                last: 2
            }
        );
    }

    // ------------------------------------------------------------------
    // Optional error propagation
    // ------------------------------------------------------------------

    #[test]
    fn error_is_none_when_not_requested() {
        let image: Array2<f64> = ndarray::array![[1.0, 2.0], [3.0, 4.0]];
        let corners = identity_corners(2, 2);
        let output = reproject_exact(image.view(), corners.view(), None).unwrap();
        assert!(output.error.is_none());
    }

    #[test]
    fn identity_error_passes_through() {
        // Identity reprojection: output pixel (i, j) is exactly input
        // pixel (i, j), so A = 1, denominator = 1, and the propagated
        // 1-sigma equals the input 1-sigma unchanged.
        let image: Array2<f64> = ndarray::array![[1.0, 2.0, 3.0], [4.0, 5.0, 6.0], [7.0, 8.0, 9.0]];
        let error: Array2<f64> = ndarray::array![[0.1, 0.2, 0.3], [0.4, 0.5, 0.6], [0.7, 0.8, 0.9]];
        let corners = identity_corners(3, 3);
        let output = reproject_exact(image.view(), corners.view(), Some(error.view())).unwrap();
        let output_error = output.error.expect("error requested");
        for i in 0..3 {
            for j in 0..3 {
                assert!(approx_eq(output_error[(i, j)], error[(i, j)], TOL));
            }
        }
    }

    #[test]
    fn half_pixel_shift_averages_error_in_quadrature() {
        // Each output pixel mixes two adjacent input pixels 50/50. With
        // equal input sigma s, sigma_out = sqrt(0.5^2 s^2 + 0.5^2 s^2)
        // = s / sqrt(2) -- the classic sigma/sqrt(N) for N = 2 equal
        // contributors. The last output pixel sees only one input pixel
        // (half area), so denom = 0.5 and sigma_out = s.
        let image: Array2<f64> = ndarray::array![[1.0, 2.0, 3.0, 4.0]];
        let sigma = 0.3_f64;
        let error: Array2<f64> = Array2::from_elem((1, 4), sigma);
        let mut corners = identity_corners(1, 4);
        for i_node in 0..corners.shape()[0] {
            for j_node in 0..corners.shape()[1] {
                corners[(i_node, j_node, 0)] += 0.5;
            }
        }
        let output = reproject_exact(image.view(), corners.view(), Some(error.view())).unwrap();
        let output_error = output.error.expect("error requested");
        let expected = sigma / 2.0_f64.sqrt();
        assert!(approx_eq(output_error[(0, 0)], expected, TOL));
        assert!(approx_eq(output_error[(0, 1)], expected, TOL));
        assert!(approx_eq(output_error[(0, 2)], expected, TOL));
        // Output pixel 3: single input pixel, A = 0.5, denom = 0.5 ->
        // sigma_out = sqrt(0.25 s^2) / 0.5 = s.
        assert!(approx_eq(output_error[(0, 3)], sigma, TOL));
    }

    #[test]
    fn nonfinite_sigma_poisons_error_but_not_image() {
        // A valid value paired with a NaN sigma: the value still enters
        // the mean, so the output image is untouched, but the output
        // error is poisoned to NaN (dropping the unknown variance would
        // bias it low).
        let image: Array2<f64> = ndarray::array![[1.0, 2.0, 3.0, 4.0]];
        let mut error: Array2<f64> = Array2::from_elem((1, 4), 0.3_f64);
        error[(0, 1)] = f64::NAN;
        let mut corners = identity_corners(1, 4);
        for i_node in 0..corners.shape()[0] {
            for j_node in 0..corners.shape()[1] {
                corners[(i_node, j_node, 0)] += 0.5;
            }
        }
        let without_error = reproject_exact(image.view(), corners.view(), None).unwrap();
        let with_error = reproject_exact(image.view(), corners.view(), Some(error.view())).unwrap();
        // Image is bitwise identical with and without the error image.
        for j in 0..4 {
            assert_eq!(
                without_error.image[(0, j)].to_bits(),
                with_error.image[(0, j)].to_bits()
            );
        }
        let output_error = with_error.error.expect("error requested");
        // Output 0 mixes input 0 (sigma 0.3) and input 1 (sigma NaN):
        // poisoned. Output 1 mixes input 1 (NaN) and input 2: poisoned.
        assert!(output_error[(0, 0)].is_nan());
        assert!(output_error[(0, 1)].is_nan());
        // Output 2 mixes input 2 and 3, both finite sigma: clean.
        assert!(output_error[(0, 2)].is_finite());
        // Output 3 sees only input 3 (finite sigma): clean.
        assert!(output_error[(0, 3)].is_finite());
    }

    #[test]
    fn negative_sigma_poisons_error() {
        // A negative sigma is just as unusable as a non-finite one.
        let image: Array2<f64> = ndarray::array![[1.0, 2.0]];
        let error: Array2<f64> = ndarray::array![[-1.0, 0.2]];
        let corners = identity_corners(1, 2);
        let output = reproject_exact(image.view(), corners.view(), Some(error.view())).unwrap();
        let output_error = output.error.expect("error requested");
        // Output 0 is exactly input 0, whose sigma is negative -> NaN.
        assert!(output_error[(0, 0)].is_nan());
        // Output 1 is exactly input 1, sigma 0.2 -> passes through.
        assert!(approx_eq(output_error[(0, 1)], 0.2, TOL));
    }

    #[test]
    fn value_nan_input_never_consults_its_sigma() {
        // The middle input value is NaN, so it is dropped before its
        // sigma is read. Even a NaN sigma there must not poison the
        // output: the output error comes only from the valid neighbour.
        let image: Array2<f64> = ndarray::array![[1.0, f64::NAN, 3.0]];
        let error: Array2<f64> = ndarray::array![[2.0, f64::NAN, 4.0]];
        let mut corners = identity_corners(1, 3);
        for i_node in 0..corners.shape()[0] {
            for j_node in 0..corners.shape()[1] {
                corners[(i_node, j_node, 0)] += 0.5;
            }
        }
        let output = reproject_exact(image.view(), corners.view(), Some(error.view())).unwrap();
        let output_error = output.error.expect("error requested");
        // Output 0: input 0 (value 1.0, A = 0.5, denom = 0.5) + input 1
        // (value NaN, dropped). sigma_out = sqrt(0.25 * 2^2) / 0.5 = 2.0.
        assert!(approx_eq(output_error[(0, 0)], 2.0, TOL));
        // Output 1: input 1 (NaN, dropped) + input 2 (value 3.0).
        // sigma_out = sqrt(0.25 * 4^2) / 0.5 = 4.0.
        assert!(approx_eq(output_error[(0, 1)], 4.0, TOL));
    }

    #[test]
    fn nan_corner_and_out_of_bounds_give_nan_error() {
        // NaN corner: the touched output pixels are NaN/0; their error
        // must be NaN too. Out-of-bounds output pixels (no valid
        // contribution) likewise carry NaN error.
        let image: Array2<f64> = ndarray::array![[1.0, 2.0], [3.0, 4.0]];
        let error: Array2<f64> = Array2::from_elem((2, 2), 0.5_f64);
        let mut corners = identity_corners(2, 2);
        corners[(1, 1, 0)] = f64::NAN;
        let output = reproject_exact(image.view(), corners.view(), Some(error.view())).unwrap();
        let output_error = output.error.expect("error requested");
        // The NaN corner (1, 1) is shared by all four output pixels.
        for i in 0..2 {
            for j in 0..2 {
                assert!(output_error[(i, j)].is_nan());
            }
        }

        // Whole output grid placed far off the input image -> no valid
        // contribution anywhere -> all-NaN error.
        let mut far_corners = Array3::<f64>::zeros((3, 3, 2));
        for i_node in 0..3 {
            for j_node in 0..3 {
                far_corners[(i_node, j_node, 0)] = 100.0 + j_node as f64;
                far_corners[(i_node, j_node, 1)] = 100.0 + i_node as f64;
            }
        }
        let far_error: Array2<f64> = Array2::from_elem((2, 2), 0.5_f64);
        let far =
            reproject_exact(image.view(), far_corners.view(), Some(far_error.view())).unwrap();
        let far_output_error = far.error.expect("error requested");
        for i in 0..2 {
            for j in 0..2 {
                assert!(far_output_error[(i, j)].is_nan());
            }
        }
    }

    #[test]
    fn error_shape_mismatch_is_rejected() {
        let image: Array2<f64> = Array2::zeros((3, 4));
        let error: Array2<f64> = Array2::zeros((3, 5));
        let corners = identity_corners(3, 4);
        let err = reproject_exact(image.view(), corners.view(), Some(error.view())).unwrap_err();
        assert_eq!(
            err,
            ReprojectError::ErrorShapeMismatch {
                image_rows: 3,
                image_cols: 4,
                error_rows: 3,
                error_cols: 5,
            }
        );
    }

    #[test]
    fn error_f32_and_f64_inputs_agree() {
        // The error path upcasts f32 -> f64 at the entry just like the
        // image path; identity reprojection makes the propagated error
        // equal the (upcast) input error.
        let image_f64: Array2<f64> = ndarray::array![[1.0, 2.0], [3.0, 4.0]];
        let error_f64: Array2<f64> = ndarray::array![[0.1, 0.2], [0.3, 0.4]];
        let image_f32: Array2<f32> = image_f64.mapv(|value| value as f32);
        let error_f32: Array2<f32> = error_f64.mapv(|value| value as f32);
        let corners = identity_corners(2, 2);
        let output_f64 =
            reproject_exact(image_f64.view(), corners.view(), Some(error_f64.view())).unwrap();
        let output_f32 =
            reproject_exact(image_f32.view(), corners.view(), Some(error_f32.view())).unwrap();
        let error_out_f64 = output_f64.error.expect("error requested");
        let error_out_f32 = output_f32.error.expect("error requested");
        for i in 0..2 {
            for j in 0..2 {
                assert!(approx_eq(
                    error_out_f64[(i, j)],
                    error_out_f32[(i, j)],
                    1e-6
                ));
            }
        }
    }

    #[test]
    fn parallel_matches_sequential_bitwise_with_error() {
        // Same moderate setup as the no-error bitwise test, now driving
        // the error path: per-output-pixel accumulation order is fixed
        // and lives inside one row task, so parallel and sequential must
        // agree bit-for-bit on the error array too.
        let height = 32usize;
        let width = 40usize;
        let mut image: Array2<f64> = Array2::zeros((height, width));
        let mut error: Array2<f64> = Array2::zeros((height, width));
        for i in 0..height {
            for j in 0..width {
                image[(i, j)] = ((i * 13 + j * 7) % 97) as f64 / 11.0 + (i as f64).sin();
                error[(i, j)] = ((i * 5 + j * 3) % 13) as f64 / 7.0 + 0.05;
            }
        }
        let mut corners = Array3::<f64>::zeros((height + 1, width + 1, 2));
        for i_node in 0..=height {
            for j_node in 0..=width {
                corners[(i_node, j_node, 0)] = j_node as f64 - 0.5 + 0.37;
                corners[(i_node, j_node, 1)] = i_node as f64 - 0.5 - 0.21;
            }
        }
        let parallel = reproject_exact(image.view(), corners.view(), Some(error.view())).unwrap();
        let sequential =
            reproject_exact_sequential(image.view(), corners.view(), Some(error.view())).unwrap();
        let parallel_error = parallel.error.expect("error requested");
        let sequential_error = sequential.error.expect("error requested");
        for i in 0..height {
            for j in 0..width {
                let p = parallel_error[(i, j)];
                let s = sequential_error[(i, j)];
                if p.is_nan() {
                    assert!(s.is_nan());
                } else {
                    assert_eq!(p.to_bits(), s.to_bits());
                }
            }
        }
    }
}