//! Exact Euclidean distance from a mask's boundary, and the shaping built on
//! it.
//!
//! # One measurement, four controls
//!
//! Feathering, growing, shrinking, closing and opening are the same number
//! read differently. Given the **signed** distance from the boundary —
//! positive inside, negative outside — dilation by `r` is the set where
//! `d >= -r`, erosion is `d >= +r`, and a feather of any shape is a function
//! of `d`. So the field is computed once and the controls are arithmetic on
//! it.
//!
//! That is also why the *field* is what gets uploaded to the GPU rather than a
//! finished alpha: growing a mask or changing its falloff then costs a uniform
//! upload and no recomputation at all, which is what makes those live
//! controls rather than ones that stall on every drag.
//!
//! Closing and opening are the exception. After the first threshold the shape
//! has changed, so the old distances describe the old one and a second field
//! is needed — [`Shaped::needs_recompute`] says so, and it is the only edge
//! control that is not free.
//!
//! # Why this is on the CPU
//!
//! ARCH §5.4 says masks rasterise on the GPU and never exist in CPU memory,
//! and the reason it says so is brush lag: a stroke rasterised per-frame on
//! the CPU is what makes darktable's drawn masks unusable. This is a different
//! operation with different economics.
//!
//! - It runs **once per mask edit**, not once per frame.
//! - Its input is already CPU-side — the model's coverage was produced here.
//! - Its output is a field the GPU then samples for free, forever after.
//!
//! And doing it here buys two things a shader could not. It is **exactly
//! deterministic**, which matters because masks reach the sidecar as indices
//! and a field that varied by vendor would mean a mask meaning one thing on
//! the desktop and another on the phone (docs/segmentation.md §6, M5). And it
//! is testable against hand-computed distances with no adapter present.
//!
//! # The transform
//!
//! Felzenszwalb and Huttenlocher's separable exact transform: the lower
//! envelope of parabolas along every row, then along every column. Linear in
//! the number of pixels, exactly Euclidean — not the chamfer approximation
//! that leaves a mask visibly octagonal when grown by more than a few pixels.

/// Larger than any squared distance in an image anyone will render.
const FAR: f32 = 1e20;

/// Signed Euclidean distance from the mask's boundary, in pixels.
///
/// Positive inside, negative outside. `coverage` is one byte per pixel;
/// `threshold` is the value at or above which a pixel counts as inside.
///
/// The sign convention is the one that makes the controls read naturally: a
/// positive offset grows the mask, matching "dilate by 3 pixels".
///
/// # The half-pixel, which is not a detail
///
/// The transform measures to the nearest pixel *centre* of the other class, so
/// the closest an inside pixel can be to an outside one is exactly 1. Reported
/// raw, that puts the boundary nowhere — no pixel is at zero, the smallest
/// magnitude either side is 1, and **eroding by anything under a pixel removes
/// nothing at all**. A control whose first notch does nothing is a broken
/// control.
///
/// So half a pixel comes off each side, which puts the boundary where it
/// physically is: between the last inside pixel and the first outside one.
/// A pixel against the edge then reads `+0.5` inside and `-0.5` outside,
/// symmetric, and eroding by 1 takes exactly the outermost ring.
pub fn signed_distance(coverage: &[u8], width: usize, height: usize, threshold: u8) -> Vec<f32> {
    assert_eq!(
        coverage.len(),
        width * height,
        "coverage must cover every pixel"
    );

    let inside: Vec<bool> = coverage.iter().map(|&c| c >= threshold).collect();

    // Two transforms, because a pixel's distance is to the nearest pixel of
    // the *opposite* class and that is a different seed set on each side.
    let to_outside = euclidean(&inside, width, height, false);
    let to_inside = euclidean(&inside, width, height, true);

    inside
        .iter()
        .enumerate()
        .map(|(i, &is_in)| {
            if is_in {
                to_outside[i] - 0.5
            } else {
                -(to_inside[i] - 0.5)
            }
        })
        .collect()
}

/// Distance from every pixel to the nearest pixel whose class is `seed`.
fn euclidean(inside: &[bool], width: usize, height: usize, seed: bool) -> Vec<f32> {
    let mut f: Vec<f32> = inside
        .iter()
        .map(|&v| if v == seed { 0.0 } else { FAR })
        .collect();

    // Scratch, allocated once and reused by every line: the parabola vertices
    // still on the lower envelope, and the crossings between consecutive ones.
    let longest = width.max(height);
    let mut v = vec![0usize; longest];
    let mut z = vec![0.0f32; longest + 1];
    let mut out = vec![0.0f32; longest];

    for y in 0..height {
        let row: Vec<f32> = f[y * width..(y + 1) * width].to_vec();
        transform(&row, &mut out[..width], &mut v, &mut z);
        f[y * width..(y + 1) * width].copy_from_slice(&out[..width]);
    }

    let mut column = vec![0.0f32; height];
    for x in 0..width {
        for y in 0..height {
            column[y] = f[y * width + x];
        }
        transform(&column, &mut out[..height], &mut v, &mut z);
        for y in 0..height {
            f[y * width + x] = out[y];
        }
    }

    // Squared until here — the transform works in squares because that is what
    // makes the parabolas parabolas.
    f.iter().map(|d| d.max(0.0).sqrt()).collect()
}

/// One-dimensional squared distance transform.
///
/// `out[q] = min over p of ( f[p] + (q - p)^2 )`, computed by walking the
/// lower envelope of those parabolas.
fn transform(f: &[f32], out: &mut [f32], v: &mut [usize], z: &mut [f32]) {
    let n = f.len();
    if n == 0 {
        return;
    }

    let mut k = 0usize;
    v[0] = 0;
    z[0] = -FAR;
    z[1] = FAR;

    for q in 1..n {
        if f[q] >= FAR {
            // An infinite parabola is never the lowest anywhere, and
            // including it would divide one infinity by another.
            continue;
        }

        loop {
            let p = v[k];
            // Where this parabola crosses the one currently on top.
            let s = ((f[q] + sq(q)) - (f[p] + sq(p))) / (2.0 * q as f32 - 2.0 * p as f32);
            if s > z[k] {
                k += 1;
                v[k] = q;
                z[k] = s;
                z[k + 1] = FAR;
                break;
            }
            if k == 0 {
                // It dominates everything before it: start the envelope again
                // from this vertex.
                v[0] = q;
                z[0] = -FAR;
                z[1] = FAR;
                break;
            }
            k -= 1;
        }
    }

    // Every parabola was infinite, so every distance is.
    if f.iter().all(|&x| x >= FAR) {
        out[..n].fill(FAR);
        return;
    }

    k = 0;
    for q in 0..n {
        while z[k + 1] < q as f32 {
            k += 1;
        }
        let p = v[k];
        out[q] = (q as f32 - p as f32).powi(2) + f[p];
    }
}

fn sq(x: usize) -> f32 {
    (x as f32) * (x as f32)
}

/// How coverage falls away from the boundary.
///
/// Mirrors `dr_pipeline::mask::Falloff`, which this crate cannot see — the
/// pipeline depends on nothing here and inverting that to share one enum would
/// be a dependency edge for five variants.
#[derive(Debug, Clone, Copy, PartialEq, Eq, Default)]
pub enum Falloff {
    Hard,
    Linear,
    #[default]
    Smooth,
    Gaussian,
    Exponential,
}

impl Falloff {
    /// Coverage at a signed distance, given a feather half-width.
    ///
    /// `t` is the distance normalised to the feather: `-1` is a feather-width
    /// outside, `+1` a feather-width inside. Every curve returns `0.5` at the
    /// boundary, which is what keeps the *edge* where the mask says it is
    /// whichever curve is chosen — changing the falloff should change how the
    /// transition looks, never where it sits.
    pub fn coverage(self, t: f32) -> f32 {
        match self {
            Self::Hard => {
                if t >= 0.0 {
                    1.0
                } else {
                    0.0
                }
            }
            Self::Linear => (t * 0.5 + 0.5).clamp(0.0, 1.0),
            Self::Smooth => {
                let x = (t * 0.5 + 0.5).clamp(0.0, 1.0);
                x * x * (3.0 - 2.0 * x)
            }
            // A logistic curve rather than a true Gaussian integral: it is the
            // same shape to the eye, has a closed form, and is exactly 0.5 at
            // the boundary by construction.
            Self::Gaussian => 1.0 / (1.0 + (-3.0 * t).exp()),
            // Reaches full coverage quickly inside and trails off slowly
            // outside, for blending an adjustment away without moving its edge.
            Self::Exponential => {
                if t >= 0.0 {
                    1.0 - 0.5 * (-3.0 * t).exp()
                } else {
                    0.5 * (3.0 * t).exp()
                }
            }
        }
    }
}

/// Growing, shrinking and tidying. Mirrors `dr_pipeline::mask::Morphology`.
#[derive(Debug, Clone, Copy, PartialEq, Eq, Default)]
pub enum Morphology {
    #[default]
    None,
    Dilate,
    Erode,
    Close,
    Open,
}

impl Morphology {
    /// Whether this needs the distance field rebuilt after a threshold.
    pub fn needs_recompute(self) -> bool {
        matches!(self, Self::Close | Self::Open)
    }

    /// The offset applied to the distance before the falloff, in pixels.
    ///
    /// Zero for the compound pair: they are applied by [`apply_morphology`]
    /// rather than by shifting the field, because their second half operates
    /// on a shape the field does not describe.
    pub fn offset(self, radius: f32) -> f32 {
        match self {
            Self::Dilate => radius,
            Self::Erode => -radius,
            Self::None | Self::Close | Self::Open => 0.0,
        }
    }
}

/// Apply a compound morphology, returning a new distance field.
///
/// Closing is a dilation followed by an erosion, opening the reverse. Each
/// half is a threshold of a field, and the second half needs a field of the
/// *thresholded* shape — so this recomputes once in the middle and is the one
/// edge control that is not free.
///
/// Returns `None` for the operations that need no rebuild, so a caller can use
/// the field it already has.
pub fn apply_morphology(
    distance: &[f32],
    width: usize,
    height: usize,
    morphology: Morphology,
    radius: f32,
) -> Option<Vec<f32>> {
    if !morphology.needs_recompute() || radius <= 0.0 {
        return None;
    }

    // First half: threshold the existing field.
    let first = match morphology {
        Morphology::Close => -radius, // dilate
        Morphology::Open => radius,   // erode
        _ => return None,
    };
    let intermediate: Vec<u8> = distance
        .iter()
        .map(|&d| if d >= first { 255 } else { 0 })
        .collect();

    // Second half: measure the new shape, and shift so the caller's threshold
    // at zero performs the opposite operation.
    let rebuilt = signed_distance(&intermediate, width, height, 128);
    let second = match morphology {
        Morphology::Close => radius, // erode
        Morphology::Open => -radius, // dilate
        _ => unreachable!("guarded above"),
    };
    Some(rebuilt.iter().map(|&d| d - second).collect())
}

/// A distance field ready to be sampled, with the offset already folded in.
#[derive(Debug, Clone, PartialEq)]
pub struct Shaped {
    pub distance: Vec<f32>,
    pub width: usize,
    pub height: usize,
}

impl Shaped {
    /// Build from coverage, applying whatever morphology needs a rebuild.
    ///
    /// The simple operations are *not* folded in here: they are an offset the
    /// shader adds when it samples, so changing "grow by 4px" to "grow by 6px"
    /// costs a uniform rather than a transform.
    pub fn build(
        coverage: &[u8],
        width: usize,
        height: usize,
        threshold: u8,
        morphology: Morphology,
        radius_px: f32,
    ) -> Self {
        let distance = signed_distance(coverage, width, height, threshold);
        let distance =
            apply_morphology(&distance, width, height, morphology, radius_px).unwrap_or(distance);
        Self {
            distance,
            width,
            height,
        }
    }

    /// Coverage at one pixel, for tests and for the CPU export path.
    pub fn coverage_at(&self, index: usize, offset: f32, feather: f32, falloff: Falloff) -> f32 {
        let d = self.distance[index] + offset;
        if feather <= 0.0 {
            return if d >= 0.0 { 1.0 } else { 0.0 };
        }
        falloff.coverage(d / feather)
    }
}

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

    /// A 2px-wide vertical bar down the middle of a 9x1 strip.
    fn bar() -> (Vec<u8>, usize, usize) {
        let mut m = vec![0u8; 9];
        m[4] = 255;
        (m, 9, 1)
    }

    #[test]
    fn distance_is_exact_along_a_line() {
        let (m, w, h) = bar();
        let d = signed_distance(&m, w, h, 128);

        // The lone inside pixel sits half a pixel from the boundary on each
        // side, and the pixels either side of it are half a pixel out.
        assert_eq!(d[4], 0.5);
        assert_eq!(d[3], -0.5);
        assert_eq!(d[5], -0.5);
        // Then one per pixel from there.
        assert_eq!(d[2], -1.5);
        assert_eq!(d[0], -3.5);
        assert_eq!(d[8], -3.5);
    }

    /// The property that separates an exact transform from a chamfer one: a
    /// diagonal neighbour is √2 away, not 1 and not 2.
    #[test]
    fn diagonals_are_euclidean_not_chamfer() {
        let mut m = vec![0u8; 25];
        m[12] = 255; // centre of 5x5
        let d = signed_distance(&m, 5, 5, 128);

        // Distance from the boundary, so the half-pixel comes back off to
        // compare against the centre-to-centre figures.
        let at = |x: usize, y: usize| -d[y * 5 + x] + 0.5;
        assert!(
            (at(1, 1) - std::f32::consts::SQRT_2).abs() < 1e-4,
            "{}",
            at(1, 1)
        );
        assert!((at(0, 0) - (8.0f32).sqrt()).abs() < 1e-4, "{}", at(0, 0));
        assert_eq!(at(2, 0), 2.0, "straight up is exactly two");
    }

    #[test]
    fn an_empty_mask_is_everywhere_outside() {
        let d = signed_distance(&vec![0u8; 16], 4, 4, 128);
        assert!(d.iter().all(|&v| v < 0.0), "no pixel can be inside");
    }

    #[test]
    fn a_full_mask_is_everywhere_inside() {
        let d = signed_distance(&vec![255u8; 16], 4, 4, 128);
        assert!(d.iter().all(|&v| v > 0.0), "no pixel can be outside");
    }

    /// Every curve must cross at the boundary, or changing the falloff would
    /// move the edge rather than soften it.
    #[test]
    fn every_falloff_is_half_at_the_boundary() {
        for f in [
            Falloff::Linear,
            Falloff::Smooth,
            Falloff::Gaussian,
            Falloff::Exponential,
        ] {
            let v = f.coverage(0.0);
            assert!((v - 0.5).abs() < 1e-5, "{f:?} gave {v} at the boundary");
        }
    }

    #[test]
    fn every_falloff_is_monotone_and_bounded() {
        for f in Falloff::ALL_FOR_TEST {
            let mut previous = -1.0;
            for i in -20..=20 {
                let v = f.coverage(i as f32 / 10.0);
                assert!((0.0..=1.0).contains(&v), "{f:?} left the range: {v}");
                assert!(v >= previous - 1e-6, "{f:?} went backwards at {i}");
                previous = v;
            }
        }
    }

    #[test]
    fn a_hard_falloff_has_no_transition() {
        assert_eq!(Falloff::Hard.coverage(-0.01), 0.0);
        assert_eq!(Falloff::Hard.coverage(0.0), 1.0);
    }

    #[test]
    fn dilating_grows_and_eroding_shrinks() {
        let mut m = vec![0u8; 81];
        for y in 3..6 {
            for x in 3..6 {
                m[y * 9 + x] = 255;
            }
        }
        let d = signed_distance(&m, 9, 9, 128);

        let area = |offset: f32| d.iter().filter(|&&v| v + offset >= 0.0).count();
        let plain = area(0.0);
        assert_eq!(plain, 9, "the 3x3 block itself");
        assert!(area(Morphology::Dilate.offset(1.0)) > plain, "dilate grows");
        assert!(area(Morphology::Erode.offset(1.0)) < plain, "erode shrinks");
    }

    /// Closing fills a hole without growing the outline — the thing it is for.
    #[test]
    fn closing_fills_a_pinhole() {
        // Generous margin on purpose. Outside the image is not "outside the
        // mask", so a dilation that reaches the border cannot erode back from
        // it, and a tight frame measures that artefact instead of the
        // operation.
        let (w, h) = (21, 21);
        let mut m = vec![0u8; w * h];
        for y in 7..14 {
            for x in 7..14 {
                m[y * w + x] = 255;
            }
        }
        // One missing pixel inside, as a soft mask commonly has.
        m[10 * w + 10] = 0;

        let before = signed_distance(&m, w, h, 128);
        assert!(
            before[10 * w + 10] < 0.0,
            "the hole starts outside the mask"
        );

        let after = apply_morphology(&before, w, h, Morphology::Close, 2.0).expect("recomputed");
        assert!(after[10 * w + 10] >= 0.0, "closing should have filled it");

        // The outline may round at a corner — closing does that, and it is the
        // price of filling — but it must not march outward across the frame.
        let grew = after
            .iter()
            .zip(&before)
            .filter(|(a, b)| **a >= 0.0 && **b < 0.0)
            .count();
        assert!(grew <= 5, "closing should fill the hole, not grow: {grew}");
    }

    /// Opening removes a speck without shrinking the body.
    #[test]
    fn opening_removes_a_speck() {
        let (w, h) = (13, 13);
        let mut m = vec![0u8; w * h];
        for y in 3..10 {
            for x in 3..10 {
                m[y * w + x] = 255;
            }
        }
        m[0] = 255; // an isolated speck in the corner

        let before = signed_distance(&m, w, h, 128);
        assert!(before[0] >= 0.0, "the speck starts inside the mask");

        let after = apply_morphology(&before, w, h, Morphology::Open, 1.5).expect("recomputed");
        assert!(after[0] < 0.0, "opening should have removed it");
        assert!(
            after[6 * w + 6] >= 0.0,
            "and left the body of the mask alone"
        );
    }

    #[test]
    fn the_simple_operations_need_no_rebuild() {
        assert!(!Morphology::None.needs_recompute());
        assert!(!Morphology::Dilate.needs_recompute());
        assert!(!Morphology::Erode.needs_recompute());
        assert!(Morphology::Close.needs_recompute());
        assert!(Morphology::Open.needs_recompute());
    }

    #[test]
    fn the_same_mask_gives_the_same_field_twice() {
        // M5: masks reach the sidecar as indices, so the field they are shaped
        // by has to be reproducible.
        let mut m = vec![0u8; 64];
        for i in 20..30 {
            m[i] = 255;
        }
        let a = signed_distance(&m, 8, 8, 128);
        let b = signed_distance(&m, 8, 8, 128);
        assert_eq!(a, b);
    }

    #[test]
    fn zero_feather_is_a_step() {
        let shaped = Shaped::build(&[0, 255, 255, 0], 4, 1, 128, Morphology::None, 0.0);
        assert_eq!(shaped.coverage_at(1, 0.0, 0.0, Falloff::Smooth), 1.0);
        assert_eq!(shaped.coverage_at(3, 0.0, 0.0, Falloff::Smooth), 0.0);
    }

    impl Falloff {
        const ALL_FOR_TEST: [Falloff; 5] = [
            Falloff::Hard,
            Falloff::Linear,
            Falloff::Smooth,
            Falloff::Gaussian,
            Falloff::Exponential,
        ];
    }
}