rust/vectorhd-core/src/window.rs

//! The **windowed** kernel — a port of `optimize._windowed_data`, which is what production
//! actually calls (hundreds of times per image, from inside the L-BFGS-B objective).
//!
//! The spike ported `jacobian.geometry_gradient`, the *full-graph oracle*, which has no call site
//! in `src/`. This module is the production-shaped contract, and answers spike risk R1.
//!
//! # Contract
//!
//! Given a window (already reduced to its own local coordinate frame by the caller):
//!
//! * `target` — the `(h, w, 3)` image patch under the window,
//! * face loops — flattened polylines **already shifted** by the window origin, in the caller's
//!   face order (this preserves the composite's reduction order),
//! * per-gradient-cubic flattened vertices, also shifted, with the labels of the regions the edge
//!   separates,
//!
//! return `E_data` for the window and `∂E_data/∂v` for every vertex of every gradient cubic.
//!
//! The Bernstein chain rule back to control points stays in Python, exactly as in the spike: it is
//! a tiny dense matmul, it is not the hot loop, and keeping the boundary at *vertex* level means
//! the crate never performs the `bernstein(linspace(...)) @ controls` BLAS product, which is not
//! bit-reproducible outside numpy.
//!
//! # Conventions inherited and one added
//!
//! Everything `coverage`/`gradient` already pin (suffix-sum direction, `left[width]` retention,
//! truncate-not-floor, frame breaks, recomputed break `y`) applies unchanged. The window adds one:
//!
//! **Clipping is by truncation into the window frame, not by rejection.** Shifted geometry
//! routinely lands at negative coordinates, and `_segment_tau_forces` locates a sub-piece with
//! `int(a + t·d)` — truncation *toward zero* — then bounds-checks `0 <= px < w`. So a sub-piece at
//! `x = −0.5` truncates to column `0` and is **accepted**, while one at `x = −1.5` truncates to
//! `−1` and is rejected. `floor()` would reject both and silently drop gradient contributions along
//! the window's left/top edges. Rust's `as i64` truncates toward zero, matching Python's `int()`.

use crate::coverage::polygon_coverage;
use crate::energy::{sum_compensated, Neumaier};
use crate::gradient::segment_tau_forces;
use crate::model::Pt;
use std::collections::HashMap;

/// A region's colour model — flat, or C-5's quadratic in a centered/scaled basis.
///
/// `c(x, y) = C₀ + C₁·u + C₂·v + C₃·u² + C₄·u·v + C₅·v²`, with `u = (x − cx)/s`, `v = (y − cy)/s`.
///
/// **The transform `(cx, cy, s)` is computed in Python and marshalled — never recomputed here.**
/// It exists to keep the `u²/uv/v²` terms O(1) instead of O(1e4–1e6) in raw pixel coordinates
/// (C-5's Trap 2, conditioning). Recomputing it on this side would create a second source of truth
/// that could silently drift from the one VarPro solved against.
#[derive(Clone, Debug)]
pub enum RegionColor {
    Flat([f64; 3]),
    Quad {
        /// Six basis terms × RGB, row-major: `coeffs[k][ch]`.
        coeffs: [[f64; 3]; 6],
        cx: f64,
        cy: f64,
        s: f64,
    },
}

impl RegionColor {
    /// Colour at a **global** pixel centre, clipped to [0,1] exactly as `color_field` does.
    #[inline]
    pub fn eval(&self, x: f64, y: f64) -> [f64; 3] {
        match self {
            RegionColor::Flat(c) => *c,
            RegionColor::Quad { coeffs, cx, cy, s } => {
                let u = (x - cx) / s;
                let v = (y - cy) / s;
                let b = [1.0, u, v, u * u, u * v, v * v];
                let mut out = [0.0_f64; 3];
                for (ch, o) in out.iter_mut().enumerate() {
                    let mut acc = 0.0;
                    for k in 0..6 {
                        acc += b[k] * coeffs[k][ch];
                    }
                    // Python clips the evaluated field before it is used, in BOTH the forward
                    // model and the colour jump — so the clip is part of the model, not cosmetic.
                    *o = acc.clamp(0.0, 1.0);
                }
                out
            }
        }
    }

    #[inline]
    pub fn is_quad(&self) -> bool {
        matches!(self, RegionColor::Quad { .. })
    }
}

/// One face's contribution to the window: its label index and its shifted, flattened loops.
pub struct Face {
    pub label: usize,
    pub loops: Vec<Vec<Pt>>,
}

/// One cubic whose vertices we differentiate, with the regions its edge separates.
pub struct GradCubic {
    pub left: usize,
    pub right: usize,
    /// Flattened, window-shifted vertices (`m` of them).
    pub pts: Vec<Pt>,
}

/// `E_data` over the window plus per-vertex gradients, mirroring `_windowed_data`.
///
/// `x0`/`y0` are the window origin in image coordinates. They matter only for quadratic regions,
/// whose colour is evaluated at **global** pixel centres (`x0 + col + 0.5`, `y0 + row + 0.5`) —
/// the transform is global, so a window-local evaluation would silently shift every gradient
/// region's colour field.
#[allow(clippy::too_many_arguments)]
pub fn windowed_data(
    width: usize,
    height: usize,
    x0: f64,
    y0: f64,
    target: &[f64],
    faces: &[Face],
    colors: &[RegionColor],
    background: f64,
    l0: f64,
    grads: &[GradCubic],
    weights: Option<&[f64]>,
) -> (f64, Vec<Vec<[f64; 2]>>) {
    let n = width * height;

    // Evaluated colour field per quadratic label, built once and shared by the composite and every
    // edge's colour jump. Python re-evaluates it per edge; caching changes no value, only work.
    let mut quad_fields: HashMap<usize, Vec<[f64; 3]>> = HashMap::new();
    let field_for = |label: usize, cache: &mut HashMap<usize, Vec<[f64; 3]>>| {
        if colors[label].is_quad() && !cache.contains_key(&label) {
            let mut f = Vec::with_capacity(n);
            for row in 0..height {
                let gy = y0 + row as f64 + 0.5;
                for col in 0..width {
                    f.push(colors[label].eval(x0 + col as f64 + 0.5, gy));
                }
            }
            cache.insert(label, f);
        }
    };

    // --- render_window: coverage per face, then composite over the background ----------------
    let mut img = vec![background; n * 3];
    for f in faces {
        if f.loops.is_empty() {
            continue;
        }
        let cov = polygon_coverage(&f.loops, width, height);
        match &colors[f.label] {
            RegionColor::Flat(c) => {
                let (d0, d1, d2) = (c[0] - background, c[1] - background, c[2] - background);
                for p in 0..n {
                    let a = cov[p].abs();
                    img[p * 3] += a * d0;
                    img[p * 3 + 1] += a * d1;
                    img[p * 3 + 2] += a * d2;
                }
            }
            RegionColor::Quad { .. } => {
                field_for(f.label, &mut quad_fields);
                let cf = &quad_fields[&f.label];
                for p in 0..n {
                    let a = cov[p].abs();
                    img[p * 3] += a * (cf[p][0] - background);
                    img[p * 3 + 1] += a * (cf[p][1] - background);
                    img[p * 3 + 2] += a * (cf[p][2] - background);
                }
            }
        }
    }
    for x in img.iter_mut() {
        *x = x.clamp(0.0, 1.0);
    }

    // --- residual + E_data --------------------------------------------------------------------
    let mut resid = vec![0.0_f64; n * 3];
    for i in 0..n * 3 {
        resid[i] = img[i] - target[i];
    }
    let mut acc = Neumaier::default();
    for p in 0..n {
        // exact 3-term per-pixel group, matching numpy's (resid*resid).sum(axis=2)
        let s = resid[p * 3] * resid[p * 3]
            + resid[p * 3 + 1] * resid[p * 3 + 1]
            + resid[p * 3 + 2] * resid[p * 3 + 2];
        acc.add(match weights {
            Some(w) => s * w[p],
            None => s,
        });
    }
    let e_data = acc.total() / l0;

    // --- boundary-integral gradient per requested cubic ---------------------------------------
    let scale = 2.0 / l0;
    let mut field = vec![0.0_f64; n];
    let mut out: Vec<Vec<[f64; 2]>> = Vec::with_capacity(grads.len());

    for g in grads {
        // C-5 Trap 1: with a quadratic region on either side the colour jump is POSITION-DEPENDENT
        // and must be evaluated per pixel. Collapsing it to a per-edge constant is the classic way
        // to get a plausible-looking but wrong gradient.
        let mixed = colors[g.left].is_quad() || colors[g.right].is_quad();
        if mixed {
            field_for(g.left, &mut quad_fields);
            field_for(g.right, &mut quad_fields);
            let fl = quad_fields.get(&g.left);
            let fr = quad_fields.get(&g.right);
            let flat_l = if let RegionColor::Flat(c) = &colors[g.left] { *c } else { [0.0; 3] };
            let flat_r = if let RegionColor::Flat(c) = &colors[g.right] { *c } else { [0.0; 3] };
            for p in 0..n {
                let cl = match fl { Some(f) => f[p], None => flat_l };
                let cr = match fr { Some(f) => f[p], None => flat_r };
                let v = resid[p * 3] * (cl[0] - cr[0])
                    + resid[p * 3 + 1] * (cl[1] - cr[1])
                    + resid[p * 3 + 2] * (cl[2] - cr[2]);
                field[p] = match weights {
                    Some(w) => v * w[p],
                    None => v,
                };
            }
        } else {
            let cl = if let RegionColor::Flat(c) = &colors[g.left] { *c } else { [0.0; 3] };
            let cr = if let RegionColor::Flat(c) = &colors[g.right] { *c } else { [0.0; 3] };
            let cd = [cl[0] - cr[0], cl[1] - cr[1], cl[2] - cr[2]];
            for p in 0..n {
                let v = resid[p * 3] * cd[0] + resid[p * 3 + 1] * cd[1] + resid[p * 3 + 2] * cd[2];
                field[p] = match weights {
                    Some(w) => v * w[p],
                    None => v,
                };
            }
        }

        let m = g.pts.len();
        let mut vgrad = vec![[0.0_f64; 2]; m];
        for i in 0..m.saturating_sub(1) {
            let sx = g.pts[i + 1][0] - g.pts[i][0];
            let sy = g.pts[i + 1][1] - g.pts[i][1];
            let length = sx.hypot(sy);
            if length < 1e-12 {
                continue;
            }
            let n0 = -sy / length; // LEFT normal
            let n1 = sx / length;
            let (i0, i1) = segment_tau_forces(g.pts[i], g.pts[i + 1], &field, width, height);
            let f0 = scale * i0;
            let f1 = scale * i1;
            vgrad[i][0] += f0 * n0;
            vgrad[i][1] += f0 * n1;
            vgrad[i + 1][0] += f1 * n0;
            vgrad[i + 1][1] += f1 * n1;
        }
        out.push(vgrad);
    }

    (e_data, out)
}

/// Sum of the window's coverage per face — diagnostic only, used by the golden-window tests to
/// localise a mismatch to the rasterizer rather than the energy.
pub fn window_coverage_sums(width: usize, height: usize, faces: &[Face]) -> Vec<f64> {
    faces
        .iter()
        .map(|f| {
            if f.loops.is_empty() {
                return 0.0;
            }
            let cov = polygon_coverage(&f.loops, width, height);
            sum_compensated(cov.iter().map(|x| x.abs()))
        })
        .collect()
}

#[cfg(feature = "python")]
pub use bindings::register;

#[cfg(feature = "python")]
mod bindings {
    use super::*;
    use numpy::{PyReadonlyArray2, PyReadonlyArray3, ToPyArray};
    use pyo3::prelude::*;
    use pyo3::types::PyList;

    /// Rebuild `Face`s from the flat (loop → face index) encoding the Python side sends.
    fn build_faces(
        n_faces: usize,
        face_labels: &[usize],
        loop_face: &[usize],
        loops: &Bound<'_, PyList>,
    ) -> PyResult<Vec<Face>> {
        let mut faces: Vec<Face> = (0..n_faces)
            .map(|i| Face {
                label: face_labels[i],
                loops: Vec::new(),
            })
            .collect();
        for (i, item) in loops.iter().enumerate() {
            let arr: PyReadonlyArray2<f64> = item.extract()?;
            let v = arr.as_array();
            let poly: Vec<Pt> = v.rows().into_iter().map(|r| [r[0], r[1]]).collect();
            faces[loop_face[i]].loops.push(poly);
        }
        Ok(faces)
    }

    /// `_windowed_data`, in Rust. See the module docs for the contract.
    ///
    /// Arrays arrive as read-only numpy views (no copy on the way in). The gradient result is the
    /// one allocation, shaped `(n_grad_cubics, m, 2)`.
    #[pyfunction]
    #[pyo3(signature = (width, height, x0, y0, target, face_labels, loop_face, loops, colors01,
                        color_kind, quad_coeffs, quad_transform, background, l0,
                        grad_left, grad_right, grad_pts, weights=None))]
    #[allow(clippy::too_many_arguments)]
    fn windowed_data<'py>(
        py: Python<'py>,
        width: usize,
        height: usize,
        x0: f64,
        y0: f64,
        target: PyReadonlyArray3<'py, f64>,
        face_labels: Vec<usize>,
        loop_face: Vec<usize>,
        loops: &Bound<'py, PyList>,
        colors01: PyReadonlyArray2<'py, f64>,
        // Per label: 0 = flat, 1 = quadratic. The coeff/transform arrays are full-length with
        // unused rows for flat labels — simpler and cheaper than a ragged structure, since the
        // label count is small.
        color_kind: Vec<u8>,
        quad_coeffs: PyReadonlyArray3<'py, f64>,
        quad_transform: PyReadonlyArray2<'py, f64>,
        background: f64,
        l0: f64,
        grad_left: Vec<usize>,
        grad_right: Vec<usize>,
        grad_pts: &Bound<'py, PyList>,
        weights: Option<PyReadonlyArray2<'py, f64>>,
    ) -> PyResult<(f64, Py<numpy::PyArray3<f64>>)> {
        let faces = build_faces(face_labels.len(), &face_labels, &loop_face, loops)?;

        let colors_view = colors01.as_array();
        let qc = quad_coeffs.as_array();
        let qt = quad_transform.as_array();
        let colors: Vec<RegionColor> = colors_view
            .rows()
            .into_iter()
            .enumerate()
            .map(|(i, r)| {
                if color_kind.get(i).copied().unwrap_or(0) == 1 {
                    let mut coeffs = [[0.0_f64; 3]; 6];
                    for (k, row) in coeffs.iter_mut().enumerate() {
                        for (ch, c) in row.iter_mut().enumerate() {
                            *c = qc[[i, k, ch]];
                        }
                    }
                    RegionColor::Quad { coeffs, cx: qt[[i, 0]], cy: qt[[i, 1]], s: qt[[i, 2]] }
                } else {
                    RegionColor::Flat([r[0], r[1], r[2]])
                }
            })
            .collect();

        let mut grads: Vec<GradCubic> = Vec::with_capacity(grad_left.len());
        for (i, item) in grad_pts.iter().enumerate() {
            let arr: PyReadonlyArray2<f64> = item.extract()?;
            let v = arr.as_array();
            grads.push(GradCubic {
                left: grad_left[i],
                right: grad_right[i],
                pts: v.rows().into_iter().map(|r| [r[0], r[1]]).collect(),
            });
        }

        let t = target.as_array();
        let t_slice: Vec<f64> = t.iter().copied().collect();
        let w_slice: Option<Vec<f64>> = weights.map(|w| w.as_array().iter().copied().collect());

        let (e, g) = super::windowed_data(
            width,
            height,
            x0,
            y0,
            &t_slice,
            &faces,
            &colors,
            background,
            l0,
            &grads,
            w_slice.as_deref(),
        );

        let m = g.first().map(|v| v.len()).unwrap_or(0);
        let mut flat = Vec::with_capacity(g.len() * m * 2);
        for vg in &g {
            for p in vg {
                flat.push(p[0]);
                flat.push(p[1]);
            }
        }
        let arr = numpy::ndarray::Array3::from_shape_vec((g.len(), m, 2), flat)
            .map_err(|e| pyo3::exceptions::PyValueError::new_err(e.to_string()))?;
        Ok((e, arr.to_pyarray(py).unbind()))
    }

    /// Per-face coverage sums for one window — lets a failing golden-window test say whether the
    /// rasterizer or the energy diverged.
    #[pyfunction]
    fn window_coverage_sums<'py>(
        py: Python<'py>,
        width: usize,
        height: usize,
        face_labels: Vec<usize>,
        loop_face: Vec<usize>,
        loops: &Bound<'py, PyList>,
    ) -> PyResult<Py<numpy::PyArray1<f64>>> {
        let faces = build_faces(face_labels.len(), &face_labels, &loop_face, loops)?;
        let sums = super::window_coverage_sums(width, height, &faces);
        Ok(numpy::ndarray::Array1::from_vec(sums)
            .to_pyarray(py)
            .unbind())
    }

    pub fn register(m: &Bound<'_, PyModule>) -> PyResult<()> {
        m.add_function(wrap_pyfunction!(self::windowed_data, m)?)?;
        m.add_function(wrap_pyfunction!(self::window_coverage_sums, m)?)?;
        Ok(())
    }
}

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

    /// A window fully covered by one face must render that face's colour exactly, so the residual
    /// against a target painted the same colour is zero.
    #[test]
    fn full_coverage_window_has_zero_energy_against_matching_target() {
        let (w, h) = (4_usize, 4_usize);
        let face = Face {
            label: 1,
            loops: vec![vec![
                [0.0, 0.0],
                [w as f64, 0.0],
                [w as f64, h as f64],
                [0.0, h as f64],
            ]],
        };
        let colors = [
            RegionColor::Flat([1.0, 1.0, 1.0]),
            RegionColor::Flat([0.25, 0.5, 0.75]),
        ];
        let mut target = vec![0.0; w * h * 3];
        for p in 0..w * h {
            target[p * 3] = 0.25;
            target[p * 3 + 1] = 0.5;
            target[p * 3 + 2] = 0.75;
        }
        let (e, g) = windowed_data(w, h, 0.0, 0.0, &target, &[face], &colors, 1.0, 1.0, &[], None);
        assert!(e.abs() < 1e-15, "energy {e} should vanish");
        assert!(g.is_empty());
    }

    /// An empty window (no faces) renders as background; energy is the background-vs-target
    /// residual, divided by l0.
    #[test]
    fn empty_window_is_background_over_l0() {
        let (w, h) = (2_usize, 3_usize);
        let target = vec![0.0; w * h * 3];
        let flat = [RegionColor::Flat([0.0, 0.0, 0.0])];
        let (e, _) = windowed_data(w, h, 0.0, 0.0, &target, &[], &flat, 1.0, 2.0, &[], None);
        // background 1.0 vs target 0.0 over 6 px x 3 channels = 18, over l0 = 2 -> 9
        assert!((e - 9.0).abs() < 1e-12, "energy {e}");
    }

    /// Per-pixel weights scale the residual exactly.
    #[test]
    fn weights_scale_the_energy() {
        let (w, h) = (2_usize, 2_usize);
        let target = vec![0.0; w * h * 3];
        let weights = vec![0.5; w * h];
        let flat = [RegionColor::Flat([0.0, 0.0, 0.0])];
        let (unweighted, _) =
            windowed_data(w, h, 0.0, 0.0, &target, &[], &flat, 1.0, 1.0, &[], None);
        let (weighted, _) =
            windowed_data(w, h, 0.0, 0.0, &target, &[], &flat, 1.0, 1.0, &[], Some(&weights));
        assert!((weighted - 0.5 * unweighted).abs() < 1e-12);
    }

    /// The added window convention: a sub-piece at x = −0.5 truncates to column 0 and is KEPT,
    /// while one at x = −1.5 truncates to −1 and is dropped. `floor()` would drop both.
    #[test]
    fn negative_coordinates_truncate_toward_zero_not_floor() {
        let field = vec![1.0; 4]; // 2x2 window
        // Segment sitting at x in (-1, 0): midpoints truncate to 0 -> accepted.
        let (a0, a1) = segment_tau_forces([-0.5, 0.5], [-0.5, 1.5], &field, 2, 2);
        assert!(
            a0.abs() + a1.abs() > 0.0,
            "x=-0.5 truncates to column 0 and must contribute"
        );
        // Segment at x in (-2, -1): midpoints truncate to -1 -> rejected.
        let (b0, b1) = segment_tau_forces([-1.5, 0.5], [-1.5, 1.5], &field, 2, 2);
        assert_eq!((b0, b1), (0.0, 0.0), "x=-1.5 truncates to -1 and is outside");
    }

    #[test]
    fn coverage_sums_report_per_face_area() {
        let (w, h) = (6_usize, 6_usize);
        let f = Face {
            label: 0,
            loops: vec![vec![[1.0, 1.0], [4.0, 1.0], [4.0, 3.0], [1.0, 3.0]]],
        };
        let sums = window_coverage_sums(w, h, &[f]);
        assert!((sums[0] - 6.0).abs() < 1e-12, "area {}", sums[0]);
    }
}

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

    fn quad(coeffs: [[f64; 3]; 6], cx: f64, cy: f64, s: f64) -> RegionColor {
        RegionColor::Quad { coeffs, cx, cy, s }
    }

    /// The evaluated field must equal a directly-computed reference of
    /// `C₀ + C₁u + C₂v + C₃u² + C₄uv + C₅v²` on the centered/scaled basis.
    #[test]
    fn quadratic_eval_matches_a_direct_reference() {
        let coeffs = [
            [0.50, 0.40, 0.30],
            [0.10, -0.05, 0.02],
            [-0.03, 0.07, 0.01],
            [0.01, 0.02, -0.01],
            [0.02, -0.01, 0.03],
            [-0.02, 0.01, 0.02],
        ];
        let (cx, cy, s) = (37.5, 21.25, 18.0);
        let c = quad(coeffs, cx, cy, s);
        for &(x, y) in &[(30.0, 20.0), (37.5, 21.25), (48.0, 33.0), (12.5, 5.5)] {
            let u = (x - cx) / s;
            let v = (y - cy) / s;
            let b = [1.0, u, v, u * u, u * v, v * v];
            let got = c.eval(x, y);
            for ch in 0..3 {
                let want: f64 = (0..6).map(|k| b[k] * coeffs[k][ch]).sum::<f64>().clamp(0.0, 1.0);
                assert!(
                    (got[ch] - want).abs() < 1e-15,
                    "ch{ch} at ({x},{y}): {} vs {want}",
                    got[ch]
                );
            }
        }
    }

    /// The evaluated field is clipped to [0,1], as `color_field` does — the clip is part of the
    /// model, and both the forward model and the colour jump see the clipped value.
    #[test]
    fn quadratic_eval_is_clipped() {
        let mut coeffs = [[0.0_f64; 3]; 6];
        coeffs[0] = [5.0, -5.0, 0.5]; // way out of gamut both ways
        let c = quad(coeffs, 0.0, 0.0, 1.0);
        assert_eq!(c.eval(1.0, 1.0), [1.0, 0.0, 0.5]);
    }

    /// Conditioning (C-5 Trap 2): the centered/scaled basis makes a region far from the origin a
    /// non-event. Same shape, same coefficients, translated by 10,000 px — identical colours.
    /// Evaluating in RAW pixel coordinates instead would blow the u² term up by ~1e8.
    #[test]
    fn far_from_origin_is_a_non_event_thanks_to_the_transform() {
        let coeffs = [
            [0.5, 0.5, 0.5],
            [0.2, 0.1, 0.0],
            [0.0, 0.1, 0.2],
            [0.05, 0.0, 0.0],
            [0.0, 0.05, 0.0],
            [0.0, 0.0, 0.05],
        ];
        let near = quad(coeffs, 50.0, 50.0, 25.0);
        let far = quad(coeffs, 10_050.0, 10_050.0, 25.0);
        for &(dx, dy) in &[(-20.0, -20.0), (0.0, 0.0), (17.0, -8.0), (25.0, 25.0)] {
            let a = near.eval(50.0 + dx, 50.0 + dy);
            let b = far.eval(10_050.0 + dx, 10_050.0 + dy);
            for ch in 0..3 {
                assert!(
                    (a[ch] - b[ch]).abs() < 1e-15,
                    "translation changed the colour: {a:?} vs {b:?}"
                );
            }
        }
        // And the basis really is O(1) there — the point of the transform.
        let u: f64 = (10_075.0 - 10_050.0) / 25.0;
        assert!(u.abs() <= 1.0 + 1e-12, "u should be O(1), got {u}");
    }

    /// Mixed windows are first class: flat on one side, quadratic on the other. The colour jump
    /// must be position-dependent (Trap 1) — if it were collapsed to a constant, moving the
    /// quadratic region's transform would not change the gradient. It must.
    #[test]
    fn mixed_flat_quadratic_jump_is_position_dependent() {
        const W: usize = 10;
        const H: usize = 10;
        let coeffs = [
            [0.9, 0.2, 0.2],
            [0.4, 0.0, 0.0], // strong horizontal ramp -> jump varies across the window
            [0.0, 0.0, 0.0],
            [0.0, 0.0, 0.0],
            [0.0, 0.0, 0.0],
            [0.0, 0.0, 0.0],
        ];
        let colors = vec![
            RegionColor::Flat([0.1, 0.1, 0.1]),
            quad(coeffs, 5.0, 5.0, 5.0),
        ];
        let left = vec![[0.0, 0.0], [5.0, 0.0], [5.0, H as f64], [0.0, H as f64]];
        let right = vec![
            [5.0, 0.0],
            [W as f64, 0.0],
            [W as f64, H as f64],
            [5.0, H as f64],
        ];
        let faces = vec![
            Face { label: 0, loops: vec![left] },
            Face { label: 1, loops: vec![right] },
        ];
        let target = vec![0.5; W * H * 3];
        let grads = vec![GradCubic {
            left: 1,
            right: 0,
            pts: vec![[5.0, 1.0], [5.0, 5.0], [5.0, 9.0]],
        }];

        let (_, g_a) = windowed_data(
            W, H, 0.0, 0.0, &target, &faces, &colors, 1.0, 1.0, &grads, None,
        );
        // Shift only the quadratic region's transform: the ramp now sits differently over the same
        // pixels, so a position-dependent jump changes the gradient.
        let colors_b = vec![
            RegionColor::Flat([0.1, 0.1, 0.1]),
            quad(coeffs, 1.0, 5.0, 5.0),
        ];
        let (_, g_b) = windowed_data(
            W, H, 0.0, 0.0, &target, &faces, &colors_b, 1.0, 1.0, &grads, None,
        );
        let diff: f64 = g_a[0]
            .iter()
            .zip(&g_b[0])
            .map(|(a, b)| (a[0] - b[0]).abs() + (a[1] - b[1]).abs())
            .sum();
        assert!(
            diff > 1e-9,
            "gradient did not respond to the quadratic transform — the colour jump was collapsed \
             to a per-edge constant (Trap 1), diff {diff:e}"
        );
    }

    /// Internal finite-difference check with a quadratic region: the analytic vertex gradient must
    /// match a central difference of E_data when one side carries a spatial colour field.
    #[test]
    fn quadratic_gradient_matches_central_difference() {
        const W: usize = 14;
        const H: usize = 14;
        let coeffs = [
            [0.30, 0.55, 0.70],
            [0.15, -0.10, 0.05],
            [-0.08, 0.12, 0.03],
            [0.02, 0.01, -0.02],
            [0.03, -0.02, 0.01],
            [-0.01, 0.02, 0.02],
        ];
        let colors = vec![
            quad(coeffs, 7.0, 7.0, 7.0), // label 0 = quadratic (left of the +y edge)
            RegionColor::Flat([0.85, 0.80, 0.75]),
        ];

        let scene = |x: f64| -> Vec<Face> {
            vec![
                Face {
                    label: 0,
                    loops: vec![vec![[0.0, 0.0], [x, 0.0], [x, H as f64], [0.0, H as f64]]],
                },
                Face {
                    label: 1,
                    loops: vec![vec![
                        [x, 0.0],
                        [W as f64, 0.0],
                        [W as f64, H as f64],
                        [x, H as f64],
                    ]],
                },
            ]
        };
        let energy = |x: f64, target: &[f64]| -> f64 {
            windowed_data(W, H, 0.0, 0.0, target, &scene(x), &colors, 1.0, 1.0, &[], None).0
        };

        // target = the same scene with the boundary elsewhere, so the minimum is not at x0
        let target = {
            let mut img = vec![1.0_f64; W * H * 3];
            let faces = scene(9.3);
            for f in &faces {
                let cov = crate::coverage::polygon_coverage(&f.loops, W, H);
                for p in 0..W * H {
                    let row = p / W;
                    let col = p % W;
                    let c = colors[f.label].eval(col as f64 + 0.5, row as f64 + 0.5);
                    for ch in 0..3 {
                        img[p * 3 + ch] += cov[p].abs() * (c[ch] - 1.0);
                    }
                }
            }
            for v in img.iter_mut() {
                *v = v.clamp(0.0, 1.0);
            }
            img
        };

        let x0 = 5.37_f64; // off the integer grid (the one-sided-derivative kink)
        // Orientation: traversed +y, the left normal (−dy, dx)/‖d‖ = (−1, 0) points toward the
        // region at smaller x — label 0. The normal points at the RIGHT label, so R = 0, L = 1.
        // Getting this backwards flips the gradient's sign while preserving its magnitude, which
        // is exactly what a first run of this test showed (+7.086 vs −7.086).
        let grads = vec![GradCubic {
            left: 1,
            right: 0,
            pts: vec![[x0, 0.0], [x0, H as f64]],
        }];
        let (_, g) = windowed_data(
            W, H, 0.0, 0.0, &target, &scene(x0), &colors, 1.0, 1.0, &grads, None,
        );
        let analytic = g[0][0][0] + g[0][1][0]; // both endpoints slide together in x

        let h = 1e-6;
        let fd = (energy(x0 + h, &target) - energy(x0 - h, &target)) / (2.0 * h);
        let rel = (analytic - fd).abs() / fd.abs().max(1e-12);
        assert!(
            rel < 1e-4,
            "quadratic analytic {analytic:.9e} vs finite-difference {fd:.9e} (rel {rel:.2e})"
        );
    }
}