radsym 0.4.1

//! Gaussian blur for response maps.
//!
//! For small sigma (≤ 2.0) uses direct separable convolution.
//! For larger sigma uses a 3-pass stacked box blur approximation,
//! which is O(1) per pixel regardless of sigma.
//!
//! Reference for stacked box blur:
//! - Wells, W.M. (1986). *Efficient Synthesis of Gaussian Filters by
//!   Cascaded Uniform Filters.* IEEE TPAMI 8(2).
//! - W3C Filter Effects Module Level 1, §12.2.

use super::image_view::OwnedImage;
use super::scalar::Scalar;

/// Gaussian blur on an `OwnedImage<f32>`, in-place.
///
/// Uses a 1D Gaussian kernel with `radius = ceil(3 * sigma)` and mirror-clamp
/// boundary handling for small sigma. For sigma > 2.0, switches to a 3-pass
/// stacked box blur approximation (O(1) per pixel).
///
/// No-op if `sigma <= 0.5` (kernel radius would be zero).
///
/// Allocates a `w * h` scratch buffer per call. Hot paths that blur repeatedly
/// (e.g. per-radius FRST) should call [`gaussian_blur_inplace_buf`] with a
/// reused buffer instead.
pub(crate) fn gaussian_blur_inplace(image: &mut OwnedImage<Scalar>, sigma: Scalar) {
    let mut scratch = Vec::new();
    gaussian_blur_inplace_buf(image, sigma, &mut scratch);
}

/// Gaussian blur, in-place, reusing a caller-owned scratch buffer.
///
/// Identical result to [`gaussian_blur_inplace`]; the only difference is that
/// the `w * h` temporary is taken from `scratch` (resized as needed) instead of
/// freshly allocated, so repeated calls at the same resolution allocate and
/// fault those pages once rather than every call. `scratch` contents are fully
/// overwritten before being read, so its incoming value is irrelevant.
pub(crate) fn gaussian_blur_inplace_buf(
    image: &mut OwnedImage<Scalar>,
    sigma: Scalar,
    scratch: &mut Vec<Scalar>,
) {
    if sigma <= 0.5 {
        return;
    }
    if sigma <= 2.0 {
        direct_gaussian_blur_inplace(image, sigma, scratch);
    } else {
        stacked_box_blur_inplace(image, sigma, scratch);
    }
}

/// Resize `buf` to exactly `len` elements (the fill value is irrelevant because
/// every blur pass overwrites the buffer before reading it).
#[inline]
fn ensure_buf(buf: &mut Vec<Scalar>, len: usize) {
    if buf.len() != len {
        buf.clear();
        buf.resize(len, 0.0);
    }
}

/// Direct separable Gaussian convolution (original implementation).
///
/// Cost: O(w * h * ceil(3*sigma)) per pass.
fn direct_gaussian_blur_inplace(
    image: &mut OwnedImage<Scalar>,
    sigma: Scalar,
    buf: &mut Vec<Scalar>,
) {
    let w = image.width();
    let h = image.height();

    let krad = (3.0 * sigma).ceil() as usize;
    if krad == 0 {
        return;
    }
    let ksize = 2 * krad + 1;

    // Build 1D Gaussian kernel
    let mut kernel = vec![0.0f32; ksize];
    let s2 = 2.0 * sigma * sigma;
    let mut sum = 0.0f32;
    for (i, k) in kernel.iter_mut().enumerate() {
        let d = i as Scalar - krad as Scalar;
        *k = (-d * d / s2).exp();
        sum += *k;
    }
    for k in &mut kernel {
        *k /= sum;
    }

    // Horizontal pass
    ensure_buf(buf, w * h);
    let data = image.data();
    for y in 0..h {
        for x in 0..w {
            let mut acc = 0.0f32;
            for (ki, &kv) in kernel.iter().enumerate() {
                let sx = (x as i32 + ki as i32 - krad as i32).clamp(0, w as i32 - 1) as usize;
                acc += data[y * w + sx] * kv;
            }
            buf[y * w + x] = acc;
        }
    }

    // Vertical pass
    let out = image.data_mut();
    for y in 0..h {
        for x in 0..w {
            let mut acc = 0.0f32;
            for (ki, &kv) in kernel.iter().enumerate() {
                let sy = (y as i32 + ki as i32 - krad as i32).clamp(0, h as i32 - 1) as usize;
                acc += buf[sy * w + x] * kv;
            }
            out[y * w + x] = acc;
        }
    }
}

// ---------------------------------------------------------------------------
// Stacked box blur: 3-pass O(1)-per-pixel Gaussian approximation
// ---------------------------------------------------------------------------

/// Compute box radii for a 3-pass stacked box blur approximating a Gaussian
/// with standard deviation `sigma`.
///
/// Returns three box half-widths (radii). A box of radius `r` has width `2r+1`.
fn box_radii_for_sigma(sigma: Scalar) -> [usize; 3] {
    // Ideal box width: w_ideal = sqrt(12 * sigma^2 / N + 1), N = 3 passes
    let w_ideal = (12.0 * sigma * sigma / 3.0 + 1.0).sqrt();
    let wl_raw = w_ideal.floor() as usize;
    let wl = if wl_raw.is_multiple_of(2) {
        wl_raw - 1
    } else {
        wl_raw
    }; // largest odd <= w_ideal
    let wu = wl + 2;

    // How many passes use wl vs wu:
    // variance of N box blurs = N * (w^2 - 1) / 12
    // We need: m * (wl^2-1)/12 + (3-m) * (wu^2-1)/12 = sigma^2
    // Solving for m:
    let target_var = 12.0 * sigma * sigma;
    let wl2 = (wl * wl) as Scalar;
    let wu2 = (wu * wu) as Scalar;
    let m_ideal = (3.0 * wu2 - 3.0 - target_var) / (wu2 - wl2);
    let m = m_ideal.round().clamp(0.0, 3.0) as usize;

    let rl = wl / 2;
    let ru = wu / 2;
    match m {
        0 => [ru, ru, ru],
        1 => [rl, ru, ru],
        2 => [rl, rl, ru],
        _ => [rl, rl, rl],
    }
}

/// 3-pass stacked box blur approximating Gaussian with standard deviation `sigma`.
///
/// Each box blur pass is O(1) per pixel using running sums, making the total
/// cost independent of sigma. Mirror-clamp boundary handling.
fn stacked_box_blur_inplace(image: &mut OwnedImage<Scalar>, sigma: Scalar, buf: &mut Vec<Scalar>) {
    let w = image.width();
    let h = image.height();
    ensure_buf(buf, w * h);

    let radii = box_radii_for_sigma(sigma);
    for &r in &radii {
        if r == 0 {
            continue;
        }
        box_blur_horizontal(image.data(), buf, w, h, r);
        box_blur_vertical(buf, image.data_mut(), w, h, r);
    }
}

/// Horizontal box blur pass using running sums. O(1) per pixel.
///
/// Reads from `src`, writes to `dst`. Box radius `r` gives window width `2r+1`.
/// Uses mirror-clamp boundary handling.
fn box_blur_horizontal(src: &[Scalar], dst: &mut [Scalar], w: usize, h: usize, r: usize) {
    let diameter = 2 * r + 1;
    let inv = 1.0 / diameter as Scalar;
    let w_i32 = w as i32;

    for y in 0..h {
        let row = y * w;

        // Initialize running sum with the first window (centered at x=0)
        let mut sum = 0.0f32;
        for i in 0..diameter {
            let sx = (i as i32 - r as i32).clamp(0, w_i32 - 1) as usize;
            sum += src[row + sx];
        }
        dst[row] = sum * inv;

        // Slide the window across the row
        for x in 1..w {
            // Add entering pixel (right edge of new window)
            let enter = (x as i32 + r as i32).clamp(0, w_i32 - 1) as usize;
            // Remove leaving pixel (left edge of old window)
            let leave = (x as i32 - r as i32 - 1).clamp(0, w_i32 - 1) as usize;
            sum += src[row + enter] - src[row + leave];
            dst[row + x] = sum * inv;
        }
    }
}

/// Number of adjacent columns processed together in the vertical pass.
///
/// A vertical box blur walks *down* columns; on a row-major buffer that is a
/// stride-`w` access that touches a fresh cache line almost every step. Handling
/// a strip of `VSTRIP` contiguous columns at once turns each row access into one
/// fully-used 64-byte cache line and exposes `VSTRIP` independent running sums
/// that the autovectorizer maps onto SIMD lanes. `16` = one cache line of `f32`
/// and a clean multiple of both NEON (4) and AVX (8) widths.
const VSTRIP: usize = 16;

/// Vertical box blur pass using running sums. O(1) per pixel.
///
/// Reads from `src`, writes to `dst`. Box radius `r` gives window width `2r+1`.
/// Uses mirror-clamp boundary handling.
///
/// Columns are processed in [`VSTRIP`]-wide strips so the memory access is
/// cache-line-contiguous and vectorizable. Each column's running sum is updated
/// in exactly the same order as a naive per-column sweep, so the output is
/// bit-identical to the scalar reference (see `vertical_strip_matches_naive`).
fn box_blur_vertical(src: &[Scalar], dst: &mut [Scalar], w: usize, h: usize, r: usize) {
    let diameter = 2 * r + 1;
    let inv = 1.0 / diameter as Scalar;
    let h_i32 = h as i32;

    let mut x0 = 0;
    while x0 < w {
        let cw = VSTRIP.min(w - x0);

        // One running sum per column in the strip.
        // One running sum per column in the strip.
        let mut sums = [0.0f32; VSTRIP];
        for i in 0..diameter {
            let sy = (i as i32 - r as i32).clamp(0, h_i32 - 1) as usize;
            let base = sy * w + x0;
            for c in 0..cw {
                sums[c] += src[base + c];
            }
        }
        for c in 0..cw {
            dst[x0 + c] = sums[c] * inv;
        }

        // Slide all columns of the strip down together.
        for y in 1..h {
            let enter = (y as i32 + r as i32).clamp(0, h_i32 - 1) as usize;
            let leave = (y as i32 - r as i32 - 1).clamp(0, h_i32 - 1) as usize;
            let ebase = enter * w + x0;
            let lbase = leave * w + x0;
            let obase = y * w + x0;
            for c in 0..cw {
                sums[c] += src[ebase + c] - src[lbase + c];
                dst[obase + c] = sums[c] * inv;
            }
        }

        x0 += VSTRIP;
    }
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::core::image_view::OwnedImage;

    /// Stacked box blur should approximate a Gaussian within tolerance.
    ///
    /// Applies both methods to a centered impulse and checks that the
    /// peak-normalized maximum absolute difference is small.
    #[test]
    fn box_blur_approximates_gaussian() {
        let size = 64;
        let sigma = 10.0;

        // Create impulse image
        let mut data = vec![0.0f32; size * size];
        data[size / 2 * size + size / 2] = 1.0;

        let mut img_direct = OwnedImage::from_vec(data.clone(), size, size).unwrap();
        let mut img_box = OwnedImage::from_vec(data, size, size).unwrap();

        direct_gaussian_blur_inplace(&mut img_direct, sigma, &mut Vec::new());
        stacked_box_blur_inplace(&mut img_box, sigma, &mut Vec::new());

        let d = img_direct.data();
        let b = img_box.data();

        let peak = d.iter().copied().fold(0.0f32, Scalar::max);
        assert!(peak > 0.0, "direct Gaussian peak should be positive");

        let max_err = d
            .iter()
            .zip(b.iter())
            .map(|(&dv, &bv)| (dv - bv).abs())
            .fold(0.0f32, Scalar::max);

        let relative_err = max_err / peak;
        assert!(
            relative_err < 0.10,
            "box blur should approximate Gaussian within 10%, got {:.1}%",
            relative_err * 100.0
        );
    }

    /// The cache-blocked vertical pass must be bit-identical to a naive
    /// per-column sweep, including partial (non-VSTRIP-multiple) tail strips.
    #[test]
    fn vertical_strip_matches_naive() {
        fn naive_vertical(src: &[Scalar], dst: &mut [Scalar], w: usize, h: usize, r: usize) {
            let diameter = 2 * r + 1;
            let inv = 1.0 / diameter as Scalar;
            let h_i32 = h as i32;
            for x in 0..w {
                let mut sum = 0.0f32;
                for i in 0..diameter {
                    let sy = (i as i32 - r as i32).clamp(0, h_i32 - 1) as usize;
                    sum += src[sy * w + x];
                }
                dst[x] = sum * inv;
                for y in 1..h {
                    let enter = (y as i32 + r as i32).clamp(0, h_i32 - 1) as usize;
                    let leave = (y as i32 - r as i32 - 1).clamp(0, h_i32 - 1) as usize;
                    sum += src[enter * w + x] - src[leave * w + x];
                    dst[y * w + x] = sum * inv;
                }
            }
        }

        // Widths chosen to exercise full strips, a single strip, and tails of
        // various widths (17 = 16 + 1, 30 = 16 + 14, 7 < 16).
        for &(w, h) in &[(7usize, 9usize), (16, 16), (17, 5), (30, 13), (64, 40)] {
            let mut src = vec![0.0f32; w * h];
            for (i, v) in src.iter_mut().enumerate() {
                // Deterministic non-trivial pattern.
                *v = ((i * 37 % 251) as f32) - 125.0;
            }
            for &r in &[1usize, 3, 6] {
                let mut a = vec![0.0f32; w * h];
                let mut b = vec![0.0f32; w * h];
                box_blur_vertical(&src, &mut a, w, h, r);
                naive_vertical(&src, &mut b, w, h, r);
                assert_eq!(a, b, "strip != naive for w={w} h={h} r={r}");
            }
        }
    }

    /// Box radii computation produces reasonable values.
    #[test]
    fn box_radii_sanity() {
        // For large sigma, all radii should be > 0
        let radii = box_radii_for_sigma(10.0);
        for &r in &radii {
            assert!(r > 0, "radius should be positive for sigma=10");
        }

        // Variance of 3 box blurs should approximate sigma^2
        let sigma = 10.0;
        let radii = box_radii_for_sigma(sigma);
        let total_var: f32 = radii
            .iter()
            .map(|&r| {
                let w = (2 * r + 1) as f32;
                (w * w - 1.0) / 12.0
            })
            .sum();
        let target_var = sigma * sigma;
        let var_err = (total_var - target_var).abs() / target_var;
        assert!(
            var_err < 0.1,
            "variance should approximate sigma^2 within 10%, got err={:.1}%",
            var_err * 100.0
        );
    }

    /// Blur preserves total energy (sum of pixels).
    #[test]
    fn blur_preserves_energy() {
        let size = 32;
        let mut data = vec![0.0f32; size * size];
        // Place some energy away from boundaries
        data[size / 2 * size + size / 2] = 100.0;
        data[size / 3 * size + size / 3] = 50.0;

        let sum_before: f32 = data.iter().sum();
        let mut img = OwnedImage::from_vec(data, size, size).unwrap();
        gaussian_blur_inplace(&mut img, 5.0);
        let sum_after: f32 = img.data().iter().sum();

        let energy_err = (sum_after - sum_before).abs() / sum_before;
        assert!(
            energy_err < 0.01,
            "blur should preserve energy within 1%, got err={:.2}%",
            energy_err * 100.0
        );
    }

    /// Dispatch: small sigma uses direct, large sigma uses box blur.
    /// Both should produce non-zero output for an impulse input.
    #[test]
    fn dispatch_both_paths() {
        for &sigma in &[1.0, 5.0, 15.0] {
            let size = 32;
            let mut data = vec![0.0f32; size * size];
            data[size / 2 * size + size / 2] = 1.0;
            let mut img = OwnedImage::from_vec(data, size, size).unwrap();
            gaussian_blur_inplace(&mut img, sigma);
            let peak = img.data().iter().copied().fold(0.0f32, Scalar::max);
            assert!(peak > 0.0, "sigma={sigma}: peak should be positive");
        }
    }
}