AI Flow: Reverse Engineering

//! `check fidelity`: what KiCad draws against what the photos show, layer by layer, at 40 px/mm. It uses the photos and
//! the board and nothing else. Silk, mask openings, exposed metal, holes and outer copper as seen under the mask each
//! get IoU and boundary distance over the part of the board the photo can show that layer on. The rest is
//! excluded, its size is reported, and it counts for nothing either way.
//!
//! Three samples: tile set A (5 mm checkerboard, used to set thresholds), tile set B (held out), and the second photo
//! of each face on tile set B. The class map is the same reader the tracer uses, so a good IoU on the primary photo
//! alone would be circular: the second photo and the threshold-free edge test exist because of that. Every layer is
//! also scored shifted by 0.25 mm, which must score worse or the check fails.
use crate::gerber::{self, Frame};
use crate::photolayers;
use crate::photomat::{self, Face};
use crate::raster::{self, Mask, Plane};
use crate::{board, geom};
use serde_json::{json, Map, Value};
use std::path::Path;

const TILE: usize = 200; // 5 mm
fn tile_b(u: usize, v: usize) -> bool { (u / TILE + v / TILE) % 2 == 1 }
fn r3(x: f64) -> f64 { (x * 1000.0).round() / 1000.0 }
fn r2(x: f64) -> f64 { (x * 100.0).round() / 100.0 }

struct Scored { iou: f64, precision: f64, recall: f64, truth_mm2: f64, board_mm2: f64, scored_mm2: f64 }

fn score(render: &Mask, truth: &Mask, scored: &Mask, sel: &dyn Fn(usize, usize) -> bool) -> Scored {
    let (w, h) = (render.w, render.h);
    let (mut i, mut u, mut t, mut r, mut s) = (0u64, 0u64, 0u64, 0u64, 0u64);
    for y in 0..h { for x in 0..w { if !sel(x, y) { continue; } let k = y * w + x; if !scored.d[k] { continue; }
        s += 1; let (a, b) = (render.d[k], truth.d[k]); if a && b { i += 1; } if a || b { u += 1; } if b { t += 1; } if a { r += 1; } } }
    Scored { iou: if u > 0 { i as f64 / u as f64 } else { 1.0 }, precision: if r > 0 { i as f64 / r as f64 } else { 0.0 }, recall: if t > 0 { i as f64 / t as f64 } else { 0.0 }, truth_mm2: t as f64 / photomat::MM2, board_mm2: r as f64 / photomat::MM2, scored_mm2: s as f64 / photomat::MM2 }
}

fn shifted_iou(render: &Mask, truth: &Mask, scored: &Mask, dx: i64, dy: i64) -> f64 {
    let (w, h) = (render.w as i64, render.h as i64);
    let (mut i, mut u) = (0u64, 0u64);
    for y in 0..h { for x in 0..w { let k = (y * w + x) as usize; if !scored.d[k] || !tile_b(x as usize, y as usize) { continue; }
        let (sx, sy) = (x - dx, y - dy);
        let a = sx >= 0 && sy >= 0 && sx < w && sy < h && render.d[(sy * w + sx) as usize];
        let b = truth.d[k]; if a && b { i += 1; } if a || b { u += 1; } } }
    if u > 0 { i as f64 / u as f64 } else { 1.0 }
}

fn boundary(m: &Mask) -> Mask {
    let (w, h) = (m.w, m.h);
    let mut o = Mask::new(w, h);
    for y in 1..h - 1 { for x in 1..w - 1 { let k = y * w + x; if m.d[k] && (!m.d[k - 1] || !m.d[k + 1] || !m.d[k - w] || !m.d[k + w]) { o.d[k] = true; } } }
    o
}

/// Symmetric boundary distance over the scored tiles: median and 95th percentile in mm.
fn boundary_dist(render: &Mask, truth: &Mask, scored: &Mask) -> (f64, f64, usize) {
    let (br, bt) = (boundary(render), boundary(truth));
    if br.count() == 0 || bt.count() == 0 { return (f64::NAN, f64::NAN, 0); }
    let (dr, dt) = (br.distance_to_set(), bt.distance_to_set());
    let inner = scored.erode(2.0);
    let mut v = vec![];
    for y in 0..render.h { for x in 0..render.w { if !tile_b(x, y) { continue; } let k = y * render.w + x; if !inner.d[k] { continue; } if br.d[k] { v.push(dt.d[k]); } if bt.d[k] { v.push(dr.d[k]); } } }
    if v.is_empty() { return (0.0, 0.0, 0); }
    v.sort_by(|a, b| a.partial_cmp(b).unwrap());
    (r3(v[v.len() / 2] as f64 / board::RES), r3(v[(v.len() as f64 * 0.95) as usize] as f64 / board::RES), v.len())
}

/// Threshold-free: from each boundary point of the drawn layer, walk along its normal to the steepest luminance
/// edge within 0.2 mm. The offset says whether the drawn edge sits on the photographed edge, whatever threshold made it.
fn edge_offset(render: &Mask, fc: &Face, scored: &Mask, inside_bright: bool) -> Value {
    let b = boundary(render);
    let bl = raster::blur(&Plane { w: render.w, h: render.h, d: render.d.iter().map(|x| if *x { 1.0 } else { 0.0 }).collect() }, 2, 1);
    let inner = scored.erode(10.0);
    let (w, h) = (render.w, render.h);
    let mut v = vec![];
    for y in (10..h - 10).step_by(2) { for x in (10..w - 10).step_by(2) {
        let k = y * w + x; if !b.d[k] || !inner.d[k] || !tile_b(x, y) { continue; }
        // the normal points outward: mask falls that way
        let (gx, gy) = (bl.d[k - 1] - bl.d[k + 1], bl.d[k - w] - bl.d[k + w]);
        let n = (gx * gx + gy * gy).sqrt(); if n < 1e-3 { continue; }
        let (nx, ny) = (gx / n, gy / n);
        let prof = |s: f32| fc.y.sample(x as f32 + nx * s, y as f32 + ny * s);
        let (mut best, mut at) = (0f32, 0f32);
        let mut s = -8.0f32;
        while s <= 8.0 { let d = prof(s - 1.0) - prof(s + 1.0); let d = if inside_bright { d } else { -d }; if d > best { best = d; at = s; } s += 0.5; }
        if best > 12.0 { v.push(at as f64 / board::RES); }
    } }
    if v.is_empty() { return json!({"samples": 0}); }
    v.sort_by(|a, b| a.partial_cmp(b).unwrap());
    let med = v[v.len() / 2];
    let mut abs: Vec<f64> = v.iter().map(|x| x.abs()).collect(); abs.sort_by(|a, b| a.partial_cmp(b).unwrap());
    json!({"samples": v.len(), "medianSignedMm": r3(med), "medianAbsMm": r3(abs[abs.len() / 2]), "p95AbsMm": r3(abs[(abs.len() as f64 * 0.95) as usize]), "meaning": "positive: the photographed edge lies outside the drawn shape"})
}

fn class_mask(fc: &Face, f: impl Fn(u8) -> bool) -> Mask { Mask { w: fc.w, h: fc.h, d: fc.cls.iter().map(|c| f(*c)).collect() } }

/// Pixels this close to anything but mask do not show the surface under the mask: within it the camera's sharpening dips and glare
/// lifts. On the Driver top the median green of mask pixels sits 3 to 7 percent off its plateau out to 8 px from a bright part, and laminate
/// reads 10 percent high out to 12 px, with 6 px the point of the dip's recovery. From these photos, tile set A.
pub const BESIDE_PX: f32 = 6.0;
const CLOSE_PX: usize = 24;

/// Copper under mask, read independently of the islands the board was drawn from. Mask over copper is the lighter of the
/// two surfaces. The light across the board is taken from the upper envelope of the green level (85th percentile of plain mask
/// pixels in 6 mm windows), the ratio to it is split in two by one Otsu threshold for the whole face, and specks under
/// 0.1 mm go. Nothing here looks at gap lines, so a pour the tracer missed shows up as photo-only. The envelope and the split use only
/// pixels `BESIDE_PX` clear of anything else: with them in, the glare beside silk blocks and pads set the envelope at 107 where the copper
/// read 64 (Driver top, MOTOR_TEMP box), and whole pours read as laminate.
pub fn copper_read(img: &image::RgbImage, fc: &Face) -> Mask {
    let (w, h) = (fc.w, fc.h);
    // a narrow gap's dip must survive this blur: at radius 3 the 5 px gaps between tracks read part copper, at radius 2 the read edges sit closer to the steepest fall of brightness
    // (share within 2 px on tile set A, Driver top 0.625 to 0.641, bottom 0.482 to 0.565) for a small loss of agreement between the two shots (top 0.949 to 0.937, bottom 0.954 to 0.953)
    let g = raster::blur(&Plane { w, h, d: img.pixels().map(|p| p[1] as f32).collect() }, 2, 2);
    let g_open = g.clone();
    // A dome's shadow is a small dark blob on copper, not a gap. A gap is long. Closing along a line fills a dark feature shorter than the line in that direction, so the lowest of four
    // directional closings fills blobs and keeps long lines. 0.6 mm is the longest dome shadow on the Driver bottom, tile set A.
    let g = { let closed: Vec<Plane> = [(1i64, 0i64), (0, 1), (1, 1), (-1, 1)].iter().map(|(dx, dy)| { let up = crate::copper::line_filter(&g, *dx, *dy, CLOSE_PX, true); crate::copper::line_filter(&up, *dx, *dy, CLOSE_PX, false) }).collect(); Plane { w, h, d: (0..w * h).map(|i| closed.iter().map(|c| c.d[i]).fold(f32::MAX, f32::min)).collect() } };
    let plain = Mask { w, h, d: fc.cls.iter().map(|c| *c == photomat::MASK).collect() }.erode(BESIDE_PX);
    let step = 120usize; let win = 240usize;
    let (cols, rows) = (w / step + 2, h / step + 2);
    let mut env = vec![f32::NAN; cols * rows];
    for r in 0..rows { for c in 0..cols {
        let (cx, cy) = (c * step, r * step);
        let (x0, x1, y0, y1) = (cx.saturating_sub(win / 2), (cx + win / 2).min(w), cy.saturating_sub(win / 2), (cy + win / 2).min(h));
        let mut v = vec![];
        for y in y0..y1 { for x in x0..x1 { let k = y * w + x; if plain.d[k] { v.push(g.d[k]); } } }
        if v.len() < 500 { continue; }
        let k = (v.len() as f64 * 0.85) as usize; let (_, e, _) = v.select_nth_unstable_by(k, |a, b| a.partial_cmp(b).unwrap());
        env[r * cols + c] = *e;
    } }
    let mut known: Vec<f32> = env.iter().copied().filter(|v| !v.is_nan()).collect();
    if known.is_empty() { return Mask::new(w, h); }
    known.sort_by(|a, b| a.partial_cmp(b).unwrap());
    let med = known[known.len() / 2];
    for v in env.iter_mut() { if v.is_nan() { *v = med; } }
    let mut full = Plane::new(w, h, 0.0);
    for y in 0..h { for x in 0..w {
        let (fx, fy) = (x as f32 / step as f32, y as f32 / step as f32);
        let (c0, r0) = (fx as usize, fy as usize); let (tx, ty) = (fx - c0 as f32, fy - r0 as f32);
        let at = |r: usize, c: usize| env[r.min(rows - 1) * cols + c.min(cols - 1)];
        full.d[y * w + x] = (at(r0, c0) * (1.0 - tx) + at(r0, c0 + 1) * tx) * (1.0 - ty) + (at(r0 + 1, c0) * (1.0 - tx) + at(r0 + 1, c0 + 1) * tx) * ty;
    } }
    let full = raster::blur(&full, 40, 2);
    let ratio: Vec<f32> = (0..w * h).map(|i| g.d[i] / full.d[i].max(1.0)).collect();
    let mut hist = [0u64; 256]; let mut tot = 0u64; let mut sum = 0f64;
    let bin = |v: f32| ((v / 1.6).clamp(0.0, 0.999) * 256.0) as usize;
    // above 1.15 is glare off domes and copper edges, not a surface: left out so its long tail cannot pull the split
    for i in 0..w * h { if plain.d[i] && ratio[i] <= 1.15 { let b = bin(ratio[i]); hist[b] += 1; tot += 1; sum += b as f64; } }
    let (mut wb, mut sb, mut best, mut bt) = (0f64, 0f64, 0f64, 128usize);
    for t in 0..256 { wb += hist[t] as f64; if wb == 0.0 { continue; } let wf = tot as f64 - wb; if wf == 0.0 { break; } sb += t as f64 * hist[t] as f64; let (mb, mf) = (sb / wb, (sum - sb) / wf); let bv = wb * wf * (mb - mf).powi(2); if bv > best { best = bv; bt = t; } }
    let thr = (bt as f32 + 1.0) / 256.0 * 1.6;
    // Otsu needs dark pixels to be laminate. On a bright sharpened photo a long tail of shadow and sharpening bands below the laminate pulls the split down under it, and then everything reads as copper (CM5 bottom:
    // photo 540 mm2 where the board draws 228). There the plain mask has two clear modes, and the split is half way between them. The Driver photos do not take this path.
    let tone = crate::prep::mask_tone(img);
    // (Read before the directional closing, which fills the narrow laminate between tracks and blurs the two modes together.)
    let modes = if tone.bright() { let lr: Vec<f32> = (0..w * h).step_by(3).filter(|i| plain.d[*i]).map(|i| g_open.d[i].max(1.0).ln()).collect(); crate::trace::two_modes(&lr) } else { None };
    let thr = modes.map(|m| ((m.0 + m.1) / 2.0).exp() / med.max(1.0)).unwrap_or(thr);
    // and with two modes over the whole face the light is flat, so the face's level is the reference: a window that is mostly laminate has its envelope on the laminate, and then the laminate reads as copper
    let ratio: Vec<f32> = if modes.is_some() { g.d.iter().map(|v| v / med.max(1.0)).collect() } else { ratio };
    if std::env::var("TRACE_DEBUG").is_ok() { eprintln!("copper_read: otsu {:.3}, split {thr:.3}, modes {:?}, env median {med:.0}, plain share {:.3}", (bt as f32 + 1.0) / 256.0 * 1.6, modes.map(|m| (m.0.exp(), m.1.exp())), plain.count() as f64 / (w * h) as f64); }
    Mask { w, h, d: ratio.iter().map(|r| *r > thr).collect() }.close(2.0).open(2.0)
}

fn overlay(img: &image::RgbImage, render: &Mask, truth: &Mask, scored: &Mask, path: &Path) -> Result<(), String> {
    let mut o = img.clone();
    for (k, p) in o.pixels_mut().enumerate() {
        if !scored.d[k] { *p = image::Rgb([p[0] / 4 + 60, p[1] / 4, p[2] / 4 + 60]); continue; }
        *p = match (render.d[k], truth.d[k]) { (true, true) => image::Rgb([30, 190, 60]), (true, false) => image::Rgb([255, 40, 40]), (false, true) => image::Rgb([60, 200, 255]), _ => image::Rgb([p[0] / 3, p[1] / 3, p[2] / 3]) };
    }
    o.save(path).map_err(|e| e.to_string())
}

/// One layer: every figure, on each sample, and the shift test.
struct Sample<'a> { name: &'a str, render: Mask, truth: Mask, scored: Mask, truth2: Option<(Mask, Mask)>, edge: Option<bool>, img: &'a image::RgbImage, img2: Option<&'a image::RgbImage>, face: &'a Face, excluded: Value }

fn measure(s: &Sample, board_area: f64, ev: &Path, problems: &mut Vec<String>) -> Result<Value, String> {
    let a = score(&s.render, &s.truth, &s.scored, &|u, v| !tile_b(u, v));
    let b = score(&s.render, &s.truth, &s.scored, &tile_b);
    let (med, p95, n) = boundary_dist(&s.render, &s.truth, &s.scored);
    let mut d = json!({
        "iouTuningTiles": r3(a.iou),
        "heldOutTiles": {"iou": r3(b.iou), "precision": r3(b.precision), "recall": r3(b.recall), "photoMm2": r2(b.truth_mm2), "boardMm2": r2(b.board_mm2), "boundaryMedianMm": med, "boundaryP95Mm": p95, "boundarySamples": n},
    });
    let scored_total = a.scored_mm2 + b.scored_mm2;
    d["scoredMm2"] = json!(r2(scored_total));
    d["excludedMm2"] = json!(r2(board_area - scored_total));
    d["scoredShareOfBoard"] = json!(r3(scored_total / board_area));
    d["excludedBecause"] = s.excluded.clone();
    if let Some((t2, sc2)) = &s.truth2 {
        let q = score(&s.render, t2, sc2, &tile_b);
        let (m2, p2, _) = boundary_dist(&s.render, t2, sc2);
        d["secondPhotoHeldOutTiles"] = json!({"iou": r3(q.iou), "precision": r3(q.precision), "recall": r3(q.recall), "boundaryMedianMm": m2, "boundaryP95Mm": p2});
        // how well the photo agrees with itself: no drawing can be checked more finely than this
        let both = s.scored.and(sc2);
        d["photoAgainstItsSecondShot"] = json!({"iou": r3(score(&s.truth, t2, &both, &tile_b).iou)});
    }
    // deliberately wrong: the same layer moved by 0.25 mm must score worse
    let px = 10i64;
    let worse: Vec<f64> = [(px, 0), (-px, 0), (0, px), (0, -px)].iter().map(|(x, y)| shifted_iou(&s.render, &s.truth, &s.scored, *x, *y)).collect();
    let moved = worse.iter().cloned().fold(0.0, f64::max);
    d["shiftTest"] = json!({"shiftMm": 0.25, "iouAtBestOfFourShifts": r3(moved), "iouAtZero": r3(b.iou), "scoresWorse": moved < b.iou});
    if b.iou >= 0.3 && moved >= b.iou { problems.push(format!("{}: moving the layer 0.25 mm did not lower its IoU ({:.3} against {:.3}): this measure cannot see it", s.name, moved, b.iou)); }
    // the systematic offset: which small shift this layer would rather have
    let scan = |vert: bool, fix: i64| -> (i64, f64) { let (mut bi, mut bv) = (0i64, -1.0); let mut k = -16i64; while k <= 16 { let v = if vert { shifted_iou(&s.render, &s.truth, &s.scored, fix, k) } else { shifted_iou(&s.render, &s.truth, &s.scored, k, fix) }; if v > bv { bv = v; bi = k; } k += 2; } (bi, bv) };
    let (bx, _) = scan(false, 0); let (by, bv) = scan(true, bx);
    d["bestShiftMm"] = json!({"x": r3(bx as f64 / board::RES), "y": r3(by as f64 / board::RES), "iouThere": r3(bv)});
    if let (Some((t2, sc2)), Some(i2)) = (&s.truth2, s.img2) { overlay(i2, &s.render, t2, sc2, &ev.join(format!("fidelity-{}-second-photo.jpg", s.name.replace('.', "_"))))?; }
    if let Some(inside_bright) = s.edge { d["edgeAgainstLuminance"] = edge_offset(&s.render, s.face, &s.scored, inside_bright); }
    overlay(s.img, &s.render, &s.truth, &s.scored, &ev.join(format!("fidelity-{}.jpg", s.name.replace('.', "_"))))?;
    Ok(d)
}

fn big_disk_holes(holes: &[gerber::Hole], fr: &Frame) -> Mask { let mut m = Mask::new(fr.w, fr.h); for h in holes { gerber::fill_disk(&mut m, fr, h.x, h.y, h.d / 2.0); } m }

/// Photographed holes against the drilled ones: per hole.
fn hole_report(fc: &Face, drilled: &[gerber::Hole]) -> Value {
    let m = class_mask(fc, |c| c == photomat::HOLE);
    let (_, cs) = raster::components(&m);
    let (mut d_err, mut c_err, mut miss) = (vec![], vec![], 0);
    let mut used = vec![false; drilled.len()];
    for c in &cs {
        let (x, y) = board::to_mm(c.cx + 0.5, c.cy + 0.5);
        let dia = 2.0 * (c.area as f64 / std::f64::consts::PI).sqrt() / board::RES;
        let best = (0..drilled.len()).min_by(|i, j| geom::dist((drilled[*i].x, drilled[*i].y), (x, y)).partial_cmp(&geom::dist((drilled[*j].x, drilled[*j].y), (x, y))).unwrap());
        match best.filter(|i| geom::dist((drilled[*i].x, drilled[*i].y), (x, y)) < 0.8) {
            Some(i) => { used[i] = true; c_err.push(geom::dist((drilled[i].x, drilled[i].y), (x, y))); d_err.push((drilled[i].d - dia).abs()); }
            None => miss += 1,
        }
    }
    let q = |v: &mut Vec<f64>, p: f64| { if v.is_empty() { return 0.0; } v.sort_by(|a, b| a.partial_cmp(b).unwrap()); r3(v[((v.len() - 1) as f64 * p) as usize]) };
    json!({"holesOnPhoto": cs.len(), "withADrill": c_err.len(), "withNoDrill": miss, "drillsWithNoPhotoHole": used.iter().filter(|u| !**u).count(), "centreErrorMedianMm": q(&mut c_err, 0.5), "centreErrorP95Mm": q(&mut c_err, 0.95), "diameterErrorMedianMm": q(&mut d_err, 0.5), "diameterErrorP95Mm": q(&mut d_err, 0.95)})
}

/// The copper of one layer as read: the runs of each pour and each track's centre line at its measured widths, from copper.json.
fn read_layer(flow: &Path, layer: &str, fr: &Frame) -> Mask {
    let mut m = Mask::new(fr.w, fr.h);
    let Ok(c) = crate::items::read_json(&flow.join("copper.json")) else { return m };
    for z in c["zones"].as_array().into_iter().flatten().filter(|z| z["layer"].as_str() == Some(layer)) { for r in z["rle_px"].as_array().into_iter().flatten() { let Some(a) = r.as_array() else { continue }; let y = a[0].as_i64().unwrap_or(-1); if y < 0 || y >= fr.h as i64 { continue; } for s in a[1..].chunks(2) { if s.len() != 2 { continue; } for x in s[0].as_i64().unwrap_or(0).max(0)..s[1].as_i64().unwrap_or(0).min(fr.w as i64) { m.d[y as usize * fr.w + x as usize] = true; } } } }
    for t in c["traces"].as_array().into_iter().flatten().filter(|t| t["layer"].as_str() == Some(layer)) { let path: Vec<(f64, f64)> = t["path_mm"].as_array().map(|a| a.iter().map(|q| (q[0].as_f64().unwrap_or(0.0), q[1].as_f64().unwrap_or(0.0))).collect()).unwrap_or_default(); for (j, s) in path.windows(2).enumerate() { let wd = crate::kicad_out::seg_width(t, j); gerber::fill_seg(&mut m, fr, s[0], s[1], wd); gerber::fill_disk(&mut m, fr, s[0].0, s[0].1, wd / 2.0); gerber::fill_disk(&mut m, fr, s[1].0, s[1].1, wd / 2.0); } }
    m
}

/// Pads: each blob of metal the photo shows, against what the board exposes in a window around it.
fn pad_report(fc: &Face, exposed: &Mask, scored: &Mask, sel: &dyn Fn(usize, usize) -> bool) -> Value {
    let truth = class_mask(fc, |c| c == photomat::GOLD || c == photomat::SILVER);
    let (_, cs) = raster::components(&truth);
    let mut v = vec![];
    for c in cs.iter().filter(|c| (c.area as f64) > 0.05 * photomat::MM2) {
        let (mx, my) = ((c.x0 + c.x1) / 2, (c.y0 + c.y1) / 2); if !sel(mx, my) { continue; }
        let (x0, x1, y0, y1) = (c.x0.saturating_sub(8), (c.x1 + 8).min(fc.w - 1), c.y0.saturating_sub(8), (c.y1 + 8).min(fc.h - 1));
        let (mut i, mut u) = (0u32, 0u32);
        for y in y0..=y1 { for x in x0..=x1 { let k = y * fc.w + x; if !scored.d[k] { continue; } let (a, b) = (exposed.d[k], truth.d[k]); if a && b { i += 1; } if a || b { u += 1; } } }
        if u > 0 { v.push(i as f64 / u as f64); }
    }
    v.sort_by(|a, b| a.partial_cmp(b).unwrap());
    if v.is_empty() { return json!({"blobs": 0}); }
    json!({"blobs": v.len(), "iouMedian": r3(v[v.len() / 2]), "iouP10": r3(v[v.len() / 10]), "under05": v.iter().filter(|x| **x < 0.5).count()})
}

/// A crude, independent look at the via layer: small bright round domes on the mask of the bottom face. Under frontal light a via has no bright dome but a dark pit, so on a bright mask the same plain test looks for small dark spots.
fn via_report(fc: &Face, drilled: &[gerber::Hole], dark: bool) -> Value {
    let (a, b) = if dark { let (s, l) = (raster::blur(&fc.y, 2, 2), raster::blur(&fc.y, 8, 2)); (Plane { w: s.w, h: s.h, d: s.d.iter().map(|v| -*v).collect() }, Plane { w: l.w, h: l.h, d: l.d.iter().map(|v| -*v).collect() }) } else { (raster::blur(&fc.y, 4, 2), raster::blur(&fc.y, 11, 2)) };
    let (w, h) = (fc.w, fc.h);
    let mut peaks: Vec<(f64, f64, f32)> = vec![];
    for y in (12..h - 12).step_by(1) { for x in (12..w - 12).step_by(1) {
        let k = y * w + x; if fc.cls[k] != photomat::MASK || !tile_b(x, y) { continue; }
        let d = a.d[k] - b.d[k];
        if d < 4.0 { continue; }
        let mut top = true; 'o: for dy in -7i64..=7 { for dx in -7i64..=7 { let j = ((y as i64 + dy) * w as i64 + x as i64 + dx) as usize; if a.d[j] - b.d[j] > d { top = false; break 'o; } } }
        if top && !peaks.iter().rev().take(300).any(|p| geom::dist((p.0, p.1), (x as f64, y as f64)) < 12.0) { peaks.push((x as f64 + 0.5, y as f64 + 0.5, d)); }
    } }
    // the strengths split in two, grain and domes, by Otsu: no threshold chosen by looking at the vias
    let peaks: Vec<(f64, f64, f32)> = { let mut hist = [0u32; 128]; for p in &peaks { hist[(p.2 as usize).min(127)] += 1; } let tot: f64 = peaks.len() as f64; let sum: f64 = hist.iter().enumerate().map(|(i, v)| i as f64 * *v as f64).sum(); let (mut wb, mut sb, mut bv, mut bt) = (0f64, 0f64, 0f64, 0usize); for t in 0..128 { wb += hist[t] as f64; if wb == 0.0 { continue; } let wf = tot - wb; if wf == 0.0 { break; } sb += t as f64 * hist[t] as f64; let v = wb * wf * (sb / wb - (sum - sb) / wf).powi(2); if v > bv { bv = v; bt = t; } } peaks.into_iter().filter(|p| p.2 > bt as f32).collect() };
    let vias: Vec<(f64, f64)> = drilled.iter().filter(|d| d.via).map(|d| board::to_px(d.x, d.y)).filter(|p| { let (u, v) = (p.0 as i64, p.1 as i64); u >= 0 && v >= 0 && (u as usize) < w && (v as usize) < h && tile_b(u as usize, v as usize) && fc.cls[v as usize * w + u as usize] == photomat::MASK }).collect();
    let near = |p: (f64, f64), set: &[(f64, f64)]| set.iter().any(|q| geom::dist(p, *q) < 10.0);
    let pk: Vec<(f64, f64)> = peaks.iter().map(|p| (p.0, p.1)).collect();
    let hit_v = vias.iter().filter(|v| near(**v, &pk)).count();
    let hit_p = pk.iter().filter(|p| near(**p, &vias)).count();
    let (pr, rc) = (hit_v as f64 / vias.len().max(1) as f64, hit_p as f64 / pk.len().max(1) as f64);
    let f1 = if pr + rc > 0.0 { 2.0 * pr * rc / (pr + rc) } else { 0.0 };
    let on_silk = drilled.iter().filter(|d| d.via).filter(|d| { let p = board::to_px(d.x, d.y); let (u, v) = (p.0 as i64, p.1 as i64); u >= 0 && v >= 0 && (u as usize) < w && (v as usize) < h && fc.cls[v as usize * w + u as usize] == photomat::SILK }).count();
    json!({"viasInHeldOutTiles": vias.len(), "domesInHeldOutTiles": pk.len(), "viasOnADome": hit_v, "domesWithAVia": hit_p, "precision": r3(pr), "recall": r3(rc), "f1": r3(f1), "viasOnSilk": on_silk, "note": if dark { "domes here are small dark spots found by a plain blob test on the bottom face, not the flow's via detector: an independent floor, not a truth" } else { "domes found by a plain bright-blob test on the bottom face, not the flow's via detector: an independent floor, not a truth" }})
}

pub fn run(flow: &Path, pcb: &Path, ev: &Path) -> Result<(Value, Vec<String>), String> {
    let open = |n: &str| image::open(flow.join(n)).map(|i| i.to_rgb8()).map_err(|e| format!("{n}: {e}"));
    let (ti, bi) = (open("board/top.jpg")?, open("board/bottom.jpg")?);
    let (mut top, mut bot) = (photomat::classify(&ti, photomat::carries_parts(&ti)), photomat::classify(&bi, photomat::carries_parts(&bi)));
    photomat::link_holes(&mut top, &mut bot);
    let (w, h) = (top.w, top.h);
    let fr = Frame { w, h, res: board::RES, ox: board::ORIGIN, oy: board::ORIGIN };
    let plot_dir = flow.join(".fidelity-plot");
    let plot = gerber::plot(pcb, &plot_dir)?;
    let layer = |k: &str| -> Result<Mask, String> { gerber::rasterize(plot.layers.get(k).ok_or(format!("the plot has no {k}"))?, &fr) };
    let (fcu, bcu, fmask, bmask, fsilk, bsilk) = (layer("F_Cu")?, layer("B_Cu")?, layer("F_Mask")?, layer("B_Mask")?, layer("F_Silkscreen")?, layer("B_Silkscreen")?);
    let pad_holes: Vec<gerber::Hole> = plot.drill.iter().filter(|d| !d.via).cloned().collect();
    let all_holes = big_disk_holes(&plot.drill, &fr);
    let drilled = big_disk_holes(&pad_holes, &fr);
    let _ = std::fs::remove_dir_all(&plot_dir);
    let edge_text = plot.layers.get("Edge_Cuts").cloned().unwrap_or_default();
    // second photos, registered through their own holes
    let prep = crate::items::read_json(&flow.join("prep.json"))?; let scale = crate::items::read_json(&flow.join("scale.json"))?;
    let others = |side: &str| -> Option<(image::RgbImage, Value)> {
        let prim = prep["primary"][side].as_str()?;
        let other = prep["photos"].as_array()?.iter().find(|p| p["side"].as_str() == Some(side) && p["file"].as_str() != Some(prim))?;
        let file = other["file"].as_str()?;
        let (i, _, rep) = photolayers::photo_quad(flow, file, side, &prep, &scale, Some(if side == "top" { &ti } else { &bi }))?;
        Some((i, json!({"file": file, "registration": rep})))
    };
    let (t2, b2) = (others("top"), others("bottom"));
    let mut second_doc = Map::new();
    let (ft2, fb2) = (t2.as_ref().map(|x| photomat::classify(&x.0, photomat::carries_parts(&x.0))), b2.as_ref().map(|x| photomat::classify(&x.0, photomat::carries_parts(&x.0))));
    if let Some(x) = &t2 { second_doc.insert("top".into(), x.1.clone()); }
    if let Some(x) = &b2 { second_doc.insert("bottom".into(), x.1.clone()); }
    let round = {
        let mut m = Mask { w, h, d: vec![true; w * h] };
        let r = 2.5 * board::RES;
        for y in 0..h { for x in 0..w { let dx = if (x as f64) < r { r - x as f64 } else if x as f64 > w as f64 - r { x as f64 - (w as f64 - r) } else { 0.0 }; let dy = if (y as f64) < r { r - y as f64 } else if y as f64 > h as f64 - r { y as f64 - (h as f64 - r) } else { 0.0 }; if dx > 0.0 && dy > 0.0 && dx * dx + dy * dy > r * r { m.d[y * w + x] = false; } } }
        m
    };
    let board_area = round.count() as f64 / photomat::MM2;
    let mut problems = vec![];
    let mut layers_doc = Map::new();
    let mut held: Vec<(String, f64)> = vec![];
    let faces: [(&str, &Face, &image::RgbImage, &Mask, &Mask, &Mask, Option<&Face>, Option<&image::RgbImage>); 2] = [
        ("F", &top, &ti, &fcu, &fmask, &fsilk, ft2.as_ref(), t2.as_ref().map(|x| &x.0)),
        ("B", &bot, &bi, &bcu, &bmask, &bsilk, fb2.as_ref(), b2.as_ref().map(|x| &x.0)),
    ];
    let area_of = |fc: &Face, f: &dyn Fn(u8) -> bool| -> f64 { r2((0..w * h).filter(|i| round.d[*i] && f(fc.cls[*i])).count() as f64 / photomat::MM2) };
    for (fname, fc, img, cu, mk, sk, f2, _i2) in faces {
        let exposed = cu.and(mk).andnot(&all_holes);
        let objects = area_of(fc, &|c| c == photomat::OBJECT);
        let (metal_a, hole_a, silk_a) = (area_of(fc, &|c| c == photomat::GOLD || c == photomat::SILVER), area_of(fc, &|c| c == photomat::HOLE), area_of(fc, &|c| c == photomat::SILK));
        let is = |c: u8, set: &[u8]| set.contains(&c);
        let scored_of = |g: &Face, set: &'static [u8]| -> Mask { Mask { w, h, d: (0..w * h).map(|i| round.d[i] && is(g.cls[i], set)).collect() } };
        const SILK_S: &[u8] = &[photomat::MASK, photomat::SILK, photomat::GOLD, photomat::SILVER];
        const OPEN_S: &[u8] = &[photomat::MASK, photomat::SILK, photomat::GOLD, photomat::SILVER, photomat::HOLE];
        const CU_S: &[u8] = &[photomat::MASK];
        const HOLE_S: &[u8] = &[photomat::MASK, photomat::SILK, photomat::GOLD, photomat::HOLE];
        let silk_t = class_mask(fc, |c| c == photomat::SILK);
        let open_t = class_mask(fc, |c| c == photomat::GOLD || c == photomat::SILVER || c == photomat::HOLE);
        let metal_t = class_mask(fc, |c| c == photomat::GOLD || c == photomat::SILVER);
        let hole_t = class_mask(fc, |c| c == photomat::HOLE);
        // the two silks, openings, exposed metal and holes
        let sets: Vec<(String, Mask, Mask, Mask, Option<bool>, Value)> = vec![
            (format!("{fname}.SilkS"), sk.clone(), silk_t, scored_of(fc, SILK_S), Some(true), json!({"partBodiesMm2": objects, "holesMm2": hole_a})),
            (format!("{fname}.Mask opening"), mk.clone(), open_t, scored_of(fc, OPEN_S), None, json!({"partBodiesMm2": objects})),
            (format!("{fname} exposed metal"), exposed.clone(), metal_t, scored_of(fc, OPEN_S), None, json!({"partBodiesMm2": objects})),
            (format!("{fname} holes"), drilled.clone(), hole_t, scored_of(fc, HOLE_S), None, json!({"partBodiesMm2": objects, "solderedMetalMm2": area_of(fc, &|c| c == photomat::SILVER), "note": "a hole with a pin and solder in it does not show"})),
        ];
        for (name, render, truth, scored, edge, excl) in sets {
            let truth2 = f2.map(|g| {
                let t = if name.ends_with("SilkS") { class_mask(g, |c| c == photomat::SILK) } else if name.ends_with("opening") { class_mask(g, |c| c == photomat::GOLD || c == photomat::SILVER || c == photomat::HOLE) } else if name.ends_with("metal") { class_mask(g, |c| c == photomat::GOLD || c == photomat::SILVER) } else { class_mask(g, |c| c == photomat::HOLE) };
                (t, scored_of(g, if name.ends_with("SilkS") { SILK_S } else if name.ends_with("holes") { HOLE_S } else { OPEN_S }))
            });
            let s = Sample { name: &name, render, truth, scored, truth2, edge, img, img2: _i2, face: fc, excluded: excl };
            let d = measure(&s, board_area, ev, &mut problems)?;
            let so = d["secondPhotoHeldOutTiles"]["iou"].as_f64().unwrap_or(1.0).min(d["heldOutTiles"]["iou"].as_f64().unwrap_or(0.0));
            held.push((name.clone(), so));
            layers_doc.insert(name, d);
        }
        // outer copper as seen under the mask: only where the photo shows plain mask, and clear of glare at its edges
        let cu_name = format!("{fname}.Cu");
        let plain = |g: &Face| -> Mask { let m = scored_of(g, CU_S); m.erode(BESIDE_PX) };
        let truth = copper_read(img, fc).and(&plain(fc));
        let truth2 = f2.and_then(|g| _i2.map(|i| (copper_read(i, g).and(&plain(g)), plain(g))));
        if let Ok(dir) = std::env::var("FIDELITY_DUMP") { let save = |m: &Mask, n: &str| { let _ = image::GrayImage::from_fn(w as u32, h as u32, |x, y| image::Luma([if m.at(x as usize, y as usize) { 255 } else { 0 }])).save(format!("{dir}/{fname}-{n}.png")); }; save(cu, "plot"); save(&truth, "truth"); save(&plain(fc), "plain"); if let Some((t, _)) = &truth2 { save(t, "truth2"); } }
        let s = Sample { name: &cu_name, render: cu.clone(), truth, scored: plain(fc), truth2, edge: None, img, img2: _i2, face: fc, excluded: json!({"partBodiesMm2": objects, "silkMm2": silk_a, "metalMm2": metal_a, "holesMm2": hole_a, "note": "copper can be read only under bare mask; a margin of 6 px around anything else goes too"}) };
        let d = measure(&s, board_area, ev, &mut problems)?;
        let so = d["secondPhotoHeldOutTiles"]["iou"].as_f64().unwrap_or(1.0).min(d["heldOutTiles"]["iou"].as_f64().unwrap_or(0.0));
        held.push((cu_name.clone(), so));
        layers_doc.insert(cu_name.clone(), d);
        // what KiCad plots against the copper as read, item by item; and how much of the face the photo can show copper on
        let read = read_layer(flow, &format!("{fname}.Cu"), &fr);
        let visible = plain(fc);
        let iou_over = |over: Option<&Mask>| -> f64 { let (mut i, mut u) = (0u64, 0u64); for k in 0..w * h { if over.map(|o| !o.d[k]).unwrap_or(!round.d[k]) { continue; } if cu.d[k] && read.d[k] { i += 1; } if cu.d[k] || read.d[k] { u += 1; } } if u > 0 { i as f64 / u as f64 } else { 1.0 } };
        let only = |a: &Mask, b: &Mask, o: &Mask| -> f64 { r2((0..w * h).filter(|k| o.d[*k] && a.d[*k] && !b.d[*k]).count() as f64 / photomat::MM2) };
        let margin = scored_of(fc, CU_S).andnot(&visible).count() as f64 / photomat::MM2;
        let vis_mm2 = visible.count() as f64 / photomat::MM2;
        layers_doc.get_mut(&cu_name).unwrap()["plotAgainstRead"] = json!({"iouWholeBoard": r3(iou_over(None)), "iouVisibleOnly": r3(iou_over(Some(&visible))), "plotOnlyVisibleMm2": only(cu, &read, &visible), "readOnlyVisibleMm2": only(&read, cu, &visible), "meaning": "KiCad's plotted copper against the pours and tracks of copper.json drawn directly. The whole-board figure also counts pads and vias, which the read layer does not draw"});
        layers_doc.get_mut(&cu_name).unwrap()["visibility"] = json!({"visibleMm2": r2(vis_mm2), "visibleShareOfBoard": r3(vis_mm2 / board_area), "hiddenBy": {"partBodiesMm2": objects, "silkMm2": silk_a, "metalMm2": metal_a, "holesMm2": hole_a, "marginBesideThemMm2": r2(margin)}, "agreementVisibleOnly": layers_doc[&cu_name]["heldOutTiles"]["iou"]});
        let hr = hole_report(fc, &pad_holes);
        layers_doc.get_mut(&format!("{fname} holes")).unwrap()["perHole"] = hr;
        let sc = class_mask(fc, |c| c != photomat::OBJECT);
        layers_doc.get_mut(&format!("{fname} exposed metal")).unwrap()["perBlob"] = pad_report(fc, &exposed, &sc, &tile_b);
    }
    let outline = photolayers::outline_check(flow, &edge_text, &prep, &scale).unwrap_or_else(|e| json!({"error": e}));
    // the drawn outline against the stated board: a reader that loses an edge must not draw a board of another size
    let off_stated = pho_outline_off(flow, &scale);
    let lowest = |names: &[&str]| -> f64 { held.iter().filter(|(n, _)| names.contains(&n.as_str())).map(|(_, v)| *v).fold(1.0, f64::min) };
    let vias = via_report(&bot, &plot.drill, crate::prep::mask_tone(&bi).bright());
    // the blob test looks for lit domes; on a face lit evenly there are none to find (the CM5 close-up set), only pits
    let domes_lit = crate::items::read_json(&flow.join("detect/domes.json")).ok().and_then(|d| d["light"]["bottom"]["oneLight"].as_bool()).unwrap_or(true);
    let pho = photolayers::load(flow);
    let lib_share = pho.as_ref().and_then(|d| d["libraryPadsAgainstPhoto"]["share"].as_f64());
    let cu_score = 100.0 * lowest(&["F.Cu", "B.Cu"]);
    let metal = lowest(&["F exposed metal", "F holes"]);
    let place = 100.0 * lib_share.map(|s| s.min(metal)).unwrap_or(metal);
    let via_score = 100.0 * vias["f1"].as_f64().unwrap_or(0.0);
    let all = 100.0 * held.iter().map(|x| x.1).fold(1.0, f64::min);
    let pass = held.iter().all(|x| x.1 >= TARGET);
    let below: Vec<String> = held.iter().filter(|x| x.1 < TARGET).map(|x| format!("{} {:.3}", x.0, x.1)).collect();
    let inputs: Map<String, Value> = ["copper.json", "parts.json", photolayers::FILE].iter().filter(|f| flow.join(f).is_file()).map(|f| (f.to_string(), json!(if *f == "copper.json" { crate::attrib::geometry_hash(&flow.join(f)) } else { crate::sha256(&flow.join(f)) }))).collect();
    let doc = json!({
        "generatedBy": "check fidelity", "at": crate::ledger::now(), "board": pcb.display().to_string(), "pxPerMm": board::RES, "outlineOffStatedMm": off_stated,
        "inputs": inputs,
        "categories": {
            "copper": {"score": r2(cu_score), "formula": "100 x the lower IoU of F.Cu and B.Cu, each the lower of its held-out tiles and the second photo. Only bare mask counts: pixels under parts, silk or metal are in neither the numerator nor the denominator, so they cannot raise it. See layers.*.scoredShareOfBoard for how much of the board this covers."},
            "vias": if domes_lit { json!({"score": r2(via_score), "formula": "100 x F1 of written vias against bright domes on the bottom mask in the held-out tiles, from a plain blob test"}) } else { json!({"score": Value::Null, "notMeasured": "the bottom photo shows no lit domes (one light across the face is not seen), so a blob test for domes measures nothing here; the category stays capped until another check covers it", "blobTestF1": r2(via_score)}) },
            "placement": {"score": r2(place), "formula": "100 x the lowest of: F exposed metal IoU, F holes IoU (held-out, second photo), and the share of parts whose library pads sit within 0.25 mm of the pads the photo shows"}
        },
        "layersWorst": r2(all), "targetIoU": TARGET, "belowTarget": below,
        "boardAreaMm2": r2(board_area),
        "samples": {"tuning": "tile set A of a 5 mm checkerboard, primary photos: thresholds come from here", "heldOut": "tile set B, primary photos", "secondPhoto": "tile set B, the other photo of each face, registered through its holes; same camera position and light, so it tests noise and registration, not a choice of threshold", "edgeAgainstLuminance": "threshold-free: the nearest steepest edge in the luminance, for bias in where the silk edge was put"},
        "secondPhotos": second_doc,
        "layers": layers_doc,
        "viaLayer": vias,
        "outline": outline,
        "scaleCheck": {"pxPerMm": scale["pxPerMm"], "boardMm": scale["board_mm"], "rulers": pho.as_ref().map(|d| d["scaleFromPads"].clone()), "outlineMm": pho.as_ref().map(|d| d["outline"]["points_mm"].clone()), "edges": pho.as_ref().map(|d| d["outline"]["edges"].clone()), "bottomRegistration": scale["alignment"]["registration"]},
    });
    let _ = pass;
    Ok((doc, problems))
}

/// A layer under this is reported as not matching; 0.8 is where the crops start to look like the same board drawn thin.
pub const TARGET: f64 = 0.8;

pub fn check(flow: &Path, pcb: &Path) -> Result<crate::checks::Pass, Vec<String>> {
    let ev = flow.join("evidence");
    let (doc, problems) = run(flow, pcb, &ev).map_err(|e| vec![e])?;
    std::fs::write(flow.join("fidelity.json"), serde_json::to_string_pretty(&doc).unwrap() + "\n").map_err(|e| vec![e.to_string()])?;
    let mut hints: Vec<String> = doc["layers"].as_object().into_iter().flatten().map(|(k, v)| format!("{k}: IoU {} held out, {} on the second photo; boundary {} / {} mm; scored {} mm2, excluded {} mm2", v["heldOutTiles"]["iou"], v["secondPhotoHeldOutTiles"]["iou"], v["heldOutTiles"]["boundaryMedianMm"], v["heldOutTiles"]["boundaryP95Mm"], v["scoredMm2"], v["excludedMm2"])).collect();
    hints.push(format!("categories: copper {}, vias {}, placement {}", doc["categories"]["copper"]["score"], doc["categories"]["vias"]["score"], doc["categories"]["placement"]["score"]));
    hints.push("evidence/fidelity-<layer>.jpg: green both, red board only, cyan photo only, purple excluded".into());
    let mut errs = problems;
    for b in doc["belowTarget"].as_array().into_iter().flatten() { errs.push(format!("below IoU {TARGET}: {}", b.as_str().unwrap_or(""))); }
    if let Some(o) = doc["outlineOffStatedMm"].as_f64() { if o > OUTLINE_TOL_MM { errs.push(format!("the drawn outline is {o:.2} mm from the stated {} x {} mm board at its worst corner: look at evidence/outline-*.jpg and photo-layers.json outline.edges", doc["scaleCheck"]["boardMm"][0], doc["scaleCheck"]["boardMm"][1])); } }
    if !errs.is_empty() { errs.extend(hints); return Err(errs); }
    Ok(crate::checks::Pass { summary: json!({"worstLayer": doc["layersWorst"], "target": TARGET}), hints, files: vec!["fidelity.json".into()] })
}

#[cfg(test)]
mod tests {
    use super::*;
    fn blob(x0: usize, y0: usize, s: usize) -> Mask { let mut m = Mask::new(800, 800); for y in y0..y0 + s { for x in x0..x0 + s { m.set(x, y, true); } } m }
    #[test]
    fn a_shifted_layer_scores_worse_and_excluded_area_counts_for_nothing() {
        let truth = blob(100, 100, 600);
        let scored = Mask { w: 800, h: 800, d: vec![true; 640000] };
        let same = score(&truth, &truth, &scored, &tile_b);
        assert!((same.iou - 1.0).abs() < 1e-9);
        assert!(shifted_iou(&truth, &truth, &scored, 10, 0) < 0.99);
        // the same wrong board, with its error excluded or not: excluding the error raises nothing but the error's own pixels
        let wrong = blob(100, 100, 300);
        let a = score(&wrong, &truth, &scored, &|_, _| true);
        assert!((a.iou - 0.25).abs() < 0.01, "{}", a.iou);
        let mut none = Mask::new(800, 800); for i in 0..640000 { none.d[i] = false; }
        let b = score(&wrong, &truth, &none, &|_, _| true);
        assert_eq!(b.scored_mm2, 0.0);
        assert_eq!(b.truth_mm2, 0.0);
    }
    #[test]
    fn distance_is_zero_on_itself() { let m = blob(300, 300, 200); let s = Mask { w: 800, h: 800, d: vec![true; 640000] }; let (a, b, n) = boundary_dist(&m, &m, &s); assert_eq!((a, b), (0.0, 0.0)); assert!(n > 0); }
}

/// Corners of the board drawn from the photos may sit this far from the stated board before `check fidelity` fails.
const OUTLINE_TOL_MM: f64 = 0.5;

/// The worst distance from a drawn outline corner to the stated board's corner, in mm.
fn pho_outline_off(flow: &std::path::Path, scale: &Value) -> Option<f64> {
    let d = photolayers::load(flow)?; let pts = d["outline"]["points_mm"].as_array()?;
    let (w, h) = (scale["board_mm"][0].as_f64()?, scale["board_mm"][1].as_f64()?);
    let o = board::ORIGIN; let want = [(o, o), (o + w, o), (o + w, o + h), (o, o + h)];
    let got: Vec<(f64, f64)> = pts.iter().filter_map(|p| Some((p[0].as_f64()?, p[1].as_f64()?))).collect();
    if got.len() != 4 { return None; }
    Some((0..4).map(|i| ((got[i].0 - want[i].0).powi(2) + (got[i].1 - want[i].1).powi(2)).sqrt()).fold(0.0, f64::max))
}