rust/vectorhd-core/src/coverage.rs

//! Exact signed-area cell-coverage rasterizer — a line-for-line port of
//! `src/vectorhd/analytic/raster.py` (`polygon_coverage` and the flattening around it).
//!
//! The identity being evaluated, per the Python module docstring:
//!
//! ```text
//! coverage[r, c] = Σ_edges  w · ∫_{y∈[r,r+1]} clamp( x_e(y) − c, 0, 1 ) dy
//! ```
//!
//! **Equivalence discipline.** Every accumulation order here is deliberate, because the S-3 gate
//! compares against Python at 1e-9. Three places matter:
//!
//! 1. `add_piece` accumulates into `partial`/`left` in polygon → edge → scanline → sub-piece
//!    order, matching the Python loop nesting exactly.
//! 2. The final suffix sum runs from column `width` downward, which is precisely what
//!    `np.cumsum(left[:, ::-1], axis=1)[:, ::-1]` does. A left-to-right accumulation would be
//!    algebraically identical and numerically different.
//! 3. `int(xm)` / `int(a + t·d)` in Python truncate toward zero; Rust's `as i64` does the same.
//!    This is load-bearing for negative coordinates, where truncation and floor disagree.

use crate::model::{Edge, Pt, Vectors};
use std::collections::HashMap;

/// Push only if not already present under exact float equality — mirrors Python's `set` of floats.
/// (`-0.0 == 0.0` in both languages, so the two agree on that edge case too.)
fn push_unique(v: &mut Vec<f64>, x: f64) {
    if !v.iter().any(|&y| y == x) {
        v.push(x);
    }
}

/// One monotone edge piece contained in a single pixel column (x already split at integers).
#[inline]
fn add_piece(prow: &mut [f64], lrow: &mut [f64], xa: f64, xb: f64, dyp: f64, w: f64, width: usize) {
    let xm = 0.5 * (xa + xb);
    if xm <= 0.0 {
        // entirely left of the image → no column lies to its left
        return;
    }
    if xm >= width as f64 {
        // entirely right of the image → every column lies to its left
        lrow[width] += w * dyp;
        return;
    }
    let c = xm as usize; // Python `int(xm)`; xm > 0 here, so truncation == floor
    prow[c] += w * dyp * (xm - c as f64); // partial area inside column c
    lrow[c] += w * dyp; // full coverage carried to every column strictly left of c
}

#[allow(clippy::too_many_arguments)]
fn add_subrow(
    prow: &mut [f64],
    lrow: &mut [f64],
    xlo: f64,
    ylo: f64,
    xhi: f64,
    yhi: f64,
    w: f64,
    width: usize,
) {
    let dy_total = yhi - ylo;
    if dy_total <= 0.0 {
        return;
    }
    if (xhi - xlo).abs() < 1e-12 {
        add_piece(prow, lrow, xlo, xhi, dy_total, w, width);
        return;
    }
    let (lo_x, hi_x) = if xlo < xhi { (xlo, xhi) } else { (xhi, xlo) };

    let mut breaks: Vec<f64> = Vec::with_capacity(8);
    push_unique(&mut breaks, xlo);
    push_unique(&mut breaks, xhi);
    let k0 = lo_x.ceil() as i64;
    let k1 = hi_x.floor() as i64;
    for k in k0..=k1 {
        let kf = k as f64;
        if lo_x < kf && kf < hi_x && k > 0 && k < width as i64 {
            push_unique(&mut breaks, kf);
        }
    }
    for b in [0.0_f64, width as f64] {
        if lo_x < b && b < hi_x {
            push_unique(&mut breaks, b);
        }
    }

    let inv = (yhi - ylo) / (xhi - xlo);
    // Python: sorted((y, x) for x in breaks) — ordered by y, x breaking ties.
    let mut pieces: Vec<(f64, f64)> = breaks.iter().map(|&x| (ylo + (x - xlo) * inv, x)).collect();
    pieces.sort_by(|a, b| a.partial_cmp(b).expect("no NaN in break coordinates"));

    for pair in pieces.windows(2) {
        let (ya_, xa_) = pair[0];
        let (yb_, xb_) = pair[1];
        if yb_ > ya_ {
            add_piece(prow, lrow, xa_, xb_, yb_ - ya_, w, width);
        }
    }
}

#[allow(clippy::too_many_arguments)]
fn add_edge(
    partial: &mut [f64],
    left: &mut [f64],
    x0: f64,
    y0: f64,
    x1: f64,
    y1: f64,
    width: usize,
    height: usize,
) {
    if y0 == y1 {
        return; // horizontal edges contribute nothing to a y-integral
    }
    let (w, xa, ya, xb, yb) = if y0 < y1 {
        (1.0, x0, y0, x1, y1)
    } else {
        (-1.0, x1, y1, x0, y0)
    };
    let dxdy = (xb - xa) / (yb - ya);
    let y_top = ya.max(0.0);
    let y_bot = yb.min(height as f64);
    if y_bot <= y_top {
        return;
    }
    let mut r = y_top.floor() as i64;
    while (r as f64) < y_bot {
        let ylo = y_top.max(r as f64);
        let yhi = y_bot.min(r as f64 + 1.0);
        if yhi > ylo {
            let xr0 = xa + (ylo - ya) * dxdy;
            let xr1 = xa + (yhi - ya) * dxdy;
            let ri = r as usize;
            let prow = &mut partial[ri * width..(ri + 1) * width];
            let lrow = &mut left[ri * (width + 1)..(ri + 1) * (width + 1)];
            add_subrow(prow, lrow, xr0, ylo, xr1, yhi, w, width);
        }
        r += 1;
    }
}

/// Exact per-pixel coverage of the region bounded by `loops` (closed polylines), row-major
/// `height × width`. The sign follows loop orientation, as in Python.
pub fn polygon_coverage(loops: &[Vec<Pt>], width: usize, height: usize) -> Vec<f64> {
    let mut partial = vec![0.0_f64; height * width];
    let mut left = vec![0.0_f64; height * (width + 1)];

    for poly in loops {
        let n = poly.len();
        if n < 2 {
            continue;
        }
        for i in 0..n {
            let p0 = poly[i];
            let p1 = poly[(i + 1) % n];
            add_edge(
                &mut partial,
                &mut left,
                p0[0],
                p0[1],
                p1[0],
                p1[1],
                width,
                height,
            );
        }
    }

    // coverage[r, q] = partial[r, q] + Σ_{c>q} left[r, c], accumulated high→low so the float
    // reduction order matches np.cumsum on the reversed axis.
    let mut out = vec![0.0_f64; height * width];
    for r in 0..height {
        let lrow = &left[r * (width + 1)..(r + 1) * (width + 1)];
        let prow = &partial[r * width..(r + 1) * width];
        let orow = &mut out[r * width..(r + 1) * width];
        let mut acc = 0.0_f64;
        for q in (0..width).rev() {
            acc += lrow[q + 1];
            orow[q] = prow[q] + acc;
        }
    }
    out
}

/// numpy's `allclose` defaults: `|a − b| <= atol + rtol·|b|`, elementwise.
fn allclose(a: Pt, b: Pt) -> bool {
    const RTOL: f64 = 1e-5;
    const ATOL: f64 = 1e-8;
    (0..2).all(|i| (a[i] - b[i]).abs() <= ATOL + RTOL * b[i].abs())
}

/// Dense polyline for one edge: first cubic whole, later cubics minus their shared joint —
/// exactly `flatten_edge`. The samples are already fixed-t evaluated by the exporter.
pub fn flatten_edge(e: &Edge) -> Vec<Pt> {
    let mut out: Vec<Pt> = Vec::new();
    for (i, c) in e.cubics.iter().enumerate() {
        if i == 0 {
            out.extend_from_slice(c);
        } else {
            out.extend_from_slice(&c[1..]);
        }
    }
    out
}

/// Closed polyline for one face loop — exactly `flatten_loop`.
pub fn flatten_loop(darts: &[[i64; 2]], by_id: &HashMap<i64, &Edge>) -> Vec<Pt> {
    let mut out: Vec<Pt> = Vec::new();
    for d in darts {
        let e = by_id
            .get(&d[0])
            .unwrap_or_else(|| panic!("loop references edge {} absent from the export", d[0]));
        let mut p = flatten_edge(e);
        if d[1] < 0 {
            p.reverse();
        }
        if out.is_empty() {
            out.extend_from_slice(&p);
        } else {
            out.extend_from_slice(&p[1..]);
        }
    }
    if out.len() > 1 && allclose(out[0], out[out.len() - 1]) {
        out.pop();
    }
    out
}

/// Per-interior-region absolute coverage, in **face order** — the same order Python's dict
/// preserves, which fixes the reduction order of the composite in `energy::compose`.
pub fn region_coverages(v: &Vectors) -> Vec<(i64, Vec<f64>)> {
    let by_id = v.edge_by_id();
    let mut out = Vec::with_capacity(v.faces.len());
    for f in &v.faces {
        let loops: Vec<Vec<Pt>> = f
            .loops
            .iter()
            .map(|l| flatten_loop(l, &by_id))
            .filter(|p| p.len() >= 3)
            .collect();
        if loops.is_empty() {
            continue;
        }
        let mut cov = polygon_coverage(&loops, v.width, v.height);
        for x in cov.iter_mut() {
            *x = x.abs();
        }
        out.push((f.label, cov));
    }
    out
}

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

    /// Mirrors the Python R-04 gate: coverage exact to machine precision vs closed-form area.
    #[test]
    fn axis_aligned_rectangle_is_exact() {
        // Rectangle [2.25, 7.75] x [1.5, 6.5] — area 5.5 * 5.0 = 27.5
        let poly = vec![[2.25, 1.5], [7.75, 1.5], [7.75, 6.5], [2.25, 6.5]];
        let cov = polygon_coverage(&[poly], 10, 10);
        let total: f64 = cov.iter().map(|x| x.abs()).sum();
        assert!(
            (total - 27.5).abs() < 1e-12,
            "rectangle area {total} != 27.5"
        );
    }

    #[test]
    fn rectangle_interior_pixels_are_fully_covered() {
        let poly = vec![[2.0, 2.0], [8.0, 2.0], [8.0, 8.0], [2.0, 8.0]];
        let cov = polygon_coverage(&[poly], 10, 10);
        for r in 2..8 {
            for c in 2..8 {
                assert!(
                    (cov[r * 10 + c].abs() - 1.0).abs() < 1e-12,
                    "interior pixel ({r},{c}) = {}",
                    cov[r * 10 + c]
                );
            }
        }
        // and pixels outside are empty
        assert!(cov[0].abs() < 1e-12);
    }

    /// Half-open pixel: a rectangle covering exactly half of one column.
    #[test]
    fn partial_column_is_exact() {
        let poly = vec![[0.0, 0.0], [1.5, 0.0], [1.5, 1.0], [0.0, 1.0]];
        let cov = polygon_coverage(&[poly], 4, 1);
        assert!((cov[0].abs() - 1.0).abs() < 1e-12, "col0 {}", cov[0]);
        assert!((cov[1].abs() - 0.5).abs() < 1e-12, "col1 {}", cov[1]);
        assert!(cov[2].abs() < 1e-12);
    }

    /// Disk coverage vs the analytic circle area — the R-1 disk gate.
    #[test]
    fn disk_area_matches_closed_form() {
        let (cx, cy, rad) = (32.0_f64, 32.0_f64, 20.0_f64);
        let n = 4096;
        let poly: Vec<Pt> = (0..n)
            .map(|i| {
                let t = 2.0 * std::f64::consts::PI * (i as f64) / (n as f64);
                [cx + rad * t.cos(), cy + rad * t.sin()]
            })
            .collect();
        let cov = polygon_coverage(&[poly], 64, 64);
        let total: f64 = cov.iter().map(|x| x.abs()).sum();
        let exact = std::f64::consts::PI * rad * rad;
        // The polygon is inscribed, so it under-covers by the sagitta area; at n=4096 that is
        // ~1e-6 relative. The gate is that the rasterizer adds no error of its own beyond it.
        let rel = (total - exact).abs() / exact;
        assert!(rel < 1e-5, "disk area {total} vs {exact} (rel {rel:.2e})");
    }

    /// The partition property: coverages of complementary regions sum to 1 per pixel.
    #[test]
    fn complementary_regions_partition_to_one() {
        // Left half and right half of a 6x4 raster, split at x = 2.4.
        let left = vec![[0.0, 0.0], [2.4, 0.0], [2.4, 4.0], [0.0, 4.0]];
        let right = vec![[2.4, 0.0], [6.0, 0.0], [6.0, 4.0], [2.4, 4.0]];
        let a = polygon_coverage(&[left], 6, 4);
        let b = polygon_coverage(&[right], 6, 4);
        for i in 0..a.len() {
            let s = a[i].abs() + b[i].abs();
            assert!((s - 1.0).abs() < 1e-12, "pixel {i} partition = {s}");
        }
    }

    #[test]
    fn winding_sign_flips_with_orientation() {
        let ccw = vec![[1.0, 1.0], [3.0, 1.0], [3.0, 3.0], [1.0, 3.0]];
        let cw: Vec<Pt> = ccw.iter().rev().copied().collect();
        let a = polygon_coverage(&[ccw], 5, 5);
        let b = polygon_coverage(&[cw], 5, 5);
        for i in 0..a.len() {
            assert!((a[i] + b[i]).abs() < 1e-12, "signs did not mirror at {i}");
        }
    }

    #[test]
    fn geometry_outside_the_raster_is_clipped_not_wrapped() {
        // A box hanging off every side; only the on-raster part counts.
        let poly = vec![[-5.0, -5.0], [3.0, -5.0], [3.0, 3.0], [-5.0, 3.0]];
        let cov = polygon_coverage(&[poly], 4, 4);
        let total: f64 = cov.iter().map(|x| x.abs()).sum();
        assert!((total - 9.0).abs() < 1e-12, "clipped area {total} != 9");
    }
}