AI Flow: Reverse Engineering

//! The board frame. `prep` measures every photo (outline, straightened image, sharpness, holes, rows of
//! holes, part bodies) and decides how the bottom maps onto the top. `scale` solves pixels per mm from
//! the hole-row pitches (which standard pitches explain every row at one scale) and a part body whose
//! size a datasheet gives, then writes both sides into one KiCad board frame: x right, y down, mm,
//! origin offset (100, 100), `RES` pixels per mm, bottom mirrored into the top view.
use crate::geom::{self, P};
use crate::prep;
use image::RgbImage;
use serde_json::{json, Value};
use std::path::Path;

pub const RES: f64 = 40.0;
pub const ORIGIN: f64 = 100.0;
pub const PITCHES: [f64; 11] = [1.0, 1.25, 1.27, 1.5, 2.0, 2.5, 2.54, 3.5, 3.81, 5.0, 5.08];

fn f(v: &Value, k: &str) -> f64 { v[k].as_f64().unwrap_or(0.0) }

/// Board-frame mm for a board-image pixel, and back.
pub fn to_mm(u: f64, v: f64) -> P { (ORIGIN + u / RES, ORIGIN + v / RES) }
pub fn to_px(x: f64, y: f64) -> P { ((x - ORIGIN) * RES, (y - ORIGIN) * RES) }

/// Measure one photo. Returns its record and its straightened image (photo-native scale, not mirrored).
pub fn prep_photo(path: &Path, side: &str, ev_dir: &Path) -> Result<Value, String> {
    let img = image::open(path).map_err(|e| format!("{}: {e}", path.display()))?.to_rgb8();
    prep_image(&img, path, side, ev_dir, None)
}

/// The same on a photo already read, straightened through `given` corners when both faces share one outline (`facereg::unify`): (corners, where they came from, the corners fitted to this photo alone).
pub fn prep_image(img: &RgbImage, path: &Path, side: &str, ev_dir: &Path, given: Option<([P; 4], String, Value, Value)>) -> Result<Value, String> {
    let img = img.clone();
    // a mosaic from `stitch` carries its board outline: used when the mask is not green and the colour test finds no board
    let fitted = given.as_ref().map(|g| g.2.clone());
    let own_c: Option<[P; 4]> = fitted.as_ref().and_then(|v| v.as_array()).filter(|a| a.len() == 4).map(|a| [0, 1, 2, 3].map(|i| (a[i][0].as_f64().unwrap_or(0.0), a[i][1].as_f64().unwrap_or(0.0))));
    let beyond_given = given.as_ref().map(|g| g.3.clone());
    let (c, outline_from): ([P; 4], String) = match given { Some(g) => (g.0, g.1), None => match crate::outline::corners(&img).or_else(|_| prep::board_corners(&img)) { Ok(c) => (c, "the photo's green mask".to_string()), Err(e) => (crate::intake::board_box(path, side).ok_or(e)?, "stitch.json (the board box of the mosaic: the mask is not green)".to_string()) } };
    // a tight outline that leaves solder mask beyond an edge has an edge along parts, not along the board. When the wide green test closes every edge, its outline
    // is the better start for a face that may end up resting on its own outline (test board 6: a header hid the top edge and the tight outline lost 13 mm of board)
    let (c, outline_from) = if beyond_given.is_none() && crate::outline::mask_beyond(&img, &c).iter().any(|v| *v > crate::outline::OPEN_EDGE_SHARE) {
        match prep::board_corners(&img) { Ok(wide) if crate::outline::mask_beyond(&img, &wide).iter().all(|v| *v <= crate::outline::OPEN_EDGE_SHARE) => (wide, "the photo's green mask, wide test: the tight outline left board beyond an edge".to_string()), _ => (c, outline_from) }
    } else { (c, outline_from) };
    // sides where the photo stops before the board does, by this photo's own outline or the one it shares with the other face
    let (pw, ph) = (img.width() as f64, img.height() as f64);
    let mut cut = crate::outline::cut_sides(&c, pw, ph);
    if let Some(o) = &own_c { for s in crate::outline::cut_sides(o, pw, ph) { if !cut.contains(&s) { cut.push(s); } } }
    let w = ((geom::dist(c[0], c[1]) + geom::dist(c[3], c[2])) / 2.0).round() as u32;
    let h = ((geom::dist(c[0], c[3]) + geom::dist(c[1], c[2])) / 2.0).round() as u32;
    let (r, _) = prep::rectify(&img, &c, w, h, false)?;
    let holes = prep::holes(&r, 20.0, 200.0);
    let gold_holes = prep::gold_holes(&r, 20.0, 200.0);
    let rows = prep::hole_rows(&holes);
    // what goes through the board, in the photo's own pixels, for laying the faces on each other
    let through = geom::homography(&[(0.0, 0.0), (w as f64, 0.0), (w as f64, h as f64), (0.0, h as f64)], &c).map(|hm| crate::facereg::through(&crate::facereg::pads(&img, 20.0, 520.0), &holes, &hm)).unwrap_or_default();
    let bodies = if side == "top" { prep::bodies(&r, 4.0, 18.0) } else { vec![] };
    let stem = path.file_stem().unwrap().to_string_lossy().to_string();
    // evidence: the outline on the photo, and the straightened photo with holes, rows and bodies numbered
    let mut ov = img.clone(); crate::draw::poly(&mut ov, &c, [255, 0, 255], 7);
    image::imageops::resize(&ov, img.width() / 3, img.height() / 3, image::imageops::FilterType::Triangle).save(ev_dir.join(format!("outline-{stem}.jpg"))).map_err(|e| e.to_string())?;
    let mut rv = r.clone();
    for (k, b) in bodies.iter().enumerate() { draw_rect(&mut rv, b, [0, 255, 255], 3); crate::draw::label(&mut rv, (f(b, "x") - 10.0, f(b, "y") - 8.0), &format!("{k}"), [255, 90, 0], 4); }
    for (k, hh) in holes.iter().enumerate() { crate::draw::circle(&mut rv, (f(hh, "x"), f(hh, "y")), f(hh, "hole") / 2.0 + 6.0, [255, 0, 255], 4); crate::draw::label(&mut rv, (f(hh, "x") + 22.0, f(hh, "y") - 10.0), &format!("{k}"), [255, 255, 0], 5); }
    for rw in &rows { let ids: Vec<usize> = rw["holes"].as_array().unwrap().iter().map(|x| x.as_u64().unwrap() as usize).collect();
        crate::draw::line(&mut rv, (f(&holes[ids[0]], "x"), f(&holes[ids[0]], "y")), (f(&holes[ids[ids.len() - 1]], "x"), f(&holes[ids[ids.len() - 1]], "y")), [255, 255, 0], 3); }
    rv.save(ev_dir.join(format!("prep-{stem}.jpg"))).map_err(|e| e.to_string())?;
    Ok(json!({"file": path.file_name().unwrap().to_string_lossy(), "side": side, "outlineFrom": outline_from, "corners": c.iter().map(|p| json!([p.0, p.1])).collect::<Vec<_>>(),
              "rectPx": [w, h], "sharpness": prep::sharpness(&r), "holes": holes, "goldHoles": gold_holes, "through": through, "photoPx": [img.width(), img.height()], "cornersOwn": fitted, "notBoardEdge": cut, "maskBeyondOwnEdges": beyond_given.unwrap_or_else(|| json!(crate::outline::mask_beyond(&img, &c))), "rows": rows, "bodies": bodies,
              "evidence": [format!("evidence/outline-{stem}.jpg"), format!("evidence/prep-{stem}.jpg")]}))
}

pub fn draw_rect(img: &mut RgbImage, b: &Value, c: [u8; 3], th: i64) {
    let (cx, cy, l, s, a) = (f(b, "x"), f(b, "y"), f(b, "long"), f(b, "short"), f(b, "angle").to_radians());
    let (co, sn) = (a.cos(), a.sin());
    let pts: Vec<P> = [(-1.0, -1.0), (1.0, -1.0), (1.0, 1.0), (-1.0, 1.0)].iter().map(|(u, v)| (cx + u * l / 2.0 * co - v * s / 2.0 * sn, cy + u * l / 2.0 * sn + v * s / 2.0 * co)).collect();
    crate::draw::poly(img, &pts, c, th);
}

/// Scale candidates: for each row and standard pitch, the px/mm it implies, scored by how many rows
/// (weighted by hole count) land on some standard pitch at that scale. Close candidates merge.
pub fn scale_candidates(rows: &[Value]) -> Vec<Value> {
    let mut cands: Vec<(f64, f64, Vec<Value>)> = vec![];
    for r in rows {
        for s in PITCHES {
            let k = f(r, "pitch") / s;
            if cands.iter().any(|c| (c.0 - k).abs() / k < 0.01) { continue; }
            let mut score = 0.0; let mut explain = vec![];
            for q in rows {
                let mm = f(q, "pitch") / k;
                if let Some(std) = PITCHES.iter().find(|p| ((mm - **p) / **p).abs() < 0.012) {
                    score += f(q, "count"); explain.push(json!({"holes": q["holes"], "pitchPx": f(q, "pitch"), "mm": std}));
                }
            }
            cands.push((k, score, explain));
        }
    }
    cands.sort_by(|a, b| b.1.partial_cmp(&a.1).unwrap().then(a.0.partial_cmp(&b.0).unwrap()));
    // every candidate is kept: a board whose rows all share one pitch ties every standard, and cutting the list
    // at six dropped the true one (2.0 mm) on the first such board (USB hub, 2026-11-10)
    cands.into_iter().map(|(k, s, e)| json!({"pxPerMm": (k * 1000.0).round() / 1000.0, "score": s, "rows": e})).collect()
}

/// Corners of a photo mapped into the board frame (pixels at RES), mirrored for the bottom.
pub fn frame_corners(wmm: f64, hmm: f64, side: &str, axis: &str) -> [P; 4] {
    let (w, h) = (wmm * RES, hmm * RES);
    match (side, axis) {
        ("bottom", "y") => [(0.0, h), (w, h), (w, 0.0), (0.0, 0.0)],
        ("bottom", _) => [(w, 0.0), (0.0, 0.0), (0.0, h), (w, h)],
        _ => [(0.0, 0.0), (w, 0.0), (w, h), (0.0, h)],
    }
}

/// Warp a photo into the board frame through its own corners.
pub fn to_frame(path: &Path, corners: &[P; 4], dst: &[P; 4], wmm: f64, hmm: f64, shift: P) -> Result<RgbImage, String> {
    let img = image::open(path).map_err(|e| e.to_string())?.to_rgb8();
    let d: [P; 4] = dst.map(|p| (p.0 + shift.0, p.1 + shift.1));
    let hm = geom::homography(corners, &d).ok_or("degenerate corners")?;
    let inv = geom::invert(&hm).ok_or("not invertible")?;
    let (w, h) = ((wmm * RES).round() as u32, (hmm * RES).round() as u32);
    let planes: Vec<crate::raster::Plane> = (0..3).map(|c| crate::raster::Plane { w: img.width() as usize, h: img.height() as usize, d: img.pixels().map(|p| p[c] as f32).collect() }).collect();
    Ok(RgbImage::from_fn(w, h, |u, v| {
        let q = geom::apply(&inv, (u as f64 + 0.5, v as f64 + 0.5));
        image::Rgb([0, 1, 2].map(|c| planes[c].sample(q.0 as f32 - 0.5, q.1 as f32 - 0.5).clamp(0.0, 255.0) as u8))
    }))
}

/// Translation that best lays the bottom holes onto the top holes (board-frame pixels), with the
/// number of pairs and the residual after the shift.
pub fn align(top: &[Value], bot: &[Value]) -> (P, usize, f64) {
    let pairs = |shift: P| -> Vec<(P, P)> {
        bot.iter().filter_map(|b| {
            let q = (f(b, "x") + shift.0, f(b, "y") + shift.1);
            top.iter().map(|t| (f(t, "x"), f(t, "y"))).min_by(|a, c| geom::dist(*a, q).partial_cmp(&geom::dist(*c, q)).unwrap())
                .filter(|t| geom::dist(*t, q) < 1.5 * RES).map(|t| (t, (f(b, "x"), f(b, "y"))))
        }).collect()
    };
    let mut shift = (0.0, 0.0);
    for _ in 0..3 {
        let p = pairs(shift);
        if p.is_empty() { break; }
        let n = p.len() as f64;
        shift = (p.iter().map(|(t, b)| t.0 - b.0).sum::<f64>() / n, p.iter().map(|(t, b)| t.1 - b.1).sum::<f64>() / n);
    }
    let p = pairs(shift);
    let res = if p.is_empty() { f64::NAN } else { (p.iter().map(|(t, b)| geom::dist(*t, (b.0 + shift.0, b.1 + shift.1)).powi(2)).sum::<f64>() / p.len() as f64).sqrt() };
    (shift, p.len(), res)
}

/// A photo framed through four mounting holes at a spacing a datasheet gives.
pub struct HoleFrame { /// where the board's corners lie in the photo, by those holes
    pub corners: [P; 4], pub holes: [P; 4], /// px/mm along the photo's x and y
    pub px_mm: (f64, f64), /// how far the corners fitted to the board's edges lie from these, mm, largest
    pub fitted_off_mm: f64, /// how far the opening's centre lies from the ring's, mm, largest: what an angled view moves
    pub opening_off_mm: f64, /// the RMS error of the ring circles, photo pixels, largest
    pub ring_fit_px: f64 }

/// The frame of one photo from its four mounting holes. A board with no rows of holes has only these for scale, and four points also give the perspective outright, so the frame need not
/// lean on corners fitted to edges. Holes are found by `prep` on the straightened picture; each centre is then read off the ring on the photo itself.
pub fn hole_frame(path: &Path, rec: &Value, spacing: (f64, f64), board: (f64, f64)) -> Result<HoleFrame, String> {
    let img = image::open(path).map_err(|e| format!("{}: {e}", path.display()))?.to_rgb8();
    // the holes were found on the picture straightened through the fitted corners, whatever an earlier `scale` has since written as corners
    let c = rec["cornersFitted"].as_array().or(rec["corners"].as_array()).ok_or("no corners in prep.json")?;
    let fitted: [P; 4] = [0, 1, 2, 3].map(|i| (c[i][0].as_f64().unwrap_or(0.0), c[i][1].as_f64().unwrap_or(0.0)));
    let (rw, rh) = (rec["rectPx"][0].as_f64().unwrap_or(0.0), rec["rectPx"][1].as_f64().unwrap_or(0.0));
    let rect: [P; 4] = [(0.0, 0.0), (rw, 0.0), (rw, rh), (0.0, rh)];
    let to_photo = geom::homography(&rect, &fitted).ok_or("degenerate corners")?;
    let holes = rec["goldHoles"].as_array().ok_or("no gold holes in prep.json: run prep again")?;
    let biggest = holes.iter().map(|h| f(h, "hole")).fold(0.0, f64::max);
    if biggest <= 0.0 { return Err("prep found no holes in this photo".into()); }
    // the hole of mounting size nearest each corner of the board
    let big: Vec<&Value> = holes.iter().filter(|h| f(h, "hole") >= 0.6 * biggest).collect();
    // where the stated spacing puts each hole, in the straightened picture
    let (sx, sy) = (rw / board.0, rh / board.1);
    let (ix, iy) = ((board.0 - spacing.0) / 2.0, (board.1 - spacing.1) / 2.0);
    let expected: [P; 4] = [(ix * sx, iy * sy), ((board.0 - ix) * sx, iy * sy), ((board.0 - ix) * sx, (board.1 - iy) * sy), (ix * sx, (board.1 - iy) * sy)];
    let reach = HOLE_REACH_MM * sx;
    let (mut pts, mut opening, mut worst_fit) = ([(0.0, 0.0); 4], [(0.0, 0.0); 4], 0f64);
    for (k, e) in expected.iter().enumerate() {
        let near = big.iter().filter(|h| geom::dist((f(h, "x"), f(h, "y")), *e) < reach).min_by(|a, b| geom::dist((f(a, "x"), f(a, "y")), *e).partial_cmp(&geom::dist((f(b, "x"), f(b, "y")), *e)).unwrap());
        match near {
            Some(h) => {
                let q = geom::apply(&to_photo, (f(h, "x"), f(h, "y")));
                let (centre, fit) = fit_ring(&img, q, f(h, "hole") / 2.0, f(h, "outer") / 2.0);
                if fit.is_finite() { worst_fit = worst_fit.max(fit); }
                pts[k] = centre; opening[k] = q;
            }
            // a ring joined to something beside it does not read as a ring in prep: the opening is still where the spacing says, and its ring is fitted from there.
            // The opening alone will not do: the camera sees the wall of the hole, which moves its middle toward the photo's centre, about 0.15 mm on the CM5 set.
            None => {
                let c = dark_opening(&img, &to_photo, *e, reach, biggest).ok_or_else(|| format!("no hole of mounting size within {HOLE_REACH_MM} mm of where the stated spacing puts hole {} ({} gold holes of that size found)", k + 1, big.len()))?;
                let q = geom::apply(&to_photo, c);
                let (centre, fit) = fit_ring(&img, q, biggest / 2.0, biggest * 1.15);
                if !fit.is_finite() { return Err(format!("hole {} has an opening but no ring that can be fitted", k + 1)); }
                worst_fit = worst_fit.max(fit);
                pts[k] = centre; opening[k] = q;
            }
        }
    }
    let (dx, dy) = ((geom::dist(pts[0], pts[1]) + geom::dist(pts[3], pts[2])) / 2.0, (geom::dist(pts[0], pts[3]) + geom::dist(pts[1], pts[2])) / 2.0);
    let (kx, ky) = (dx / spacing.0, dy / spacing.1);
    // an angled photo leaves some foreshortening, which the map takes out; beyond that these four are not the stated pattern
    if (kx - ky).abs() / kx.max(ky) > 0.05 { return Err(format!("these holes do not fit the stated spacing: {dx:.1} x {dy:.1} px gives {kx:.2} and {ky:.2} px/mm")); }
    let inset = ((board.0 - spacing.0) / 2.0, (board.1 - spacing.1) / 2.0);
    let at = [(inset.0, inset.1), (inset.0 + spacing.0, inset.1), (inset.0 + spacing.0, inset.1 + spacing.1), (inset.0, inset.1 + spacing.1)];
    let to_mm = geom::homography(&pts, &at).ok_or("degenerate holes")?;
    let back = geom::invert(&to_mm).ok_or("not invertible")?;
    let corners = [(0.0, 0.0), (board.0, 0.0), (board.0, board.1), (0.0, board.1)].map(|p| geom::apply(&back, p));
    let k = (kx + ky) / 2.0;
    let fitted_off_mm = (0..4).map(|i| geom::dist(corners[i], fitted[i]) / k).fold(0.0, f64::max);
    let opening_off_mm = (0..4).map(|i| geom::dist(pts[i], opening[i]) / k).fold(0.0, f64::max);
    Ok(HoleFrame { corners, holes: pts, px_mm: (kx, ky), fitted_off_mm, opening_off_mm, ring_fit_px: worst_fit })
}

/// The ring around an opening, by the ring's own colour first since all four must be measured the same way, then by the metal test when that finds no ring.
fn fit_ring(img: &RgbImage, q: P, hole_r: f64, ring_r: f64) -> (P, f64) {
    let warm = prep::ring_centre_by(img, q, hole_r, ring_r, prep::is_ring_gold);
    if warm.1.is_finite() { warm } else { prep::ring_centre(img, q, hole_r, ring_r) }
}

/// How far from its stated place a mounting hole may be found, in mm: the fitted corners a photo starts from are off by about a millimetre on rounded boards.
const HOLE_REACH_MM: f64 = 3.0;
/// An opening darker than this share of the window's brightness range, over what shows through a hole against a dark table. CM5 photos, one close-up set: openings under 15 percent.
const OPENING_DARK: f64 = 0.15;

/// The round dark opening nearest a stated hole position, searched in the straightened picture's coordinates. `None` when the table shows lighter than the board there, or nothing round and of mounting size is near.
fn dark_opening(img: &RgbImage, to_photo: &[f64; 9], at: P, reach: f64, hole_d: f64) -> Option<P> {
    let r = reach.ceil() as i64; let n = (2 * r + 1) as usize;
    let mut luma = vec![0f64; n * n];
    for v in 0..n { for u in 0..n {
        let q = geom::apply(to_photo, (at.0 + u as f64 - r as f64, at.1 + v as f64 - r as f64));
        let (x, y) = (q.0 as i64, q.1 as i64);
        luma[v * n + u] = if x < 0 || y < 0 || x >= img.width() as i64 || y >= img.height() as i64 { 255.0 } else { let p = img.get_pixel(x as u32, y as u32); 0.299 * p[0] as f64 + 0.587 * p[1] as f64 + 0.114 * p[2] as f64 };
    } }
    let lo = luma.iter().cloned().fold(f64::MAX, f64::min); let hi = luma.iter().cloned().fold(0.0, f64::max);
    let lim = lo + OPENING_DARK * (hi - lo);
    let mut seen = vec![false; n * n]; let mut best: Option<(f64, P)> = None;
    for s in 0..n * n {
        if seen[s] || luma[s] > lim { continue; }
        let (mut st, mut px) = (vec![s], vec![]); seen[s] = true;
        while let Some(i) = st.pop() {
            px.push(i); let (u, v) = (i % n, i / n);
            if u > 0 && !seen[i - 1] && luma[i - 1] <= lim { seen[i - 1] = true; st.push(i - 1); }
            if u + 1 < n && !seen[i + 1] && luma[i + 1] <= lim { seen[i + 1] = true; st.push(i + 1); }
            if v > 0 && !seen[i - n] && luma[i - n] <= lim { seen[i - n] = true; st.push(i - n); }
            if v + 1 < n && !seen[i + n] && luma[i + n] <= lim { seen[i + n] = true; st.push(i + n); }
        }
        let (us, vs): (Vec<f64>, Vec<f64>) = px.iter().map(|i| ((i % n) as f64, (i / n) as f64)).unzip();
        let (w, h) = (us.iter().cloned().fold(0.0, f64::max) - us.iter().cloned().fold(f64::MAX, f64::min) + 1.0, vs.iter().cloned().fold(0.0, f64::max) - vs.iter().cloned().fold(f64::MAX, f64::min) + 1.0);
        let a = px.len() as f64; let d = (w + h) / 2.0;
        // round, full, and about the size of the other mounting holes
        if !(0.8..1.25).contains(&(w / h)) || a < 0.68 * w * h || !(0.7 * hole_d..1.4 * hole_d).contains(&d) { continue; }
        let c = (at.0 - r as f64 + us.iter().sum::<f64>() / a, at.1 - r as f64 + vs.iter().sum::<f64>() / a);
        let off = geom::dist(c, at);
        if best.map_or(true, |b| off < b.0) { best = Some((off, c)); }
    }
    best.map(|b| b.1)
}

/// High-pass luminance and where it can be compared: plain mask, away from anything else. Via rims and the copper's edges survive it, parts and print do not.
fn plain_detail(img: &RgbImage) -> (crate::raster::Plane, crate::raster::Mask) {
    let tone = prep::mask_tone(img);
    let (w, h) = (img.width() as usize, img.height() as usize);
    let l = crate::raster::luma(img);
    let (a, b) = (crate::raster::blur(&l, 1, 3), crate::raster::blur(&l, 7, 3));
    // a mask that is not green is told by its own hue (`maskcolour`); the green rule stands as it was
    let mc = crate::maskcolour::of(img);
    let m = crate::raster::Mask { w, h, d: img.pixels().map(|p| if mc.green { (p[1] as f32 - p[0] as f32) > 0.3 * tone.gr50 && p[1] > p[2] } else { mc.is_mask(p[0] as f32, p[1] as f32, p[2] as f32) }).collect() }.erode(6.0);
    (crate::raster::Plane { w, h, d: a.d.iter().zip(&b.d).map(|(x, y)| x - y).collect() }, m)
}

/// Normalised correlation of `a` with `b` moved by (dx, dy), over pixels both show, on every second pixel. None with too few pixels.
fn ncc_at(a: &(crate::raster::Plane, crate::raster::Mask), b: &(crate::raster::Plane, crate::raster::Mask), area: (usize, usize, usize, usize), dx: i64, dy: i64) -> Option<f64> {
    let (w, h) = (a.0.w as i64, a.0.h as i64);
    let (mut sa, mut sb, mut saa, mut sbb, mut sab, mut n) = (0f64, 0f64, 0f64, 0f64, 0f64, 0f64);
    let mut y = area.1 as i64;
    while y < area.3 as i64 { let mut x = area.0 as i64; while x < area.2 as i64 {
        let (u, v) = (x + dx, y + dy);
        if u >= 0 && v >= 0 && u < w && v < h { let (i, j) = ((y * w + x) as usize, (v * w + u) as usize);
            if a.1.d[i] && b.1.d[j] { let (p, q) = (a.0.d[i] as f64, b.0.d[j] as f64); sa += p; sb += q; saa += p * p; sbb += q * q; sab += p * q; n += 1.0; } }
        x += 2; } y += 2; }
    if n < 300.0 { return None; }
    let d = (saa - sa * sa / n) * (sbb - sb * sb / n);
    if d <= 0.0 { None } else { Some((sab - sa * sb / n) / d.sqrt()) }
}

/// The best correlation over a small search of shifts, whole frame.
fn best_ncc(a: &(crate::raster::Plane, crate::raster::Mask), b: &(crate::raster::Plane, crate::raster::Mask), area: (usize, usize, usize, usize), reach: i64, step: i64) -> Option<(f64, i64, i64)> {
    let mut best: Option<(f64, i64, i64)> = None;
    let mut dy = -reach; while dy <= reach { let mut dx = -reach; while dx <= reach {
        if let Some(c) = ncc_at(a, b, area, dx, dy) { if best.map(|q| c > q.0).unwrap_or(true) { best = Some((c, dx, dy)); } }
        dx += step; } dy += step; }
    best
}

/// Which way the board was turned over, when the holes cannot say because their pattern is symmetric. Through features show on both faces, so the right mirror correlates and the
/// wrong one does not. Returns (axis, correlation of x, correlation of y), axis None when neither stands out: at least 0.1 and twice the other.
pub fn flip_by_detail(top: &RgbImage, bottom_x: &RgbImage, bottom_y: &RgbImage) -> (Option<String>, f64, f64) {
    let t = plain_detail(top); let area = (0, 0, t.0.w, t.0.h);
    let cx = best_ncc(&t, &plain_detail(bottom_x), area, 16, 2).map(|c| c.0).unwrap_or(0.0);
    let cy = best_ncc(&t, &plain_detail(bottom_y), area, 16, 2).map(|c| c.0).unwrap_or(0.0);
    let axis = if cx >= 0.1 && cx >= 2.0 * cy.max(0.0) { Some("x".to_string()) } else if cy >= 0.1 && cy >= 2.0 * cx.max(0.0) { Some("y".to_string()) } else { None };
    (axis, cx, cy)
}

/// Where each 5 mm tile of `top` lies in `bottom`, from the same through features: (tile centre, its place in `bottom`), tiles that correlate at `TILE_NCC` or better. For boards with too few shared holes to fit a warp.
pub fn tile_pairs(top: &RgbImage, bottom: &RgbImage) -> Vec<(P, P)> {
    let (a, b) = (plain_detail(top), plain_detail(bottom));
    let (w, h) = (a.0.w, a.0.h); let t = 200usize;
    let mut out = vec![];
    let mut y = 0; while y + t <= h { let mut x = 0; while x + t <= w {
        if let Some((c, dx, dy)) = best_ncc(&a, &b, (x, y, x + t, y + t), 16, 1) { if c >= TILE_NCC { let m = ((x + t / 2) as f64, (y + t / 2) as f64); out.push((m, (m.0 + dx as f64, m.1 + dy as f64))); } }
        x += t; } y += t; }
    out
}
/// Vias that show on both faces, paired top to bottom: the same positions on both sides of the board, so they hold the faces together more tightly than correlating tiles.
/// Pairs start wide for what four holes leave, then the fit pulls them in. None unless the pits line up well beyond chance. Pairs are (top, bottom) in board px.
pub fn pit_pairs(top: &RgbImage, bottom: &RgbImage) -> Option<Vec<(P, P)>> {
    let kept = |i: &RgbImage| -> Vec<P> { crate::domes::read_pits(i).pits.iter().filter(|p| p.keep).map(|p| (p.x as f64, p.y as f64)).collect() };
    let (a, b) = (kept(top), kept(bottom));
    let (_, _, decisive) = crate::domes::pits_confirmed(&a, &b, (top.width() as f64, top.height() as f64));
    if !decisive { return None; }
    let mut q: Option<crate::regframe::Quad> = None;
    let mut out: Vec<(P, P)> = vec![];
    for reach_mm in [0.3, 0.2, 0.15, 0.1, 0.1] {
        let moved = |t: P| q.as_ref().map_or(t, |f| f.apply(t));
        out = b.iter().filter_map(|o| a.iter().min_by(|x, y| geom::dist(moved(**x), *o).partial_cmp(&geom::dist(moved(**y), *o)).unwrap()).filter(|t| geom::dist(moved(**t), *o) < reach_mm * RES).map(|t| (*t, *o))).collect();
        q = crate::regframe::Quad::fit(&out).or(q);
    }
    (out.len() >= PIT_PAIRS).then_some(out)
}
/// Fewer through vias than this and a quadratic over the whole board is not well held. The CM5 close-up set gives 154.
const PIT_PAIRS: usize = 40;

/// On the CM5 photos, tiles of plain mask correlate between faces at 0.3 to 0.75 where vias go through and under 0.2 where they do not.
const TILE_NCC: f64 = 0.4;

#[cfg(test)]
mod tests {
    use super::*;
    #[test]
    fn four_holes_frame_a_tilted_photo() {
        // a 55 x 40 board at 24 px/mm with the top edge shrunk by perspective, holes as gold rings; the fitted corners are off on purpose
        let truth = [(100.0, 80.0), (1400.0, 90.0), (1420.0, 1060.0), (90.0, 1050.0)];
        let h = geom::homography(&[(0.0, 0.0), (55.0, 0.0), (55.0, 40.0), (0.0, 40.0)], &truth).unwrap();
        let hi = geom::invert(&h).unwrap();
        let img = RgbImage::from_fn(1500, 1150, |x, y| { let q = geom::apply(&hi, (x as f64 + 0.5, y as f64 + 0.5));
            if !(0.0..55.0).contains(&q.0) || !(0.0..40.0).contains(&q.1) { return image::Rgb([90, 58, 40]); }
            let d = [(3.5, 3.5), (51.5, 3.5), (51.5, 36.5), (3.5, 36.5)].iter().map(|c: &(f64, f64)| geom::dist(*c, q)).fold(f64::MAX, f64::min);
            if d < 1.25 { image::Rgb([60, 40, 30]) } else if d < 2.8 { image::Rgb([230, 210, 60]) } else { image::Rgb([50, 186, 140]) } });
        let dir = std::env::temp_dir().join(format!("hole-frame-{}", std::process::id())); std::fs::create_dir_all(&dir).unwrap();
        let path = dir.join("p.png"); img.save(&path).unwrap();
        let off = [(truth[0].0 + 6.0, truth[0].1 + 5.0), (truth[1].0 - 4.0, truth[1].1 + 7.0), (truth[2].0 + 3.0, truth[2].1 - 6.0), (truth[3].0 - 2.0, truth[3].1 - 5.0)];
        let (rw, rh) = (1320u32, 970u32);
        let (r, _) = prep::rectify(&img, &off, rw, rh, false).unwrap();
        let holes = prep::gold_holes(&r, 20.0, 200.0);
        assert_eq!(holes.len(), 4, "{holes:?}");
        let rec = json!({"corners": off.iter().map(|p| json!([p.0, p.1])).collect::<Vec<_>>(), "rectPx": [rw, rh], "goldHoles": holes});
        let hf = hole_frame(&path, &rec, (48.0, 33.0), (55.0, 40.0)).unwrap();
        for i in 0..4 { assert!(geom::dist(hf.corners[i], truth[i]) < 2.0, "{i}: {:?} vs {:?}", hf.corners[i], truth[i]); }
        assert!(hf.fitted_off_mm > 0.1 && hf.fitted_off_mm < 1.0, "{}", hf.fitted_off_mm);
        std::fs::remove_dir_all(&dir).ok();
    }

    #[test]
    fn a_ring_joined_to_its_neighbour_is_found_by_its_opening() {
        // the same board with a gold bar running into the upper-right ring: that shape is too long to read as a ring, so three gold holes, against a dark table
        let truth = [(100.0, 80.0), (1400.0, 90.0), (1420.0, 1060.0), (90.0, 1050.0)];
        let h = geom::homography(&[(0.0, 0.0), (55.0, 0.0), (55.0, 40.0), (0.0, 40.0)], &truth).unwrap();
        let hi = geom::invert(&h).unwrap();
        let img = RgbImage::from_fn(1500, 1150, |x, y| { let q = geom::apply(&hi, (x as f64 + 0.5, y as f64 + 0.5));
            if !(0.0..55.0).contains(&q.0) || !(0.0..40.0).contains(&q.1) { return image::Rgb([20, 24, 22]); }
            let c = [(3.5, 3.5), (51.5, 3.5), (51.5, 36.5), (3.5, 36.5)];
            let (mut d, mut k) = (f64::MAX, 0); for (i, p) in c.iter().enumerate() { let e = geom::dist(*p, q); if e < d { d = e; k = i; } }
            let _ = k; if d < 1.25 { image::Rgb([12, 12, 12]) } else if d < 2.8 || ((3.0..4.0).contains(&q.1) && (25.0..50.0).contains(&q.0)) { image::Rgb([230, 170, 60]) } else { image::Rgb([50, 186, 140]) } });
        let dir = std::env::temp_dir().join(format!("hole-frame-joined-{}", std::process::id())); std::fs::create_dir_all(&dir).unwrap();
        let path = dir.join("p.png"); img.save(&path).unwrap();
        let (rw, rh) = (1320u32, 970u32);
        let (r, _) = prep::rectify(&img, &truth, rw, rh, false).unwrap();
        let holes = prep::gold_holes(&r, 20.0, 200.0);
        assert_eq!(holes.len(), 3, "{holes:?}");
        let rec = json!({"corners": truth.iter().map(|p| json!([p.0, p.1])).collect::<Vec<_>>(), "rectPx": [rw, rh], "goldHoles": holes});
        let hf = hole_frame(&path, &rec, (48.0, 33.0), (55.0, 40.0)).unwrap();
        for i in 0..4 { assert!(geom::dist(hf.corners[i], truth[i]) < 2.0, "{i}: {:?} vs {:?}", hf.corners[i], truth[i]); }
        std::fs::remove_dir_all(&dir).ok();
    }

    #[test]
    fn vias_on_both_faces_recover_a_warp_between_them() {
        // dark pits with a pale ring of plain mask on both faces, the bottom's moved by a shift and a slow stretch
        let (w, h) = (1200u32, 800u32);
        let mut seed = 777u64; let mut next = move || { seed = seed.wrapping_mul(6364136223846793005).wrapping_add(1442695040888963407); (seed >> 33) as f64 };
        let dots: Vec<P> = (0..160).map(|_| (40.0 + next() % 1120.0, 40.0 + next() % 720.0)).collect();
        let moved = |p: P| (p.0 + 3.0 + 0.002 * p.0, p.1 - 2.0 + 0.001 * p.1);
        let draw = |pts: &[P]| { let mut grain = 99u64; RgbImage::from_fn(w, h, |x, y| { grain = grain.wrapping_mul(25214903917).wrapping_add(11); let n = ((grain >> 40) % 5) as u8;
            if pts.iter().any(|d| geom::dist(*d, (x as f64, y as f64)) < 3.5) { image::Rgb([24, 60, 34]) } else { image::Rgb([70 + n, 150 + n, 90 + n]) } }) };
        let (top, bottom) = (draw(&dots), draw(&dots.iter().map(|d| moved(*d)).collect::<Vec<_>>()));
        let prs = pit_pairs(&top, &bottom).expect("the pits should pair");
        assert!(prs.len() >= 100, "{}", prs.len());
        let q = crate::regframe::Quad::fit(&prs).unwrap();
        for d in dots.iter().step_by(20) { assert!(geom::dist(q.apply(*d), moved(*d)) < 1.0, "{d:?}"); }
    }

    #[test]
    fn the_mirror_that_correlates_is_the_flip() {
        // via rims scattered over plain mask, mirrored left to right on the bottom photo; the wrong mirror shares nothing
        let (w, h) = (1200u32, 800u32);
        let mut seed = 12345u64; let mut next = move || { seed = seed.wrapping_mul(6364136223846793005).wrapping_add(1442695040888963407); (seed >> 33) as i64 };
        let dots: Vec<(i64, i64)> = (0..400).map(|_| (next() % 1100 + 50, next() % 700 + 50)).collect();
        let paint = |flip: bool| RgbImage::from_fn(w, h, |x, y| { let u = if flip { w as i64 - 1 - x as i64 } else { x as i64 }; if dots.iter().any(|d| (d.0 - u).pow(2) + (d.1 - y as i64).pow(2) < 16) { image::Rgb([30, 130, 90]) } else { image::Rgb([50, 186, 140]) } });
        let (top, bx) = (paint(false), paint(true));
        let by = image::imageops::flip_vertical(&bx);
        let bm = image::imageops::flip_horizontal(&bx); // the bottom photo laid in the top view, right mirror
        let (axis, cx, _) = flip_by_detail(&top, &bm, &by);
        assert_eq!(axis.as_deref(), Some("x"), "{cx}");
    }

    #[test]
    fn driver_pitches_solve_to_one_scale() {
        // the Driver top photo: rows at 93.3 px (six of them) and one at 118.5 px
        let rows = vec![json!({"pitch": 93.2, "count": 6}), json!({"pitch": 93.2, "count": 7}), json!({"pitch": 118.5, "count": 4}), json!({"pitch": 93.5, "count": 4}), json!({"pitch": 745.4, "count": 3})];
        let c = scale_candidates(&rows);
        let k = c[0]["pxPerMm"].as_f64().unwrap();
        assert!((k - 46.6).abs() < 0.3, "{c:?}");
        assert_eq!(c[0]["rows"].as_array().unwrap().len(), 4);
    }
}