fovea 0.5.1

//! Pyramid construction primitives and strategies.
//!
//! [`pyr_down`] and [`pyr_up`] are the two fundamental resolution-halving /
//! -doubling operations, useful on their own and composed by the
//! [`PyramidMethod`] strategies ([`Gaussian`]) that build a
//! [`Pyramid`](crate::image::Pyramid).
//!
//! Both operations are **named standard operations with a pinned contract**:
//! the smoothing filter is the binomial 5-tap `[1, 4, 6, 4, 1] / 16` per
//! axis (effective σ exactly 1.0) and out-of-bounds access reflects at the
//! edge without duplicating the edge pixel — matching OpenCV's `pyrDown` /
//! `pyrUp` defaults for cross-library comparability. Callers who need a
//! different anti-aliasing filter or border treatment build their own
//! reducer from the parameterized
//! [`gaussian_blur`](crate::transform::gaussian_blur) /
//! [`convolve_separable`](crate::transform::convolve_separable); that path
//! stays fully available.

use crate::border::Mirror;
use crate::error::Error;
use crate::image::{
    Dyadic, Image, ImageView, ImageViewMut, LevelChain, OriginOffset, PlacedImage, PlacedPyramid,
    Pyramid, PyramidLevel, RasterImage, ScaledImage, ScaledPyramid, SeparableKernel,
};
use crate::pixel::{FromLinear, LinearPixel, LinearSpace, ZeroablePixel};
use crate::transform::convolve_separable::convolve_separable;
use crate::{PixelDistance, Sigma, Size};

/// The `pyr_up` interpolation kernel: the binomial `[1, 4, 6, 4, 1] / 8`
/// per axis — the `pyr_down` kernel with weights ×2 per axis (×4 combined),
/// compensating for the zero-inserted samples so brightness is preserved.
const PYR_UP_WEIGHTS: [f32; 5] = [0.125, 0.5, 0.75, 0.5, 0.125];

/// The same kernel normalized to sum 1, for an axis that received **no**
/// zero-inserted samples: a length-1 axis whose target is also 1 (the only
/// case where `target == source` passes validation). The ×2 compensation
/// above assumes half of each axis's samples are inserted zeros; with none,
/// it would double the brightness of every pixel along that axis.
const PYR_UP_WEIGHTS_UNDOUBLED: [f32; 5] = [0.0625, 0.25, 0.375, 0.25, 0.0625];

// ─── pyr_down / pyr_up ──────────────────────────────────────────────────────

/// Blurs and decimates the image by a factor of 2.
///
/// Applies the binomial 5×5 Gaussian (`[1, 4, 6, 4, 1] / 16` per axis:
/// [`SeparableKernel::gaussian_5`](crate::image::SeparableKernel::gaussian_5),
/// the same weights [`gaussian_blur_5x5`](crate::transform::gaussian_blur_5x5)
/// applies, effective σ exactly 1.0) followed by 2× downsampling that keeps the
/// even-indexed samples (pixels 0, 2, 4, …).
///
/// The output dimensions are `((width + 1) / 2, (height + 1) / 2)`
/// (ceiling division): an `n`-wide row has `ceil(n / 2)` even samples, so
/// the last column/row of an odd-sized image stays represented, and the
/// sizes match OpenCV's `pyrDown`.
///
/// The kernel and the border treatment (reflection without edge
/// duplication, OpenCV's `BORDER_REFLECT_101`) are part of this function's
/// contract — a named standard operation, not a moving definition. For a
/// different anti-aliasing filter, compose your own reducer from
/// [`gaussian_blur`](crate::transform::gaussian_blur).
///
/// Because the smoothing blends neighboring samples, the pixel type must
/// live in a linear space ([`LinearSpace`]) — linearize sRGB first.
///
/// # Example
///
/// ```
/// use fovea::Size;
/// use fovea::image::{Image, ImageView};
/// use fovea::pixel::MonoF32;
/// use fovea::transform::pyr_down;
///
/// let src = Image::fill(9, 6, MonoF32::new(0.5));
/// let half: Image<MonoF32> = pyr_down(&src);
///
/// // Ceiling division: 9 → 5, 6 → 3.
/// assert_eq!(half.size(), Size::new(5, 3));
/// // A flat image stays flat — the kernel preserves brightness.
/// assert!((half.pixel_at(2, 1).0 - 0.5).abs() < 1e-6);
/// ```
#[must_use]
pub fn pyr_down<I, P, Acc>(image: &I) -> Image<P>
where
    I: RasterImage<Pixel = P>,
    P: LinearPixel<f32, Accumulator = Acc> + LinearSpace + ZeroablePixel + FromLinear<Acc>,
    Acc: Copy
        + Default
        + ZeroablePixel
        + LinearPixel<f32, Accumulator = Acc>
        + std::ops::Add<Output = Acc>,
{
    let blurred: Image<P> = convolve_separable(image, &SeparableKernel::gaussian_5(), &Mirror);
    let out_width = image.width().div_ceil(2);
    let out_height = image.height().div_ceil(2);
    Image::generate(out_width, out_height, |x, y| blurred.pixel_at(2 * x, 2 * y))
}

/// Upsamples the image by a factor of 2 to an explicit target size.
///
/// Inserts zero rows/columns (input pixel `(x, y)` lands at output
/// `(2x, 2y)`) and interpolates the missing values with the same binomial
/// kernel as [`pyr_down`], weights ×4 so brightness is preserved after
/// zero-insertion. The border treatment matches `pyr_down` (reflection
/// without edge duplication).
///
/// The explicit `target` is deliberate: [`pyr_down`] maps both an odd and
/// an even dimension onto the same output size, so an upsampler that always
/// doubles would reconstruct the wrong size for odd parents. Naming the
/// parent size removes the ambiguity — `target` must be a size whose
/// `pyr_down` result is this image's size.
///
/// Because the interpolation blends neighboring samples, the pixel type
/// must live in a linear space ([`LinearSpace`]) — linearize sRGB first.
///
/// # Errors
///
/// Returns [`Error::InvalidPyrUpTarget`] if `target.width ∉ {2·w − 1, 2·w}`
/// or `target.height ∉ {2·h − 1, 2·h}` — the caller named a size this
/// image cannot be the `pyr_down` of. Because the valid target is a
/// relation between two runtime sizes (often originating from camera or
/// file dimensions), this is a recoverable error, not a panic.
///
/// # Example
///
/// ```
/// use fovea::Size;
/// use fovea::image::{Image, ImageView};
/// use fovea::pixel::MonoF32;
/// use fovea::transform::{pyr_down, pyr_up};
///
/// let src = Image::fill(9, 7, MonoF32::new(0.25));
/// let half: Image<MonoF32> = pyr_down(&src);
/// assert_eq!(half.size(), Size::new(5, 4));
///
/// // The explicit target restores the odd parent size exactly.
/// let restored: Image<MonoF32> = pyr_up(&half, src.size())?;
/// assert_eq!(restored.size(), Size::new(9, 7));
/// assert!((restored.pixel_at(4, 3).0 - 0.25).abs() < 1e-6);
/// # Ok::<(), fovea::Error>(())
/// ```
pub fn pyr_up<I, P, Acc>(image: &I, target: Size) -> Result<Image<P>, Error>
where
    I: RasterImage<Pixel = P>,
    P: LinearPixel<f32, Accumulator = Acc> + LinearSpace + ZeroablePixel + FromLinear<Acc>,
    Acc: Copy
        + Default
        + ZeroablePixel
        + LinearPixel<f32, Accumulator = Acc>
        + std::ops::Add<Output = Acc>,
{
    /// `target` doubles `dim` (even parent) or doubles it minus one (odd
    /// parent). Checked: a zero-area image can carry a dimension past
    /// `usize::MAX / 2` and any `Size` can be named as the target, and both
    /// must be rejected rather than wrapped over (or aborted on, in debug).
    fn doubles(target: usize, dim: usize) -> bool {
        match dim.checked_mul(2) {
            Some(two) => target == two || target.checked_add(1) == Some(two),
            None => false,
        }
    }
    let (w, h) = (image.width(), image.height());
    if !doubles(target.width, w) || !doubles(target.height, h) {
        return Err(Error::InvalidPyrUpTarget {
            source: image.size(),
            target,
        });
    }

    // Zero-insertion: every input sample keeps its even-even position; the
    // in-between positions start at zero and are filled by the smoothing.
    let mut upsampled = Image::<P>::zero(target.width, target.height);
    for y in 0..h {
        let row = image.row(y);
        for (x, &pixel) in row.iter().enumerate() {
            *upsampled.pixel_at_mut(2 * x, 2 * y) = pixel;
        }
    }

    // Per axis: the doubled weights compensate for the interleaved zeros;
    // an axis that stayed at length 1 has none, so it takes the normalized
    // weights instead (a flat field must stay flat either way).
    let h_weights = if target.width == w {
        PYR_UP_WEIGHTS_UNDOUBLED
    } else {
        PYR_UP_WEIGHTS
    };
    let v_weights = if target.height == h {
        PYR_UP_WEIGHTS_UNDOUBLED
    } else {
        PYR_UP_WEIGHTS
    };
    let kernel = SeparableKernel::new(h_weights, v_weights);
    Ok(convolve_separable(&upsampled, &kernel, &Mirror))
}

// ─── expand ─────────────────────────────────────────────────────────────────

impl<C: Pyramid> Dyadic<C> {
    /// Lifts the level at `child` back to the size of its parent.
    ///
    /// This is [`pyr_up`] with the target taken from the neighbouring
    /// level instead of from the caller. The free function keeps existing
    /// and stays the right call for an image of unknown provenance: a
    /// 51-wide image could be the reduction of 101 or of 102, so its
    /// parent size is genuinely input. Inside a [`Dyadic`] pyramid it is
    /// not input, it is already stored, and the ambiguity `pyr_up`'s
    /// `target` argument exists to resolve cannot arise.
    ///
    /// The result is a plain [`Image<P>`], not a level: a raised child is
    /// an approximation of its parent, not the parent, and giving it the
    /// parent's geometry would be a claim nothing backs.
    ///
    /// Returns `None` exactly when there is no parent to lift to: for
    /// `child == 0`, which is the finest level, and for any `child` at or
    /// past [`depth`](Pyramid::depth). That is the same question
    /// [`get`](Pyramid::get) answers, and it is the only one left: the
    /// halving guarantee makes every index that does have a parent
    /// succeed.
    ///
    /// # Example
    ///
    /// ```
    /// use fovea::image::{Image, ImageView, Pyramid};
    /// use fovea::pixel::MonoF32;
    /// use fovea::transform::{Gaussian, PyramidMethod};
    ///
    /// let img = Image::fill(20, 12, MonoF32::new(0.5));
    /// let pyramid = Gaussian.build(&img, 3);
    ///
    /// let raised: Image<MonoF32> = pyramid.expand(2).expect("level 2 has a parent");
    /// assert_eq!(raised.size(), pyramid.level(1).size());
    ///
    /// assert!(pyramid.expand(0).is_none()); // the finest level has no parent
    /// assert!(pyramid.expand(3).is_none()); // no such level
    /// ```
    pub fn expand<P, Acc>(&self, child: usize) -> Option<Image<P>>
    where
        C::Level: PyramidLevel<Pixel = P>,
        P: LinearPixel<f32, Accumulator = Acc> + LinearSpace + ZeroablePixel + FromLinear<Acc>,
        Acc: Copy
            + Default
            + ZeroablePixel
            + LinearPixel<f32, Accumulator = Acc>
            + std::ops::Add<Output = Acc>,
    {
        let source = self.get(child)?;
        let parent = self.get(child.checked_sub(1)?)?;
        Some(
            pyr_up(source.as_image(), parent.as_image().size())
                .expect("Dyadic guarantees the parent is a valid pyr_up target"),
        )
    }
}

// ─── PyramidMethod strategy ─────────────────────────────────────────────────

/// Strategy trait for pyramid construction.
///
/// A `PyramidMethod` produces a [`Pyramid`] and is then discarded — the
/// result does not remember how it was built, following the same pattern as
/// [`ResizeMethod`](crate::transform::ResizeMethod) and
/// [`ConvertPixel`](crate::transform::ConvertPixel). Implement this trait
/// for custom decomposition schemes; assemble the result with
/// [`LevelChain::try_from_levels`](crate::image::LevelChain::try_from_levels).
///
/// The container is named by [`Output`](Self::Output) rather than fixed, so
/// a method that halves can return a [`Dyadic`](crate::image::Dyadic) one
/// and hand its callers [`expand`](crate::image::Dyadic::expand), while a
/// method that does not returns the plain chain. The level type is reached
/// through it as `<Self::Output as Pyramid>::Level`.
///
/// # Example
///
/// ```
/// use fovea::image::{Image, ImageView, Pyramid};
/// use fovea::pixel::MonoF32;
/// use fovea::transform::{Gaussian, PyramidMethod};
///
/// let img = Image::fill(32, 32, MonoF32::new(1.0));
/// let pyramid = Gaussian.build(&img, 4);
///
/// assert_eq!(pyramid.depth(), 4);
/// assert_eq!(pyramid.coarsest().size().width, 4);
/// ```
pub trait PyramidMethod<P: Copy> {
    /// The pyramid this method produces, container and level type at once.
    type Output: Pyramid;

    /// Builds a pyramid from the given image.
    ///
    /// The input is any [`RasterImage`] view — an owned [`Image`], a
    /// borrowed buffer, or an ROI — so building a pyramid of a camera
    /// frame's sub-rectangle needs no intermediate copy beyond the base
    /// level the pyramid owns anyway.
    ///
    /// `max_depth` is an **upper bound, not a promise**. If the image is
    /// too small to support the requested depth, `build` clamps at the
    /// method's minimum usable level size — it never panics, never errors,
    /// and never mutates the caller's parameters. The resolved depth is
    /// whatever [`Pyramid::depth`] reports afterwards. The result always
    /// contains at least one level.
    fn build<I>(&self, image: &I, max_depth: usize) -> Self::Output
    where
        I: RasterImage<Pixel = P>;
}

/// Gaussian pyramid construction: repeated [`pyr_down`].
///
/// Level 0 is a copy of the input image; each further level is the
/// [`pyr_down`] of the previous one, halving the resolution (ceiling
/// division) with the pinned binomial smoothing.
///
/// Builds a [`PlacedPyramid<P>`](crate::image::PlacedPyramid): every level
/// carries the sampling geometry the builder computed anyway, a
/// `pixel_distance` of `2^k` and an unshifted origin, so a caller never
/// hand-writes `x * 2^level` again. It carries **no** σ, and that is the
/// point: this method is told nothing about how sharp its input was, so any
/// σ it attached would be a claim it cannot back. Name the assumption with
/// [`assuming_input_sigma`](Self::assuming_input_sigma) and the levels come
/// back as [`ScaledImage`](crate::image::ScaledImage) instead.
///
/// Because every step halves, the result is a
/// [`Dyadic`](crate::image::Dyadic) pyramid and therefore carries
/// [`expand`](crate::image::Dyadic::expand): the way back up needs no
/// target size. The halving is established by construction here, so nothing
/// re-checks it.
///
/// The build stops early once a level cannot shrink further (1×1), so the
/// resolved depth may be smaller than requested; a `max_depth` of 0 is
/// treated as 1, because a pyramid always contains at least its base level.
///
/// # Example
///
/// ```
/// use fovea::Size;
/// use fovea::image::{Decimated, Image, ImageView, PlacedPyramid, Pyramid};
/// use fovea::pixel::MonoF32;
/// use fovea::transform::{Gaussian, PyramidMethod};
///
/// let img = Image::fill(20, 12, MonoF32::new(0.5));
/// let pyramid: PlacedPyramid<MonoF32> = Gaussian.build(&img, 3);
///
/// let sizes: Vec<Size> = pyramid.iter().map(|l| l.size()).collect();
/// assert_eq!(sizes, [Size::new(20, 12), Size::new(10, 6), Size::new(5, 3)]);
///
/// // The geometry is answered, not remembered by the caller.
/// assert_eq!(pyramid.level(2).pixel_distance().get(), 4.0);
/// ```
#[derive(Debug, Clone, Copy, Default, PartialEq, Eq)]
pub struct Gaussian;

impl Gaussian {
    /// Names the σ the input image is assumed to already carry, so the
    /// levels can state an absolute σ of their own.
    ///
    /// A pyramid's σ ladder is only defined relative to the blur its input
    /// already had. A photograph off a sensor is never perfectly sharp;
    /// Lowe's SIFT assumes σ = 0.5 for exactly this reason. Since the
    /// library cannot know that number, it refuses to invent one: without
    /// this call the levels carry no σ at all.
    ///
    /// The ladder is `σ_k² = σ_in² + (4^k − 1)/3`, in base-image pixels.
    /// With `σ_in = 0.5`, levels 1 to 3 come out at 1.118, 2.291 and 4.610,
    /// not at 1, 2 and 4: a caller extrapolating from `pixel_distance` is
    /// 13% wrong by level 3, silently. That gap is the reason this method
    /// exists.
    ///
    /// **The ladder is nominal, and it is worth being precise about which
    /// part.** The arithmetic is exact: [`pyr_down`]'s pinned binomial 5-tap
    /// has a kernel variance of exactly 1, convolution adds variances
    /// exactly, and the composite kernel from the base image to level `k`
    /// therefore has variance `(4^k − 1)/3` away from the borders. What is
    /// not exact is calling the result a Gaussian. The 5-tap is not one, and
    /// repeated convolution converges to one without ever being one; the
    /// `σ_in` term likewise assumes the input's own blur is Gaussian. A
    /// level reports a **nominal effective** scale under Gaussian
    /// composition, the conventional number every comparable library uses,
    /// not a measurement.
    ///
    /// # Example
    ///
    /// ```
    /// use fovea::image::{Image, ScaleLevel, ScaledPyramid, Pyramid};
    /// use fovea::pixel::MonoF32;
    /// use fovea::sigma;
    /// use fovea::transform::{Gaussian, PyramidMethod};
    ///
    /// let img = Image::fill(32, 32, MonoF32::new(0.5));
    /// let pyramid: ScaledPyramid<MonoF32> =
    ///     Gaussian.assuming_input_sigma(sigma!(0.5)).build(&img, 4);
    ///
    /// assert_eq!(pyramid.level(0).sigma().get(), 0.5);
    /// assert!((pyramid.level(1).sigma().get() - 1.118).abs() < 1e-3);
    /// ```
    #[must_use]
    pub fn assuming_input_sigma(self, input_sigma: Sigma) -> ScaledGaussian {
        ScaledGaussian { input_sigma }
    }
}

impl<P, Acc> PyramidMethod<P> for Gaussian
where
    P: LinearPixel<f32, Accumulator = Acc> + LinearSpace + ZeroablePixel + FromLinear<Acc>,
    Acc: Copy
        + Default
        + ZeroablePixel
        + LinearPixel<f32, Accumulator = Acc>
        + std::ops::Add<Output = Acc>,
{
    type Output = PlacedPyramid<P>;

    fn build<I>(&self, image: &I, max_depth: usize) -> Self::Output
    where
        I: RasterImage<Pixel = P>,
    {
        let levels = gaussian_levels(image, max_depth)
            .into_iter()
            .enumerate()
            .map(|(index, image)| {
                PlacedImage::new(image, level_pixel_distance(index), OriginOffset::ZERO)
            })
            .collect();
        finish(levels)
    }
}

/// Gaussian pyramid construction with a named input σ: [`Gaussian`] plus the
/// one fact it cannot derive.
///
/// Reached through
/// [`Gaussian::assuming_input_sigma`](Gaussian::assuming_input_sigma), which
/// is where the ladder and its caveats are documented. Builds a
/// [`ScaledPyramid<P>`](crate::image::ScaledPyramid): the same levels as
/// [`Gaussian`], each additionally carrying its finished absolute σ.
///
/// The levels store the **result**, not the schedule. Asked afterwards, a
/// level says what σ it has, not which assumption produced it; σ_in is
/// consumed here and retained nowhere. That keeps the level type usable for
/// a level some other reducer produced, and it is the crate's standing
/// pattern: results do not remember their strategy.
///
/// # Example
///
/// ```
/// use fovea::image::{Decimated, Image, ScaleLevel, Pyramid};
/// use fovea::pixel::MonoF32;
/// use fovea::sigma;
/// use fovea::transform::{Gaussian, PyramidMethod};
///
/// let img = Image::fill(16, 16, MonoF32::new(1.0));
/// let pyramid = Gaussian.assuming_input_sigma(sigma!(0.5)).build(&img, 3);
///
/// let level = pyramid.level(1);
/// assert_eq!(level.pixel_distance().get(), 2.0);
/// assert!((level.sigma().get() - 1.118).abs() < 1e-3);
/// ```
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct ScaledGaussian {
    input_sigma: Sigma,
}

impl ScaledGaussian {
    /// Returns the assumed input σ this method was given.
    #[must_use]
    pub fn input_sigma(&self) -> Sigma {
        self.input_sigma
    }
}

impl<P, Acc> PyramidMethod<P> for ScaledGaussian
where
    P: LinearPixel<f32, Accumulator = Acc> + LinearSpace + ZeroablePixel + FromLinear<Acc>,
    Acc: Copy
        + Default
        + ZeroablePixel
        + LinearPixel<f32, Accumulator = Acc>
        + std::ops::Add<Output = Acc>,
{
    type Output = ScaledPyramid<P>;

    fn build<I>(&self, image: &I, max_depth: usize) -> Self::Output
    where
        I: RasterImage<Pixel = P>,
    {
        // Each `pyr_down` adds the variance of a σ = 1 blur measured in its
        // *parent's* pixels, which is 4^(k-1) base pixels squared. Summed,
        // that is the (4^k - 1)/3 of the ladder, accumulated here rather
        // than re-derived per level so the exponent cannot overflow.
        let mut variance = f64::from(self.input_sigma.get()).powi(2);
        let mut added = 1.0_f64;
        let levels = gaussian_levels(image, max_depth)
            .into_iter()
            .enumerate()
            .map(|(index, image)| {
                if index > 0 {
                    variance += added;
                    added *= 4.0;
                }
                let sigma = Sigma::try_new(variance.sqrt() as f32)
                    .expect("a positive input sigma stays positive and finite along the ladder");
                ScaledImage::new(
                    image,
                    level_pixel_distance(index),
                    OriginOffset::ZERO,
                    sigma,
                )
            })
            .collect();
        finish(levels)
    }
}

/// The base-pixel spacing of level `index`: `2^index`, since every step
/// keeps the even samples of the one before it.
fn level_pixel_distance(index: usize) -> PixelDistance {
    let spacing = (1_u64 << index.min(63)) as f64;
    PixelDistance::try_new(spacing).expect("a power of two is finite and strictly positive")
}

/// Assembles built levels into the dyadic chain both builders return.
///
/// Every level is the `pyr_down` of its predecessor, so the halving relation
/// holds by construction and needs no second check here.
fn finish<L: PyramidLevel>(levels: Vec<L>) -> Dyadic<LevelChain<L>> {
    let chain = LevelChain::try_from_levels(levels)
        .expect("Gaussian construction produces non-empty, non-growing levels");
    Dyadic::new_unchecked(chain)
}

/// The shared level cascade: level 0 is an owned copy of the input, each
/// further level the [`pyr_down`] of the one before it.
fn gaussian_levels<I, P, Acc>(image: &I, max_depth: usize) -> Vec<Image<P>>
where
    I: RasterImage<Pixel = P>,
    P: LinearPixel<f32, Accumulator = Acc> + LinearSpace + ZeroablePixel + FromLinear<Acc>,
    Acc: Copy
        + Default
        + ZeroablePixel
        + LinearPixel<f32, Accumulator = Acc>
        + std::ops::Add<Output = Acc>,
{
    let resolved = max_depth.max(1);
    // Level 0 is an owned copy of whatever view came in, row by row.
    let base = {
        let mut data = Vec::with_capacity(image.width() * image.height());
        for y in 0..image.height() {
            data.extend_from_slice(image.row(y));
        }
        Image::from_vec(image.width(), image.height(), data)
            .expect("rows fill width * height exactly")
    };
    let mut levels = vec![base];
    while levels.len() < resolved {
        let prev = levels.last().expect("levels start non-empty");
        let Size { width, height } = prev.size();
        // Minimum usable level size: a level that cannot shrink
        // further (or has no pixels at all) ends the chain.
        if width <= 1 && height <= 1 || width == 0 || height == 0 {
            break;
        }
        let next = pyr_down(prev);
        levels.push(next);
    }
    levels
}

// ─── Tests ──────────────────────────────────────────────────────────────────

#[cfg(test)]
mod tests {
    use super::*;
    use crate::CoordinateF64;
    use crate::image::{Decimated, PyramidLevel, ScaleLevel};
    use crate::pixel::{Mono8, MonoF32};
    use crate::{pixel_distance, sigma};

    // ── pyr_down: size contract ─────────────────────────────────────────

    #[test]
    fn pyr_down_even_dimensions_halve() {
        let src = Image::fill(8, 6, MonoF32::new(0.0));
        let out: Image<MonoF32> = pyr_down(&src);
        assert_eq!(out.size(), Size::new(4, 3));
    }

    #[test]
    fn pyr_down_odd_dimensions_use_ceiling_division() {
        let src = Image::fill(7, 5, MonoF32::new(0.0));
        let out: Image<MonoF32> = pyr_down(&src);
        assert_eq!(out.size(), Size::new(4, 3));
    }

    #[test]
    fn pyr_down_one_pixel_image_stays_one_pixel() {
        let src = Image::fill(1, 1, MonoF32::new(0.3));
        let out: Image<MonoF32> = pyr_down(&src);
        assert_eq!(out.size(), Size::new(1, 1));
        assert!((out.pixel_at(0, 0).0 - 0.3).abs() < 1e-6);
    }

    #[test]
    fn pyr_down_single_row_and_column() {
        let row = Image::fill(9, 1, MonoF32::new(0.5));
        let out: Image<MonoF32> = pyr_down(&row);
        assert_eq!(out.size(), Size::new(5, 1));

        let col = Image::fill(1, 8, MonoF32::new(0.5));
        let out: Image<MonoF32> = pyr_down(&col);
        assert_eq!(out.size(), Size::new(1, 4));
    }

    // ── pyr_down: value contract ────────────────────────────────────────

    #[test]
    fn pyr_down_flat_image_preserves_brightness() {
        let src = Image::fill(10, 10, MonoF32::new(0.7));
        let out: Image<MonoF32> = pyr_down(&src);
        for y in 0..out.height() {
            for x in 0..out.width() {
                assert!(
                    (out.pixel_at(x, y).0 - 0.7).abs() < 1e-6,
                    "flat value drifted at ({x}, {y}): {}",
                    out.pixel_at(x, y).0
                );
            }
        }
    }

    #[test]
    fn pyr_down_flat_mono8_preserves_brightness() {
        let src = Image::fill(12, 8, Mono8::new(100));
        let out: Image<Mono8> = pyr_down(&src);
        for y in 0..out.height() {
            for x in 0..out.width() {
                assert_eq!(out.pixel_at(x, y), Mono8::new(100));
            }
        }
    }

    #[test]
    fn pyr_down_keeps_even_samples() {
        // The symmetric binomial kernel preserves a linear ramp in the
        // interior, so out(x) must equal ramp(2x) — the even-sample
        // convention (origin offset (0, 0)), not 2x + 0.5 (area average).
        let src = Image::generate(16, 16, |x, _| MonoF32::new(x as f32));
        let out: Image<MonoF32> = pyr_down(&src);
        for y in 2..out.height() - 2 {
            for x in 2..out.width() - 2 {
                assert!(
                    (out.pixel_at(x, y).0 - 2.0 * x as f32).abs() < 1e-4,
                    "expected even sample 2·{x} at ({x}, {y}), got {}",
                    out.pixel_at(x, y).0
                );
            }
        }
    }

    #[test]
    fn pyr_down_impulse_center_weight() {
        // A unit impulse picks out the kernel's center weight: the 2D
        // binomial center is (6/16)² = 0.140625.
        let src = Image::generate(9, 9, |x, y| {
            if x == 4 && y == 4 {
                MonoF32::new(1.0)
            } else {
                MonoF32::new(0.0)
            }
        });
        let out: Image<MonoF32> = pyr_down(&src);
        assert!((out.pixel_at(2, 2).0 - 0.140625).abs() < 1e-6);
    }

    // ── pyr_up: size contract ───────────────────────────────────────────

    #[test]
    fn pyr_up_accepts_both_valid_widths() {
        let src = Image::fill(4, 4, MonoF32::new(0.5));
        let a: Image<MonoF32> = pyr_up(&src, Size::new(8, 8)).unwrap();
        assert_eq!(a.size(), Size::new(8, 8));
        let b: Image<MonoF32> = pyr_up(&src, Size::new(7, 7)).unwrap();
        assert_eq!(b.size(), Size::new(7, 7));
    }

    #[test]
    fn pyr_up_handles_one_pixel_sources() {
        // 1x1 and 1xN sources: the zero-insertion and reflection paths must
        // survive the degenerate shapes, and a flat field must stay inside
        // its own range.
        let dot = Image::fill(1, 1, MonoF32::new(0.5));
        let up: Image<MonoF32> = pyr_up(&dot, Size::new(2, 2)).unwrap();
        assert_eq!(up.size(), Size::new(2, 2));
        let same: Image<MonoF32> = pyr_up(&dot, Size::new(1, 1)).unwrap();
        assert_eq!(same.size(), Size::new(1, 1));

        let bar = Image::fill(1, 4, MonoF32::new(0.5));
        let up: Image<MonoF32> = pyr_up(&bar, Size::new(2, 8)).unwrap();
        assert_eq!(up.size(), Size::new(2, 8));
        let odd: Image<MonoF32> = pyr_up(&bar, Size::new(1, 7)).unwrap();
        assert_eq!(odd.size(), Size::new(1, 7));
        for y in 0..odd.height() {
            let v = odd.pixel_at(0, y).0;
            assert!((0.0..=0.5 + 1e-6).contains(&v), "({y}) = {v}");
        }
    }

    #[test]
    fn pyr_up_rejects_extreme_sizes_without_overflowing() {
        // A zero-area image can legally carry a dimension past
        // usize::MAX / 2, and any Size can be named as the target. Both
        // used to wrap in the validation arithmetic (aborting in debug
        // builds); they must simply be rejected.
        let wide: Image<MonoF32> = Image::zero(usize::MAX, 0);
        let result: Result<Image<MonoF32>, Error> = pyr_up(&wide, Size::new(4, 4));
        assert!(result.is_err());

        let src = Image::fill(4, 4, MonoF32::new(0.5));
        let result: Result<Image<MonoF32>, Error> = pyr_up(&src, Size::new(usize::MAX, usize::MAX));
        assert!(result.is_err());
    }

    #[test]
    fn pyr_up_rejects_invalid_targets() {
        // Too-large width, too-small width, invalid height: each must
        // report the rejected target and the source size.
        let src = Image::fill(4, 4, MonoF32::new(0.5));
        for target in [Size::new(9, 8), Size::new(6, 8), Size::new(8, 10)] {
            let result: Result<Image<MonoF32>, Error> = pyr_up(&src, target);
            assert_eq!(
                result.unwrap_err(),
                Error::InvalidPyrUpTarget {
                    source: Size::new(4, 4),
                    target,
                },
                "target {target:?} must be rejected"
            );
        }
    }

    #[test]
    fn pyr_up_round_trips_odd_sizes() {
        // The reason target is explicit: odd parents reconstruct exactly.
        let src = Image::fill(9, 7, MonoF32::new(0.25));
        let half: Image<MonoF32> = pyr_down(&src);
        assert_eq!(half.size(), Size::new(5, 4));
        let restored: Image<MonoF32> = pyr_up(&half, src.size()).unwrap();
        assert_eq!(restored.size(), src.size());
    }

    // ── pyr_up: value contract ──────────────────────────────────────────

    #[test]
    fn pyr_up_flat_image_preserves_brightness() {
        // The ×4 weight compensation must hold at every position parity
        // (even/odd × even/odd) and at the borders, for both target
        // parities.
        let src = Image::fill(5, 4, MonoF32::new(0.6));
        for target in [Size::new(10, 8), Size::new(9, 7)] {
            let out: Image<MonoF32> = pyr_up(&src, target).unwrap();
            for y in 0..out.height() {
                for x in 0..out.width() {
                    assert!(
                        (out.pixel_at(x, y).0 - 0.6).abs() < 1e-6,
                        "flat value drifted at ({x}, {y}) for target {target:?}: {}",
                        out.pixel_at(x, y).0
                    );
                }
            }
        }
    }

    #[test]
    fn pyr_up_impulse_spreads_interpolation_weights() {
        // Input sample (1, 1) lands at output (2, 2); the separable
        // interpolation weights around it are the per-axis
        // [1, 4, 6, 4, 1] / 8 taps that hit non-zero samples.
        let src = Image::generate(3, 3, |x, y| {
            if x == 1 && y == 1 {
                MonoF32::new(1.0)
            } else {
                MonoF32::new(0.0)
            }
        });
        let out: Image<MonoF32> = pyr_up(&src, Size::new(6, 6)).unwrap();
        // Even-even: center weight 0.75².
        assert!((out.pixel_at(2, 2).0 - 0.5625).abs() < 1e-6);
        // Odd-even: 0.5 · 0.75.
        assert!((out.pixel_at(3, 2).0 - 0.375).abs() < 1e-6);
        // Odd-odd: 0.5 · 0.5.
        assert!((out.pixel_at(3, 3).0 - 0.25).abs() < 1e-6);
    }

    #[test]
    fn pyr_up_mono8_flat() {
        let src = Image::fill(6, 6, Mono8::new(80));
        let out: Image<Mono8> = pyr_up(&src, Size::new(12, 12)).unwrap();
        for y in 0..out.height() {
            for x in 0..out.width() {
                assert_eq!(out.pixel_at(x, y), Mono8::new(80));
            }
        }
    }

    // ── Gaussian PyramidMethod ──────────────────────────────────────────

    #[test]
    fn gaussian_build_level_zero_is_the_input() {
        let src = Image::generate(8, 8, |x, y| MonoF32::new((x + y) as f32));
        let pyramid = Gaussian.build(&src, 3);
        let level0 = pyramid.finest();
        for y in 0..src.height() {
            for x in 0..src.width() {
                assert_eq!(level0.pixel_at(x, y), src.pixel_at(x, y));
            }
        }
    }

    #[test]
    fn gaussian_build_halves_each_level() {
        let src = Image::fill(20, 12, MonoF32::new(0.5));
        let pyramid = Gaussian.build(&src, 3);
        let sizes: Vec<Size> = pyramid.iter().map(|l| l.size()).collect();
        assert_eq!(
            sizes,
            [Size::new(20, 12), Size::new(10, 6), Size::new(5, 3)]
        );
    }

    #[test]
    fn gaussian_build_clamps_depth_on_small_images() {
        // Resolved ≠ requested: 4×4 supports 4, 2, 1 — then 1×1 stops.
        let src = Image::fill(4, 4, MonoF32::new(0.5));
        let pyramid = Gaussian.build(&src, 100);
        assert_eq!(pyramid.depth(), 3);
        assert_eq!(pyramid.coarsest().size(), Size::new(1, 1));
    }

    #[test]
    fn gaussian_build_max_depth_zero_yields_base_level() {
        let src = Image::fill(8, 8, MonoF32::new(0.5));
        let pyramid = Gaussian.build(&src, 0);
        assert_eq!(pyramid.depth(), 1);
        assert_eq!(pyramid.finest().size(), Size::new(8, 8));
    }

    #[test]
    fn gaussian_build_respects_requested_depth() {
        let src = Image::fill(64, 64, MonoF32::new(0.5));
        let pyramid = Gaussian.build(&src, 3);
        assert_eq!(pyramid.depth(), 3);
    }

    #[test]
    fn gaussian_build_flat_stays_flat_at_every_level() {
        let src = Image::fill(16, 16, MonoF32::new(0.4));
        let pyramid = Gaussian.build(&src, 5);
        for (i, level) in pyramid.iter().enumerate() {
            for y in 0..level.height() {
                for x in 0..level.width() {
                    assert!(
                        (level.pixel_at(x, y).0 - 0.4).abs() < 1e-5,
                        "level {i} drifted at ({x}, {y})"
                    );
                }
            }
        }
    }

    #[test]
    fn gaussian_build_mono8() {
        let src = Image::fill(16, 12, Mono8::new(200));
        let pyramid = Gaussian.build(&src, 3);
        assert_eq!(pyramid.depth(), 3);
        assert_eq!(pyramid.coarsest().size(), Size::new(4, 3));
        assert_eq!(pyramid.coarsest().pixel_at(0, 0), Mono8::new(200));
    }

    // ── expand ──────────────────────────────────────────────────────────

    #[test]
    fn expand_recovers_every_parent_size() {
        // Odd on both axes, so the doubling `pyr_up` cannot guess is
        // exercised at every step rather than only at the last.
        let src = Image::fill(21, 13, MonoF32::new(0.5));
        let pyramid = Gaussian.build(&src, 4);
        for child in 1..pyramid.depth() {
            let raised: Image<MonoF32> = pyramid.expand(child).expect("child has a parent");
            assert_eq!(raised.size(), pyramid.level(child - 1).size());
        }
    }

    #[test]
    fn expand_matches_pyr_up_with_the_parent_size() {
        let src = Image::fill(20, 12, MonoF32::new(0.25));
        let pyramid = Gaussian.build(&src, 3);
        let by_hand: Image<MonoF32> =
            pyr_up(pyramid.level(2), pyramid.level(1).size()).expect("valid target");
        let raised: Image<MonoF32> = pyramid.expand(2).expect("level 2 has a parent");
        assert_eq!(raised.size(), by_hand.size());
        for y in 0..raised.height() {
            for x in 0..raised.width() {
                assert_eq!(raised.pixel_at(x, y), by_hand.pixel_at(x, y));
            }
        }
    }

    #[test]
    fn expand_has_no_parent_for_the_finest_level() {
        let src = Image::fill(16, 16, MonoF32::new(0.5));
        let pyramid = Gaussian.build(&src, 3);
        assert!(pyramid.expand::<MonoF32, _>(0).is_none());
    }

    #[test]
    fn expand_has_no_level_past_the_depth() {
        let src = Image::fill(16, 16, MonoF32::new(0.5));
        let pyramid = Gaussian.build(&src, 3);
        assert!(pyramid.expand::<MonoF32, _>(3).is_none());
        assert!(pyramid.expand::<MonoF32, _>(99).is_none());
    }

    #[test]
    fn expand_works_on_an_imported_chain() {
        // Not builder-made: the validating constructor is what unlocks it.
        let chain = LevelChain::try_from_levels(vec![
            Image::fill(9, 7, MonoF32::new(0.5)),
            Image::fill(5, 4, MonoF32::new(0.5)),
        ])
        .unwrap();
        let pyramid = Dyadic::try_new(chain).unwrap();
        let raised: Image<MonoF32> = pyramid.expand(1).expect("level 1 has a parent");
        assert_eq!(raised.size(), Size::new(9, 7));
    }

    // ── Level→base coordinate lift property ─────────────────────────────

    #[test]
    fn decimated_lift_recovers_base_position() {
        // A bright Gaussian-ish blob at base position (12, 8): after two
        // pyr_down steps its maximum sits at level coordinates that must
        // lift back to (12, 8) via the even-sample convention.
        let (cx, cy) = (12.0f32, 8.0f32);
        let src = Image::generate(33, 25, |x, y| {
            let dx = x as f32 - cx;
            let dy = y as f32 - cy;
            MonoF32::new((-(dx * dx + dy * dy) / 18.0).exp())
        });

        // Two pyr_down steps, from the builder rather than by hand, so the
        // sampling distance under test is the one the build derived and not
        // a number this fixture restated. The lift reads geometry only, and
        // a `PlacedImage` level has nothing else to offer it.
        let pyramid: PlacedPyramid<MonoF32> = Gaussian.build(&src, 3);
        let scaled = pyramid.level(2);
        assert_eq!(scaled.pixel_distance(), pixel_distance!(4.0));

        // Find the argmax on the coarse level.
        let img = scaled.as_image();
        let mut best = (0usize, 0usize, f32::MIN);
        for y in 0..img.height() {
            for x in 0..img.width() {
                let v = img.pixel_at(x, y).0;
                if v > best.2 {
                    best = (x, y, v);
                }
            }
        }

        let lifted = scaled.to_base(CoordinateF64::new(best.0 as f64, best.1 as f64));
        assert_eq!(lifted, CoordinateF64::new(f64::from(cx), f64::from(cy)));
    }

    // ── Gaussian: the geometry the builder no longer discards ───────────

    #[test]
    fn gaussian_levels_carry_their_sampling_grid() {
        let src = Image::fill(32, 20, MonoF32::new(0.5));
        let pyramid: PlacedPyramid<MonoF32> = Gaussian.build(&src, 4);
        let distances: Vec<f64> = pyramid.iter().map(|l| l.pixel_distance().get()).collect();
        assert_eq!(distances, [1.0, 2.0, 4.0, 8.0]);
        for level in pyramid.iter() {
            assert_eq!(level.origin_offset(), CoordinateF64::new(0.0, 0.0));
        }
    }

    #[test]
    fn a_coarse_position_lifts_without_the_caller_naming_the_grid() {
        let src = Image::fill(32, 32, MonoF32::new(0.5));
        let pyramid: PlacedPyramid<MonoF32> = Gaussian.build(&src, 3);
        // pyr_down keeps even samples, so level-2 pixel 3 is base pixel 12.
        assert_eq!(
            pyramid.level(2).to_base(CoordinateF64::new(3.0, 1.0)),
            CoordinateF64::new(12.0, 4.0)
        );
    }

    // ── ScaledGaussian: the sigma ladder ────────────────────────────────

    #[test]
    fn the_sigma_ladder_matches_the_published_numbers_for_lowes_assumption() {
        let src = Image::fill(64, 64, MonoF32::new(0.5));
        let pyramid: ScaledPyramid<MonoF32> =
            Gaussian.assuming_input_sigma(sigma!(0.5)).build(&src, 4);
        let sigmas: Vec<f32> = pyramid.iter().map(|l| l.sigma().get()).collect();
        for (actual, expected) in sigmas.iter().zip([0.5, 1.118, 2.291, 4.610]) {
            assert!(
                (actual - expected).abs() < 1e-3,
                "sigma ladder: got {sigmas:?}"
            );
        }
    }

    #[test]
    fn the_sigma_ladder_matches_the_published_numbers_for_a_sharp_input() {
        let src = Image::fill(64, 64, MonoF32::new(0.5));
        // A perfectly sharp input is not expressible: Sigma is strictly
        // positive, and that is the point. This is as close as it gets.
        let pyramid: ScaledPyramid<MonoF32> =
            Gaussian.assuming_input_sigma(sigma!(0.001)).build(&src, 4);
        let sigmas: Vec<f32> = pyramid.iter().skip(1).map(|l| l.sigma().get()).collect();
        for (actual, expected) in sigmas.iter().zip([1.000, 2.236, 4.583]) {
            assert!(
                (actual - expected).abs() < 1e-3,
                "sigma ladder: got {sigmas:?}"
            );
        }
    }

    #[test]
    fn the_naive_extrapolation_is_the_error_the_ladder_prevents() {
        // The tempting guess is that sigma scales with the sampling
        // distance: sigma_in * 2^k, so 0.5, 1, 2, 4. It does not, because
        // variances add rather than sigmas. By level 3 the guess is 13% short
        // of the truth, silently. This test pins the gap so it cannot be
        // closed by accident.
        let src = Image::fill(64, 64, MonoF32::new(0.5));
        let pyramid: ScaledPyramid<MonoF32> =
            Gaussian.assuming_input_sigma(sigma!(0.5)).build(&src, 4);
        let level = pyramid.level(3);
        let naive = 0.5 * level.pixel_distance().get() as f32;
        let truth = level.sigma().get();
        assert!((naive - 4.0).abs() < 1e-6, "naive guess {naive}");
        assert!(
            ((truth - naive) / truth - 0.132).abs() < 5e-3,
            "truth {truth}, naive {naive}"
        );
    }

    #[test]
    fn scaled_levels_keep_the_same_grid_as_placed_ones() {
        let src = Image::fill(32, 20, MonoF32::new(0.5));
        let placed: PlacedPyramid<MonoF32> = Gaussian.build(&src, 4);
        let scaled: ScaledPyramid<MonoF32> =
            Gaussian.assuming_input_sigma(sigma!(0.5)).build(&src, 4);
        assert_eq!(placed.depth(), scaled.depth());
        for index in 0..placed.depth() {
            assert_eq!(placed.level(index).size(), scaled.level(index).size());
            assert_eq!(
                placed.level(index).pixel_distance(),
                scaled.level(index).pixel_distance()
            );
        }
    }

    #[test]
    fn the_input_sigma_is_readable_on_the_method_and_nowhere_else() {
        // The method keeps it because the caller named it; the levels store
        // their finished sigma, not the assumption that produced it.
        let method = Gaussian.assuming_input_sigma(sigma!(0.5));
        assert_eq!(method.input_sigma(), sigma!(0.5));
    }

    #[test]
    fn both_builders_produce_dyadic_pyramids_that_expand() {
        let src = Image::fill(21, 13, MonoF32::new(0.5));
        let placed: PlacedPyramid<MonoF32> = Gaussian.build(&src, 3);
        let scaled: ScaledPyramid<MonoF32> =
            Gaussian.assuming_input_sigma(sigma!(0.5)).build(&src, 3);
        let from_placed: Image<MonoF32> = placed.expand(2).expect("level 2 has a parent");
        let from_scaled: Image<MonoF32> = scaled.expand(2).expect("level 2 has a parent");
        assert_eq!(from_placed.size(), placed.level(1).size());
        assert_eq!(from_scaled.size(), scaled.level(1).size());
    }
}