main
Caleb Version 0.3.0: new copper reader, dome via filter, every layer read from the photos, net evidence, measured confidence;… e874d7d
//! Copper from the photos. Copper under solder mask cannot be told from bare laminate by colour in a
//! single photo (measured on the Driver bottom, 2026-11-03: the two overlap in hue, saturation and
//! lightness once the lighting varies). What the photo does show reliably is the edges: a pour's
//! outline is a narrow dark clearance gap with a light ridge beside it. So the AI sketches each pour
//! or trace roughly, in board millimetres, on the gridded board image, and this module fits the sketch
//! to the photo: every edge moves to the strongest gap line within reach, corners are re-solved as edge
//! intersections, and the share of the outline that sits on a measured gap is reported as the evidence.
use crate::board::{to_mm, to_px, RES};
use crate::geom::{self, P};
use crate::raster::{self, Plane};
use image::{Rgb, RgbImage};
use serde_json::{json, Value};

/// Grey-level max filter with a square window of half-size r (separable).
fn max_filter(p: &Plane, r: usize, min: bool) -> Plane {
    let (w, h) = (p.w, p.h);
    let pick = |a: f32, b: f32| if min { a.min(b) } else { a.max(b) };
    let mut t = vec![0f32; w * h];
    for y in 0..h { for x in 0..w {
        let (a, b) = (x.saturating_sub(r), (x + r).min(w - 1));
        let mut m = p.d[y * w + a]; for k in a + 1..=b { m = pick(m, p.d[y * w + k]); } t[y * w + x] = m;
    } }
    let mut o = vec![0f32; w * h];
    for y in 0..h { for x in 0..w {
        let (a, b) = (y.saturating_sub(r), (y + r).min(h - 1));
        let mut m = t[a * w + x]; for k in a + 1..=b { m = pick(m, t[k * w + x]); } o[y * w + x] = m;
    } }
    Plane { w, h, d: o }
}

/// The gap map: a black-hat (closing minus image) of the blurred lightness, about 0.4 mm wide, so
/// narrow dark lines (clearance gaps, trace edges) stand out and broad areas do not.
pub fn gap_map(img: &RgbImage) -> Plane {
    let l = raster::blur(&raster::luma(img), 2, 3);
    let r = (0.2 * RES) as usize;
    let closed = max_filter(&max_filter(&l, r, false), r, true);
    Plane { w: l.w, h: l.h, d: closed.d.iter().zip(&l.d).map(|(c, v)| c - v).collect() }
}

/// The board image with a labelled millimetre grid every `step` mm (minor lines every mm), for sketching.
pub fn grid_image(img: &RgbImage, step: f64) -> RgbImage {
    let mut g = img.clone();
    let (w, h) = (g.width() as f64, g.height() as f64);
    let (x1, y1) = to_mm(w, h);
    let mut x = (crate::board::ORIGIN).ceil();
    while x <= x1 {
        let u = to_px(x, 0.0).0; let major = ((x - crate::board::ORIGIN) % step).abs() < 1e-6;
        if major { crate::draw::line(&mut g, (u, 0.0), (u, h), [255, 255, 255], 1); crate::draw::label(&mut g, (u + 3.0, 4.0), &format!("{}", x as i64), [255, 255, 0], 3); }
        else { for v in (0..h as i64).step_by(6) { crate::draw::px(&mut g, u as i64, v, [200, 200, 200]); } }
        x += 1.0;
    }
    let mut y = (crate::board::ORIGIN).ceil();
    while y <= y1 {
        let v = to_px(0.0, y).1; let major = ((y - crate::board::ORIGIN) % step).abs() < 1e-6;
        if major { crate::draw::line(&mut g, (0.0, v), (w, v), [255, 255, 255], 1); crate::draw::label(&mut g, (4.0, v + 3.0), &format!("{}", y as i64), [255, 255, 0], 3); }
        else { for u in (0..w as i64).step_by(6) { crate::draw::px(&mut g, u, v as i64, [200, 200, 200]); } }
        y += 1.0;
    }
    g
}

/// Fit one sketched polygon (board mm) to the gap map. Returns the fitted polygon (mm), the share of
/// sampled edge points that found a gap (0..1) and the largest corner move (mm).
pub fn snap_polygon(gap: &Plane, poly_mm: &[P], reach_mm: f64, strength: f32) -> (Vec<P>, f64, f64) { snap_polygon_in(gap, poly_mm, reach_mm, strength, None) }

/// As `snap_polygon`, with the board box: an edge lying along the board outline has no gap line to
/// find, so it is kept as sketched (clipped to the board later) and left out of the share.
pub fn snap_polygon_in(gap: &Plane, poly_mm: &[P], reach_mm: f64, strength: f32, board: Option<(f64, f64, f64, f64)>) -> (Vec<P>, f64, f64) {
    let n = poly_mm.len();
    let pts: Vec<P> = poly_mm.iter().map(|p| to_px(p.0, p.1)).collect();
    let reach = reach_mm * RES;
    let mut lines: Vec<Option<(P, P)>> = vec![]; // fitted edge as (point, unit direction)
    // which side of each edge is inside: the left normal points in for a counter-clockwise ring (y down)
    let signed: f64 = (0..n).map(|i| { let (a, b) = (pts[i], pts[(i + 1) % n]); a.0 * b.1 - b.0 * a.1 }).sum();
    let inward = if signed > 0.0 { 1.0 } else { -1.0 };
    let (mut hit, mut tried) = (0usize, 0usize);
    for i in 0..n {
        let (a, b) = (pts[i], pts[(i + 1) % n]);
        let len = geom::dist(a, b); if len < 1.0 { lines.push(None); continue; }
        if let Some((x0, y0, x1, y1)) = board {
            let on = |p: P| { let q = to_mm(p.0, p.1); (q.0 - x0).abs() < 0.6 || (q.0 - x1).abs() < 0.6 || (q.1 - y0).abs() < 0.6 || (q.1 - y1).abs() < 0.6 };
            let (am, bm) = (to_mm(a.0, a.1), to_mm(b.0, b.1));
            let same_side = ((am.0 - bm.0).abs() < 0.6 && ((am.0 - x0).abs() < 0.6 || (am.0 - x1).abs() < 0.6)) || ((am.1 - bm.1).abs() < 0.6 && ((am.1 - y0).abs() < 0.6 || (am.1 - y1).abs() < 0.6));
            if on(a) && on(b) && same_side { lines.push(None); continue; }
        }
        let d = ((b.0 - a.0) / len, (b.1 - a.1) / len); let nrm = (-d.1, d.0);
        let mut found: Vec<(f64, f64)> = vec![]; // (along, offset)
        let samples = ((len / (0.25 * RES)) as usize).max(4);
        for k in 1..samples { // skip the corners themselves
            let t = k as f64 / samples as f64 * len;
            let base = (a.0 + d.0 * t, a.1 + d.1 * t);
            // score each offset by the gap averaged along the edge (+-0.4 mm): a continuous line beats the
            // dimples of a via field; a small penalty per mm keeps the fit near the sketch on a tie
            let along = (0.4 * RES) as i64;
            let (mut best, mut bo, mut braw) = (f32::MIN, 0f64, 0f32);
            let mut o = -reach;
            while o <= reach {
                let (mut s, mut n) = (0f32, 0f32);
                for k in (-along..=along).step_by(2) {
                    let q = (base.0 + nrm.0 * o + d.0 * k as f64, base.1 + nrm.1 * o + d.1 * k as f64);
                    if q.0 >= 0.0 && q.1 >= 0.0 && q.0 < gap.w as f64 && q.1 < gap.h as f64 { s += gap.sample(q.0 as f32, q.1 as f32); n += 1.0; }
                }
                if n > 0.0 {
                    let g = s / n;
                    // a channel between two pours has two edges: the pour's own is on its inside, so an
                    // outward move costs three times an inward one
                    let outward = o * inward < 0.0;
                    let score = g - strength * (if outward { 0.45 } else { 0.15 }) * (o.abs() / RES) as f32;
                    if score > best { best = score; bo = o; braw = g; }
                }
                o += 1.0;
            }
            tried += 1;
            if braw >= strength * 0.7 { hit += 1; found.push((t, bo)); }
        }
        // a robust straight fit of offset against position along the edge
        let line = if found.len() >= 3 { geom::robust_line(&found, 0.06 * RES) } else { None };
        lines.push(line.map(|(m, c)| {
            let p0 = (a.0 + nrm.0 * c, a.1 + nrm.1 * c);
            let dir = (d.0 + nrm.0 * m, d.1 + nrm.1 * m); let l = (dir.0 * dir.0 + dir.1 * dir.1).sqrt();
            (p0, (dir.0 / l, dir.1 / l))
        }));
    }
    // corners: intersect neighbouring fitted edges; keep the sketch where an edge found nothing
    let mut out = vec![];
    let mut moved: f64 = 0.0;
    for i in 0..n {
        let prev = &lines[(i + n - 1) % n]; let cur = &lines[i];
        let p = match (prev, cur) {
            (Some((p1, d1)), Some((p2, d2))) => {
                let den = d1.0 * d2.1 - d1.1 * d2.0;
                if den.abs() < 0.05 { pts[i] } else { let t = ((p2.0 - p1.0) * d2.1 - (p2.1 - p1.1) * d2.0) / den; (p1.0 + d1.0 * t, p1.1 + d1.1 * t) }
            }
            (Some((p0, d0)), None) | (None, Some((p0, d0))) => { // project onto the one fitted edge
                let v = (pts[i].0 - p0.0, pts[i].1 - p0.1); let t = v.0 * d0.0 + v.1 * d0.1; (p0.0 + d0.0 * t, p0.1 + d0.1 * t)
            }
            _ => pts[i],
        };
        let p = if geom::dist(p, pts[i]) > 2.0 * reach { pts[i] } else { p };
        moved = moved.max(geom::dist(p, pts[i]) / RES);
        out.push(to_mm(p.0, p.1));
    }
    (out.iter().map(|p| ((p.0 * 1000.0).round() / 1000.0, (p.1 * 1000.0).round() / 1000.0)).collect(), if tried > 0 { hit as f64 / tried as f64 } else { 0.0 }, moved)
}

/// The gap strength that counts as an edge: a high percentile of the gap map over the board.
pub fn edge_strength(gap: &Plane) -> f32 {
    let mut v: Vec<f32> = gap.d.iter().step_by(7).copied().collect();
    v.sort_by(|a, b| a.partial_cmp(b).unwrap());
    v[(v.len() as f64 * 0.93) as usize]
}

pub fn draw_poly_mm(img: &mut RgbImage, p: &[P], c: [u8; 3], th: i64) {
    let px: Vec<P> = p.iter().map(|q| to_px(q.0, q.1)).collect();
    crate::draw::poly(img, &px, c, th);
}

pub fn gap_image(g: &Plane) -> RgbImage {
    let mut v: Vec<f32> = g.d.iter().step_by(7).copied().collect(); v.sort_by(|a, b| a.partial_cmp(b).unwrap());
    let hi = v[(v.len() as f64 * 0.99) as usize].max(1.0);
    RgbImage::from_fn(g.w as u32, g.h as u32, |x, y| { let t = (g.at(x as usize, y as usize) / hi * 255.0).clamp(0.0, 255.0) as u8; Rgb([t, t, t]) })
}

#[cfg(test)]
mod tests {
    use super::*;
    #[test]
    fn a_rough_sketch_snaps_to_the_drawn_pour() {
        // a light pour 20 x 12 mm with a 0.3 mm dark clearance gap around it, on a mid background
        let (w, h) = ((40.0 * RES) as u32, (30.0 * RES) as u32);
        let pour = [(105.0, 105.0), (125.0, 105.0), (125.0, 117.0), (105.0, 117.0)];
        let img = RgbImage::from_fn(w, h, |x, y| {
            let (mx, my) = to_mm(x as f64 + 0.5, y as f64 + 0.5);
            let inside = mx > 105.0 && mx < 125.0 && my > 105.0 && my < 117.0;
            let near = mx > 104.7 && mx < 125.3 && my > 104.7 && my < 117.3;
            if inside { Rgb([60, 120, 90]) } else if near { Rgb([20, 50, 35]) } else { Rgb([55, 110, 85]) }
        });
        let g = gap_map(&img);
        let sketch = [(105.6, 104.5), (124.4, 105.5), (125.5, 116.4), (104.6, 117.6)];
        let (fit, share, moved) = snap_polygon(&g, &sketch, 1.0, edge_strength(&g));
        for i in 0..4 { assert!(geom::dist(fit[i], (pour[i].0 + if i == 1 || i == 2 { 0.15 } else { -0.15 }, pour[i].1 + if i >= 2 { 0.15 } else { -0.15 })) < 0.2, "{i}: {:?}", fit[i]); }
        assert!(share > 0.8, "{share}"); assert!(moved < 1.5);
    }
}

/// Running min or max over a window of `len` pixels along one direction (dx, dy in {-1,0,1}), van Herk /
/// Gil-Werman: constant work per pixel whatever the window length.
pub fn line_filter(p: &Plane, dx: i64, dy: i64, len: usize, max: bool) -> Plane {
    let (w, h) = (p.w as i64, p.h as i64);
    let mut out = p.clone();
    let pick = |a: f32, b: f32| if max { a.max(b) } else { a.min(b) };
    let k = len.max(1); let half = (k / 2) as i64;
    // every line in this direction starts on the image border
    let mut starts = vec![];
    if dy == 0 { for y in 0..h { starts.push((0, y)); } }
    else if dx == 0 { for x in 0..w { starts.push((x, 0)); } }
    else if dx == 1 { for x in 0..w { starts.push((x, 0)); } for y in 1..h { starts.push((0, y)); } }
    else { for x in 0..w { starts.push((x, 0)); } for y in 1..h { starts.push((w - 1, y)); } }
    let mut vals: Vec<f32> = vec![]; let mut idx: Vec<usize> = vec![];
    for (sx, sy) in starts {
        vals.clear(); idx.clear();
        let (mut x, mut y) = (sx, sy);
        while x >= 0 && y >= 0 && x < w && y < h { let i = (y * w + x) as usize; vals.push(p.d[i]); idx.push(i); x += dx; y += dy; }
        let n = vals.len(); if n == 0 { continue; }
        let mut g = vec![0f32; n]; let mut hh = vec![0f32; n];
        for i in 0..n { g[i] = if i % k == 0 { vals[i] } else { pick(g[i - 1], vals[i]) }; }
        for i in (0..n).rev() { hh[i] = if i == n - 1 || (i + 1) % k == 0 { vals[i] } else { pick(hh[i + 1], vals[i]) }; }
        for i in 0..n {
            let lo = (i as i64 - half).max(0) as usize; let hi = ((i as i64 - half) + k as i64 - 1).min(n as i64 - 1) as usize;
            // window [lo, hi]: combine the suffix of lo's block and the prefix of hi's block
            let v = if lo / k == hi / k { let mut m = vals[lo]; for j in lo + 1..=hi { m = pick(m, vals[j]); } m } else { pick(hh[lo], g[hi]) };
            out.d[idx[i]] = v;
        }
    }
    out
}

fn bin_line(m: &crate::raster::Mask, dx: i64, dy: i64, len: usize, close: bool) -> crate::raster::Mask {
    let p = Plane { w: m.w, h: m.h, d: m.d.iter().map(|b| if *b { 1.0 } else { 0.0 }).collect() };
    let (a, b) = if close { (true, false) } else { (false, true) };
    let q = line_filter(&line_filter(&p, dx, dy, len, a), dx, dy, len, b);
    crate::raster::Mask { w: m.w, h: m.h, d: q.d.iter().map(|v| *v > 0.5).collect() }
}

/// Pour boundaries as a mask: the black-hat's long straight dark lines (0, 45, 90, 135 degrees: pours on
/// boards are drawn that way; round via dimples are not lines), each closed along its own direction to
/// bridge small breaks.
pub fn line_barrier(img: &RgbImage) -> crate::raster::Mask {
    // two features: thin dark gaps (black-hat) and the step from a pour to a wide bare channel (the
    // brightness gradient across each direction); each kept only where it runs straight for 1.6 mm,
    // then closed along its own direction over 1.2 mm breaks
    let g = gap_map(img);
    let l = raster::blur(&raster::luma(img), 2, 3);
    let (w, h) = (l.w, l.h);
    let at = |x: usize, y: usize| l.d[y.min(h - 1) * w + x.min(w - 1)];
    let mut gx = Plane::new(w, h, 0.0); let mut gy = Plane::new(w, h, 0.0);
    for y in 1..h - 1 { for x in 1..w - 1 {
        gx.d[y * w + x] = (at(x + 1, y - 1) + 2.0 * at(x + 1, y) + at(x + 1, y + 1)) - (at(x - 1, y - 1) + 2.0 * at(x - 1, y) + at(x - 1, y + 1));
        gy.d[y * w + x] = (at(x - 1, y + 1) + 2.0 * at(x, y + 1) + at(x + 1, y + 1)) - (at(x - 1, y - 1) + 2.0 * at(x, y - 1) + at(x + 1, y - 1));
    } }
    // (direction of the line, response across it)
    let resp = |k: usize| -> Plane { Plane { w, h, d: (0..w * h).map(|i| match k { 0 => gy.d[i].abs(), 1 => gx.d[i].abs(), 2 => (gx.d[i] - gy.d[i]).abs() / 1.414, _ => (gx.d[i] + gy.d[i]).abs() / 1.414 }).collect() } };
    let dirs = [(1i64, 0i64), (0, 1), (1, 1), (-1, 1)];
    let open_len = (1.6 * RES) as usize; let close_len = (1.2 * RES) as usize;
    let mut bar = crate::raster::Mask::new(w, h);
    for (k, (dx, dy)) in dirs.iter().enumerate() {
        for src in [resp(k), g.clone()] {
            let o = line_filter(&line_filter(&src, *dx, *dy, open_len, false), *dx, *dy, open_len, true);
            let mut s: Vec<f32> = o.d.iter().step_by(11).copied().collect(); s.sort_by(|a, b| a.partial_cmp(b).unwrap());
            let thr = s[(s.len() as f64 * 0.93) as usize].max(4.0);
            let b = crate::raster::Mask { w, h, d: o.d.iter().map(|v| *v > thr).collect() };
            bar = bar.or(&bin_line(&b, *dx, *dy, close_len, true));
        }
    }
    bar.dilate(1.5)
}

/// Tented vias on this face, as board-mm positions. See `domes`.
pub fn domes(img: &RgbImage) -> Vec<(f64, f64)> { crate::domes::finish(&crate::domes::measure(img)).iter().map(|d| to_mm(d.x, d.y)).collect() }

/// Candidate pours: the regions the line barrier encloses, with outline, area and via density.
pub fn regions(img: &RgbImage, side: &str) -> (Vec<Value>, Vec<u32>, usize, usize) {
    let bar = line_barrier(img);
    let free = bar.not();
    // regions joined only through a neck under ~1 mm (a break in a boundary line) are split there
    let (lab, cs) = crate::prep::split_components(&free, (0.5 * RES) as f32);
    let ds = domes(img);
    let w = free.w;
    let mut dens = std::collections::HashMap::<u32, usize>::new();
    for (x, y) in &ds { let (u, v) = to_px(*x, *y); let (u, v) = (u as usize, v as usize); if u < w && v < free.h { *dens.entry(lab[v * w + u]).or_default() += 1; } }
    let l = raster::luma(img);
    let mut out = vec![];
    for c in cs.iter().filter(|c| c.area as f64 > 2.0 * RES * RES && (c.area as f64) < 0.6 * (free.w * free.h) as f64) {
        let ring = geom::simplify(&raster::trace_boundary(&lab, w, free.h, c), 2.5);
        if ring.len() < 3 { continue; }
        let mut sum = 0f64; for y in c.y0..=c.y1 { for x in c.x0..=c.x1 { if lab[y * w + x] == c.label { sum += l.at(x, y) as f64; } } }
        let area_mm = c.area as f64 / (RES * RES);
        let nd = dens.get(&c.label).copied().unwrap_or(0);
        let cm = to_mm(c.cx, c.cy);
        let poly: Vec<Value> = ring.iter().map(|p| { let m = to_mm(p.0, p.1); json!([(m.0 * 1000.0).round() / 1000.0, (m.1 * 1000.0).round() / 1000.0]) }).collect();
        out.push(json!({"id": format!("{side}:{}", c.label), "side": side, "area_mm2": (area_mm * 10.0).round() / 10.0, "centre_mm": [(cm.0 * 100.0).round() / 100.0, (cm.1 * 100.0).round() / 100.0],
                        "domes": nd, "domesPerMm2": ((nd as f64 / area_mm) * 1000.0).round() / 1000.0, "luma": (sum / c.area as f64 * 10.0).round() / 10.0,
                        "polygon_mm": poly}));
    }
    (out, lab, w, ds.len())
}