AI Flow: Reverse Engineering

//! Copper read from one photo, pixel by pixel. Rewritten for the goal that the drawn layers match the photo.
//!
//! 1. Lighting. Copper under mask is brighter than laminate, but only relative to the same part of the photo. A
//!    two-class model is fitted to the whole face: a smooth copper level and a smooth laminate level (a polynomial,
//!    then a local correction). A dark corner lowers both levels together, so the threshold follows the light, and a
//!    wide area of one class keeps a threshold from the other class's level instead of averaging itself away.
//! 2. Material. Silk, tin, domes and glints are pale and are transparent: what is under them takes the class of the
//!    nearest visible pixels. A thin glint is the lit edge of copper, so it belongs to the copper beside it.
//!    Holes and (on a face with parts) part bodies are opaque.
//! 3. Class. Each pixel is averaged along the direction it varies least in, so grain cancels without blurring across
//!    an edge or a gap, then compared to the threshold between the two levels.
//! 4. Barrier. Laminate, plus thin dark lines from an oriented filter for gaps too thin to reach laminate level.
//!    Barrier ends are walked on to close small breaks. Gaps are real: nothing crosses a barrier by itself.
//! 5. Islands. What is left between barriers is one island. Long thin necks are cut from wide copper, because a thin
//!    tentacle on a pour is a track whose gap the photo lost at one end, and a short circuit is worse than a missing
//!    link. Fragments rejoin a neighbour only across a boundary with no evidence of a gap.
use crate::raster::{blur, Mask, Plane};
use image::{Rgb, RgbImage};
use serde_json::{json, Value};
use std::collections::HashMap;

/// Mid-point between the laminate and copper levels, as a share of the distance: from the two-class model on these photos.
pub const CLASS_AT: f32 = 0.5;
/// A pour's half width, in pixels at 40 px/mm: 0.35 mm, measured on the Driver photos' narrowest stitched pours.
pub const POUR_HALF: f32 = 14.0;
/// The two levels are never closer than this share of the face's median separation.
const SEPARATION_FLOOR: f32 = 0.8;
/// A window of one class is copper unless its brightness is this close to the laminate level, as a share of the way up to copper: the light
/// falls off faster than the coarser scale can follow towards the board's edge, and that must not turn copper into laminate.
const UNIFORM_AT: f32 = 0.3;
/// Where a window is all copper, its mean is the copper level less half its grain's spread.
const GRAIN_LIFT: f32 = 0.25;
/// A line is a directional structure only when its variance is this far below the variance all around: otherwise it is grain.
const STRUCTURE: f32 = 0.55;
/// A gap line's own pixels read below this class value; brighter is a ridge of grain on copper.
const LINE_VALUE: f32 = 0.8;
/// A gap line's mean response in each shot must reach this share of the lower ridge threshold. Grain makes ridges strong in one shot and near zero in the other.
const BOTH_SHOTS: f32 = 0.7;
/// A dark blob smaller than this beside something pale, not enclosing it, is a dome's shadow or the ink of a letter: 1000 pixels is
/// 0.6 mm2, from the Driver bottom photo's domes in stitched pours.
pub const SHADOW_PX: usize = 1000;
/// A piece is really darker or lighter than what surrounds it only if the step is this share of the separation of the levels.
const REAL_STEP: f32 = 0.4;
/// A thin piece of an island this long or longer is its own conductor: 0.8 mm at 40 px/mm.
const NECK_LEN: f64 = 32.0;
/// How far round a printed line its blurred rim reaches, frame pixels: 2, the width over which a white line's edge falls to the mask's level at 40 px/mm on the hub board's photos (76.8 px/mm
/// halved) and the Driver's.
const PRINT_RIM: f32 = 2.0;
/// See `padjoin` in `trace`: a pixel this far below the laminate level, as a share of the separation of the levels, is in shadow.
const PAD_SHADOW: f32 = 0.3;
/// Half width at which a connection between two wide parts counts as a pinch: 0.125 mm at 40 px/mm.
const PINCH_HALF: f32 = 5.0;

/// `hi` and `lo` are fractions of the photo's 99th percentile line response. `bodies`: this face carries parts.
/// `board`: x0, y0, x1, y1, corner radius of the board in pixels, when known.
pub struct Params<'a> { pub hi: f32, pub lo: f32, pub minlen: usize, pub ext: i64, pub min_area: usize, pub bodies: bool, pub board: Option<(f64, f64, f64, f64, f64)>, /// the other photo of this face, laid on the first: its grain is independent, so the mean of the two has less
    pub second: Option<&'a RgbImage>,
    /// domes read on this face, in pixels: the dark blob a dome casts is no gap in the copper
    pub domes: &'a [(f64, f64)],
    /// towards the face's one light, a unit vector in pixels, when its bumps agree on one
    pub light: Option<(f64, f64)>,
    /// the pads of the parts placed on this face, in frame pixels: each pad's copper is started from the pad (`padjoin`)
    pub pads: &'a [crate::padjoin::Seed] }
impl Default for Params<'_> { fn default() -> Self { Params { hi: 0.75, lo: 0.45, minlen: 30, ext: 26, min_area: 30, bodies: false, board: None, second: None, domes: &[], light: None, pads: &[] } } }

/// What became of the labels as first read.
#[derive(Clone, Debug, Default)]
pub struct Steps { /// how the levels sit against the face's own step edges, and whether the levels had to be taken from them (`edgelevels`)
    pub edge_check: Option<crate::edgelevels::Agreement>, pub levels_from_edges: bool, /// the share of the face whose levels were taken from the edges near it (`locallevels`)
    pub levels_local: f32, /// copper as a share of what can be seen
    pub copper_share: f32, /// pads that reached copper (`padjoin`)
    pub pads_joined: usize, /// track pieces joined in line (`trackbridge`)
    pub bridged: usize,
    pub first: Vec<u32>, pub crumbs: usize, pub silk_only: Vec<u32>, pub off_board: Vec<u32>, pub necks_cut: usize, pub pinches_cut: usize, pub joined: usize, pub refused_joins: usize }

pub struct Traced { pub w: usize, pub h: usize, pub resp: Plane, pub gap: Mask, pub labels: Vec<u32>, pub regions: Vec<Region>, pub other: Mask, pub ratio: Plane, pub dark: Mask, pub silk: Mask, pub pale: Mask, pub steps: Steps }

#[derive(Clone, Debug)]
pub struct Region { pub id: u32, pub area: usize, pub cx: f64, pub cy: f64, pub x0: usize, pub y0: usize, pub x1: usize, pub y1: usize, pub level: f32, pub half_width: f32, pub other_share: f32 }

/// Least squares for a polynomial of `deg` in x, y scaled to the image.
fn poly_terms(x: f64, y: f64, deg: usize, out: &mut Vec<f64>) { out.clear(); let (mut xp, n) = (1.0, deg + 1); for i in 0..n { let mut t = xp; for _ in 0..n - i { out.push(t); t *= y; } xp *= x; } }
fn poly_fit(xs: &[(f64, f64, f64)], deg: usize) -> Vec<f64> {
    let n = (deg + 1) * (deg + 2) / 2; let mut a = vec![0f64; n * n]; let mut b = vec![0f64; n]; let mut t = Vec::with_capacity(n);
    for (x, y, z) in xs { poly_terms(*x, *y, deg, &mut t); for i in 0..n { b[i] += t[i] * z; for j in 0..n { a[i * n + j] += t[i] * t[j]; } } }
    for i in 0..n { a[i * n + i] += 1e-6 * (xs.len() as f64 + 1.0); }
    for c in 0..n { let mut p = c; for r in c + 1..n { if a[r * n + c].abs() > a[p * n + c].abs() { p = r; } }
        if p != c { for k in 0..n { a.swap(c * n + k, p * n + k); } b.swap(c, p); }
        let d = a[c * n + c]; if d.abs() < 1e-12 { continue; }
        for r in c + 1..n { let f = a[r * n + c] / d; if f == 0.0 { continue; } for k in c..n { a[r * n + k] -= f * a[c * n + k]; } b[r] -= f * b[c]; } }
    let mut x = vec![0f64; n];
    for c in (0..n).rev() { let mut v = b[c]; for k in c + 1..n { v -= a[c * n + k] * x[k]; } x[c] = if a[c * n + c].abs() > 1e-12 { v / a[c * n + c] } else { 0.0 }; }
    x
}
fn poly_eval(c: &[f64], x: f64, y: f64, deg: usize, t: &mut Vec<f64>) -> f64 { poly_terms(x, y, deg, t); c.iter().zip(t.iter()).map(|(a, b)| a * b).sum() }

/// The two levels over the face, in log green, from coarse to fine. In each window the pixels are split in two by a local
/// two-means. Where the split is real, those are the levels. Where it is only grain, the window is one class: it takes the class
/// its brightness is nearer, keeps the separation from the coarser scale, and so the levels follow the light across a wide
/// uniform area. The coarsest scale falls back on a polynomial fitted to the whole face.
pub fn two_levels(lg: &Plane, excl: &Mask) -> (Plane, Plane) { two_levels_at(lg, excl, &[360, 180, 90], 4) }

/// How much the light changes over a face: the spread of the copper level (85th percentile of log green, seen pixels) over 8 mm tiles, 90th less 10th percentile.
/// Reads 0.32 and 0.56 on the Driver top and bottom, whose light falls by up to half across the board, and 0.18 and 0.11 on the studio-lit CM5 everypin photos (`TRACE_DEBUG` prints it).
pub fn light_spread(lg: &Plane, excl: &Mask) -> f32 {
    let (w, h) = (lg.w, lg.h); let t = 320usize; let mut lv = vec![];
    let mut y = 0; while y < h { let mut x = 0; while x < w {
        let mut v = vec![];
        for yy in (y..(y + t).min(h)).step_by(3) { for xx in (x..(x + t).min(w)).step_by(3) { if !excl.at(xx, yy) { v.push(lg.at(xx, yy)); } } }
        if v.len() > 800 { let k = v.len() * 85 / 100; let (_, e, _) = v.select_nth_unstable_by(k, |a, b| a.partial_cmp(b).unwrap()); lv.push(*e); }
        x += t; } y += t; }
    if lv.len() < 4 { return 9.0; }
    lv.sort_by(|a, b| a.partial_cmp(b).unwrap());
    lv[lv.len() * 9 / 10] - lv[lv.len() / 10]
}
/// The two modes of log green over plain mask, when it has two: the highest bin of a 0.01 wide smoothed histogram, and the highest on the other side at least `MODES_APART` away with a valley between
/// no higher than 80 percent of the lower peak, and at least 5 percent of the higher, so a stray pile-up in a dark tail is not a mode. The CM5 faces read 172 and 202, 177 and 207 (separation 0.16).
pub fn two_modes(v: &[f32]) -> Option<(f32, f32)> {
    if v.len() < 1000 { return None; }
    let (lo, hi) = v.iter().fold((f32::MAX, f32::MIN), |a, x| (a.0.min(*x), a.1.max(*x)));
    let n = (((hi - lo) / 0.01) as usize + 1).min(800);
    let mut h = vec![0f32; n]; for x in v { h[(((x - lo) / 0.01) as usize).min(n - 1)] += 1.0; }
    let s: Vec<f32> = (0..n).map(|i| { let (a, b) = (i.saturating_sub(2), (i + 3).min(n)); h[a..b].iter().sum::<f32>() / (b - a) as f32 }).collect();
    let p1 = (0..n).max_by(|a, b| s[*a].partial_cmp(&s[*b]).unwrap())?;
    let apart = (MODES_APART / 0.01) as usize;
    let cands: Vec<usize> = (0..n).filter(|i| i.abs_diff(p1) >= apart && (*i == 0 || s[*i] >= s[*i - 1]) && (*i + 1 == n || s[*i] >= s[*i + 1])).collect();
    let p2 = cands.into_iter().filter(|i| { let (a, b) = (p1.min(*i), p1.max(*i)); let valley = s[a..=b].iter().copied().fold(f32::MAX, f32::min); s[*i] >= 0.05 * s[p1] && valley <= 0.8 * s[*i].min(s[p1]) }).max_by(|a, b| s[*a].partial_cmp(&s[*b]).unwrap())?;
    let val = |i: usize| lo + (i as f32 + 0.5) * 0.01;
    Some((val(p1.min(p2)), val(p1.max(p2))))
}
const MODES_APART: f32 = 0.08;

/// Half the width at which pale and neutral stops being print, frame pixels: 0.6 mm. The strokes of the largest legend on the CM5 top are 0.5 mm wide, its metal parts 3 mm and up.
const PALE_BODY_HALF: f32 = 24.0;

/// How far a sharpened photo's dark band round bright things reaches, frame pixels, 0 when it has none to speak of. Median green of seen pixels against distance from what is hidden. A band counts at a dip of 10 percent
/// below the plateau that comes back to within 4 percent by `SHARP_MAX_PX`: a filter's band is narrow, while a cast shadow, such as those round the Driver's tall parts, runs on for 13 px and more and is not this.
pub fn sharpening_halo(green: &Plane, hidden: &Mask, pale: &Mask) -> usize {
    let dist = hidden.distance_to_set();
    let mut bins: Vec<Vec<f32>> = vec![vec![]; 25];
    for i in 0..green.d.len() { let d = dist.d[i]; if d >= 1.0 && d < 25.0 && !pale.d[i] { bins[d as usize].push(green.d[i]); } }
    let med: Vec<f32> = bins.iter_mut().map(|v| median(v)).collect();
    let mut tail: Vec<f32> = med[18..25].to_vec(); let plateau = median(&mut tail);
    if plateau <= 0.0 { return 0; }
    let low = med[1..12].iter().copied().filter(|v| *v > 0.0).fold(f32::MAX, f32::min);
    if low > 0.90 * plateau { return 0; }
    let reach = (1..18).rev().find(|d| med[*d] > 0.0 && med[*d] < 0.96 * plateau).unwrap_or(0);
    if reach > SHARP_MAX_PX { 0 } else { reach }
}
/// CM5 bottom, everypin: back within 4 percent at 5 px, 3 photo pixels at 24 px/mm.
const SHARP_MAX_PX: usize = 8;

/// Below this spread the light is flat enough that a level needs no following at fine scales. Local levels exist to follow the light, and where the light does not move they
/// follow the copper instead: a wide field of laminate with thin tracks through it splits into its own two classes and reads as copper (CM5 bottom, at a 2.25 mm window).
pub const FLAT_LIGHT: f32 = 0.2;

/// `deg`: the polynomial's degree. A flat light needs a plane at most, and a higher degree swings about where parts leave no plain mask to hold it.
pub fn two_levels_at(lg: &Plane, excl: &Mask, scales: &[usize], deg: usize) -> (Plane, Plane) {
    #[allow(non_snake_case)] let DEG = deg; let (w, h) = (lg.w, lg.h); let sc = 3000.0;
    let mut sm: Vec<(f64, f64, f64)> = vec![];
    for y in (0..h).step_by(6) { for x in (0..w).step_by(6) { if !excl.at(x, y) { sm.push((x as f64 / sc, y as f64 / sc, lg.at(x, y) as f64)); } } }
    let mut t = vec![];
    if sm.len() < 100 { let m = lg.d.iter().sum::<f32>() / lg.d.len().max(1) as f32; return (Plane::new(w, h, m + 0.1), Plane::new(w, h, m - 0.1)); }
    let c0 = poly_fit(&sm, DEG);
    let res: Vec<f64> = sm.iter().map(|s| s.2 - poly_eval(&c0, s.0, s.1, DEG, &mut t)).collect();
    let mut sorted = res.clone(); sorted.sort_by(|a, b| a.partial_cmp(b).unwrap()); let mid = sorted[sorted.len() / 2];
    let mut hi: Vec<bool> = res.iter().map(|r| *r > mid).collect();
    let (mut ch, mut cl) = (c0.clone(), c0.clone());
    for _ in 0..12 {
        let a: Vec<_> = sm.iter().zip(&hi).filter(|x| *x.1).map(|x| *x.0).collect(); let b: Vec<_> = sm.iter().zip(&hi).filter(|x| !*x.1).map(|x| *x.0).collect();
        if a.len() < 30 || b.len() < 30 { break; }
        ch = poly_fit(&a, DEG); cl = poly_fit(&b, DEG);
        let next: Vec<bool> = sm.iter().map(|s| { let (u, v) = (poly_eval(&ch, s.0, s.1, DEG, &mut t), poly_eval(&cl, s.0, s.1, DEG, &mut t)); s.2 > (u + v) / 2.0 }).collect();
        let same = next.iter().zip(&hi).filter(|x| x.0 == x.1).count(); hi = next; if same * 2000 > hi.len() * 1999 { break; }
    }
    let (mut cop, mut lam) = (Plane::new(w, h, 0.0), Plane::new(w, h, 0.0));
    for y in 0..h { for x in 0..w { let (a, b) = (x as f64 / sc, y as f64 / sc); cop.d[y * w + x] = poly_eval(&ch, a, b, DEG, &mut t) as f32; lam.d[y * w + x] = poly_eval(&cl, a, b, DEG, &mut t) as f32; } }
    // the face's own separation sets how far apart the levels must be for a split to be real
    let mut sep: Vec<f32> = (0..w * h).step_by(97).filter(|i| !excl.d[*i]).map(|i| cop.d[i] - lam.d[i]).collect(); sep.sort_by(|a, b| a.partial_cmp(b).unwrap());
    let floor = (sep.get(sep.len() / 2).copied().unwrap_or(0.2) * SEPARATION_FLOOR).max(0.08);
    for i in 0..w * h { if cop.d[i] - lam.d[i] < floor { let m = (cop.d[i] + lam.d[i]) / 2.0; cop.d[i] = m + floor / 2.0; lam.d[i] = m - floor / 2.0; } }
    let wt = Plane { w, h, d: excl.d.iter().map(|e| if *e { 0.0 } else { 1.0 }).collect() };
    for &r in scales {
        // two-means in the window, started from the levels so far
        let (mut u, mut d) = (cop.clone(), lam.clone());
        let bw = blur(&wt, r, 1);
        for _ in 0..3 {
            let (mut su, mut wu, mut sd, mut wd) = (Plane::new(w, h, 0.0), Plane::new(w, h, 0.0), Plane::new(w, h, 0.0), Plane::new(w, h, 0.0));
            for i in 0..w * h { if excl.d[i] { continue; } if lg.d[i] > (u.d[i] + d.d[i]) / 2.0 { su.d[i] = lg.d[i]; wu.d[i] = 1.0; } else { sd.d[i] = lg.d[i]; wd.d[i] = 1.0; } }
            let (bsu, bwu, bsd, bwd) = (blur(&su, r, 1), blur(&wu, r, 1), blur(&sd, r, 1), blur(&wd, r, 1));
            for i in 0..w * h {
                // a class with almost no pixels here keeps its previous level, moved with the other
                let (fu, fd) = (bwu.d[i] / bw.d[i].max(1e-3), bwd.d[i] / bw.d[i].max(1e-3));
                let (mut nu, mut nd) = (u.d[i], d.d[i]);
                if fu > 0.03 { nu = bsu.d[i] / bwu.d[i]; } if fd > 0.03 { nd = bsd.d[i] / bwd.d[i]; }
                if fu <= 0.03 && fd > 0.03 { nu = nd + (u.d[i] - d.d[i]); } if fd <= 0.03 && fu > 0.03 { nd = nu - (u.d[i] - d.d[i]); }
                u.d[i] = nu; d.d[i] = nd;
            }
        }
        let (mut c2, mut l2) = (cop.clone(), lam.clone());
        for i in 0..w * h {
            if bw.d[i] < 0.05 { continue; }
            let (cs, ls) = (cop.d[i], lam.d[i]); let was = (cs - ls).max(floor); let s = u.d[i] - d.d[i];
            let real = ((s / was - 0.45) / 0.3).clamp(0.0, 1.0);
            // one class here: which one, by nearness to the coarser levels; the local mean of both halves is the class's level
            let m = (u.d[i] + d.d[i]) / 2.0;
            let (uc, ul) = if m > ls + UNIFORM_AT * (cs - ls) { (m + GRAIN_LIFT * s, m + GRAIN_LIFT * s - was) } else { (m + was, m) };
            c2.d[i] = real * u.d[i] + (1.0 - real) * uc; l2.d[i] = real * d.d[i] + (1.0 - real) * ul;
        }
        cop = blur(&c2, r / 3, 2); lam = blur(&l2, r / 3, 2);
        for i in 0..w * h { if cop.d[i] - lam.d[i] < floor { lam.d[i] = cop.d[i] - floor; } }
    }
    (cop, lam)
}

/// Each pixel's mean along the line through it with the least variance, over the pixels that can be seen. With two shots
/// the direction comes from one and the mean from the other, both ways round. Grain differs between the shots (correlation
/// 0.08 at pixel level on the Driver top, 0.06 on the bottom), so a line chosen because grain happens to run along it in one
/// shot has no such run in the other, and the choice cannot draw an edge that is not there.
fn along_least_variance(p: &Plane, q: &Plane, skip: &Mask, half: i64) -> Plane {
    let (w, h) = (p.w, p.h); const NA: usize = 16;
    // what the whole neighbourhood says, over the seen pixels: the answer wherever no line stands out
    const ISO: usize = 6;
    let stats = |p: &Plane| { let (mut a, mut a2, mut c) = (p.clone(), p.clone(), Plane::new(w, h, 1.0));
        for i in 0..w * h { if skip.d[i] { a.d[i] = 0.0; a2.d[i] = 0.0; c.d[i] = 0.0; } else { a2.d[i] = p.d[i] * p.d[i]; } }
        (blur(&a, ISO, 1), blur(&a2, ISO, 1), blur(&c, ISO, 1)) };
    let (pa, pa2, pc) = stats(p); let (qa, qa2, _) = stats(q);
    let lines: Vec<Vec<(i64, i64)>> = (0..NA).map(|k| { let a = k as f64 * std::f64::consts::PI / NA as f64; (-half..=half).map(|t| ((t as f64 * a.cos()).round() as i64, (t as f64 * a.sin()).round() as i64)).collect() }).collect();
    let flat: Vec<Vec<isize>> = lines.iter().map(|l| l.iter().map(|q| q.0 as isize + q.1 as isize * w as isize).collect()).collect();
    let threads = std::thread::available_parallelism().map(|v| v.get()).unwrap_or(4).clamp(1, 8);
    let mut out = Plane::new(w, h, 0.0);
    let m = half as usize + 1;
    std::thread::scope(|sc| {
        for (ci, chunk) in out.d.chunks_mut((h / threads + 1) * w).enumerate() {
            let flat = &flat; let lines = &lines; let (pa, pa2, pc, qa, qa2) = (&pa, &pa2, &pc, &qa, &qa2);
            sc.spawn(move || {
                let y0 = ci * (h / threads + 1);
                for (off, o) in chunk.iter_mut().enumerate() {
                    let (x, y) = (off % w, y0 + off / w); let i = y * w + x;
                    // (variance, mean) along the best line so far, in each shot; the mean of the other shot along that line
                    let (mut bp, mut bq) = (f32::MAX, f32::MAX); let (mut p_on_q, mut q_on_p) = (p.d[i], q.d[i]);
                    let inner = x >= m && y >= m && x + m < w && y + m < h;
                    for a in 0..NA {
                        let (mut s, mut s2, mut t, mut t2, mut n) = (0f32, 0f32, 0f32, 0f32, 0f32);
                        if inner { for d in &flat[a] { let j = (i as isize + d) as usize; if !skip.d[j] { let (u, v) = (p.d[j], q.d[j]); s += u; s2 += u * u; t += v; t2 += v * v; n += 1.0; } } }
                        else { for l in &lines[a] { let (px, py) = (x as i64 + l.0, y as i64 + l.1); if px < 0 || py < 0 || px >= w as i64 || py >= h as i64 { continue; } let j = py as usize * w + px as usize; if !skip.d[j] { let (u, v) = (p.d[j], q.d[j]); s += u; s2 += u * u; t += v; t2 += v * v; n += 1.0; } } }
                        if n < 5.0 { continue; }
                        let (vp, vq) = (s2 / n - (s / n) * (s / n), t2 / n - (t / n) * (t / n));
                        if vp < bp { bp = vp; p_on_q = t / n; }
                        if vq < bq { bq = vq; q_on_p = s / n; }
                    }
                    let n = pc.d[i];
                    if n > 0.2 {
                        let (mp, mq) = (pa.d[i] / n, qa.d[i] / n);
                        let (vp, vq) = (pa2.d[i] / n - mp * mp, qa2.d[i] / n - mq * mq);
                        if bp >= STRUCTURE * vp { p_on_q = mq; }
                        if bq >= STRUCTURE * vq { q_on_p = mp; }
                    }
                    *o = (p_on_q + q_on_p) / 2.0;
                }
            });
        }
    });
    out
}

/// Dark-line response at each pixel and the best direction. `n` is clipped brightness; `skip` pixels read as neutral.
fn line_filter(n: &Plane, skip: &Mask) -> (Plane, Vec<u8>) {
    let (w, h) = (n.w, n.h);
    const NA: usize = 16;
    let threads = std::thread::available_parallelism().map(|v| v.get()).unwrap_or(4).clamp(1, 8);
    let run = |k: usize| -> Plane {
        let a = k as f64 * std::f64::consts::PI / NA as f64;
        let (ca, sa) = (a.cos(), a.sin());
        // a dark line: a Gaussian along it, the second derivative of a Gaussian across it, zero mean
        let (sig_across, sig_along) = (2.6f64, 9.0f64);
        let rad = 27i64;
        let mut kern: Vec<(i64, i64, f32)> = vec![];
        for dy in -rad..=rad { for dx in -rad..=rad {
            let (u, v) = (dx as f64 * ca + dy as f64 * sa, -(dx as f64) * sa + dy as f64 * ca);
            let g = (-u * u / (2.0 * sig_along * sig_along)).exp() * (1.0 - v * v / (sig_across * sig_across)) * (-v * v / (2.0 * sig_across * sig_across)).exp();
            if g.abs() > 0.012 { kern.push((dx, dy, g as f32)); }
        } }
        let mean = kern.iter().map(|k| k.2).sum::<f32>() / kern.len() as f32;
        let norm = kern.iter().map(|k| (k.2 - mean).abs()).sum::<f32>() / 2.0;
        for k in kern.iter_mut() { k.2 = -(k.2 - mean) / norm; }
        let mut r = Plane::new(w, h, 0.0);
        let flat: Vec<(isize, f32)> = kern.iter().map(|k| (k.0 as isize + k.1 as isize * w as isize, k.2)).collect();
        let m = rad as usize;
        for y in m..h.saturating_sub(m) { for x in m..w.saturating_sub(m) { let i = y * w + x; if skip.d[i] { continue; }
            let ii = i as isize; let mut v = 0f32; for (o, g) in &flat { v += *g * n.d[(ii + *o) as usize]; }
            r.d[i] = v; } }
        r
    };
    let mut all: Vec<Option<Plane>> = (0..NA).map(|_| None).collect();
    for chunk in (0..NA).collect::<Vec<_>>().chunks(threads) {
        let res: Vec<(usize, Plane)> = std::thread::scope(|sc| { let hs: Vec<_> = chunk.iter().map(|k| { let k = *k; let run = &run; sc.spawn(move || (k, run(k))) }).collect(); hs.into_iter().map(|h| h.join().unwrap()).collect() });
        for (k, p) in res { all[k] = Some(p); }
    }
    let mut best = Plane::new(w, h, 0.0); let mut ang = vec![0u8; w * h];
    for (k, p) in all.iter().enumerate() { let p = p.as_ref().unwrap(); for i in 0..w * h { if p.d[i] > best.d[i] { best.d[i] = p.d[i]; ang[i] = k as u8; } } }
    (best, ang)
}

/// A class is only as good as the step at its edge. A wide patch read as laminate with no real step to what surrounds it, and a
/// small patch read as copper with none, were put there by the threshold wandering through grain, not by anything in the photo.
fn settle(seed: &mut [i8], lg: &Plane, cop: &Plane, lam: &Plane, lines: &Mask, hidden: &Mask, unseen: &Mask, flat_light: bool, w: usize, h: usize) {
    let step_of = |m: &Mask, sign: f32| -> (Vec<u32>, Vec<f32>, Vec<usize>) {
        let (lab, cs) = crate::raster::components(m); let mut out = vec![f32::MAX; cs.len() + 1]; let mut areas = vec![0usize; cs.len() + 1];
        let near = m.dilate(3.0); let far = m.dilate(9.0);
        let (mut a, mut ca, mut b, mut cb, mut sp) = (vec![0f32; cs.len() + 1], vec![0f32; cs.len() + 1], vec![0f32; cs.len() + 1], vec![0f32; cs.len() + 1], vec![0f32; cs.len() + 1]);
        // the ring around a piece belongs to the nearest piece: spread labels outward a pixel at a time, 9 rounds
        let mut spread = lab.clone();
        for _ in 0..9 { let prev = spread.clone(); for y in 1..h - 1 { for x in 1..w - 1 { let i = y * w + x; if prev[i] != 0 || !far.d[i] { continue; } for j in [i - 1, i + 1, i - w, i + w] { if prev[j] != 0 { spread[i] = prev[j]; break; } } } } }
        for i in 0..w * h { let l = spread[i] as usize; if l == 0 || unseen.d[i] { continue; } if lab[i] != 0 { a[l] += lg.d[i]; ca[l] += 1.0; sp[l] += cop.d[i] - lam.d[i]; } else if !near.d[i] { b[l] += lg.d[i]; cb[l] += 1.0; } }
        for c in &cs { let l = c.label as usize; areas[l] = c.area; if ca[l] >= 8.0 && cb[l] >= 8.0 { out[l] = sign * (b[l] / cb[l] - a[l] / ca[l]) / (sp[l] / ca[l]).max(0.05); } }
        (lab, out, areas)
    };
    // wide laminate: closed over grain, opened to keep off thin gaps
    let dark = Mask { w, h, d: seed.iter().zip(&hidden.d).map(|(s, o)| *s == -1 && !*o).collect() };
    let wide = dark.close(3.0).open(6.0);
    let (lab, _, _) = step_of(&wide, 1.0);
    // Laminate has copper on some side of it, brighter by the separation of the levels. A dim patch with nothing much brighter around it, however dark the gap on its other
    // side, is copper in poor light: the 90th percentile of its surroundings is the side that counts, so a dark gap beside it cannot lower its step (Driver bottom, PROG box:
    // a 10 percent drop at a silk line and a real gap inside the box). Gap lines along most of its edge still mark a lane between tracks. Under flat light there is no poor light to excuse a dim patch.
    if !flat_light { let n = lab.iter().copied().max().unwrap_or(0) as usize + 1;
      let mut ring: Vec<Vec<f32>> = vec![vec![]; n]; let (mut a, mut ca, mut sp) = (vec![0f32; n], vec![0f32; n], vec![0f32; n]);
      let far = wide.dilate(10.0); let near = wide.dilate(3.0);
      let mut spread = lab.clone();
      for _ in 0..10 { let prev = spread.clone(); for y in 1..h - 1 { for x in 1..w - 1 { let i = y * w + x; if prev[i] != 0 || !far.d[i] { continue; } for j in [i - 1, i + 1, i - w, i + w] { if prev[j] != 0 { spread[i] = prev[j]; break; } } } } }
      for i in 0..w * h { let l = spread[i] as usize; if l == 0 || unseen.d[i] { continue; } if lab[i] != 0 { a[l] += lg.d[i]; ca[l] += 1.0; sp[l] += cop.d[i] - lam.d[i]; } else if !near.d[i] { ring[l].push(lg.d[i]); } }
      let (near_a, edge) = (lines.dilate(4.0), wide.andnot(&wide.erode(1.0)));
      let (mut on, mut all) = (vec![0u32; n], vec![0u32; n]);
      for i in 0..w * h { if edge.d[i] { let l = lab[i] as usize; all[l] += 1; if near_a.d[i] { on[l] += 1; } } }
      let up: Vec<f32> = (0..n).map(|l| { if ca[l] < 8.0 || ring[l].len() < 8 { return f32::MAX; } let k = ring[l].len() * 9 / 10; let (_, top, _) = ring[l].select_nth_unstable_by(k, |p, q| p.partial_cmp(q).unwrap()); (*top - a[l] / ca[l]) / (sp[l] / ca[l]).max(0.05) }).collect();
      for i in 0..w * h { let l = lab[i] as usize; if l != 0 && up[l] < REAL_STEP && seed[i] == -1 && !lines.d[i] && on[l] * 2 < all[l] { seed[i] = 1; } } }
    // small copper: not touching a gap line, which a real conductor always has along it
    let bright = Mask { w, h, d: seed.iter().map(|s| *s == 1).collect() };
    let (lab, step, areas) = step_of(&bright, -1.0);
    let near_line = lines.dilate(3.0);
    let mut lined = vec![false; areas.len()]; for i in 0..w * h { if lab[i] != 0 && near_line.d[i] { lined[lab[i] as usize] = true; } }
    for i in 0..w * h { let l = lab[i] as usize; if l != 0 && areas[l] < 1600 && !lined[l] && step[l] < REAL_STEP { seed[i] = -1; } }
}

/// Every unseen pixel takes the class of its nearest seen pixel (two-pass chamfer sweep). `seed`: 1 copper, -1 laminate, 0 unseen.
fn nearest_class(seed: &[i8], w: usize, h: usize) -> Vec<i8> {
    let mut d: Vec<f32> = seed.iter().map(|s| if *s != 0 { 0.0 } else { 1e9 }).collect(); let mut c = seed.to_vec();
    let mut relax = |i: usize, j: usize, cost: f32, d: &mut Vec<f32>, c: &mut Vec<i8>| { if d[j] + cost < d[i] { d[i] = d[j] + cost; c[i] = c[j]; } };
    for y in 0..h { for x in 0..w { let i = y * w + x;
        if x > 0 { relax(i, i - 1, 3.0, &mut d, &mut c); }
        if y > 0 { relax(i, i - w, 3.0, &mut d, &mut c); if x > 0 { relax(i, i - w - 1, 4.0, &mut d, &mut c); } if x + 1 < w { relax(i, i - w + 1, 4.0, &mut d, &mut c); } } } }
    for y in (0..h).rev() { for x in (0..w).rev() { let i = y * w + x;
        if x + 1 < w { relax(i, i + 1, 3.0, &mut d, &mut c); }
        if y + 1 < h { relax(i, i + w, 3.0, &mut d, &mut c); if x + 1 < w { relax(i, i + w + 1, 4.0, &mut d, &mut c); } if x > 0 { relax(i, i + w - 1, 4.0, &mut d, &mut c); } } } }
    c
}

/// Connected pieces of `keep`, within one label: a piece never crosses into another label.
fn pieces_within(labels: &[u32], keep: &[bool], w: usize, h: usize) -> (Vec<u32>, u32) {
    let mut out = vec![0u32; w * h]; let mut n = 0u32; let mut st = vec![];
    for s in 0..w * h { if !keep[s] || labels[s] == 0 || out[s] != 0 { continue; } n += 1; out[s] = n; st.push(s); let l = labels[s];
        while let Some(i) = st.pop() { let (x, y) = (i % w, i / w);
            let mut go = |j: usize| { if keep[j] && labels[j] == l && out[j] == 0 { out[j] = n; st.push(j); } };
            if x > 0 { go(i - 1); } if x + 1 < w { go(i + 1); } if y > 0 { go(i - w); } if y + 1 < h { go(i + w); } } }
    (out, n)
}

/// Cut thin tentacles and pinches off wide copper. Returns necks cut and pinches cut.
pub fn split_necks(labels: &mut [u32], w: usize, h: usize) -> (usize, usize) {
    let on = Mask { w, h, d: labels.iter().map(|l| *l != 0).collect() };
    let dt = on.not().distance_to_set();
    let wide = Mask { w, h, d: dt.d.iter().map(|v| *v >= POUR_HALF).collect() };
    let core = wide.dilate(POUR_HALF).and(&on);
    let thin: Vec<bool> = (0..w * h).map(|i| on.d[i] && !core.d[i]).collect();
    let (tl, tn) = pieces_within(labels, &thin, w, h);
    let mut area = vec![0usize; tn as usize + 1]; let mut sum = vec![0f64; tn as usize + 1]; let mut bb = vec![(usize::MAX, usize::MAX, 0usize, 0usize); tn as usize + 1]; let mut touches = vec![false; tn as usize + 1];
    for y in 0..h { for x in 0..w { let i = y * w + x; let t = tl[i] as usize; if t == 0 { continue; } area[t] += 1; sum[t] += dt.d[i] as f64; let b = &mut bb[t]; b.0 = b.0.min(x); b.1 = b.1.min(y); b.2 = b.2.max(x); b.3 = b.3.max(y);
        for j in [i.wrapping_sub(1), i + 1, i.wrapping_sub(w), i + w] { if j < w * h && core.d[j] && labels[j] == labels[i] { touches[t] = true; } } } }
    let mut next = labels.iter().copied().max().unwrap_or(0);
    let mut newl = vec![0u32; tn as usize + 1]; let mut necks = 0;
    for t in 1..=tn as usize { if !touches[t] { continue; }
        let (b, mean) = (bb[t], sum[t] / area[t].max(1) as f64);
        let diag = (((b.2 - b.0) as f64).powi(2) + ((b.3 - b.1) as f64).powi(2)).sqrt();
        let len = (area[t] as f64 / (2.0 * mean.max(1.0))).max(diag);
        if len >= NECK_LEN && area[t] as f64 >= NECK_LEN * 3.0 { next += 1; newl[t] = next; necks += 1; } }
    for i in 0..w * h { let t = tl[i] as usize; if t != 0 && newl[t] != 0 { labels[i] = newl[t]; } }
    // what is left of a label may now be in pieces: each piece is its own island
    let (pl, pn) = pieces_within(labels, &on.d, w, h);
    let mut seen: HashMap<u32, u32> = HashMap::new(); let mut map = vec![0u32; pn as usize + 1];
    for i in 0..w * h { let p = pl[i] as usize; if p == 0 || map[p] != 0 { continue; } let l = labels[i]; if seen.insert(l, p as u32).is_none() { map[p] = l; } else { next += 1; map[p] = next; } }
    for i in 0..w * h { let p = pl[i] as usize; if p != 0 { labels[i] = map[p]; } }
    // a pinch: remove what is thinner than the pinch width; if two substantial parts fall apart, the join was a leak
    let on2 = Mask { w, h, d: labels.iter().map(|l| *l != 0).collect() };
    let dt2 = on2.not().distance_to_set();
    let body = Mask { w, h, d: dt2.d.iter().map(|v| *v >= PINCH_HALF).collect() }.dilate(PINCH_HALF).and(&on2);
    let (bl, bn) = pieces_within(labels, &body.d, w, h);
    let mut barea = vec![0usize; bn as usize + 1]; for b in &bl { barea[*b as usize] += 1; }
    let large = 400usize; // 0.25 mm2
    let mut by_label: HashMap<u32, Vec<u32>> = HashMap::new();
    for i in 0..w * h { let b = bl[i]; if b != 0 && barea[b as usize] >= large { let e = by_label.entry(labels[i]).or_default(); if !e.contains(&b) { e.push(b); } } }
    let mut pinches = 0;
    let parents: std::collections::HashSet<u32> = by_label.iter().filter(|x| x.1.len() >= 2).map(|x| *x.0).collect();
    let mut px_of: HashMap<u32, Vec<usize>> = HashMap::new();
    for i in 0..w * h { if parents.contains(&labels[i]) { px_of.entry(labels[i]).or_default().push(i); } }
    let mut seen_px = vec![false; w * h];
    let mut parent_list: Vec<u32> = parents.iter().copied().collect(); parent_list.sort();
    for l in parent_list {
        // the first part keeps the label, the others get their own; loose pixels go to whatever part they touch, or their own island when two
        let mut part_label: HashMap<u32, u32> = HashMap::new();
        for (k, b) in by_label[&l].iter().enumerate() { if k == 0 { part_label.insert(*b, l); } else { next += 1; part_label.insert(*b, next); } }
        let pix = &px_of[&l];
        let is_rest = |i: usize, labels: &[u32]| labels[i] == l && !part_label.contains_key(&bl[i]);
        let mut pieces: Vec<(Vec<usize>, Vec<u32>)> = vec![]; let mut st = vec![];
        for &s in pix { if seen_px[s] || !is_rest(s, labels) { continue; } seen_px[s] = true; st.push(s); let (mut these, mut hit) = (vec![], vec![]);
            while let Some(i) = st.pop() { these.push(i);
                for j in [i.wrapping_sub(1), i + 1, i.wrapping_sub(w), i + w] { if j >= w * h { continue; } if labels[j] != l { continue; }
                    if let Some(q) = part_label.get(&bl[j]) { if !hit.contains(q) { hit.push(*q); } } else if !seen_px[j] { seen_px[j] = true; st.push(j); } } }
            pieces.push((these, hit)); }
        for &i in pix { if let Some(q) = part_label.get(&bl[i]) { labels[i] = *q; } }
        for (these, hit) in pieces { let q = if hit.len() == 1 { hit[0] } else { next += 1; pinches += 1; next }; for i in these { labels[i] = q; } }
    }
    (necks, pinches)
}

/// Rejoin fragments across a boundary with no evidence of a gap. `value` is the class value (0 laminate, 1 copper level),
/// `strong` the line response that marks a real gap. A join never connects two substantial pieces.
pub fn join_fragments(labels: &mut [u32], value: &Plane, resp: &Plane, strong: f32, other: &Mask, w: usize, h: usize) -> (usize, usize) {
    let mut area: HashMap<u32, usize> = HashMap::new(); for l in labels.iter() { if *l != 0 { *area.entry(*l).or_default() += 1; } }
    let (mut weak, mut hard): (HashMap<(u32, u32), u32>, HashMap<(u32, u32), u32>) = (HashMap::new(), HashMap::new());
    const REACH: usize = 6;
    for y in 0..h { for x in 0..w { let i = y * w + x; let a = labels[i]; if a == 0 { continue; }
        for (dx, dy) in [(1usize, 0usize), (0, 1)] {
            if x + dx > w - 1 || y + dy > h - 1 { continue; }
            if labels[i + dx + dy * w] == a { continue; }
            let (mut k, mut clean) = (1usize, !other.d[i]); let mut b = 0;
            while k <= REACH { let (px, py) = (x + dx * k, y + dy * k); if px >= w || py >= h { break; } let j = py * w + px; if labels[j] != 0 { b = labels[j]; if other.d[j] { clean = false; } break; } if other.d[j] || value.d[j] < CLASS_AT - 0.2 || resp.d[j] > strong { clean = false; } k += 1; }
            if b == 0 || b == a { continue; }
            let key = (a.min(b), a.max(b)); if clean { *weak.entry(key).or_default() += 1; } else { *hard.entry(key).or_default() += 1; }
        } } }
    const SMALL: usize = 960; // 0.6 mm2
    let mut parent: HashMap<u32, u32> = HashMap::new();
    fn root(p: &mut HashMap<u32, u32>, x: u32) -> u32 { let mut r = x; while let Some(q) = p.get(&r) { if *q == r { break; } r = *q; } r }
    let mut groups: HashMap<u32, usize> = area.clone();
    let mut cand: Vec<((u32, u32), u32)> = weak.iter().map(|(k, v)| (*k, *v)).collect();
    // the strongest clean boundaries first; ties settle by id so a run repeats itself
    cand.sort_by(|a, b| b.1.cmp(&a.1).then(a.0.cmp(&b.0)));
    let (mut joined, mut refused) = (0, 0);
    for ((a, b), n) in cand {
        let h_ = hard.get(&(a, b)).copied().unwrap_or(0);
        if n < 5 || h_ * 2 > n { continue; }
        let (ra, rb) = (root(&mut parent, a), root(&mut parent, b)); if ra == rb { continue; }
        let (sa, sb) = (groups[&ra], groups[&rb]);
        if sa.min(sb) >= SMALL { refused += 1; continue; }
        let (keep, gone) = if sa >= sb { (ra, rb) } else { (rb, ra) };
        parent.insert(gone, keep); *groups.get_mut(&keep).unwrap() = sa + sb; joined += 1;
    }
    for l in labels.iter_mut() { if *l != 0 { *l = root(&mut parent, *l); } }
    (joined, refused)
}

/// Zhang-Suen thinning, in place.
pub fn thin(m: &mut [bool], w: usize, h: usize) {
    let mut active: Vec<usize> = (0..w * h).filter(|i| m[*i] && i % w > 0 && i % w < w - 1 && i / w > 0 && i / w < h - 1).collect();
    loop {
        let mut changed = false;
        for step in 0..2 {
            let mut rem = vec![];
            for &i in &active {
                if !m[i] { continue; }
                let p = [m[i - w], m[i - w + 1], m[i + 1], m[i + w + 1], m[i + w], m[i + w - 1], m[i - 1], m[i - w - 1]];
                let b = p.iter().filter(|v| **v).count();
                if !(2..=6).contains(&b) { continue; }
                let a = (0..8).filter(|k| !p[*k] && p[(k + 1) % 8]).count();
                if a != 1 { continue; }
                let (p2, p4, p6, p8) = (p[0], p[2], p[4], p[6]);
                let c = if step == 0 { !(p2 && p4 && p6) && !(p4 && p6 && p8) } else { !(p2 && p4 && p8) && !(p2 && p6 && p8) };
                if c { rem.push(i); }
            }
            if !rem.is_empty() { changed = true; for i in rem { m[i] = false; } }
        }
        if !changed { break; }
        active.retain(|i| m[*i]);
    }
}

fn endpoints(m: &[bool], w: usize, h: usize) -> Vec<usize> {
    let mut v = vec![];
    for y in 1..h - 1 { for x in 1..w - 1 { let i = y * w + x; if !m[i] { continue; }
        let n = [i - w - 1, i - w, i - w + 1, i - 1, i + 1, i + w - 1, i + w, i + w + 1].iter().filter(|j| m[**j]).count();
        if n == 1 { v.push(i); } } }
    v
}

fn draw_line(m: &mut [bool], w: usize, a: (i64, i64), b: (i64, i64)) {
    let n = (b.0 - a.0).abs().max((b.1 - a.1).abs()).max(1);
    for t in 0..=n { let (x, y) = (a.0 + (b.0 - a.0) * t / n, a.1 + (b.1 - a.1) * t / n); m[y as usize * w + x as usize] = true; }
}

fn label4(open: &[bool], w: usize, h: usize) -> (Vec<u32>, u32) {
    let mut lab = vec![0u32; w * h]; let mut n = 0; let mut st = vec![];
    for s in 0..w * h { if !open[s] || lab[s] != 0 { continue; } n += 1; lab[s] = n; st.push(s);
        while let Some(i) = st.pop() { let (x, y) = (i % w, i / w);
            if x > 0 && open[i - 1] && lab[i - 1] == 0 { lab[i - 1] = n; st.push(i - 1); }
            if x + 1 < w && open[i + 1] && lab[i + 1] == 0 { lab[i + 1] = n; st.push(i + 1); }
            if y > 0 && open[i - w] && lab[i - w] == 0 { lab[i - w] = n; st.push(i - w); }
            if y + 1 < h && open[i + w] && lab[i + w] == 0 { lab[i + w] = n; st.push(i + w); } } }
    (lab, n)
}

fn median(v: &mut [f32]) -> f32 { if v.is_empty() { return 0.0; } v.sort_by(|a, b| a.partial_cmp(b).unwrap()); v[v.len() / 2] }

pub fn trace(img: &RgbImage, pr: &Params) -> Traced {
    let (w, h) = (img.width() as usize, img.height() as usize);
    let mean2: Option<RgbImage> = pr.second.filter(|s| s.dimensions() == img.dimensions()).map(|s| RgbImage::from_fn(w as u32, h as u32, |x, y| { let (a, b) = (img.get_pixel(x, y), s.get_pixel(x, y)); Rgb([((a[0] as u16 + b[0] as u16 + 1) / 2) as u8, ((a[1] as u16 + b[1] as u16 + 1) / 2) as u8, ((a[2] as u16 + b[2] as u16 + 1) / 2) as u8]) }));
    let shot_a = img; let shot_b = pr.second.filter(|_| mean2.is_some());
    let img = mean2.as_ref().unwrap_or(img);
    let debug = std::env::var("TRACE_DEBUG").is_ok();
    let (mut rp, mut gp, mut bp) = (Plane::new(w, h, 0.0), Plane::new(w, h, 0.0), Plane::new(w, h, 0.0));
    for y in 0..h { for x in 0..w { let p: &Rgb<u8> = img.get_pixel(x as u32, y as u32); let i = y * w + x; rp.d[i] = p[0] as f32; gp.d[i] = p[1] as f32; bp.d[i] = p[2] as f32; } }
    let (rs, gs, bs) = (blur(&rp, 1, 2), blur(&gp, 1, 2), blur(&bp, 1, 2));
    // the levels are read from a little smoothing: enough that grain alone does not look like two classes
    let sat: Vec<f32> = (0..w * h).map(|i| { let (a, b, c) = (rs.d[i], gs.d[i], bs.d[i]); let (mx, mn) = (a.max(b).max(c), a.min(b).min(c)); if mx > 1.0 { (mx - mn) / mx } else { 0.0 } }).collect();
    let val: Vec<f32> = (0..w * h).map(|i| rs.d[i].max(gs.d[i]).max(bs.d[i]) / 255.0).collect();
    // The mask's colour is measured, not assumed green (`maskcolour`). `warm` is "red over green" on a green mask and a warm hue off the mask's own on any other; from here on `gp` and `gs` are
    // the channel copper is read in, which is green itself on a green mask.
    let mc = crate::maskread::with_channel(img, crate::maskcolour::of(img), debug);
    if debug { eprintln!("mask hue {:.0}, median RGB {:?}, read in channel {}", mc.hue, mc.rgb, mc.measure); }
    let warm_at: Option<Vec<bool>> = if mc.green { None } else { Some((0..w * h).map(|i| mc.warm(rs.d[i], gs.d[i], bs.d[i], 0.0)).collect()) };
    let g_true = gs.clone();
    let warm = |i: usize, margin: f32| -> bool { match &warm_at { None => rs.d[i] > g_true.d[i] + margin, Some(v) => v[i] } };
    // a colour that cannot be the mask's is a part, not board (`maskread`): never set on a green mask
    let off_board: Option<Vec<bool>> = if mc.green { None } else { Some((0..w * h).map(|i| crate::maskread::not_board(&mc, rs.d[i], g_true.d[i], bs.d[i])).collect()) };
    let off = |i: usize| -> bool { match &off_board { None => false, Some(v) => v[i] } };
    let (gp, gs) = match mc.measure { 0 => (rp.clone(), rs.clone()), 2 => (bp.clone(), bs.clone()), _ => (gp, gs) };
    let lg_fit = { let b = blur(&gp, 3, 2); Plane { w, h, d: b.d.iter().map(|v| v.max(1.0).ln()).collect() } };
    // the fit sees what is plainly solder mask
    let wide = blur(&gp, 60, 2);
    let excl = Mask { w, h, d: (0..w * h).map(|i| warm(i, 5.0) || gs.d[i] < 18.0 || gs.d[i] > wide.d[i] * 1.5 || (pr.bodies && sat[i] < 0.3)).collect() }.dilate(4.0);
    // Wide bright things (silk blocks, pads, parts) light the mask beside them, and a window with them in reads its copper too high: a whole pour in a dim box then sits on the mid-point
    // (Driver bottom, PROG box: value 0.5 throughout). Thin glints are left in: they are everywhere along copper edges. 8 px is where the lift has fallen off on these photos, tile set A.
    let excl = excl.or(&Mask { w, h, d: (0..w * h).map(|i| gs.d[i] > wide.d[i] * 1.2).collect() }.open(2.5).dilate(8.0));
    let spread = light_spread(&lg_fit, &excl);
    let (cop, lam) = if spread < FLAT_LIGHT {
        // Under flat light one fit over the face is all the levels need, and it must see only plain mask: dark green shadow beside parts is not laminate, and it took the lower level with it (CM5 top: separation 0.42 where copper over laminate is 0.21).
        let parts = Mask { w, h, d: (0..w * h).map(|i| { let (s, v) = (sat[i], val[i]); (s < 0.27 && v > 0.30) || (warm(i, 10.0) && v > 0.35) || (v < 0.30 && s < 0.42) || off(i) }).collect() }.open(2.0).dilate(16.0);
        let plain = excl.or(&parts);
        let samples: Vec<f32> = (0..w * h).step_by(3).filter(|i| !plain.d[*i]).map(|i| lg_fit.d[i]).collect();
        // what is darker than the laminate mode by more than half the way to copper is a shadow or a pit, so it stays out of the fit: the mode, not the mean of a long dark tail, is the laminate's level
        match two_modes(&samples) {
            Some((lo, hi)) => { if debug { eprintln!("green modes: laminate {:.0}, copper {:.0}", lo.exp(), hi.exp()); } let limit = lo - 0.5 * (hi - lo); two_levels_at(&lg_fit, &Mask { w, h, d: (0..w * h).map(|i| plain.d[i] || lg_fit.d[i] < limit).collect() }, &[], 1) }
            None => two_levels_at(&lg_fit, &plain, &[], 1),
        }
    } else { two_levels(&lg_fit, &excl) };
    if debug { let at = |p: &Plane| { let mut v: Vec<f32> = p.d.iter().step_by(97).copied().collect(); v.sort_by(|a, b| a.partial_cmp(b).unwrap()); v[v.len() / 2].exp() }; eprintln!("light spread {spread:.3}; copper green {:.0}, laminate green {:.0}", at(&cop), at(&lam)); }
    // Print, from the material reader (`photomat`), which tells it from a glint by more than brightness. A printed line's blurred rim is brighter than laminate and not pale, so it read as copper
    // and made a track of every line of print on laminate (hub board, bottom: the frame of a printed table). The rim is not read, and print is never the lit edge of copper.
    let (print, exposed) = if std::env::var("TRACE_NO_PRINT").is_ok() || !crate::photomat::reads(&mc) { (Mask::new(w, h), Mask::new(w, h)) } else { let f = crate::photomat::classify(img, pr.bodies);
        (Mask { w, h, d: f.cls.iter().map(|c| *c == crate::photomat::SILK).collect() }, Mask { w, h, d: f.cls.iter().map(|c| *c == crate::photomat::GOLD || *c == crate::photomat::SILVER).collect() }) };
    let print_rim = print.dilate(PRINT_RIM);
    // What the levels make of the face. Run again when the photo's own step edges show the first levels did not come from the copper (`edgelevels`).
    let masks_of = |cop: &Plane, lam: &Plane| {
        // pale: silk, tin, domes and glints. They are not read, only filled in from what is around them. A copper level near the top of the range leaves no room for pale things to stand
        // 0.30 above it (CM5: copper's green reads 200, so white print is 0.24 up), so the step is at most half the room left.
        let pale = Mask { w, h, d: (0..w * h).map(|i| gs.d[i].max(1.0).ln() - cop.d[i] > (0.5 * (255f32.ln() - cop.d[i]).max(0.0)).min(0.30) && sat[i] < 0.38).collect() };
        let pale_wide = pale.open(3.2);
        let halo = pale_wide.dilate(3.0);
        let gold = Mask { w, h, d: (0..w * h).map(|i| warm(i, 10.0) && val[i] > 0.45).collect() }.open(2.0);
        // holes and dark solder, and on a face with parts, part bodies and metal
        let opaque = if pr.bodies {
            Mask { w, h, d: (0..w * h).map(|i| { let (s, v) = (sat[i], val[i]); (s < 0.27 && v > 0.30) || (warm(i, 10.0) && v > 0.35) || (v < 0.30 && s < 0.42) || off(i) }).collect() }.open(2.0).dilate(2.0)
        } else {
            Mask { w, h, d: (0..w * h).map(|i| gs.d[i].max(1.0).ln() < lam.d[i] - 0.65 || (warm(i, 0.0) && val[i] <= 0.45)).collect() }.open(4.0).dilate(1.0)
        };
        // silk on a face with parts stays transparent: copper runs under it
        let silk_top = Mask { w, h, d: (0..w * h).map(|i| sat[i] < 0.2 && val[i] > 0.55).collect() }.close(2.0).dilate(4.0);
        // ...but print is a thin thing: anything pale and neutral wider than 1.2 mm solid is a metal can, a lid or a plastic body, and copper does not run through a part
        let silk_top = silk_top.andnot(&silk_top.erode(PALE_BODY_HALF + 4.0).dilate(PALE_BODY_HALF + 8.0));
        let mut hidden = Mask { w, h, d: (0..w * h).map(|i| opaque.d[i] && !(pr.bodies && silk_top.d[i])).collect() };
        if let Some((x0, y0, x1, y1, r)) = pr.board {
            let r = r.max(8.0);
            // off the board and the photo's edge band: not copper we can say anything about
            for y in 0..h { for x in 0..w { let (fx, fy) = (x as f64, y as f64); let (cx, cy) = (fx.clamp(x0 + r, x1 - r), fy.clamp(y0 + r, y1 - r));
                let d = ((fx - cx).powi(2) + (fy - cy).powi(2)).sqrt(); if d > r - 8.0 { hidden.d[y * w + x] = true; } } }
        }
        // A sharpened photo digs a dark band round anything bright: it is not a gap, so it is not read, and left in it also sets the scale of every line threshold (`sharpening_halo`).
        let sharp_px = sharpening_halo(&gp, &hidden, &pale);
        if debug { eprintln!("sharpening halo {sharp_px} px"); }
        let sharp = if sharp_px > 0 { hidden.dilate(sharp_px as f32) } else { Mask::new(w, h) };
        // pale (and hidden) pixels are not looked at by the class or line filters
        let unseen = Mask { w, h, d: (0..w * h).map(|i| hidden.d[i] || pale.d[i] || halo.d[i] || sharp.d[i] || print_rim.d[i] || (!pr.bodies && gold.d[i])).collect() };
        (pale, pale_wide, halo, gold, hidden, sharp, unseen)
    };
    let (mut cop, mut lam) = (cop, lam);
    let (mut pale, mut pale_wide, mut halo, mut gold, mut hidden, mut sharp, mut unseen) = masks_of(&cop, &lam);
    let mut edge_check = None; let mut levels_from_edges = false; let mut levels_local = 0f32;
    // Levels from the edges: under a light that moves every place gives a sample by the side of its nearest edge (`locallevels::levels_by_side`); under flat light the edges' own
    // two sides averaged over a window are enough, and stay as they were when the CM5 top was read.
    let from_edges = |e: &crate::edgelevels::StepEdges, unseen: &Mask, solid: &Mask| -> (Plane, Plane) {
        let mean = crate::edgelevels::levels(e);
        if spread < FLAT_LIGHT || std::env::var("TRACE_NO_SIDE_LEVELS").is_ok() { return mean; }
        // The side of an edge is only worth asking when the face's edges are copper edges. Where even their own window means leave many of them wholly in one class, the band holds print rims,
        // part edges and fibres as well, and a place's nearest edge is as likely one of those (`locallevels::SIDES_NEED`).
        let a = crate::edgelevels::agreement(e, &mean.0, &mean.1, CLASS_AT);
        if debug { eprintln!("edges' own window means: straddle {:.2}; levels by side {}", a.straddle, a.straddle >= crate::locallevels::SIDES_NEED); }
        if a.straddle >= crate::locallevels::SIDES_NEED { crate::locallevels::levels_by_side(e, &lg_fit, unseen, solid) } else { mean } };
    if let Some(e) = crate::edgelevels::step_edges(&gp, &unseen) {
        let a = crate::edgelevels::agreement(&e, &cop, &lam, CLASS_AT);
        if debug { eprintln!("levels against {} step edges (modal step {:.3}): straddle {:.2}, both copper {:.2}, both laminate {:.2}", a.edges, a.modal_step, a.straddle, a.both_copper, a.both_laminate); }
        edge_check = Some(a);
        if let Ok(dir) = std::env::var("TRACE_DUMP") { let mut t = String::from("x,y,hi,lo\n"); for (x, y, hi, lo) in &e.at { t.push_str(&format!("{x},{y},{hi:.3},{lo:.3}\n")); } let _ = std::fs::write(format!("{dir}/{}-edges.csv", if pr.bodies { "top" } else { "bottom" }), t); }
        let forced = std::env::var("TRACE_EDGE_SHARE").ok().and_then(|v| v.parse::<f32>().ok()).map(|v| a.edges >= 500 && a.both_copper.max(a.both_laminate) > v);
        // on a mask that is not green the levels are also checked against the size of the step itself (`maskread::separation_is_off`)
        let off_step = !mc.green && crate::maskread::separation_is_off(&e, &cop, &lam, debug);
        if (forced.unwrap_or(a.levels_are_wrong()) || off_step) && std::env::var("TRACE_NO_EDGE_LEVELS").is_err() {
            let (c2, l2) = from_edges(&e, &unseen, &hidden.or(&pale_wide).or(&gold)); cop = c2; lam = l2; levels_from_edges = true;
            (pale, pale_wide, halo, gold, hidden, sharp, unseen) = masks_of(&cop, &lam);
            if let Some(e2) = crate::edgelevels::step_edges(&gp, &unseen) { let a2 = crate::edgelevels::agreement(&e2, &cop, &lam, CLASS_AT); if debug { eprintln!("levels from the edges: straddle {:.2}, both copper {:.2}, both laminate {:.2}", a2.straddle, a2.both_copper, a2.both_laminate); } edge_check = Some(a2); }
        }
        // a face that passes as a whole can still be wrong over part of itself: each place is tested against the edges near it (`locallevels::correct`)
        else if std::env::var("TRACE_NO_LOCAL_LEVELS").is_err() && std::env::var("TRACE_NO_EDGE_LEVELS").is_err() { let solid = hidden.or(&pale_wide).or(&gold); if let Some(lc) = crate::locallevels::correct(&e, &cop, &lam, CLASS_AT, &|| from_edges(&e, &unseen, &solid)) {
            cop = lc.cop; lam = lc.lam; levels_local = lc.share;
            (pale, pale_wide, halo, gold, hidden, sharp, unseen) = masks_of(&cop, &lam);
            if let Some(e2) = crate::edgelevels::step_edges(&gp, &unseen) { let a2 = crate::edgelevels::agreement(&e2, &cop, &lam, CLASS_AT); if debug { eprintln!("levels from the edges over {:.0} percent of the face: straddle {:.2}, both copper {:.2}, both laminate {:.2}", lc.share * 100.0, a2.straddle, a2.both_copper, a2.both_laminate); } edge_check = Some(a2); }
        } }
        if std::env::var("TRACE_LEVELS_ONLY").is_ok() { std::process::exit(0); }
    }
    // brightness clipped between the two levels and a little beyond: a glint or a dome must not outweigh an edge. One plane per shot.
    let clip_of = |im: &RgbImage| -> Plane { let g = Plane { w, h, d: im.pixels().map(|p| p[mc.measure] as f32).collect() }; let b = blur(&g, 1, 1);
        Plane { w, h, d: (0..w * h).map(|i| { let (c, l) = (cop.d[i], lam.d[i]); let d = (c - l).max(0.05); b.d[i].max(1.0).ln().clamp(l - 0.5 * d, c + 0.25 * d) }).collect() } };
    let clip_a = clip_of(shot_a); let clip_b = shot_b.map(&clip_of);
    let sm = along_least_variance(&clip_a, clip_b.as_ref().unwrap_or(&clip_a), &unseen, 8);
    let mut value = Plane::new(w, h, 0.0);
    for i in 0..w * h { value.d[i] = (sm.d[i] - lam.d[i]) / (cop.d[i] - lam.d[i]); }
    // gap lines, on the clipped brightness; the filter is linear, so the mean of the two responses is the response to the mean of the shots
    let filter_on = |c: &Plane| -> Plane { let mut n = Plane::new(w, h, 80.0); for i in 0..w * h { if !unseen.d[i] { n.d[i] = (c.d[i] - cop.d[i]).exp() * 80.0; } } n };
    // a gap line's pixels are dimmed by it, so the filter need not run on pixels at or above the copper level
    let line_skip = Mask { w, h, d: (0..w * h).map(|i| hidden.d[i] || halo.d[i] || sharp.d[i] || value.d[i] > 1.0).collect() };
    let (resp, ang, resp_a, resp_b) = match &clip_b {
        Some(cb) => { let ((ra, aa), (rb, _)) = (line_filter(&filter_on(&clip_a), &line_skip), line_filter(&filter_on(cb), &line_skip));
            // the direction goes with the stronger shot; the response is the mean over the two at that pixel
            let m = Plane { w, h, d: (0..w * h).map(|i| (ra.d[i] + rb.d[i]) / 2.0).collect() }; (m, aa, Some(ra), Some(rb)) }
        None => { let (r, a) = line_filter(&filter_on(&clip_a), &line_skip); (r, a, None, None) }
    };
    // The scale of the line thresholds is the response on plain board, 0.4 mm clear of parts and print. Taken over the whole face it was set by whatever clutter the face carries, not by its
    // gaps: the weave of a strip of fabric tape on test board 1's bottom put it at 14.0 where the board itself reads 4.1, and chip markings on test board 4's top at 14.6 against 9.5, so
    // gap lines the eye follows round every pour fell under the bar. `TRACE_P99_ALL` gives the old scale.
    let p99 = { let pct = |v: &mut Vec<f32>| -> Option<f32> { v.sort_by(|a, b| a.partial_cmp(b).unwrap()); v.get(v.len() * 99 / 100).copied() };
        let mut all: Vec<f32> = resp.d.iter().step_by(5).copied().filter(|x| *x > 0.0).collect();
        let clear = hidden.or(&pale_wide).dilate(16.0);
        let mut plain: Vec<f32> = (0..w * h).step_by(5).filter(|i| !clear.d[*i] && resp.d[*i] > 0.0).map(|i| resp.d[i]).collect();
        if plain.len() >= 2000 && std::env::var("TRACE_P99_ALL").is_err() { pct(&mut plain).unwrap_or(1.0) } else { pct(&mut all).unwrap_or(1.0) } };
    if debug {
        let clear = hidden.or(&pale_wide).dilate(16.0);
        let at = |v: &mut Vec<f32>, q: f64| -> f32 { if v.is_empty() { return 0.0; } v.sort_by(|a, b| a.partial_cmp(b).unwrap()); v[((v.len() - 1) as f64 * q) as usize] };
        let mut flat: Vec<f32> = (0..w * h).filter(|i| !clear.d[*i] && value.d[*i] > 0.8 && value.d[*i] < 1.2).map(|i| resp.d[i]).collect();
        let mut plain: Vec<f32> = (0..w * h).filter(|i| !clear.d[*i]).map(|i| resp.d[i]).collect();
        eprintln!("resp p99 {p99}; flat copper n {} p50 {:.2} p90 {:.2} p99 {:.2}; clear of parts n {} p50 {:.2} p90 {:.2} p99 {:.2}; sep {:.3}", flat.len(), at(&mut flat, 0.5), at(&mut flat, 0.9), at(&mut flat, 0.99), plain.len(), at(&mut plain, 0.5), at(&mut plain, 0.9), at(&mut plain, 0.99), { let mut d: Vec<f32> = (0..w * h).step_by(97).map(|i| cop.d[i] - lam.d[i]).collect(); at(&mut d, 0.5) });
    }
    let (hi, lo) = (pr.hi * p99, pr.lo * p99);
    let mut ridge = vec![false; w * h];
    for y in 3..h - 3 { for x in 3..w - 3 { let i = y * w + x; let v = resp.d[i]; if v <= lo { continue; }
        let a = ang[i] as f64 * std::f64::consts::PI / 16.0; let (nx, ny) = (-a.sin() * 1.5, a.cos() * 1.5);
        if v >= resp.sample((x as f64 + nx) as f32, (y as f64 + ny) as f32) && v >= resp.sample((x as f64 - nx) as f32, (y as f64 - ny) as f32) { ridge[i] = true; } } }
    let rm = Mask { w, h, d: ridge };
    let (lab, comps) = crate::raster::components(&rm);
    let mut mx = vec![0f32; comps.len() + 1];
    for i in 0..w * h { if lab[i] != 0 { let l = lab[i] as usize; mx[l] = mx[l].max(resp.d[i]); } }
    let (mut vs, mut vc) = (vec![0f32; comps.len() + 1], vec![0f32; comps.len() + 1]);
    for i in 0..w * h { if lab[i] != 0 && !unseen.d[i] { vs[lab[i] as usize] += value.d[i]; vc[lab[i] as usize] += 1.0; } }
    let mut keep = vec![false; comps.len() + 1];
    // a line must be there in both shots: grain lines one up in one shot, not in the same place in the other
    let (mut sa, mut sb) = (vec![0f32; comps.len() + 1], vec![0f32; comps.len() + 1]);
    if let (Some(ra), Some(rb)) = (&resp_a, &resp_b) { for i in 0..w * h { if lab[i] != 0 { sa[lab[i] as usize] += ra.d[i]; sb[lab[i] as usize] += rb.d[i]; } } }
    let both = |l: usize, area: usize| -> bool { resp_a.is_none() || { let (a, b) = (sa[l] / area.max(1) as f32, sb[l] / area.max(1) as f32); a.min(b) >= BOTH_SHOTS * lo } };
    for c in &comps { let l = c.label as usize; keep[l] = mx[l] > hi && c.area >= pr.minlen && (vc[l] < 1.0 || vs[l] / vc[l] < LINE_VALUE) && both(l, c.area); }
    // a dome's shadow lies within 0.5 mm of it (Driver bottom: 8 to 20 px on the side away from the light)
    let mut near_dome = Mask::new(w, h); for (dx, dy) in pr.domes { let (u, v) = (*dx as i64, *dy as i64); for b in -20i64..=20 { for a in -20i64..=20 { if a * a + b * b <= 400 && u + a >= 0 && v + b >= 0 && u + a < w as i64 && v + b < h as i64 { near_dome.d[(v + b) as usize * w + (u + a) as usize] = true; } } } }
    // a line that never leaves the neighbourhood of a dome is the edge of its shadow: a gap runs on past
    { let mut out = vec![false; comps.len() + 1]; for i in 0..w * h { if lab[i] != 0 && !near_dome.d[i] { out[lab[i] as usize] = true; } } for c in &comps { if !out[c.label as usize] { keep[c.label as usize] = false; } } }
    let mut lines = Mask { w, h, d: (0..w * h).map(|i| lab[i] != 0 && keep[lab[i] as usize]).collect() }.close(2.0);
    // Under one oblique light a raised edge facing it glints, and the camera's sharpening digs a dark band 2 to 3 px wide right behind the glint, on the copper. It looks like a gap but the flat of the
    // copper comes back straight after. Driver bottom, tracks across their length: glint 5 px, band 2 px with a floor near 40, flat 3 px at the pour's level, then the real gap, 4 px with a floor near 15. Read as a gap,
    // the band cut every track in two along its length. A band is taken for copper only if brightness returns to the copper side without another glint: a narrow track with a gap and then the next glint stays apart.
    let mut undershoot = Mask::new(w, h);
    if let Some((lx, ly)) = pr.light {
        let glint = |x: i64, y: i64| x >= 0 && y >= 0 && x < w as i64 && y < h as i64 && { let j = y as usize * w + x as usize; pale.d[j] && !pale_wide.d[j] };
        for y in 8..h as i64 - 8 { for x in 8..w as i64 - 8 { let i = y as usize * w + x as usize; if hidden.d[i] || pale.d[i] || value.d[i] >= CLASS_AT { continue; }
            let behind = (1..=3).any(|k| glint(x + (lx * k as f64).round() as i64, y + (ly * k as f64).round() as i64)); if !behind { continue; }
            let back = (1..=5).find_map(|k| { let (u, v) = (x - (lx * k as f64).round() as i64, y - (ly * k as f64).round() as i64); let j = v as usize * w + u as usize; if pale.d[j] { Some(false) } else if value.d[j] >= CLASS_AT { Some(true) } else { None } }).unwrap_or(false);
            if back { undershoot.d[i] = true; } } }
        let u = undershoot.dilate(1.0); for i in 0..w * h { if u.d[i] { lines.d[i] = false; } }
    }
    // the class of every pixel: a seen pixel by its value, a pale one by its nearest seen neighbours
    let mut seed = vec![0i8; w * h];
    for i in 0..w * h { if hidden.d[i] { seed[i] = -1; } else if !pr.bodies && gold.d[i] { seed[i] = 1; } else if !pale.d[i] && !halo.d[i] && !sharp.d[i] && !print_rim.d[i] { seed[i] = if value.d[i] > CLASS_AT || undershoot.d[i] { 1 } else { -1 }; } }
    // grain: laminate specks with no gap line through them are not laminate
    { let dark = Mask { w, h, d: seed.iter().zip(&hidden.d).map(|(s, o)| *s == -1 && !*o).collect() };
      let (dl, dc) = crate::raster::components(&dark);
      let mut has_line = vec![false; dc.len() + 1]; for i in 0..w * h { if dl[i] != 0 && lines.d[i] { has_line[dl[i] as usize] = true; } }
      let beside_pale = pale.dilate(2.0).or(&near_dome);
      // wholly within reach of a dome and enclosing nothing: its shadow, whatever the line filter made of its edge
      let mut leaves = vec![false; dc.len() + 1]; for i in 0..w * h { if dl[i] != 0 && !near_dome.d[i] { leaves[dl[i] as usize] = true; } }
      let mut pale_touch = vec![false; dc.len() + 1]; for i in 0..w * h { if dl[i] != 0 && beside_pale.d[i] { pale_touch[dl[i] as usize] = true; } }
      // a clearance around a via encloses the via; a dome's shadow lies to one side of it
      let encloses = |c: &crate::raster::Component| -> bool { let (sw, sh) = (c.x1 - c.x0 + 3, c.y1 - c.y0 + 3); let mut m = Mask::new(sw, sh); for y in c.y0..=c.y1 { for x in c.x0..=c.x1 { if dl[y * w + x] == c.label { m.set(x - c.x0 + 1, y - c.y0 + 1, true); } } } m.fill_holes().count() > c.area + 20 };
      let small: Vec<bool> = { let mut s = vec![false; dc.len() + 1]; for c in &dc { let l = c.label as usize; s[l] = (c.area < 60 && !has_line[l]) || (c.area < SHADOW_PX && pale_touch[l] && !has_line[l] && !encloses(c)) || (c.area < SHADOW_PX && !leaves[l] && !encloses(c)); } s };
      for i in 0..w * h { if dl[i] != 0 && small[dl[i] as usize] { seed[i] = 1; } } }
    settle(&mut seed, &lg_fit, &cop, &lam, &lines, &hidden, &unseen, spread < FLAT_LIGHT, w, h);
    // a wide pale piece (a legend, a connector) takes the class around it as a whole, when that is clear; a mixed one goes by nearest
    { let (pl, pc) = crate::raster::components(&pale_wide);
      let ring = pale_wide.dilate(6.0).andnot(&halo);
      let (mut cu, mut la) = (vec![0u32; pc.len() + 1], vec![0u32; pc.len() + 1]);
      // the nearest piece to a ring pixel: walk a few pixels each way
      for y in 8..h.saturating_sub(8) { for x in 8..w.saturating_sub(8) { let i = y * w + x; if !ring.d[i] || seed[i] == 0 { continue; }
          for (dx, dy) in [(1i64, 0i64), (-1, 0), (0, 1), (0, -1)] { let mut k = 1i64; while k <= 8 { let j = ((y as i64 + dy * k) * w as i64 + x as i64 + dx * k) as usize; if pl[j] != 0 { if seed[i] == 1 { cu[pl[j] as usize] += 1; } else { la[pl[j] as usize] += 1; } break; } k += 1; } } } }
      for i in 0..w * h { let l = pl[i] as usize; if l == 0 || hidden.d[i] { continue; } let (a, b) = (cu[l] as f32, la[l] as f32); if a + b < 8.0 { continue; } if a >= 0.75 * (a + b) { seed[i] = 1; } else if b >= 0.75 * (a + b) { seed[i] = -1; } }
      // its halo goes with it
      for y in 3..h.saturating_sub(3) { for x in 3..w.saturating_sub(3) { let i = y * w + x; if !halo.d[i] || pale_wide.d[i] || seed[i] != 0 || hidden.d[i] { continue; } 'o: for r in 1i64..=3 { for (dx, dy) in [(1i64, 0i64), (-1, 0), (0, 1), (0, -1)] { let j = ((y as i64 + dy * r) * w as i64 + x as i64 + dx * r) as usize; if pl[j] != 0 && seed[j] != 0 { seed[i] = seed[j]; break 'o; } } } } } }
    let filled = nearest_class(&seed, w, h);
    let copper_seen = Mask { w, h, d: seed.iter().map(|s| *s == 1).collect() };
    let near_copper = copper_seen.dilate(5.0);
    let mut copper = Mask { w, h, d: filled.iter().map(|c| *c == 1).collect() };
    // a thin pale line beside copper is its lit edge
    for i in 0..w * h { if pale.d[i] && !pale_wide.d[i] && near_copper.d[i] && !print_rim.d[i] { copper.d[i] = true; } if hidden.d[i] { copper.d[i] = false; } }
    // a gap that runs into print runs on under it (`printgaps`)
    if std::env::var("TRACE_NO_PRINT_GAPS").is_err() { let g = crate::printgaps::gaps_under(&copper, &print_rim, &hidden); for i in 0..w * h { if g.d[i] { copper.d[i] = false; } } }
    let dark = Mask { w, h, d: (0..w * h).map(|i| !copper.d[i] && !hidden.d[i] && !pale.d[i]).collect() };
    // the barrier: what is not copper, and gap lines, with their loose ends walked on
    let mut skel = lines.d.clone();
    thin(&mut skel, w, h);
    // a gap line is as wide as its dip at half depth, not as wide as the filter's response: the filter's width ate 1 to 1.5 px off each side of the copper,
    // measured against the steepest fall of brightness on the mean of the shots. The centre line stays, so what a line separates stays separate.
    let mut gap_line = skel.clone();
    for y in 3..h - 3 { for x in 3..w - 3 { let i = y * w + x; if !skel[i] { continue; } let half = (value.d[i].clamp(0.0, 1.0) + 1.0) / 2.0;
        for dy in -3i64..=3 { for dx in -3i64..=3 { let j = ((y as i64 + dy) * w as i64 + x as i64 + dx) as usize; if lines.d[j] && value.d[j] < half { gap_line[j] = true; } } } } }
    let mut barrier = Mask { w, h, d: (0..w * h).map(|i| !copper.d[i] || gap_line[i]).collect() };
    let before = skel.clone();
    let mut with_edges = skel.clone(); for i in 0..w * h { if barrier.d[i] && !lines.d[i] { let (x, y) = (i % w, i / w); if x > 0 && x + 1 < w && y > 0 && y + 1 < h && (copper.d[i - 1] || copper.d[i + 1] || copper.d[i - w] || copper.d[i + w]) { with_edges[i] = true; } } }
    extend_from(&mut with_edges, &before, w, h, pr.ext);
    for i in 0..w * h { if with_edges[i] && !barrier.d[i] { barrier.d[i] = true; } }
    for x in 0..w { barrier.d[x] = true; barrier.d[(h - 1) * w + x] = true; }
    for y in 0..h { barrier.d[y * w] = true; barrier.d[y * w + w - 1] = true; }
    // exposed metal is its own island: whether a pad joins the pour around it is for the ring tests to say, not for two colours touching
    let metal = |i: usize| !pr.bodies && gold.d[i];
    let open: Vec<bool> = (0..w * h).map(|i| !barrier.d[i] && !metal(i)).collect();
    let (mut labels, n_l) = { let (a, na) = label4(&open, w, h); let open_m: Vec<bool> = (0..w * h).map(|i| !barrier.d[i] && metal(i)).collect(); let (b, nb) = label4(&open_m, w, h); let mut l = a; for i in 0..w * h { if b[i] != 0 { l[i] = b[i] + na; } } (l, na + nb) };
    let first = labels.clone();
    let mut steps = Steps::default();
    // crumbs are not copper we can say anything about
    let mut area = vec![0usize; n_l as usize + 1];
    for l in &labels { area[*l as usize] += 1; }
    let mut crumb = 0; { let mut seen = vec![false; n_l as usize + 1]; for l in labels.iter_mut() { if *l != 0 && area[*l as usize] < pr.min_area { if !seen[*l as usize] { seen[*l as usize] = true; crumb += 1; } *l = 0; } } }
    steps.crumbs = crumb; steps.edge_check = edge_check; steps.levels_from_edges = levels_from_edges; steps.levels_local = levels_local; steps.copper_share = copper.count() as f32 / (w * h - hidden.count()).max(1) as f32;
    steps.first = first;
    let clock = std::time::Instant::now();
    let (nk, pk) = split_necks(&mut labels, w, h); steps.necks_cut = nk; steps.pinches_cut = pk;
    let no_join = Mask { w, h, d: (0..w * h).map(|i| hidden.d[i] || metal(i) || print_rim.d[i]).collect() };
    let (jn, rf) = join_fragments(&mut labels, &value, &resp, hi * 0.5, &no_join, w, h); steps.joined = jn; steps.refused_joins = rf;
    // a track that stops and goes on in line within 0.5 mm, with no gap line between, is one track (`trackbridge`)
    if std::env::var("TRACE_NO_BRIDGE").is_err() { let b = crate::trackbridge::bridge(&mut labels, &lines, &value, &unseen, w, h); for i in 0..w * h { if b.painted.d[i] { barrier.d[i] = false; } } steps.bridged = b.count;
        if let Ok(dir) = std::env::var("TRACE_DUMP") { let mut t = String::from("x0,y0,x1,y1,half\n"); for (a, c, hf) in &b.segments { t.push_str(&format!("{},{},{},{},{}\n", a.0, a.1, c.0, c.1, hf)); } let _ = std::fs::write(format!("{dir}/{}-bridges.csv", if pr.bodies { "top" } else { "bottom" }), t); }
        if debug { eprintln!("track pieces bridged {}", b.count); } }
    // A pad, its solder and the copper leaving it are one conductor (`padjoin`): from the footprints' pads once parts are placed, else from the exposed metal on a face with parts.
    if std::env::var("TRACE_NO_PADJOIN").is_err() && (!pr.pads.is_empty() || pr.bodies) {
        let seeds = if !pr.pads.is_empty() { crate::padjoin::seed_map(pr.pads, w, h) } else { crate::padjoin::metal_seeds(&exposed) };
        // What stops a pad's growth is laminate the photo shows, not shadow: beside a tall part both levels fall together and everything reads below the laminate level.
        // Darker than laminate by `PAD_SHADOW` of the separation is shadow, and says nothing (hub board top: 5 percent of pixels further than 1.5 mm from any part read that dark, 19 percent of those nearer).
        let blocked = Mask { w, h, d: (0..w * h).map(|i| (seed[i] == -1 && !unseen.d[i] && value.d[i] >= -PAD_SHADOW) || lines.d[i] || gap_line[i]).collect() };
        let j = crate::padjoin::join(&mut labels, &seeds, &blocked, &hidden, &if pr.pads.is_empty() { vec![] } else { crate::padjoin::mates(pr.pads) }, w, h);
        for i in 0..w * h { if j.painted.d[i] { barrier.d[i] = false; } }
        steps.pads_joined = j.pads_joined;
        if debug { eprintln!("pads joined to copper {}, islands merged through a pad {}", j.pads_joined, j.islands_merged); }
    }
    // the local brightness over the copper level, for the ring tests around pads
    let (fine, fine_w) = { let mut a = gp.clone(); let mut c = Plane::new(w, h, 1.0); for i in 0..w * h { if unseen.d[i] { a.d[i] = 0.0; c.d[i] = 0.0; } } (blur(&a, 3, 1), blur(&c, 3, 1)) };
    let mut ratio = Plane::new(w, h, 1.0);
    for i in 0..w * h { if !unseen.d[i] && fine_w.d[i] > 0.02 { ratio.d[i] = (fine.d[i] / fine_w.d[i]) / cop.d[i].exp().max(1.0); } }
    if debug { eprintln!("necks and joins {:?}", clock.elapsed()); eprintln!("p99 {p99}, regions {n_l}, crumbs {crumb}, necks {nk}, pinches {pk}, joined {jn}, refused {rf}, copper {:.3}, pale {:.3}, hidden {:.3}", copper.count() as f64 / (w * h) as f64, pale.count() as f64 / (w * h) as f64, hidden.count() as f64 / (w * h) as f64); }
    if let Ok(dir) = std::env::var("TRACE_DUMP") { let side = if pr.bodies { "top" } else { "bottom" };
        let save = |m: &Mask, n: &str| { let im = 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 }])); let _ = im.save(format!("{dir}/{side}-{n}.png")); };
        save(&copper, "copper"); save(&pale, "pale"); { let im = image::GrayImage::from_fn(w as u32, h as u32, |x, y| image::Luma([(value.at(x as usize, y as usize).clamp(-0.5, 1.5) * 127.0 + 63.0) as u8])); let _ = im.save(format!("{dir}/{side}-value.png")); } save(&hidden, "hidden"); save(&lines, "lines"); save(&barrier, "barrier");
        for (pl, n) in [(&cop, "level-copper"), (&lam, "level-laminate")] { let im = image::GrayImage::from_fn(w as u32, h as u32, |x, y| image::Luma([pl.at(x as usize, y as usize).exp().clamp(0.0, 255.0) as u8])); let _ = im.save(format!("{dir}/{side}-{n}.png")); } }
    let acc = measure_regions(&mut labels, &ratio, &hidden, w, h);
    Traced { w, h, resp, gap: barrier, labels, regions: acc, other: hidden, ratio, dark, silk: pale_wide, pale, steps }
}

/// Walk on from each end of the gap lines until another barrier is met.
fn extend_from(all: &mut Vec<bool>, lines: &[bool], w: usize, h: usize, ext: i64) {
    let snap = all.clone();
    for e in endpoints(lines, w, h) {
        let (x, y) = ((e % w) as i64, (e / w) as i64);
        let (mut sx, mut sy, mut c) = (0f64, 0f64, 0f64);
        for dy in -6i64..=6 { for dx in -6i64..=6 { let (px, py) = (x + dx, y + dy); if px >= 0 && py >= 0 && px < w as i64 && py < h as i64 && lines[py as usize * w + px as usize] { sx += px as f64; sy += py as f64; c += 1.0; } } }
        if c < 2.0 { continue; }
        let base = (y as f64 - sy / c).atan2(x as f64 - sx / c);
        'walk: for t in 3..ext { for da in [0.0, 0.17, -0.17, 0.35, -0.35, 0.6, -0.6] {
            let a = base + da;
            let (px, py) = (x + (a.cos() * t as f64).round() as i64, y + (a.sin() * t as f64).round() as i64);
            if px < 0 || py < 0 || px >= w as i64 || py >= h as i64 { continue; }
            if snap[py as usize * w + px as usize] && (px - x).abs() + (py - y).abs() > 3 { draw_line(all, w, (x, y), (px, py)); break 'walk; }
        } }
    }
}

/// Renumber the labels 1..n and measure each island.
pub fn measure_regions(labels: &mut [u32], ratio: &Plane, other: &Mask, w: usize, h: usize) -> Vec<Region> {
    let mut remap: HashMap<u32, u32> = HashMap::new();
    for l in labels.iter_mut() { if *l != 0 { let k = remap.len() as u32 + 1; *l = *remap.entry(*l).or_insert(k); } }
    let count = remap.len();
    let inv = Mask { w, h, d: labels.iter().map(|l| *l == 0).collect() };
    let dt = inv.distance_to_set();
    let mut acc: Vec<Region> = (0..count).map(|i| Region { id: i as u32 + 1, area: 0, cx: 0.0, cy: 0.0, x0: usize::MAX, y0: usize::MAX, x1: 0, y1: 0, level: 0.0, half_width: 0.0, other_share: 0.0 }).collect();
    let mut levels: Vec<Vec<f32>> = vec![vec![]; count]; let mut oth = vec![0usize; count];
    for y in 0..h { for x in 0..w { let i = y * w + x; let l = labels[i]; if l == 0 { continue; } let r = &mut acc[l as usize - 1];
        r.area += 1; r.cx += x as f64; r.cy += y as f64; r.x0 = r.x0.min(x); r.y0 = r.y0.min(y); r.x1 = r.x1.max(x); r.y1 = r.y1.max(y); r.half_width = r.half_width.max(dt.d[i]);
        if other.d[i] { oth[l as usize - 1] += 1; } else if (x + y) % 3 == 0 { levels[l as usize - 1].push(ratio.d[i]); } } }
    for (i, r) in acc.iter_mut().enumerate() { if r.area > 0 { r.cx /= r.area as f64; r.cy /= r.area as f64; } r.level = median(&mut levels[i]); r.other_share = oth[i] as f32 / r.area.max(1) as f32; }
    acc
}

/// An overlay: every region its own colour, boundaries dark.
pub fn overlay(img: &RgbImage, t: &Traced) -> RgbImage {
    let mut o = img.clone();
    let col = |l: u32| -> [u8; 3] { let x = l.wrapping_mul(2654435761); [90 + (x >> 8 & 0x9f) as u8, 90 + (x >> 16 & 0x9f) as u8, 90 + (x >> 24 & 0x9f) as u8] };
    for y in 0..t.h { for x in 0..t.w { let i = y * t.w + x; let p = o.get_pixel_mut(x as u32, y as u32);
        if t.gap.d[i] { *p = Rgb([0, 0, 0]); } else if t.labels[i] != 0 { let c = col(t.labels[i]); for k in 0..3 { p[k] = ((p[k] as u32 * 2 + c[k] as u32 * 3) / 5) as u8; } } } }
    o
}

pub fn regions_json(t: &Traced, to_mm: &dyn Fn(f64, f64) -> (f64, f64), px_per_mm: f64) -> Value {
    let v: Vec<Value> = t.regions.iter().map(|r| { let c = to_mm(r.cx, r.cy);
        json!({"id": r.id, "areaMm2": ((r.area as f64 / (px_per_mm * px_per_mm)) * 100.0).round() / 100.0, "centreMm": [(c.0 * 100.0).round() / 100.0, (c.1 * 100.0).round() / 100.0], "halfWidthMm": ((r.half_width as f64 + 1.5) / px_per_mm * 100.0).round() / 100.0, "level": r.level, "otherShare": (r.other_share * 100.0).round() / 100.0,
               "boxPx": [r.x0, r.y0, r.x1, r.y1]}) }).collect();
    json!({"regions": v})
}

#[cfg(test)]
mod tests {
    use super::*;
    fn board(w: u32, h: u32, paint: &dyn Fn(u32, u32) -> u8) -> RgbImage {
        let mut im = RgbImage::new(w, h);
        for y in 0..h { for x in 0..w { let n = ((x.wrapping_mul(7919) ^ y.wrapping_mul(104729)) % 7) as u8; let g = paint(x, y).saturating_add(n); im.put_pixel(x, y, Rgb([g / 6, g, g * 2 / 3])); } }
        im
    }
    #[test]
    fn two_islands_split_by_a_dark_lane_and_grain_does_not_split_them() {
        let im = board(400, 240, &|x, _| if (190..210).contains(&x) { 45 } else { 95 });
        let t = trace(&im, &Params::default());
        let (a, b) = (t.labels[120 * 400 + 80], t.labels[120 * 400 + 320]);
        assert!(a != 0 && b != 0 && a != b, "{a} {b}");
        assert_eq!(t.labels[30 * 400 + 80], a, "grain must not split the left island");
    }
    #[test]
    fn a_track_between_two_lanes_is_its_own_narrow_island() {
        let im = board(400, 240, &|_, y| if (100..112).contains(&y) || (124..136).contains(&y) { 45 } else { 95 });
        let t = trace(&im, &Params::default());
        let mid = t.labels[118 * 400 + 200];
        assert!(mid != 0);
        assert_ne!(t.labels[50 * 400 + 200], mid);
        assert_ne!(t.labels[200 * 400 + 200], mid);
        let r = t.regions.iter().find(|r| r.id == mid).unwrap();
        assert!(r.half_width < 12.0, "half width {}", r.half_width);
    }
    #[test]
    fn grain_in_one_shot_draws_no_line_when_the_other_shot_lacks_it() {
        // a flat pour; one shot has a dark scratch of grain in it, the other not. Read from the first alone the scratch is a gap; read from both it is not
        let flat = board(400, 240, &|_, _| 95);
        let mut scratched = flat.clone();
        for x in 80..320 { for y in 118..122 { scratched.put_pixel(x, y, Rgb([8, 50, 33])); } }
        let lanes = |im: &RgbImage| { let mut o = im.clone(); for y in 0..240 { for x in 0..30 { o.put_pixel(x, y, Rgb([8, 45, 30])); } } o };
        let (a, b) = (lanes(&scratched), lanes(&flat));
        let both = trace(&a, &Params { second: Some(&b), ..Default::default() });
        assert_eq!(both.labels[100 * 400 + 200], both.labels[140 * 400 + 200], "the two sides of a scratch one shot lacks are one piece");
        assert!(both.labels[100 * 400 + 200] != 0);
    }
    #[test]
    fn two_modes_finds_laminate_and_copper_past_a_dark_tail() {
        // copper at ln 205, laminate at ln 175, and a long tail of shadow below that
        let mut v = vec![];
        for k in 0..6000 { v.push((205.0f32).ln() + 0.012 * ((k % 21) as f32 - 10.0) / 10.0); }
        for k in 0..4000 { v.push((175.0f32).ln() + 0.012 * ((k % 21) as f32 - 10.0) / 10.0); }
        for k in 0..2000 { v.push((100.0f32).ln() + 0.5 * (k as f32 / 2000.0)); }
        let (lo, hi) = two_modes(&v).unwrap();
        assert!((lo.exp() - 175.0).abs() < 4.0 && (hi.exp() - 205.0).abs() < 4.0, "{} {}", lo.exp(), hi.exp());
        assert!(two_modes(&v[..6000]).is_none(), "one mode is not two");
    }

    #[test]
    fn thinning_leaves_a_line() {
        let (w, h) = (20, 10); let mut m = vec![false; w * h];
        for y in 3..7 { for x in 2..18 { m[y * w + x] = true; } }
        thin(&mut m, w, h);
        let c = m.iter().filter(|v| **v).count();
        assert!(c > 8 && c < 22, "{c}");
    }
    #[test]
    fn a_dark_corner_reads_like_the_rest_of_the_face() {
        // the light falls by nearly half across the face, as on the Driver bottom photo; the same copper and laminate must read the same
        let paint = |x: u32, y: u32| { let base: f32 = if (100..120).contains(&y) { 66.0 } else { 100.0 }; let light: f32 = 1.0 - 0.45 * x as f32 / 1200.0; (base * light) as u8 };
        let im = board(1200, 240, &paint);
        let t = trace(&im, &Params::default());
        let (on_bright, on_dark) = (t.labels[50 * 1200 + 150], t.labels[50 * 1200 + 1050]);
        assert!(on_bright != 0 && on_dark != 0, "both ends read as copper: {on_bright} {on_dark}");
        assert_eq!(t.labels[110 * 1200 + 150], 0, "the lane at the bright end is a gap");
        assert_eq!(t.labels[110 * 1200 + 1050], 0, "the lane at the dark end is a gap");
    }
    #[test]
    fn a_thin_tentacle_is_cut_from_a_pour() {
        // a wide block with a thin track leaving it, in one label
        let (w, h) = (300, 200); let mut labels = vec![0u32; w * h];
        for y in 40..160 { for x in 20..140 { labels[y * w + x] = 1; } }
        for y in 96..104 { for x in 140..280 { labels[y * w + x] = 1; } }
        let (necks, _) = split_necks(&mut labels, w, h);
        assert_eq!(necks, 1);
        assert_ne!(labels[100 * w + 80], labels[100 * w + 220]);
    }
    #[test]
    fn a_pinch_between_two_pours_is_not_a_connection() {
        let (w, h) = (300, 120); let mut labels = vec![0u32; w * h];
        for y in 10..110 { for x in 10..130 { labels[y * w + x] = 1; } for x in 170..290 { labels[y * w + x] = 1; } }
        for y in 58..62 { for x in 130..170 { labels[y * w + x] = 1; } }
        split_necks(&mut labels, w, h);
        assert_ne!(labels[60 * w + 60], labels[60 * w + 230], "two pours joined only through a thin neck stay apart");
    }
    #[test]
    fn a_join_never_connects_two_large_pieces() {
        let (w, h) = (200, 60); let mut labels = vec![0u32; w * h];
        for y in 5..55 { for x in 5..97 { labels[y * w + x] = 1; } for x in 100..195 { labels[y * w + x] = 2; } }
        let value = Plane::new(w, h, 1.0); let resp = Plane::new(w, h, 0.0); let o = Mask::new(w, h);
        let (j, r) = join_fragments(&mut labels, &value, &resp, 10.0, &o, w, h);
        assert_eq!((j, r), (0, 1));
    }
}