vello_cpu/filter/

gaussian_blur.rs

// Copyright 2025 the Vello Authors
// SPDX-License-Identifier: Apache-2.0 OR MIT

//! Gaussian blur filter implementation using multi-scale separable convolution.
//!
//! This implementation uses a multi-scale approach for efficient blurring:
//! - **Small blurs** (σ ≤ 2): Direct separable convolution at full resolution
//! - **Large blurs** (σ > 2): Iterative downsample → blur → upsample pyramid
//!
//! The algorithm automatically determines the optimal number of decimation levels
//! using variance analysis. Each 2× decimation applies a \[1,3,3,1\]/8 binomial filter
//! (adding variance = 3.0), then downsamples, reducing the remaining blur work needed.
//! This exploits the variance additivity property: `σ²_total = σ²_downsample + σ²_blur`.
//!
//! **Variance Addition Reference**: The convolution of two Gaussians with variances σ₁²
//! and σ₂² produces a Gaussian with variance σ₁² + σ₂². This fundamental property comes
//! from probability theory and applies to Gaussian convolution in image processing. See:
//! - Torralba & Freeman, "Foundations of Computer Vision" (MIT Press), Section 2.2:
//!   <https://visionbook.mit.edu/blurring_2.html#properties-of-the-continuous-gaussian>

use super::FilterEffect;
use crate::filter::context::ScratchBuffer;
use vello_common::filter::gaussian_blur::{DecimationSizer, GaussianBlur};
use vello_common::filter_effects::EdgeMode;
use vello_common::peniko::color::PremulRgba8;
#[cfg(not(feature = "std"))]
use vello_common::peniko::kurbo::common::FloatFuncs as _;
use vello_common::pixmap::Pixmap;

impl FilterEffect for GaussianBlur {
    fn execute_lowp(&self, pixmap: &mut Pixmap, filter_scratch: &mut ScratchBuffer) {
        // No blur if std_deviation is zero or negative
        if self.std_deviation <= 0.0 {
            return;
        }

        let scratch = filter_scratch.get_scratch_buffer(pixmap.width(), pixmap.height());
        apply_blur(
            pixmap,
            scratch,
            self.n_decimations,
            &self.kernel[..usize::from(self.kernel_size)],
            self.edge_mode,
        );
    }

    fn execute_highp(&self, pixmap: &mut Pixmap, filter_scratch: &mut ScratchBuffer) {
        // TODO: Currently only lowp is implemented and used for highp as well.
        // This needs to be updated to use proper high-precision arithmetic.
        Self::execute_lowp(self, pixmap, filter_scratch);
    }
}

/// Apply Gaussian blur using multi-scale decimation and upsampling.
///
/// Uses a precomputed decimation plan and kernel for optimal performance.
/// Operates in-place using a single pixmap buffer with logical dimension tracking
/// to minimize memory allocations. For `n_decimations=0`, applies direct convolution.
///
/// The `scratch` buffer is used for separable convolution and must be at least as
/// large as the source pixmap.
pub(crate) fn apply_blur(
    pixmap: &mut Pixmap,
    scratch: &mut Pixmap,
    n_decimations: usize,
    kernel: &[f32],
    edge_mode: EdgeMode,
) {
    let radius = (kernel.len() / 2) as u8;
    let width = pixmap.width();
    let height = pixmap.height();

    // Small blur: apply direct convolution at full resolution
    if n_decimations == 0 {
        convolve(pixmap, scratch, width, height, kernel, radius, edge_mode);
        return;
    }

    // Track logical dimensions through decimation (physical buffer stays the same size)
    let mut sizer = DecimationSizer::new(width, height);

    // Downsample n times (each step reduces resolution by 2×)
    for _ in 0..n_decimations {
        let (w, h) = sizer.current();
        downscale(pixmap, w, h, edge_mode);
        sizer.downscale();
    }

    // Apply the reduced blur at the coarsest resolution
    let (w, h) = sizer.current();
    convolve(pixmap, scratch, w, h, kernel, radius, edge_mode);

    // Upsample back to original resolution (each step doubles resolution by 2×)
    for _ in 0..n_decimations {
        let (w, h) = sizer.current();
        upscale(pixmap, w, h, edge_mode);
        sizer.upscale();
    }

    debug_assert_eq!(
        sizer.current(),
        (width, height),
        "Final dimensions should match original"
    );
}

/// Apply separable Gaussian convolution with logical dimensions.
///
/// Performs horizontal blur followed by vertical blur. Works with a logical view
/// of the pixmap, using only the top-left region defined by width × height.
/// The `temp` buffer is provided by the caller to avoid allocations.
pub(crate) fn convolve(
    src: &mut Pixmap,
    scratch: &mut Pixmap,
    width: u16,
    height: u16,
    kernel: &[f32],
    radius: u8,
    edge_mode: EdgeMode,
) {
    convolve_x(src, scratch, width, height, kernel, radius, edge_mode);
    convolve_y(scratch, src, width, height, kernel, radius, edge_mode);
}

/// Apply horizontal blur pass (1D convolution along x-axis).
///
/// For each output pixel, computes a weighted sum of horizontally neighboring pixels
/// using the Gaussian kernel. Handles edge cases according to the specified edge mode.
/// Writes results to a destination buffer to avoid overwriting source data.
pub(crate) fn convolve_x(
    src: &Pixmap,
    dst: &mut Pixmap,
    src_width: u16,
    src_height: u16,
    kernel: &[f32],
    radius: u8,
    edge_mode: EdgeMode,
) {
    for y in 0..src_height {
        for x in 0..src_width {
            let mut rgba = [0.0_f32; 4];

            // Sum contributions from all kernel positions: output = Σ(weight[j] × pixel[x+j-radius])
            for (j, &k) in kernel.iter().enumerate() {
                #[expect(
                    clippy::cast_possible_wrap,
                    reason = "This cast never wraps because `kernel.len()` is never greater than `u8::MAX` due to the restriction on `MAX_KERNEL_SIZE`"
                )]
                let j = j as i32;
                let src_x = x as i32 + j - radius as i32;
                let p = sample_x(src, src_x, y, src_width, edge_mode);

                rgba[0] += p.r as f32 * k;
                rgba[1] += p.g as f32 * k;
                rgba[2] += p.b as f32 * k;
                rgba[3] += p.a as f32 * k;
            }

            // Convert back to u8 with rounding
            dst.set_pixel(
                x,
                y,
                PremulRgba8 {
                    r: rgba[0].round() as u8,
                    g: rgba[1].round() as u8,
                    b: rgba[2].round() as u8,
                    a: rgba[3].round() as u8,
                },
            );
        }
    }
}

/// Apply vertical blur pass (1D convolution along y-axis).
///
/// For each output pixel, computes a weighted sum of vertically neighboring pixels
/// using the Gaussian kernel. Handles edge cases according to the specified edge mode.
/// Writes results to a destination buffer to avoid overwriting source data.
pub(crate) fn convolve_y(
    src: &Pixmap,
    dst: &mut Pixmap,
    src_width: u16,
    src_height: u16,
    kernel: &[f32],
    radius: u8,
    edge_mode: EdgeMode,
) {
    for y in 0..src_height {
        for x in 0..src_width {
            let mut rgba = [0.0_f32; 4];

            // Sum contributions from all kernel positions: output = Σ(weight[j] × pixel[y+j-radius])
            for (j, &k) in kernel.iter().enumerate() {
                #[expect(
                    clippy::cast_possible_wrap,
                    reason = "This cast never wraps because `kernel.len()` is never greater than `u8::MAX` due to the restriction on `MAX_KERNEL_SIZE`"
                )]
                let j = j as i32;
                let src_y = y as i32 + j - radius as i32;
                let p = sample_y(src, x, src_y, src_height, edge_mode);

                rgba[0] += p.r as f32 * k;
                rgba[1] += p.g as f32 * k;
                rgba[2] += p.b as f32 * k;
                rgba[3] += p.a as f32 * k;
            }

            // Convert back to u8 with rounding
            dst.set_pixel(
                x,
                y,
                PremulRgba8 {
                    r: rgba[0].round() as u8,
                    g: rgba[1].round() as u8,
                    b: rgba[2].round() as u8,
                    a: rgba[3].round() as u8,
                },
            );
        }
    }
}

/// Downsample image by 2x using separable \[1,3,3,1\]/8 binomial filter.
///
/// Performs horizontal and vertical decimation in sequence. Returns the new
/// logical dimensions (ceil(width/2), ceil(height/2)).
pub(crate) fn downscale(
    src: &mut Pixmap,
    src_width: u16,
    src_height: u16,
    edge_mode: EdgeMode,
) -> (u16, u16) {
    let dst_width = src_width.div_ceil(2);
    let dst_height = src_height.div_ceil(2);
    downscale_x(src, src_width, src_height, dst_width, edge_mode);
    // We can pass `dst_width` instead of `src_width` here, since we already decimated
    // horizontally.
    downscale_y(src, dst_width, src_height, dst_height, edge_mode);
    (dst_width, dst_height)
}

/// Horizontal decimation pass using \[1,3,3,1\]/8 filter.
///
/// Reduces width by 2x while applying a binomial blur kernel. The \[1,3,3,1\] weights
/// approximate a Gaussian and contribute variance=0.75 before the 2x downsampling.
fn downscale_x(
    src: &mut Pixmap,
    src_width: u16,
    src_height: u16,
    dst_width: u16,
    edge_mode: EdgeMode,
) {
    for y in 0..src_height {
        // Sliding window: maintains 2 previous pixels (x0, x1) to form 4-tap filter
        // with current pixels (x2, x3). Start with sample at -1 for implicit left padding.
        let mut p0 = sample_x(src, -1, y, src_width, edge_mode);
        let mut p1 = sample_x(src, 0, y, src_width, edge_mode);

        for x in 0..dst_width {
            // Sample 4 horizontally adjacent pixels for [1,3,3,1] kernel
            // Pattern: [x*2-1, x*2, x*2+1, x*2+2]
            let src_x = (x * 2) as i32;
            let p2 = sample_x(src, src_x + 1, y, src_width, edge_mode);
            let p3 = sample_x(src, src_x + 2, y, src_width, edge_mode);

            // Apply [1,3,3,1]/8 weights → output = (p0 + 3×p1 + 3×p2 + p3) / 8
            src.set_pixel(x, y, decimate_weighted(p0, p1, p2, p3));

            // Advance window: previous p2,p3 become next p0,p1
            p0 = p2;
            p1 = p3;
        }
    }
}

/// Vertical decimation pass using \[1,3,3,1\]/8 filter.
///
/// Reduces logical height by 2x while applying a binomial blur kernel.
/// Operates in-place by writing to the beginning of the same buffer.
fn downscale_y(
    src: &mut Pixmap,
    src_width: u16,
    src_height: u16,
    dst_height: u16,
    edge_mode: EdgeMode,
) {
    for x in 0..src_width {
        // Sliding window: maintains 2 previous pixels (y0, y1) to form 4-tap filter
        // with current pixels (y2, y3). Start with sample at -1 for implicit top padding.
        let mut p0 = sample_y(src, x, -1, src_height, edge_mode);
        let mut p1 = sample_y(src, x, 0, src_height, edge_mode);

        for y in 0..dst_height {
            // Sample 4 vertically adjacent pixels for [1,3,3,1] kernel
            // Pattern: [y*2-1, y*2, y*2+1, y*2+2]
            let src_y = (y * 2) as i32;
            let p2 = sample_y(src, x, src_y + 1, src_height, edge_mode);
            let p3 = sample_y(src, x, src_y + 2, src_height, edge_mode);

            // Apply [1,3,3,1]/8 weights → output = (p0 + 3×p1 + 3×p2 + p3) / 8
            src.set_pixel(x, y, decimate_weighted(p0, p1, p2, p3));

            // Advance window: previous p2,p3 become next p0,p1
            p0 = p2;
            p1 = p3;
        }
    }
}

/// Upsample a pixmap by 2x using linear interpolation with [0.75, 0.25] weights.
///
/// Uses separable passes: horizontal doubling followed by vertical doubling.
///
/// ## Phase Alignment Theory
///
/// The downsampling \[1,3,3,1\]/8 filter creates a half-pixel offset. Its center of mass
/// is at position `2k + 0.5`, not `2k`. This means each downsampled pixel at discrete
/// position `k` represents a sample at continuous position `2k + 0.5`.
///
/// When upsampling, we reconstruct pixels at positions `2k` and `2k+1` from downsampled
/// pixels whose centers are at `..., 2k-1.5, 2k+0.5, 2k+2.5, ...`
///
/// **Interpolation weights** are derived from linear interpolation based on distances:
/// - Position `2k`:   distance 0.5 from center at `2k+0.5`, distance 1.5 from center at `2k-1.5`
///   → weights: 0.75×pixel\[k\] + 0.25×pixel\[k-1\]
/// - Position `2k+1`: distance 0.5 from center at `2k+0.5`, distance 1.5 from center at `2k+2.5`
///   → weights: 0.75×pixel\[k\] + 0.25×pixel\[k+1\]
pub(crate) fn upscale(
    src: &mut Pixmap,
    src_width: u16,
    src_height: u16,
    edge_mode: EdgeMode,
) -> (u16, u16) {
    let dst_width = src_width * 2;
    let dst_height = src_height * 2;
    upscale_x(src, src_width, src_height, edge_mode);
    upscale_y(src, dst_width, src_height, edge_mode);
    (dst_width, dst_height)
}

/// Horizontal upsampling pass using [0.75, 0.25] interpolation with logical dimensions.
///
/// Doubles the logical width using phase-aligned interpolation. Each input pixel
/// generates two output pixels with different weights based on their distance from
/// the downsampled pixel's center position.
/// Operates in-place by processing backwards to avoid overwriting source data.
fn upscale_x(src: &mut Pixmap, src_width: u16, src_height: u16, edge_mode: EdgeMode) {
    // Process backwards (right to left) to avoid overwriting source data
    for y in 0..src_height {
        // Maintain sliding window of three pixels: prev, current, next
        // This allows us to compute both output pixels that depend on current pixel x
        let mut p0 = sample_x(src, src_width as i32, y, src_width, edge_mode);
        let mut p1 = sample_x(src, src_width as i32 - 1, y, src_width, edge_mode);

        for x in (0..src_width).rev() {
            let src_x = x as i32;
            let p2 = sample_x(src, src_x - 1, y, src_width, edge_mode);

            // Generate two output pixels per input with phase-aligned interpolation:
            // output[2x]   = 0.25×p2 + 0.75×p1  (position 2x   is 0.5 from center at 2x+0.5)
            // output[2x+1] = 0.75×p1 + 0.25×p0  (position 2x+1 is 0.5 from center at 2x+0.5)
            let dst_x = x * 2;
            src.set_pixel(dst_x, y, interpolate_25_75(p2, p1));
            src.set_pixel(dst_x + 1, y, interpolate_75_25(p1, p0));

            // Advance sliding window for next iteration
            p0 = p1;
            p1 = p2;
        }
    }
}

/// Vertical upsampling pass using [0.75, 0.25] interpolation with logical dimensions.
///
/// Doubles the logical height using phase-aligned interpolation. Each input pixel
/// generates two output pixels with different weights based on their distance from
/// the downsampled pixel's center position.
/// Operates in-place by processing backwards to avoid overwriting source data.
fn upscale_y(src: &mut Pixmap, src_width: u16, src_height: u16, edge_mode: EdgeMode) {
    // Process backwards (bottom to top) to avoid overwriting source data
    for x in 0..src_width {
        // Maintain sliding window of three pixels: prev, current, next
        // This allows us to compute both output pixels that depend on current pixel y
        let mut p0 = sample_y(src, x, src_height as i32, src_height, edge_mode);
        let mut p1 = sample_y(src, x, src_height as i32 - 1, src_height, edge_mode);

        for y in (0..src_height).rev() {
            let src_y = y as i32;
            let p2 = sample_y(src, x, src_y - 1, src_height, edge_mode);

            // Generate two output rows per input with phase-aligned interpolation:
            // output[2y]   = 0.25×p2 + 0.75×p1  (position 2y   is 0.5 from center at 2y+0.5)
            // output[2y+1] = 0.75×p1 + 0.25×p0  (position 2y+1 is 0.5 from center at 2y+0.5)
            let dst_y = y * 2;
            src.set_pixel(x, dst_y, interpolate_25_75(p2, p1));
            src.set_pixel(x, dst_y + 1, interpolate_75_25(p1, p0));

            // Advance sliding window for next iteration
            p0 = p1;
            p1 = p2;
        }
    }
}

/// Transparent black pixel constant.
const TRANSPARENT_BLACK: PremulRgba8 = PremulRgba8 {
    r: 0,
    g: 0,
    b: 0,
    a: 0,
};

/// Sample a pixel with edge mode handling for horizontal sampling.
#[inline(always)]
fn sample_x(src: &Pixmap, x: i32, y: u16, width: u16, edge_mode: EdgeMode) -> PremulRgba8 {
    sample(x, width, edge_mode, |src_x| src.sample(src_x, y))
}

/// Sample a pixel with edge mode handling for vertical sampling.
#[inline(always)]
fn sample_y(src: &Pixmap, x: u16, y: i32, height: u16, edge_mode: EdgeMode) -> PremulRgba8 {
    sample(y, height, edge_mode, |src_y| src.sample(x, src_y))
}

/// Sample a pixel with edge mode handling (generic implementation).
///
/// Handles both horizontal and vertical sampling based on the provided closure.
/// The `sample_fn` closure receives the clamped/extended coordinate and returns the pixel.
/// For `EdgeMode::None`, returns transparent black if the coordinate is out of bounds.
#[inline(always)]
fn sample<F>(coord: i32, size: u16, edge_mode: EdgeMode, sample_fn: F) -> PremulRgba8
where
    F: FnOnce(u16) -> PremulRgba8,
{
    // For EdgeMode::None, return transparent black if out of bounds
    if edge_mode == EdgeMode::None && (coord < 0 || coord >= size as i32) {
        return TRANSPARENT_BLACK;
    }
    let extended_coord = extend(coord, size, edge_mode);
    sample_fn(extended_coord)
}

/// Extend a coordinate beyond image boundaries according to the edge mode.
///
/// Transforms out-of-bounds coordinates for sampling: clamped, wrapped, or mirrored
/// depending on the mode. For `EdgeMode::None`, the coordinate is guaranteed to be
/// in-bounds (already checked by caller, which returns transparent black for out-of-bounds).
#[inline(always)]
fn extend(coord: i32, size: u16, edge_mode: EdgeMode) -> u16 {
    match edge_mode {
        EdgeMode::Duplicate => {
            // Clamp to image bounds: pixels outside use nearest edge pixel
            coord.clamp(0, size as i32 - 1) as u16
        }
        EdgeMode::None => {
            // Coordinate is already validated as in-bounds by caller
            coord as u16
        }
        EdgeMode::Wrap => {
            // Wrap around using modulo: image tiles infinitely
            let mut c = coord % size as i32;
            if c < 0 {
                c += size as i32;
            }
            c as u16
        }
        EdgeMode::Mirror => {
            // Mirror at boundaries: image reflects across edges
            let period = size as i32 * 2;
            let mut c = coord % period;
            if c < 0 {
                c += period;
            }
            if c >= size as i32 {
                c = period - c - 1;
            }
            c as u16
        }
    }
}

/// Blend 4 RGBA pixels using \[1,3,3,1\]/8 binomial weights.
///
/// Computes `(p0 + 3×p1 + 3×p2 + p3) / 8` using efficient integer arithmetic (right shift).
/// This binomial pattern approximates a Gaussian and contributes variance=0.75 to the blur.
/// Adds 4 before the shift to implement round-to-nearest instead of floor division.
#[inline(always)]
fn decimate_weighted(
    p0: PremulRgba8,
    p1: PremulRgba8,
    p2: PremulRgba8,
    p3: PremulRgba8,
) -> PremulRgba8 {
    let r = ((p0.r as u32 + p1.r as u32 * 3 + p2.r as u32 * 3 + p3.r as u32 + 4) >> 3) as u8;
    let g = ((p0.g as u32 + p1.g as u32 * 3 + p2.g as u32 * 3 + p3.g as u32 + 4) >> 3) as u8;
    let b = ((p0.b as u32 + p1.b as u32 * 3 + p2.b as u32 * 3 + p3.b as u32 + 4) >> 3) as u8;
    let a = ((p0.a as u32 + p1.a as u32 * 3 + p2.a as u32 * 3 + p3.a as u32 + 4) >> 3) as u8;
    PremulRgba8 { r, g, b, a }
}

/// Blend 2 RGBA pixels using [0.25, 0.75] weights (right-weighted interpolation).
///
/// Used for upsampling to generate the second of two output pixels.
/// Favors the right/next pixel (p1) with 75% weight.
/// Adds 2 before the shift to implement round-to-nearest instead of floor division.
#[inline(always)]
fn interpolate_25_75(p0: PremulRgba8, p1: PremulRgba8) -> PremulRgba8 {
    let r = ((p0.r as u32 + p1.r as u32 * 3 + 2) >> 2) as u8;
    let g = ((p0.g as u32 + p1.g as u32 * 3 + 2) >> 2) as u8;
    let b = ((p0.b as u32 + p1.b as u32 * 3 + 2) >> 2) as u8;
    let a = ((p0.a as u32 + p1.a as u32 * 3 + 2) >> 2) as u8;
    PremulRgba8 { r, g, b, a }
}

/// Blend 2 RGBA pixels using [0.75, 0.25] weights.
///
/// Computes `0.75×p0 + 0.25×p1` using efficient integer arithmetic.
/// Used during upsampling for positions closer to the first pixel.
/// Adds 2 before the shift to implement round-to-nearest instead of floor division.
#[inline(always)]
fn interpolate_75_25(p0: PremulRgba8, p1: PremulRgba8) -> PremulRgba8 {
    let r = ((p0.r as u32 * 3 + p1.r as u32 + 2) >> 2) as u8;
    let g = ((p0.g as u32 * 3 + p1.g as u32 + 2) >> 2) as u8;
    let b = ((p0.b as u32 * 3 + p1.b as u32 + 2) >> 2) as u8;
    let a = ((p0.a as u32 * 3 + p1.a as u32 + 2) >> 2) as u8;
    PremulRgba8 { r, g, b, a }
}

#[cfg(test)]
mod tests {
    use super::*;
    use vello_common::filter::gaussian_blur::{
        MAX_KERNEL_SIZE, compute_gaussian_kernel, plan_decimated_blur,
    };

    /// Test edge extension with Duplicate mode.
    #[test]
    fn test_extend_duplicate() {
        let size = 10;

        // In-bounds coordinate
        assert_eq!(extend(5, size, EdgeMode::Duplicate), 5);

        // Below bounds: clamp to 0
        assert_eq!(extend(-1, size, EdgeMode::Duplicate), 0);
        assert_eq!(extend(-10, size, EdgeMode::Duplicate), 0);

        // Above bounds: clamp to size-1
        assert_eq!(extend(10, size, EdgeMode::Duplicate), 9);
        assert_eq!(extend(20, size, EdgeMode::Duplicate), 9);
    }

    /// Test edge extension with Wrap mode.
    #[test]
    fn test_extend_wrap() {
        let size = 10;

        // In-bounds: identity
        assert_eq!(extend(5, size, EdgeMode::Wrap), 5);

        // Above bounds: wrap around
        assert_eq!(extend(10, size, EdgeMode::Wrap), 0);
        assert_eq!(extend(11, size, EdgeMode::Wrap), 1);
        assert_eq!(extend(25, size, EdgeMode::Wrap), 5);

        // Below bounds: wrap from end
        assert_eq!(extend(-1, size, EdgeMode::Wrap), 9);
        assert_eq!(extend(-2, size, EdgeMode::Wrap), 8);
    }

    /// Test edge extension with Mirror mode.
    #[test]
    fn test_extend_mirror() {
        let size = 10;

        // In-bounds: identity
        assert_eq!(extend(5, size, EdgeMode::Mirror), 5);

        // Just past boundary: mirror back
        // For coord=10: period=20, c=10, c>=10 so c = 20-10-1 = 9
        assert_eq!(extend(10, size, EdgeMode::Mirror), 9);
        // For coord=11: period=20, c=11, c>=10 so c = 20-11-1 = 8
        assert_eq!(extend(11, size, EdgeMode::Mirror), 8);

        // Full reflection cycle (period = 2*size = 20)
        // For coord=19: c=19, c>=10 so c = 20-19-1 = 0
        assert_eq!(extend(19, size, EdgeMode::Mirror), 0);
        // For coord=20: c=20%20=0, c<10 so c=0
        assert_eq!(extend(20, size, EdgeMode::Mirror), 0);

        // Negative coordinates
        // For coord=-1: c=-1%20=-1, c<0 so c+=20 → c=19, c>=10 so c=20-19-1=0
        assert_eq!(extend(-1, size, EdgeMode::Mirror), 0);
        // For coord=-2: c=-2%20=-2, c<0 so c+=20 → c=18, c>=10 so c=20-18-1=1
        assert_eq!(extend(-2, size, EdgeMode::Mirror), 1);
    }

    /// Test edge extension with None mode (coordinate should pass through).
    #[test]
    fn test_extend_none() {
        let size = 10;
        // None mode just passes the coordinate as u16
        assert_eq!(extend(5, size, EdgeMode::None), 5);
        assert_eq!(extend(9, size, EdgeMode::None), 9);
        assert_eq!(extend(10, size, EdgeMode::None), 10);
        assert_eq!(extend(11, size, EdgeMode::None), 11);
    }

    /// Test [1,3,3,1]/8 decimation weights.
    #[test]
    fn test_decimate_weighted() {
        let p0 = PremulRgba8 {
            r: 0,
            g: 0,
            b: 0,
            a: 0,
        };
        let p1 = PremulRgba8 {
            r: 8,
            g: 8,
            b: 8,
            a: 8,
        };
        let p2 = PremulRgba8 {
            r: 8,
            g: 8,
            b: 8,
            a: 8,
        };
        let p3 = PremulRgba8 {
            r: 0,
            g: 0,
            b: 0,
            a: 0,
        };

        // (0 + 3*8 + 3*8 + 0 + 4) / 8 = (48 + 4) / 8 = 52 / 8 = 6.5 → 6 (rounds down)
        let result = decimate_weighted(p0, p1, p2, p3);
        assert_eq!(result.r, 6);
        assert_eq!(result.g, 6);
        assert_eq!(result.b, 6);
        assert_eq!(result.a, 6);
    }

    /// Test [1,3,3,1]/8 with all same values (should be identity).
    #[test]
    fn test_decimate_weighted_uniform() {
        let p = PremulRgba8 {
            r: 100,
            g: 150,
            b: 200,
            a: 255,
        };
        let result = decimate_weighted(p, p, p, p);
        // All same: (100 + 300 + 300 + 100 + 4) / 8 = 804/8 = 100.5 → 100 (rounds down)
        assert_eq!(result.r, 100);
        assert_eq!(result.g, 150);
        assert_eq!(result.b, 200);
        assert_eq!(result.a, 255);
    }

    /// Test [0.25, 0.75] interpolation.
    #[test]
    fn test_interpolate_25_75() {
        let p0 = PremulRgba8 {
            r: 0,
            g: 0,
            b: 0,
            a: 0,
        };
        let p1 = PremulRgba8 {
            r: 100,
            g: 100,
            b: 100,
            a: 100,
        };

        // (0 + 300 + 2) / 4 = 302 / 4 = 75.5 → 75 (rounds down)
        let result = interpolate_25_75(p0, p1);
        assert_eq!(result.r, 75);
        assert_eq!(result.g, 75);
        assert_eq!(result.b, 75);
        assert_eq!(result.a, 75);
    }

    /// Test [0.75, 0.25] interpolation.
    #[test]
    fn test_interpolate_75_25() {
        let p0 = PremulRgba8 {
            r: 100,
            g: 100,
            b: 100,
            a: 100,
        };
        let p1 = PremulRgba8 {
            r: 0,
            g: 0,
            b: 0,
            a: 0,
        };

        // (300 + 0 + 2) / 4 = 302 / 4 = 75.5 → 75 (rounds down)
        let result = interpolate_75_25(p0, p1);
        assert_eq!(result.r, 75);
        assert_eq!(result.g, 75);
        assert_eq!(result.b, 75);
        assert_eq!(result.a, 75);
    }

    /// Test interpolation symmetry.
    #[test]
    fn test_interpolation_symmetry() {
        let p0 = PremulRgba8 {
            r: 50,
            g: 100,
            b: 150,
            a: 200,
        };
        let p1 = PremulRgba8 {
            r: 200,
            g: 150,
            b: 100,
            a: 50,
        };

        let r1 = interpolate_25_75(p0, p1);
        let r2 = interpolate_75_25(p1, p0);

        // Should be symmetric
        assert!((r1.r as i32 - r2.r as i32).abs() <= 0);
        assert!((r1.g as i32 - r2.g as i32).abs() <= 0);
        assert!((r1.b as i32 - r2.b as i32).abs() <= 0);
        assert!((r1.a as i32 - r2.a as i32).abs() <= 0);
    }

    /// Test that very small image sizes don't panic.
    #[test]
    fn test_small_image_sizes() {
        let mut pixmap = Pixmap::new(1, 1);
        let (n_decimations, kernel, kernel_size) = plan_decimated_blur(2.0);

        // Should not panic
        let result = std::panic::catch_unwind(move || {
            let mut scratch = Pixmap::new(1, 1);
            apply_blur(
                &mut pixmap,
                &mut scratch,
                n_decimations,
                &kernel[..usize::from(kernel_size)],
                EdgeMode::None,
            );
        });

        assert!(result.is_ok());
    }

    /// Test downscale with odd dimensions.
    #[test]
    fn test_downscale_odd_dimensions() {
        let mut pixmap = Pixmap::new(5, 5);
        // Fill with white
        for y in 0..5 {
            for x in 0..5 {
                pixmap.set_pixel(
                    x,
                    y,
                    PremulRgba8 {
                        r: 255,
                        g: 255,
                        b: 255,
                        a: 255,
                    },
                );
            }
        }

        let (new_width, new_height) = downscale(&mut pixmap, 5, 5, EdgeMode::Duplicate);
        // 5 / 2 = 2.5 → ceil = 3
        assert_eq!(new_width, 3);
        assert_eq!(new_height, 3);
    }

    /// Test upscale dimensions.
    #[test]
    fn test_upscale_dimensions() {
        let mut pixmap = Pixmap::new(6, 6);
        let (new_width, new_height) = upscale(&mut pixmap, 3, 3, EdgeMode::Duplicate);
        // 3 * 2 = 6
        assert_eq!(new_width, 6);
        assert_eq!(new_height, 6);
    }

    /// Test that horizontal convolution preserves uniform colors.
    #[test]
    fn test_convolve_x_uniform() {
        let mut src = Pixmap::new(5, 3);
        let mut dst = Pixmap::new(5, 3);
        // Fill with uniform gray
        src.data_mut().fill(PremulRgba8 {
            r: 128,
            g: 128,
            b: 128,
            a: 255,
        });

        let (kernel, kernel_size) = compute_gaussian_kernel(1.0);
        convolve_x(
            &src,
            &mut dst,
            5,
            3,
            &kernel[..usize::from(kernel_size)],
            kernel_size / 2,
            EdgeMode::Duplicate,
        );

        // Uniform input should produce uniform output
        for y in 0..3 {
            for x in 0..5 {
                let pixel = dst.sample(x, y);
                assert_eq!(
                    pixel,
                    PremulRgba8 {
                        r: 128,
                        g: 128,
                        b: 128,
                        a: 255,
                    }
                );
            }
        }
    }

    /// Test that vertical convolution preserves uniform colors.
    #[test]
    fn test_convolve_y_uniform() {
        let mut src = Pixmap::new(3, 5);
        let mut dst = Pixmap::new(3, 5);
        // Fill with uniform gray
        src.data_mut().fill(PremulRgba8 {
            r: 128,
            g: 128,
            b: 128,
            a: 255,
        });

        let (kernel, kernel_size) = compute_gaussian_kernel(1.0);
        convolve_y(
            &src,
            &mut dst,
            3,
            5,
            &kernel[..usize::from(kernel_size)],
            kernel_size / 2,
            EdgeMode::Duplicate,
        );

        // Uniform input should produce uniform output
        for y in 0..5 {
            for x in 0..3 {
                let pixel = dst.sample(x, y);
                assert_eq!(
                    pixel,
                    PremulRgba8 {
                        r: 128,
                        g: 128,
                        b: 128,
                        a: 255,
                    }
                );
            }
        }
    }

    /// Test that convolution with identity kernel is a no-op.
    #[test]
    fn test_convolve_identity_kernel() {
        let mut src = Pixmap::new(3, 3);
        let mut dst = Pixmap::new(3, 3);
        src.set_pixel(
            1,
            1,
            PremulRgba8 {
                r: 255,
                g: 100,
                b: 50,
                a: 200,
            },
        );

        let kernel = [1.0]; // Identity kernel
        convolve_x(&src, &mut dst, 3, 3, &kernel, 0, EdgeMode::None);

        // Should be unchanged
        assert_eq!(dst.sample(1, 1), src.sample(1, 1));
    }

    /// Test that large sigma values are clamped to `MAX_KERNEL_SIZE`.
    #[test]
    fn test_large_sigma_clamped_to_max() {
        let (kernel, kernel_size) = compute_gaussian_kernel(100.0);
        // For σ=100, radius = ceil(300) = 300, size would be 601
        // But it should be clamped to MAX_KERNEL_SIZE
        assert_eq!(usize::from(kernel_size), MAX_KERNEL_SIZE);

        // The clamped kernel should still sum to 1.0 (normalized)
        let sum: f32 = kernel.iter().take(usize::from(kernel_size)).sum();
        assert!((sum - 1.0).abs() < 1e-6);
    }

    /// Test that decimation prevents kernel from exceeding `MAX_KERNEL_SIZE`.
    #[test]
    fn test_decimation_reduces_kernel_size() {
        let (n_decimations, _kernel, kernel_size) = plan_decimated_blur(100.0);
        assert_eq!(kernel_size, 11);
        assert_eq!(n_decimations, 6, "Large sigma should trigger decimation");
    }

    /// Test sampling behavior at exact boundaries for each edge mode.
    #[test]
    fn test_sample_x_at_boundaries() {
        let mut pixmap = Pixmap::new(3, 1);
        pixmap.set_pixel(
            0,
            0,
            PremulRgba8 {
                r: 10,
                g: 0,
                b: 0,
                a: 255,
            },
        );
        pixmap.set_pixel(
            1,
            0,
            PremulRgba8 {
                r: 20,
                g: 0,
                b: 0,
                a: 255,
            },
        );
        pixmap.set_pixel(
            2,
            0,
            PremulRgba8 {
                r: 30,
                g: 0,
                b: 0,
                a: 255,
            },
        );

        // Test left boundary with Duplicate
        let mut p = sample_x(&pixmap, -1, 0, 3, EdgeMode::Duplicate);
        assert_eq!(p.r, 10); // Should clamp to first pixel

        // Test right boundary with Duplicate
        p = sample_x(&pixmap, 3, 0, 3, EdgeMode::Duplicate);
        assert_eq!(p.r, 30); // Should clamp to last pixel

        // Test with None mode (should return transparent black)
        p = sample_x(&pixmap, -1, 0, 3, EdgeMode::None);
        assert_eq!(p.a, 0);
    }

    /// Test sampling behavior at exact boundaries for vertical sampling.
    #[test]
    fn test_sample_y_at_boundaries() {
        let mut pixmap = Pixmap::new(1, 3);
        pixmap.set_pixel(
            0,
            0,
            PremulRgba8 {
                r: 10,
                g: 0,
                b: 0,
                a: 255,
            },
        );
        pixmap.set_pixel(
            0,
            1,
            PremulRgba8 {
                r: 20,
                g: 0,
                b: 0,
                a: 255,
            },
        );
        pixmap.set_pixel(
            0,
            2,
            PremulRgba8 {
                r: 30,
                g: 0,
                b: 0,
                a: 255,
            },
        );

        // Test top boundary with Duplicate
        let mut p = sample_y(&pixmap, 0, -1, 3, EdgeMode::Duplicate);
        assert_eq!(p.r, 10); // Should clamp to first pixel

        // Test bottom boundary with Duplicate
        p = sample_y(&pixmap, 0, 3, 3, EdgeMode::Duplicate);
        assert_eq!(p.r, 30); // Should clamp to last pixel

        // Test with None mode (should return transparent black)
        p = sample_y(&pixmap, 0, -1, 3, EdgeMode::None);
        assert_eq!(p.a, 0);
    }

    /// Test that kernel normalization is precise for various sigma values.
    #[test]
    fn test_kernel_normalization_precision() {
        for sigma in [0.1, 0.5, 1.0, 2.0, 5.0, 10.0] {
            let (kernel, kernel_size) = compute_gaussian_kernel(sigma);
            let sum: f32 = kernel.iter().take(usize::from(kernel_size)).sum();
            assert!(
                (sum - 1.0).abs() < 1e-6,
                "Kernel for σ={} not normalized: sum={}",
                sigma,
                sum
            );
        }
    }

    /// Test that downscale → upscale preserves dimensions.
    #[test]
    fn test_downscale_upscale_roundtrip() {
        let mut pixmap = Pixmap::new(8, 8);
        // Fill with a pattern
        pixmap.data_mut().fill(PremulRgba8 {
            r: 128,
            g: 128,
            b: 128,
            a: 255,
        });

        let (w1, h1) = downscale(&mut pixmap, 8, 8, EdgeMode::Duplicate);
        assert_eq!(w1, 4);
        assert_eq!(h1, 4);

        let (w2, h2) = upscale(&mut pixmap, w1, h1, EdgeMode::Duplicate);
        assert_eq!(w2, 8);
        assert_eq!(h2, 8);
    }

    fn pixmap_from_red(width: u16, height: u16, values: &[&[u8]]) -> Pixmap {
        let mut pixmap = Pixmap::new(width, height);
        for (y, row) in values.iter().enumerate() {
            for (x, &r) in row.iter().enumerate() {
                pixmap.set_pixel(
                    x as u16,
                    y as u16,
                    PremulRgba8 {
                        r,
                        g: 0,
                        b: 0,
                        a: 255,
                    },
                );
            }
        }
        pixmap
    }

    fn assert_red_values<const W: usize, const H: usize>(pixmap: &Pixmap, expected: [[u8; W]; H]) {
        for (y, row) in expected.iter().enumerate() {
            for (x, &want) in row.iter().enumerate() {
                let got = pixmap.sample(x as u16, y as u16).r;

                assert_eq!(
                    got, want,
                    "red mismatch at ({x}, {y}): got {got}, expected {want}"
                );
            }
        }
    }

    #[test]
    fn test_downscale_x_non_uniform() {
        //   (0,0)=0    (1,0)=40   (2,0)=80   (3,0)=120
        //   (0,1)=20   (1,1)=60   (2,1)=100  (3,1)=140
        let mut pixmap = pixmap_from_red(4, 2, &[&[0, 40, 80, 120], &[20, 60, 100, 140]]);

        let dst_width = 4_u16.div_ceil(2);
        downscale_x(&mut pixmap, 4, 2, dst_width, EdgeMode::Duplicate);

        // Row 0: p[-1]=0,p[0]=0,p[1]=40,p[2]=80    → (0+0*3+40*3+80+4)>>3       = 25
        //        p[1]=40,p[2]=80,p[3]=120,p[4]=120  → (40+80*3+120*3+120+4)>>3   = 95
        // Row 1: p[-1]=20,p[0]=20,p[1]=60,p[2]=100  → (20+20*3+60*3+100+4)>>3    = 45
        //        p[1]=60,p[2]=100,p[3]=140,p[4]=140  → (60+100*3+140*3+140+4)>>3   = 115
        assert_red_values(&pixmap, [[25, 95], [45, 115]]);
    }

    #[test]
    fn test_downscale_y_non_uniform() {
        // Use a 2x4 image (output of a prior downscale_x).
        //   (0,0)=25   (1,0)=95
        //   (0,1)=45   (1,1)=115
        //   (0,2)=35   (1,2)=105
        //   (0,3)=55   (1,3)=125
        let mut pixmap = pixmap_from_red(2, 4, &[&[25, 95], &[45, 115], &[35, 105], &[55, 125]]);

        let dst_height = 4_u16.div_ceil(2);
        downscale_y(&mut pixmap, 2, 4, dst_height, EdgeMode::Duplicate);

        // Col 0: p[-1]=25,p[0]=25,p[1]=45,p[2]=35 → (25+25*3+45*3+35+4)>>3   = 34
        //        p[1]=45,p[2]=35,p[3]=55,p[4]=55   → (45+35*3+55*3+55+4)>>3  = 46
        // Col 1: p[-1]=95,p[0]=95,p[1]=115,p[2]=105 → (95+95*3+115*3+105+4)>>3 = 104
        //        p[1]=115,p[2]=105,p[3]=125,p[4]=125 → (115+105*3+125*3+125+4)>>3 = 116
        assert_red_values(&pixmap, [[34, 104], [46, 116]]);
    }

    #[test]
    fn test_upscale_x_non_uniform() {
        let mut pixmap = pixmap_from_red(4, 2, &[&[34, 104], &[46, 116]]);

        upscale_x(&mut pixmap, 2, 2, EdgeMode::Duplicate);

        // Row 0 [34, 104]:
        //   x=0: interp25_75(34,34)=34,  interp75_25(34,104)=(34*3+104+2)>>2 = 52
        //   x=1: interp25_75(34,104)=(34+104*3+2)>>2=87, interp75_25(104,104)=104
        // Row 1 [46, 116]:
        //   x=0: 46, interp75_25(46,116)=(46*3+116+2)>>2=64
        //   x=1: interp25_75(46,116)=(46+116*3+2)>>2=99, 116
        assert_red_values(&pixmap, [[34, 52, 87, 104], [46, 64, 99, 116]]);
    }

    #[test]
    fn test_upscale_y_non_uniform() {
        let mut pixmap = pixmap_from_red(4, 4, &[&[34, 52, 87, 104], &[46, 64, 99, 116]]);

        upscale_y(&mut pixmap, 4, 2, EdgeMode::Duplicate);

        #[rustfmt::skip]
        assert_red_values(&pixmap, [
            [ 34,  52,  87, 104],
            [ 37,  55,  90, 107],
            [ 43,  61,  96, 113],
            [ 46,  64,  99, 116],
        ]);
    }
}