ridge-core 0.3.2

//! `scipy.ndimage.rotate` — frozen port.
//!
//! Replicates `scipy.ndimage.rotate(input, angle, axes=(1, 0), reshape,
//! order, mode='constant', cval=0.0)` for spline orders 0 (nearest) and 1
//! (bilinear), which covers everything ridge_map uses (`interpolation` in
//! 0..=1; higher orders upstream tend to produce all-NaN maps anyway).
//!
//! Algorithm (from `scipy/ndimage/_interpolation.py::rotate`):
//!
//! ```text
//! M = [[cos, sin], [-sin, cos]]      # cosdg/sindg: exact at 90-degree steps
//! in_plane  = (nrows, ncols)
//! if reshape:
//!     bounds  = M @ [(0,0), (0,ncols), (nrows,0), (nrows,ncols)]
//!     out_plane = floor(ptp(bounds) + 0.5)
//! else:
//!     out_plane = in_plane
//! out_center = M @ ((out_plane - 1) / 2)
//! in_center  = (in_plane - 1) / 2
//! offset     = in_center - out_center
//! output[o]  = interpolate(input, M @ o + offset)   # cval outside
//! ```

use ndarray::Array2;

/// `scipy.special.cosdg/sindg` — exact values at multiples of 90 degrees.
pub(crate) fn sincosdg(angle_deg: f64) -> (f64, f64) {
    let a = angle_deg.rem_euclid(360.0);
    match a {
        0.0 => (1.0, 0.0),
        90.0 => (0.0, 1.0),
        180.0 => (-1.0, 0.0),
        270.0 => (0.0, -1.0),
        v => {
            let rad = v * std::f64::consts::PI / 180.0;
            (rad.cos(), rad.sin())
        }
    }
}

fn output_shape(shape: (usize, usize), angle_deg: f64, reshape: bool) -> (usize, usize) {
    let (c, s) = sincosdg(angle_deg);
    let (nrows, ncols) = shape;
    if !reshape {
        return shape;
    }
    // M @ the four corners (0,0), (0,ncols), (nrows,0), (nrows,ncols).
    let mut rs = [0.0f64; 4];
    let mut cs = [0.0f64; 4];
    for (k, (r, col)) in [
        (0.0, 0.0),
        (0.0, ncols as f64),
        (nrows as f64, 0.0),
        (nrows as f64, ncols as f64),
    ]
    .into_iter()
    .enumerate()
    {
        rs[k] = c * r + s * col;
        cs[k] = -s * r + c * col;
    }
    let ptp_r =
        rs.iter().cloned().fold(f64::MIN, f64::max) - rs.iter().cloned().fold(f64::MAX, f64::min);
    let ptp_c =
        cs.iter().cloned().fold(f64::MIN, f64::max) - cs.iter().cloned().fold(f64::MAX, f64::min);
    // (ptp + 0.5).astype(int): numpy casts truncate toward zero; values are positive.
    ((ptp_r + 0.5) as usize, (ptp_c + 0.5) as usize)
}

/// Rotate `input` by `angle_deg` degrees, replicating scipy's `mode='constant'`,
/// `cval=0.0` behavior (out-of-bounds samples become exactly 0.0).
pub fn rotate(input: &Array2<f64>, angle_deg: f64, reshape: bool, order: u32) -> Array2<f64> {
    assert!(order <= 1, "only spline orders 0 and 1 are supported");
    let (nrows, ncols) = (input.nrows(), input.ncols());
    let (c, s) = sincosdg(angle_deg);
    let (out_rows, out_cols) = output_shape((nrows, ncols), angle_deg, reshape);

    // Centers are at plane midpoints ((n - 1) / 2 each).
    let in_center = ((nrows as f64 - 1.0) / 2.0, (ncols as f64 - 1.0) / 2.0);
    let out_center_r = c * (out_rows as f64 - 1.0) / 2.0 + s * (out_cols as f64 - 1.0) / 2.0;
    let out_center_c = -s * (out_rows as f64 - 1.0) / 2.0 + c * (out_cols as f64 - 1.0) / 2.0;
    let offset = (in_center.0 - out_center_r, in_center.1 - out_center_c);

    let mut out = Array2::<f64>::zeros((out_rows, out_cols));
    for or in 0..out_rows {
        for oc in 0..out_cols {
            let pr = c * or as f64 + s * oc as f64 + offset.0;
            let pc = -s * or as f64 + c * oc as f64 + offset.1;
            out[(or, oc)] = match order {
                0 => sample_nearest(input, pr, pc, 0.0),
                _ => sample_bilinear(input, pr, pc, 0.0),
            };
        }
    }
    out
}

pub(crate) fn sample_nearest(input: &Array2<f64>, pr: f64, pc: f64, fill: f64) -> f64 {
    // scipy mode='constant': the *coordinate* must lie within [0, n-1]; anything
    // outside is `cval`, even when it would round to an in-range cell (e.g.
    // -0.1 or 3.1 for a length-4 axis). Checking the rounded index instead
    // samples the edge cell where scipy returns cval.
    let max_r = input.nrows() as f64 - 1.0;
    let max_c = input.ncols() as f64 - 1.0;
    if pr < 0.0 || pc < 0.0 || pr > max_r || pc > max_c {
        return fill;
    }
    let r = (pr + 0.5).floor() as usize; // scipy order-0: round half up
    let c = (pc + 0.5).floor() as usize;
    input[(r, c)]
}

pub(crate) fn sample_bilinear(input: &Array2<f64>, pr: f64, pc: f64, fill: f64) -> f64 {
    let max_r = input.nrows() as f64 - 1.0;
    let max_c = input.ncols() as f64 - 1.0;
    if pr < 0.0 || pc < 0.0 || pr > max_r || pc > max_c {
        return fill;
    }
    let r0 = pr.floor();
    let c0 = pc.floor();
    let tr = pr - r0;
    let tc = pc - c0;
    let (ri, ci) = (r0 as usize, c0 as usize);
    let r1 = (ri + 1).min(input.nrows() - 1);
    let c1 = (ci + 1).min(input.ncols() - 1);
    let top = input[(ri, ci)] * (1.0 - tc) + input[(ri, c1)] * tc;
    let bot = input[(r1, ci)] * (1.0 - tc) + input[(r1, c1)] * tc;
    top * (1.0 - tr) + bot * tr
}

#[cfg(test)]
mod tests {
    use super::*;
    use ndarray::array;

    #[test]
    fn identity_at_zero() {
        let a = array![[1.0, 2.0, 3.0], [4.0, 5.0, 6.0]];
        let out = rotate(&a, 0.0, true, 0);
        assert_eq!(out, a);
    }

    #[test]
    fn quarter_turn_matches_scipy() {
        // scipy.ndimage.rotate([[1,2,3],[4,5,6]], 90, order=0)
        // -> [[3, 6], [2, 5], [1, 4]]
        let a = array![[1.0, 2.0, 3.0], [4.0, 5.0, 6.0]];
        let out = rotate(&a, 90.0, true, 0);
        assert_eq!(out.dim(), (3, 2));
        assert_eq!(out, array![[3.0, 6.0], [2.0, 5.0], [1.0, 4.0]]);
    }

    #[test]
    fn crop_keeps_shape_and_zeroes_outside() {
        // Ground truth from scipy itself:
        //   ndimage.rotate(u, 45, reshape=False, order=0, mode='constant', cval=0.0)
        //   -> [[0, 0], [4, 0], [0, 6], [0, 0]]
        // Note (0,0) and (3,1) map to coordinates just outside [0, n-1], so
        // scipy returns cval there rather than the nearest edge cell.
        let a = array![[1.0, 2.0], [3.0, 4.0], [5.0, 6.0], [7.0, 8.0]];
        let out = rotate(&a, 45.0, false, 0);
        assert_eq!(out.dim(), (4, 2));
        assert_eq!(out, array![[0.0, 0.0], [4.0, 0.0], [0.0, 6.0], [0.0, 0.0]]);
    }

    #[test]
    fn zero_fill_mode_untouched() {
        // The scipy-parity `rotate` zero-fills out-of-bounds samples.
        let a = array![[1.0, 2.0], [3.0, 4.0], [5.0, 6.0], [7.0, 8.0]];
        let out = rotate(&a, 45.0, false, 0);
        assert_eq!(out[(3, 0)], 0.0);
    }

    #[test]
    fn reshape_grows_the_canvas() {
        let a = Array2::<f64>::from_elem((80, 300), 1.0);
        let out45 = rotate(&a, 45.0, true, 0);
        // ptp for 45 deg: (80+300) * sqrt(2)/2 = 268.7 -> 269 both ways
        assert_eq!(out45.dim(), (269, 269));
        let out11 = rotate(&a, 11.0, true, 0);
        // bounds r: 80*cos11 + 300*sin11 = 78.53 + 57.24 = 135.77 -> 136
        // bounds c: 300*cos11 + 80*sin11 = 294.49 + 15.26 = 309.75 -> 310
        assert_eq!(out11.dim(), (136, 310));
    }

    #[test]
    fn negative_and_large_angles() {
        let a = array![[1.0, 2.0, 3.0], [4.0, 5.0, 6.0]];
        let out = rotate(&a, -270.0, true, 0); // == +90
        assert_eq!(out, rotate(&a, 90.0, true, 0));
        let out2 = rotate(&a, 450.0, true, 0);
        assert_eq!(out2, rotate(&a, 90.0, true, 0));
    }
}