BREP_kernel 0.7.0

//! Region-scoped offset REGULARITY — one definition of "does the equidistant
//! surface fold over THIS trimmed region", shared by `thicken`'s sheet gate and
//! `offset_shell`'s carrier-collapse lane.
//!
//! The condition itself is the classical one and neither caller changes it:
//! moving a surface by `δ` along its parametrization normal `n = Su × Sv`
//! scales the area element, per principal direction, by `1 − δ·κ` with `κ`
//! signed against that same normal ([`NurbsSurface::principal_curvatures`]).
//! The offset folds through the evolute where that factor reaches zero.
//!
//! What this module changes is WHERE the condition is asked.  Both callers used
//! to sweep a fixed grid over the whole surface DOMAIN, which answers about the
//! carrier and not about the face: a cone's apex row rejected a sheet trimmed to
//! the cone's wide end, and a fold narrower than the grid step INSIDE the trim
//! was not seen at all.  Here the samples are laid over the trim's own parameter
//! region, and their density is derived rather than fixed:
//!
//! * **Per knot span.**  A NURBS surface is one rational polynomial patch per
//!   knot span, so curvature cannot spike between spans without the spike living
//!   inside one.  Every span the region touches therefore gets its own samples,
//!   and a span covering 0.2% of the domain gets as many as one covering half of
//!   it.  A grid uniform over the domain is exactly what steps over the narrow
//!   one.
//! * **Per offset distance.**  Inside one span the fold's own scale is the
//!   offset distance — the fold happens where the curvature RADIUS falls to
//!   `|δ|` — so the arc-length step is additionally capped at a fraction of
//!   `|δ|`, converted to parameter steps through the first fundamental form.
//! * **Plus the trim boundary itself.**  The loop pcurves are sampled directly:
//!   the worst point of a region is often on its rim, and a region small enough
//!   to fall between grid nodes still gets a population this way.
//!
//! Refinement around the worst node sharpens the location that goes into a
//! refusal.  It never votes in the sample census, so the census stays an
//! unbiased population — which is what the shell lane's "collapsed EVERYWHERE in
//! this trim" question needs, and what a refinement biased toward the minimum
//! would quietly destroy.
//!
//! # Between the nodes
//!
//! The census cannot see a fold narrower than its node spacing, and that is
//! not hypothetical: on a free-form patch it let the carve keep a piece with a
//! folded ISLAND inside it (factor −2.1e-3, 0.03 across, between nodes 0.028
//! apart), which `thicken` then refused only as a self-intersecting result.
//! [`refine_between_nodes`] looks inside every grid cell whose corners agree but
//! read close to the level against the slope their neighbourhood shows, halving
//! the riskiest cells first. What it finds is reported beside the census
//! (`between_collapsed`, `between_regular`), never in it, and it seeds the fold
//! trace so an island comes back as the closed branch it is.
//!
//! What this does NOT do: no finite sampling can PROVE a region regular, and the
//! refinement's slope is estimated from the very nodes an island hides between.
//! So the claim is a measurement, on the island family in this module's tests
//! (an island is the folded top of a bump, at every width and depth the family
//! reaches, placed on a knot halfway between nodes or at random): with the SHEET
//! budget every island at least **0.02 of the node spacing** wide is caught, where
//! the grid alone missed 14 of those at least 0.1 wide; with SHELL_FACE every
//! island at least **0.25 of the node spacing** is. A narrower island is missed
//! when the nodes round it read far from the level (0.5 and more on the misses
//! measured — nothing computed from those nodes can see it), or when it is
//! narrower than the finest cell, `2⁻⁷` of the node spacing. The refinement has
//! no evaluation budget to run out of: each risky cell costs at most
//! `5·(4⁷ − 1)/3` = 27,305 samples, so it always finishes, and nothing it did
//! not reach is reported regular. What it costs on the lib, integration and
//! case populations is in the 2026-09-17 record. Away from the between-node question
//! the density rules still hold: a fold has to hide inside a single knot span,
//! below the `|δ|`-derived arc-length step, and away from the rim, to be missed.

use crate::{KernelRefusal, KernelStage, OrRefuse};
use crate::topology::FaceRecord;
use crate::{NurbsCurve, NurbsSurface};

/// The offset family's own trace switch, the same `BREP_OS_DEBUG` the shell
/// reads — this module answers for both callers, so it cannot borrow either
/// one's macro.
fn debug_enabled() -> bool {
    std::env::var("BREP_OS_DEBUG").is_ok_and(|value| !value.is_empty() && value != "0")
}

/// The worst sample a scan saw: the fold factor, where it is, and the
/// displacement and curvature that produced it.
#[derive(Clone, Copy, Debug)]
pub(crate) struct FoldSample {
    /// `1 − δ·κ`; at or below the caller's collapse factor the offset folds.
    pub(crate) factor: f64,
    pub(crate) u: f64,
    pub(crate) v: f64,
    /// The principal curvature that produced `factor`, signed against `Su × Sv`.
    pub(crate) kappa: f64,
    /// The displacement along `Su × Sv` that produced `factor`.
    pub(crate) displacement: f64,
}

impl FoldSample {
    /// The concave curvature radius the displacement is measured against.
    pub(crate) fn radius(&self) -> f64 {
        1.0 / self.kappa.abs().max(1e-300)
    }
}

/// What a scan found over a region.
#[derive(Clone, Copy, Debug, Default)]
pub(crate) struct RegularityScan {
    /// Census population: grid and trim-boundary samples inside the region at
    /// which the curvature could be read at all.
    pub(crate) sampled: usize,
    /// How many of those collapsed.
    pub(crate) collapsed: usize,
    /// The worst factor seen, refinement included.
    pub(crate) worst: Option<FoldSample>,
    /// What the between-node refinement found that the census could not: cells
    /// whose four nodes all read regular but which hold a collapsed sample, and
    /// cells whose nodes all read collapsed but which hold a regular one. Not
    /// census members — the census stays the unbiased grid — but a region with
    /// either is not what its census says it is. See [`refine_between_nodes`].
    pub(crate) between_collapsed: usize,
    pub(crate) between_regular: usize,
    /// What that refinement cost.
    pub(crate) between_evaluations: usize,
    /// Cells still risky at the finest spacing the refinement reaches.
    pub(crate) between_unresolved: usize,
}

/// How hard a caller wants its region swept.
#[derive(Clone, Copy, Debug)]
pub(crate) struct ScanBudget {
    /// Samples per knot span per direction, before the node cap.
    pub(crate) per_span: usize,
    /// Arc-length steps per offset distance: the step is `|δ| / this`.
    pub(crate) steps_per_distance: f64,
    /// Ceiling on the grid (nodes, before the region test).
    pub(crate) max_nodes: usize,
    /// Refine around the worst node to sharpen the reported location.
    pub(crate) refine: bool,
}

impl ScanBudget {
    /// The sheet gate: one scan per thicken, so it can afford to hunt.  The
    /// per-span floor is what sees a tight bend that occupies a sliver of the
    /// parameter domain, and the refinement is what puts a usable `(u, v)` in
    /// the refusal instead of the nearest node.
    pub(crate) const SHEET: Self = Self {
        per_span: 9,
        steps_per_distance: 4.0,
        max_nodes: 40_000,
        refine: true,
    };

    /// The shell's per-face collapse probe.  It runs on every non-analytic face
    /// of a shell and its question is "is this trim collapsed EVERYWHERE", which
    /// buys coverage rather than spike-hunting — so it takes a coverage budget
    /// of the same order as the 7×7 whole-face grid it replaces, spent over the
    /// trim instead of over the carrier.
    pub(crate) const SHELL_FACE: Self = Self {
        per_span: 3,
        steps_per_distance: 1.0,
        max_nodes: 512,
        refine: false,
    };
}

/// A trimmed parameter region: the trim loops as `(u, v)` polylines, plus the
/// parameter box they span.
///
/// Winding is deliberately not consulted.  Containment is crossing parity over
/// every loop at once, which puts a point in a hole OUTSIDE the region however
/// the hole is wound, and which needs no convexity or nesting assumption.
pub(crate) struct TrimRegion {
    loops: Vec<Vec<[f64; 2]>>,
    bounds: [f64; 4],
}

impl TrimRegion {
    /// The whole carrier: what a face with no loops covers, and the region a
    /// caller means when it has no trim to name.
    pub(crate) fn whole_domain(surface: &NurbsSurface) -> Result<Self, KernelRefusal> {
        let [u0, u1] = surface.domain_u().or_refuse(KernelStage::Classify, "domain_u")?;
        let [v0, v1] = surface.domain_v().or_refuse(KernelStage::Classify, "domain_v")?;
        Ok(Self {
            loops: Vec::new(),
            bounds: [u0, u1, v0, v1],
        })
    }

    /// A region from bare parameter-space loops — `thicken`'s form, where the
    /// trim is passed in rather than stored on a face.
    ///
    /// Malformed loops are not this module's refusal to make: a loop that does
    /// not close is polylined and closed implicitly, so the region is still a
    /// region and the caller's own loop validation keeps the right message.
    pub(crate) fn from_pcurve_loops(
        surface: &NurbsSurface,
        loops: &[Vec<NurbsCurve>],
    ) -> Result<Self, KernelRefusal> {
        let mut polylines = Vec::with_capacity(loops.len());
        for loop_curves in loops {
            let mut polyline = Vec::new();
            for curve in loop_curves {
                append_pcurve_polyline(curve, &mut polyline)?;
            }
            if polyline.len() >= 3 {
                polylines.push(polyline);
            }
        }
        Self::from_polylines(surface, polylines)
    }

    /// A region from a face's own stored trim.
    pub(crate) fn from_face(face: &FaceRecord) -> Result<Self, KernelRefusal> {
        let mut polylines = Vec::with_capacity(face.loops.len());
        for loop_record in &face.loops {
            let mut polyline = Vec::new();
            for coedge in &loop_record.coedges {
                append_pcurve_polyline(&coedge.pcurve, &mut polyline)?;
            }
            if polyline.len() >= 3 {
                polylines.push(polyline);
            }
        }
        Self::from_polylines(&face.surface, polylines)
    }

    fn from_polylines(
        surface: &NurbsSurface,
        loops: Vec<Vec<[f64; 2]>>,
    ) -> Result<Self, KernelRefusal> {
        let [u0, u1] = surface.domain_u().or_refuse(KernelStage::Classify, "domain_u")?;
        let [v0, v1] = surface.domain_v().or_refuse(KernelStage::Classify, "domain_v")?;
        if loops.is_empty() {
            return Self::whole_domain(surface);
        }
        let (mut lo_u, mut hi_u) = (f64::INFINITY, f64::NEG_INFINITY);
        let (mut lo_v, mut hi_v) = (f64::INFINITY, f64::NEG_INFINITY);
        for polyline in &loops {
            for point in polyline {
                lo_u = lo_u.min(point[0]);
                hi_u = hi_u.max(point[0]);
                lo_v = lo_v.min(point[1]);
                hi_v = hi_v.max(point[1]);
            }
        }
        // A pcurve may run a rounding step outside its carrier's domain; the
        // region is what the carrier can be evaluated over.
        let bounds = [
            lo_u.max(u0).min(u1),
            hi_u.max(u0).min(u1),
            lo_v.max(v0).min(v1),
            hi_v.max(v0).min(v1),
        ];
        if !(bounds[1] > bounds[0]) || !(bounds[3] > bounds[2]) {
            // A degenerate trim box is not a region to judge: fall back to the
            // carrier so nothing is silently skipped — and SAY so, because
            // "why did this scan the whole carrier again" is otherwise a
            // bisect.
            if debug_enabled() {
                eprintln!(
                    "regularity: trim box degenerate (u[{:.6},{:.6}] v[{:.6},{:.6}]) — scanning \
                     the whole carrier instead",
                    bounds[0], bounds[1], bounds[2], bounds[3]
                );
            }
            return Self::whole_domain(surface);
        }
        Ok(Self { loops, bounds })
    }

    /// `[u0, u1, v0, v1]` — the trim's own parameter box, clamped to the
    /// carrier's domain.
    pub(crate) fn bounds(&self) -> [f64; 4] {
        self.bounds
    }

    /// Crossing parity against every loop at once.
    pub(crate) fn contains(&self, u: f64, v: f64) -> bool {
        if self.loops.is_empty() {
            return u >= self.bounds[0]
                && u <= self.bounds[1]
                && v >= self.bounds[2]
                && v <= self.bounds[3];
        }
        let mut inside = false;
        for polyline in &self.loops {
            let count = polyline.len();
            let mut previous = count - 1;
            for current in 0..count {
                let (a, b) = (polyline[previous], polyline[current]);
                if (a[1] > v) != (b[1] > v) {
                    let t = (v - a[1]) / (b[1] - a[1]);
                    if u < a[0] + t * (b[0] - a[0]) {
                        inside = !inside;
                    }
                }
                previous = current;
            }
        }
        inside
    }

    fn boundary_points(&self) -> impl Iterator<Item = [f64; 2]> + '_ {
        self.loops.iter().flat_map(|polyline| polyline.iter().copied())
    }
}

/// Sample a parameter-space curve into the polyline that stands for it.
///
/// Per knot span, for the same reason the surface grid is: a pcurve's shape is
/// one polynomial per span, so a span is the unit that cannot be stepped over.
/// The last parameter is deliberately omitted — the next pcurve of the loop
/// starts there, and a closed single-pcurve loop closes on its own.
fn append_pcurve_polyline(curve: &NurbsCurve, into: &mut Vec<[f64; 2]>) -> Result<(), KernelRefusal> {
    let [t0, t1] = curve.domain().or_refuse(KernelStage::Classify, "domain")?;
    if !(t1 > t0) {
        return Ok(());
    }
    let spans = breakpoints(&curve.knots, t0, t1).len().saturating_sub(1).max(1);
    let count = (8 * spans).clamp(16, 256);
    for index in 0..count {
        let t = t0 + (t1 - t0) * index as f64 / count as f64;
        let point = curve.evaluate(t).or_refuse(KernelStage::Classify, "evaluate")?;
        into.push([point.x, point.y]);
    }
    Ok(())
}

/// The distinct knot values bounding the spans that `[lo, hi]` touches,
/// `lo` and `hi` included.
fn breakpoints(knots: &[f64], lo: f64, hi: f64) -> Vec<f64> {
    let width = (hi - lo).abs();
    let epsilon = 1e-12 * width.max(1.0);
    let mut points = vec![lo];
    for &knot in knots {
        if knot > lo + epsilon && knot < hi - epsilon && knot > *points.last().unwrap() + epsilon {
            points.push(knot);
        }
    }
    points.push(hi);
    points
}

/// Drop interior breakpoints until the span grid fits under `max_spans`.
fn thin(points: &mut Vec<f64>, max_spans: usize) {
    while points.len() > max_spans.max(1) + 1 {
        let kept: Vec<f64> = points
            .iter()
            .enumerate()
            .filter(|(index, _)| *index == 0 || *index == points.len() - 1 || index % 2 == 0)
            .map(|(_, value)| *value)
            .collect();
        if kept.len() == points.len() {
            break;
        }
        *points = kept;
    }
}

/// Strictly interior samples, `per_span` of them inside every span.
///
/// Interior on purpose: a pole or a seam sits ON a domain boundary, where the
/// parametrization is degenerate and the curvature is unreadable, and a knot is
/// where a one-sided curvature would be read.
fn span_samples(breaks: &[f64], per_span: usize) -> Vec<f64> {
    let per_span = per_span.max(1);
    let mut samples = Vec::with_capacity(breaks.len() * per_span);
    for window in breaks.windows(2) {
        let (lo, hi) = (window[0], window[1]);
        for index in 0..per_span {
            samples.push(lo + (hi - lo) * (index as f64 + 0.5) / per_span as f64);
        }
    }
    samples
}

/// Largest `|Su|` and `|Sv|` over a coarse pass of the region box — the first
/// fundamental form's diagonal, which is what turns a 3D step into a parameter
/// step.
fn metric_scale(surface: &NurbsSurface, bounds: [f64; 4]) -> (f64, f64) {
    const PROBES: usize = 5;
    let (mut su, mut sv) = (0.0f64, 0.0f64);
    for iu in 0..PROBES {
        let u = bounds[0] + (bounds[1] - bounds[0]) * (iu as f64 + 0.5) / PROBES as f64;
        for iv in 0..PROBES {
            let v = bounds[2] + (bounds[3] - bounds[2]) * (iv as f64 + 0.5) / PROBES as f64;
            if let Ok(derivatives) = surface.derivatives(u, v, 1) {
                su = su.max(derivatives[1][0].length());
                sv = sv.max(derivatives[0][1].length());
            }
        }
    }
    (su, sv)
}

/// The sample nodes a scan lays over a region: `per_span` per knot span in each
/// direction, floored by the arc-length cap and capped by the node budget.
///
/// Extracted so the fold TRACE (`offset/fold_locus.rs`) brackets the locus on
/// exactly the nodes the scan censuses.  "The scan can already bracket the fold
/// but does not trace it" is then literally true rather than approximately: a
/// tracer with a grid of its own could seed on a crossing this scan never saw,
/// and the two would disagree about the same face.
pub(crate) struct RegionGrid {
    pub(crate) samples_u: Vec<f64>,
    pub(crate) samples_v: Vec<f64>,
}

pub(crate) fn region_grid(
    surface: &NurbsSurface,
    bounds: [f64; 4],
    reach: f64,
    budget: ScanBudget,
) -> RegionGrid {
    let mut breaks_u = breakpoints(&surface.knots_u, bounds[0], bounds[1]);
    let mut breaks_v = breakpoints(&surface.knots_v, bounds[2], bounds[3]);
    // A fitted carrier can carry hundreds of knots; the span grid alone must
    // still fit the budget before anything is spent inside a span.
    let max_spans = (budget.max_nodes as f64).sqrt().max(2.0) as usize;
    thin(&mut breaks_u, max_spans);
    thin(&mut breaks_v, max_spans);
    let (spans_u, spans_v) = (breaks_u.len() - 1, breaks_v.len() - 1);

    let mut per_u = budget.per_span.max(1);
    let mut per_v = budget.per_span.max(1);
    // Arc-length cap: the fold's own scale is the offset distance, so no step
    // may cover more 3D length than a fraction of it.
    let step = reach / budget.steps_per_distance.max(1.0);
    let (metric_u, metric_v) = metric_scale(surface, bounds);
    let needed = |metric: f64, width: f64, spans: usize| -> usize {
        if metric <= 0.0 || step <= 0.0 {
            return 1;
        }
        ((metric * width / step) / spans as f64).ceil().max(1.0).min(512.0) as usize
    };
    per_u = per_u.max(needed(metric_u, bounds[1] - bounds[0], spans_u));
    per_v = per_v.max(needed(metric_v, bounds[3] - bounds[2], spans_v));
    while spans_u * per_u * spans_v * per_v > budget.max_nodes && (per_u > 1 || per_v > 1) {
        if per_u >= per_v {
            per_u -= 1;
        } else {
            per_v -= 1;
        }
    }
    RegionGrid {
        samples_u: span_samples(&breaks_u, per_u),
        samples_v: span_samples(&breaks_v, per_v),
    }
}

/// Worst `1 − δ·κ` at one parameter, over every requested displacement and both
/// principal directions — the scalar FIELD whose zero set is the fold locus.
///
/// `min` over the two principal factors is continuous (both principal
/// curvatures are continuous functions of the fundamental forms), so the field
/// is continuous wherever the parametrization is regular, and smooth away from
/// the umbilics where the two branches cross.  At the locus itself the branch
/// that vanishes is strictly the smaller one, so the crossing is generically
/// transversal and the trace in `offset/fold_locus.rs` can march it.
///
/// `None` where the curvature cannot be read at all (a pole, a seam, a
/// degenerate parametrization) — not zero, which would be a fold.
pub(crate) fn fold_sample_at(
    surface: &NurbsSurface,
    u: f64,
    v: f64,
    displacements: &[f64],
) -> Option<FoldSample> {
    let (kappa_min, kappa_max) = surface.principal_curvatures(u, v).ok()?;
    let mut worst: Option<FoldSample> = None;
    for &displacement in displacements {
        if displacement == 0.0 {
            continue;
        }
        for kappa in [kappa_min, kappa_max] {
            let factor = 1.0 - displacement * kappa;
            if worst.is_none_or(|current| factor < current.factor) {
                worst = Some(FoldSample {
                    factor,
                    u,
                    v,
                    kappa,
                    displacement,
                });
            }
        }
    }
    worst
}

impl RegularityScan {
    fn keep_worst(&mut self, candidate: FoldSample) {
        if self.worst.is_none_or(|current| candidate.factor < current.factor) {
            self.worst = Some(candidate);
        }
    }
}

/// Sweep `region` for the worst offset fold factor.
///
/// `displacements` are signed moves ALONG the parametrization normal `Su × Sv`
/// — the one convention this module has, so each caller converts its own signed
/// distance once at the call site instead of re-deriving a sign here.
pub(crate) fn scan_offset_regularity(
    surface: &NurbsSurface,
    region: &TrimRegion,
    displacements: &[f64],
    collapse_factor: f64,
    budget: ScanBudget,
) -> Result<RegularityScan, KernelRefusal> {
    let mut scan = RegularityScan::default();
    let reach = displacements
        .iter()
        .fold(0.0f64, |worst, distance| worst.max(distance.abs()));
    // A plane has no curvature to read and no offset that can fold.
    if reach == 0.0 || surface.is_affine().or_refuse(KernelStage::Classify, "is_affine")? {
        return Ok(scan);
    }
    let bounds = region.bounds();
    let RegionGrid {
        samples_u,
        samples_v,
    } = region_grid(surface, bounds, reach, budget);

    let record = |scan: &mut RegularityScan, sample: FoldSample| {
        scan.sampled += 1;
        if sample.factor <= collapse_factor {
            scan.collapsed += 1;
        }
        scan.keep_worst(sample);
    };
    // Every node of the box is read, inside the region or not: the census counts
    // only the ones inside, and the between-node refinement below needs the
    // others as the far corners of the cells the region's own rim passes through.
    let nodes = samples_u
        .iter()
        .map(|&u| {
            samples_v
                .iter()
                .map(|&v| fold_sample_at(surface, u, v, displacements))
                .collect::<Vec<_>>()
        })
        .collect::<Vec<_>>();
    for (iu, &u) in samples_u.iter().enumerate() {
        for (iv, &v) in samples_v.iter().enumerate() {
            if let Some(sample) = nodes[iu][iv] {
                if region.contains(u, v) {
                    record(&mut scan, sample);
                }
            }
        }
    }
    // The rim is part of the region whatever parity says about a point sitting
    // exactly on it, and on a trim too small to catch a grid node it is the
    // whole population.  It is also the one part of the population a caller did
    // not size: a booleaned face carries as many polyline points as it has
    // knotty pcurves, which on the shell lane would outvote the grid the budget
    // was spent on.  Stride it to the same ceiling.
    let rim: Vec<[f64; 2]> = region.boundary_points().collect();
    let stride = (rim.len() / budget.max_nodes.max(1)).max(1);
    for point in rim.iter().step_by(stride) {
        let u = point[0].clamp(bounds[0], bounds[1]);
        let v = point[1].clamp(bounds[2], bounds[3]);
        if let Some(sample) = fold_sample_at(surface, u, v, displacements) {
            record(&mut scan, sample);
        }
    }

    let values = nodes
        .iter()
        .map(|column| {
            column
                .iter()
                .map(|sample| sample.map(|sample| sample.factor - collapse_factor))
                .collect::<Vec<_>>()
        })
        .collect::<Vec<_>>();
    let between = refine_between_nodes(
        surface,
        region,
        &samples_u,
        &samples_v,
        &values,
        displacements,
        collapse_factor,
    );
    scan.between_collapsed = between.collapsed;
    scan.between_regular = between.regular;
    scan.between_evaluations = between.evaluations;
    scan.between_unresolved = between.unresolved;
    if let Some(worst) = between.worst {
        scan.keep_worst(worst);
    }
    if debug_enabled() && between.evaluations > 0 {
        eprintln!(
            "regularity: between-node refinement read {} sample(s) in {} risky cell(s) — {} \
             collapsed and {} regular sample(s) the census's own nodes did not show, {} cell(s) \
             still risky at the finest spacing",
            between.evaluations,
            between.cells,
            between.collapsed,
            between.regular,
            between.unresolved
        );
    }

    if budget.refine {
        refine_worst(surface, region, displacements, &mut scan, bounds);
    }
    Ok(scan)
}

/// The least `1 − δ·κ` over a face's own trim, for a caller asking its own
/// question of the same field rather than the census's.
///
/// What it reports is what [`scan_offset_regularity`] read: the grid, the rim,
/// the refinement around the worst node when the budget asks for it, and the
/// between-node refinement. `bound` says how far that reading can be trusted
/// (see [`LeastFoldFactor::bound`]).
// Built for the outward-shell guard (`outward_fold_the_carve_missed`), which
// sweeps the same field on a grid of its own; it moves onto this entry point as
// a separate change, and until then only the tests call it.
#[cfg_attr(not(test), allow(dead_code))]
#[derive(Clone, Copy, Debug)]
pub(crate) struct LeastFoldFactor {
    /// The least factor found and where, or `None` when nothing was readable.
    pub(crate) least: Option<FoldSample>,
    /// The scan grid's largest node spacing, per direction.
    pub(crate) node_spacing: [f64; 2],
    /// The finest cell the between-node refinement may reach, per direction:
    /// `node_spacing / 2^REFINE_DEPTH`.
    pub(crate) finest_spacing: [f64; 2],
    /// Every cell whose nodes agreed was either cleared by its neighbourhood's
    /// slope bound or refined to a sample of the other sign. `false` when a cell
    /// was still risky at the finest spacing — near a fold locus that is the
    /// normal case, and the least factor is then the least SEEN at that spacing.
    pub(crate) bound: bool,
    pub(crate) evaluations: usize,
}

/// [`LeastFoldFactor`] over `face`'s own trim, with `displacements` along the
/// parametrization normal `Su × Sv` (this module's one convention) and `level`
/// the fold factor at which a sample counts as collapsed.
#[cfg_attr(not(test), allow(dead_code))]
pub(crate) fn least_fold_factor(
    face: &FaceRecord,
    displacements: &[f64],
    level: f64,
    budget: ScanBudget,
) -> Result<LeastFoldFactor, KernelRefusal> {
    let region = TrimRegion::from_face(face)?;
    let reach = displacements
        .iter()
        .fold(0.0f64, |worst, distance| worst.max(distance.abs()));
    let grid = region_grid(&face.surface, region.bounds(), reach, budget);
    let spacing = |samples: &[f64]| {
        samples
            .windows(2)
            .map(|pair| pair[1] - pair[0])
            .fold(0.0f64, f64::max)
    };
    let node_spacing = [spacing(&grid.samples_u), spacing(&grid.samples_v)];
    let scan = scan_offset_regularity(&face.surface, &region, displacements, level, budget)?;
    let finest = (1u64 << REFINE_DEPTH) as f64;
    Ok(LeastFoldFactor {
        least: scan.worst,
        node_spacing,
        finest_spacing: [node_spacing[0] / finest, node_spacing[1] / finest],
        bound: scan.between_unresolved == 0,
        evaluations: scan.between_evaluations,
    })
}

/// How many times the variation a cell's neighbourhood shows across half the
/// cell its own nodes have to clear before the cell is taken to hold no sign
/// change. The slope is ESTIMATED from node differences, which is exactly what
/// underestimates a feature between nodes, so the margin is doubled; the
/// island family in this module's tests measures what that buys.
const REFINE_SAFETY: f64 = 2.0;
/// How many times a risky cell may be halved: 2⁷ = 128 subdivisions per node
/// spacing.
const REFINE_DEPTH: usize = 7;

/// What [`refine_between_nodes`] found.
#[derive(Clone, Debug, Default)]
pub(crate) struct BetweenNodes {
    /// Refined samples inside the region that read at or below the level, in
    /// cells whose four nodes all read above it.
    pub(crate) collapsed: usize,
    /// Refined samples inside the region that read above the level, in cells
    /// whose four nodes all read at or below it.
    pub(crate) regular: usize,
    /// A `(u, v)` pair straddling the level for every such find — one refined
    /// sample and the corner of its cell with the other sign — so the fold
    /// TRACE can seed on a branch the grid never bracketed.
    pub(crate) brackets: Vec<([f64; 2], [f64; 2])>,
    /// The worst refined sample inside the region.
    pub(crate) worst: Option<FoldSample>,
    pub(crate) cells: usize,
    pub(crate) evaluations: usize,
    /// Cells still risky at [`REFINE_DEPTH`]: the bound was not cleared there,
    /// and nothing of the other sign was found either.
    pub(crate) unresolved: usize,
}

/// Look BETWEEN the scan's nodes for a sign change the nodes cannot show.
///
/// A cell whose four corners agree can still hold an island of the other sign
/// — a shallow fold narrower than the node spacing, which is how a carved
/// piece was once kept with a folded island inside it. A cell is RISKY when the
/// smallest `|f|` at its corners is below [`REFINE_SAFETY`] times the
/// variation its 3×3 neighbourhood of cells shows over half the cell,
/// `(L_u·Δu + L_v·Δv)/2` with `L` the largest node-difference slope per
/// direction. Risky cells are halved, most risky first, recomputing the bound
/// on each child with its own corners (never below the parent's slopes), until
/// the bound clears, a sample of the other sign turns up, or the depth reaches
/// [`REFINE_DEPTH`]. There is no evaluation budget: at most
/// `5·(4^REFINE_DEPTH − 1)/3` samples per risky cell, so it always finishes.
///
/// Near a real fold locus every cell is risky and the refinement costs up to
/// `REFINE_DEPTH` levels along it; away from one `|f|` is large against the
/// neighbourhood's variation and nothing is read.
///
/// `values[iu][iv]` is `factor − level` at `(samples_u[iu], samples_v[iv])`.
#[allow(clippy::too_many_arguments)]
pub(crate) fn refine_between_nodes(
    surface: &NurbsSurface,
    region: &TrimRegion,
    samples_u: &[f64],
    samples_v: &[f64],
    values: &[Vec<Option<f64>>],
    displacements: &[f64],
    level: f64,
) -> BetweenNodes {
    let mut found = BetweenNodes::default();
    let (nu, nv) = (samples_u.len(), samples_v.len());
    if nu < 2 || nv < 2 {
        return found;
    }
    let value = |iu: usize, iv: usize| values[iu][iv];
    // Largest node-difference slope per direction over the 3×3 cells around
    // cell (iu, iv).
    let slopes = |iu: usize, iv: usize| -> (f64, f64) {
        let (mut lu, mut lv) = (0.0f64, 0.0f64);
        let u_lo = iu.saturating_sub(1);
        let u_hi = (iu + 2).min(nu - 1);
        let v_lo = iv.saturating_sub(1);
        let v_hi = (iv + 2).min(nv - 1);
        for a in u_lo..=u_hi {
            for b in v_lo..=v_hi {
                if a < u_hi {
                    if let (Some(p), Some(q)) = (value(a, b), value(a + 1, b)) {
                        lu = lu.max((q - p).abs() / (samples_u[a + 1] - samples_u[a]));
                    }
                }
                if b < v_hi {
                    if let (Some(p), Some(q)) = (value(a, b), value(a, b + 1)) {
                        lv = lv.max((q - p).abs() / (samples_v[b + 1] - samples_v[b]));
                    }
                }
            }
        }
        (lu, lv)
    };

    struct Cell {
        u: [f64; 2],
        v: [f64; 2],
        /// Corner values in the order (u0,v0), (u1,v0), (u0,v1), (u1,v1).
        corners: [f64; 4],
        slopes: (f64, f64),
        depth: usize,
        risk: f64,
    }
    // Ordered so the heap pops the SMALLEST risk ratio — the riskiest cell — first.
    impl PartialEq for Cell {
        fn eq(&self, other: &Self) -> bool {
            self.risk.total_cmp(&other.risk).is_eq()
        }
    }
    impl Eq for Cell {}
    impl PartialOrd for Cell {
        fn partial_cmp(&self, other: &Self) -> Option<std::cmp::Ordering> {
            Some(self.cmp(other))
        }
    }
    impl Ord for Cell {
        fn cmp(&self, other: &Self) -> std::cmp::Ordering {
            other.risk.total_cmp(&self.risk)
        }
    }
    let risk_of = |corners: &[f64; 4], du: f64, dv: f64, (lu, lv): (f64, f64)| -> Option<f64> {
        let margin = REFINE_SAFETY * (lu * du + lv * dv) / 2.0;
        let smallest = corners.iter().fold(f64::INFINITY, |m, c| m.min(c.abs()));
        (margin > 0.0 && smallest < margin).then(|| smallest / margin)
    };
    let mut work = std::collections::BinaryHeap::<Cell>::new();
    for iu in 0..nu - 1 {
        for iv in 0..nv - 1 {
            let (Some(a), Some(b), Some(c), Some(d)) =
                (value(iu, iv), value(iu + 1, iv), value(iu, iv + 1), value(iu + 1, iv + 1))
            else {
                continue;
            };
            let corners = [a, b, c, d];
            let negative = corners.iter().filter(|x| **x <= 0.0).count();
            if negative != 0 && negative != 4 {
                continue; // the grid already brackets this cell
            }
            let (u0, u1) = (samples_u[iu], samples_u[iu + 1]);
            let (v0, v1) = (samples_v[iv], samples_v[iv + 1]);
            let touches = region.contains(u0, v0)
                || region.contains(u1, v0)
                || region.contains(u0, v1)
                || region.contains(u1, v1)
                || region.contains(0.5 * (u0 + u1), 0.5 * (v0 + v1));
            if !touches {
                continue;
            }
            let cell_slopes = slopes(iu, iv);
            if let Some(risk) = risk_of(&corners, u1 - u0, v1 - v0, cell_slopes) {
                work.push(Cell {
                    u: [u0, u1],
                    v: [v0, v1],
                    corners,
                    slopes: cell_slopes,
                    depth: 0,
                    risk,
                });
            }
        }
    }
    found.cells = work.len();
    // No evaluation budget: the quadtree below one risky cell has at most
    // (4^REFINE_DEPTH − 1)/3 cells to process at five samples each, so the
    // refinement always finishes, and a cap could only ever stop it early —
    // which is the island escape again with a different trigger.
    while let Some(cell) = work.pop() {
        let collapsed_side = cell.corners[0] <= 0.0;
        let (um, vm) = (0.5 * (cell.u[0] + cell.u[1]), 0.5 * (cell.v[0] + cell.v[1]));
        let points = [
            [um, vm],
            [um, cell.v[0]],
            [cell.u[0], vm],
            [cell.u[1], vm],
            [um, cell.v[1]],
        ];
        let mut read = [None; 5];
        let mut flipped = false;
        for (index, point) in points.iter().enumerate() {
            found.evaluations += 1;
            let Some(sample) = fold_sample_at(surface, point[0], point[1], displacements) else {
                continue;
            };
            let f = sample.factor - level;
            read[index] = Some(f);
            let inside = region.contains(point[0], point[1]);
            if inside && found.worst.is_none_or(|worst| sample.factor < worst.factor) {
                found.worst = Some(sample);
            }
            if (f <= 0.0) != collapsed_side {
                // A sign change inside this cell whether or not the sample is in
                // the region: the trace needs the bracket either way (its branch
                // may enter the region), and a child cell built on this sample
                // would no longer have corners of one sign. Only a sample in
                // the region counts against the region's census.
                flipped = true;
                if inside && f <= 0.0 {
                    found.collapsed += 1;
                } else if inside {
                    found.regular += 1;
                }
                // Straddle the level: this sample and the corner nearest it.
                let corner = [
                    [cell.u[0], cell.v[0]],
                    [cell.u[1], cell.v[0]],
                    [cell.u[0], cell.v[1]],
                    [cell.u[1], cell.v[1]],
                ]
                .into_iter()
                .min_by(|p, q| {
                    let dp = (p[0] - point[0]).hypot(p[1] - point[1]);
                    let dq = (q[0] - point[0]).hypot(q[1] - point[1]);
                    dp.partial_cmp(&dq).unwrap_or(std::cmp::Ordering::Equal)
                })
                .expect("four corners");
                found.brackets.push((*point, corner));
            }
        }
        if flipped {
            continue;
        }
        let [Some(centre), Some(bottom), Some(left), Some(right), Some(top)] = read else {
            continue;
        };
        let [c00, c10, c01, c11] = cell.corners;
        let children = [
            ([cell.u[0], um], [cell.v[0], vm], [c00, bottom, left, centre]),
            ([um, cell.u[1]], [cell.v[0], vm], [bottom, c10, centre, right]),
            ([cell.u[0], um], [vm, cell.v[1]], [left, centre, c01, top]),
            ([um, cell.u[1]], [vm, cell.v[1]], [centre, right, top, c11]),
        ];
        for (u, v, corners) in children {
            let (du, dv) = (u[1] - u[0], v[1] - v[0]);
            let own = (
                ((corners[1] - corners[0]).abs().max((corners[3] - corners[2]).abs())) / du,
                ((corners[2] - corners[0]).abs().max((corners[3] - corners[1]).abs())) / dv,
            );
            let child_slopes = (cell.slopes.0.max(own.0), cell.slopes.1.max(own.1));
            if let Some(risk) = risk_of(&corners, du, dv, child_slopes) {
                if cell.depth + 1 >= REFINE_DEPTH {
                    found.unresolved += 1;
                    continue;
                }
                work.push(Cell {
                    u,
                    v,
                    corners,
                    slopes: child_slopes,
                    depth: cell.depth + 1,
                    risk,
                });
            }
        }
    }
    found
}

/// Bisect toward the worst node so the refusal names where the fold IS, not the
/// nearest node to it.  Census counts are untouched: refinement hunts a minimum
/// and would bias any population it joined.
fn refine_worst(
    surface: &NurbsSurface,
    region: &TrimRegion,
    displacements: &[f64],
    scan: &mut RegularityScan,
    bounds: [f64; 4],
) {
    const LEVELS: usize = 4;
    const PROBES: usize = 4;
    let Some(start) = scan.worst else { return };
    let mut radius_u = (bounds[1] - bounds[0]) / 16.0;
    let mut radius_v = (bounds[3] - bounds[2]) / 16.0;
    let (mut centre_u, mut centre_v) = (start.u, start.v);
    for _ in 0..LEVELS {
        let mut best: Option<FoldSample> = None;
        for iu in 0..=2 * PROBES {
            let u = (centre_u + radius_u * (iu as f64 / PROBES as f64 - 1.0))
                .clamp(bounds[0], bounds[1]);
            for iv in 0..=2 * PROBES {
                let v = (centre_v + radius_v * (iv as f64 / PROBES as f64 - 1.0))
                    .clamp(bounds[2], bounds[3]);
                if !region.contains(u, v) {
                    continue;
                }
                if let Some(sample) = fold_sample_at(surface, u, v, displacements) {
                    if best.is_none_or(|current| sample.factor < current.factor) {
                        best = Some(sample);
                    }
                }
            }
        }
        if let Some(sample) = best {
            scan.keep_worst(sample);
            centre_u = sample.u;
            centre_v = sample.v;
        }
        radius_u *= 0.5;
        radius_v *= 0.5;
    }
}