geop_ops_drawing/

silhouette.rs

//! Silhouettes: where a curved face turns away from the eye, the outline it
//! shows without any edge there — the sides of a cylinder seen side on.
//!
//! For a view along `d`, a point of a surface `Scalar(u, v)` is on its silhouette
//! when its normal is perpendicular to `d`: `g(u, v) = (S_u x S_v) . d = 0`.
//! That is a curve in `(u, v)`, traced here as such:
//!
//! 1. **Seeds.** `g`'s sign is sampled along a grid of iso lines of the
//!    surface's domain. Between two samples of opposite sign a silhouette
//!    crosses the line, and bisection finds where. Sampling only decides
//!    where to look: what is traced is the exact zero set. The interior grid
//!    lines sit at irrational fractions of the domain, a free choice made so
//!    they do not run along the symmetric iso lines silhouettes so often are.
//!    A closed loop crossing no grid line lies within one cell, and encloses
//!    a critical point of `g` (an extremum): each cell where `g` could
//!    vanish is halved into boxes, those that could hold a critical point
//!    halved again, and through each left a line across the cell is sampled
//!    too ([`loop_seats`]), crossing such a loop. A trace seeded on a short
//!    line steps no further than half of it, so a small loop is traced in
//!    steps of its own size.
//! 2. **Tracing.** From a seed, predictor-corrector steps follow the curve:
//!    a step along the tangent `(-g_v, g_u)`, sized to a fixed length in 3-D,
//!    then Newton back onto `g = 0`. A trace ends where the curve leaves the
//!    domain, landing exactly on its boundary, or where it closes on itself.
//!    Every grid line a trace crosses marks the seed there as traced.
//! 3. **Trimming.** The traced curve in `(u, v)` is cut where it crosses the
//!    face's trim curves, and only the pieces inside the face are kept.
//! 4. **Curves.** Each piece becomes a 3-D curve through exact silhouette
//!    points along it, widened to enclose the silhouette between them too
//!    (as a traced intersection branch is, `NurbCurve3D::fit_enclosing`):
//!    two silhouettes that are one line on the paper — mirror images seen
//!    along the mirror's normal — are then seen to run along each other.
//!
//! The curve encloses the silhouette, so a point of it is near the
//! silhouette, within its width, not on it. Whoever needs a point exactly
//! on it — the visibility test does: a ray grazing a face at its silhouette
//! tells the truth only from the exact point — asks
//! [`Silhouette::exact_point`].

use std::convert::Infallible;

use geop_core_geometry::{
    intersection::curve_curve_overlaps_and_crossings,
    nurb_curve::{NurbCurve2D, NurbCurve3D},
    nurb_surface::NurbSurface3D,
};
use geop_core_math::{
    geop_error::{GeopError, GeopResult, WithContext},
    plain,
    scalars::Scalar,
    solvers::{
        bisect::sign_changes,
        iterate::{End, Repeat, Step, iterate},
        newton::{Arrive, newton_1d},
        repeat::{self, repeat},
    },
    vector::{Vector2, Vector3},
};
use geop_core_topology::{
    CoedgeId, FaceId, Model,
    contains::face::{PointClassification, face_contains},
};

use crate::hidden_lines::{COARSE_SUBDIVISION, cut_points};
use geop_core_geometry::search::{MAX_NODES, min_subdivision_size};

/// Interior grid lines per direction of the domain, besides its boundary.
const GRID_LINES: usize = 12;
/// Sample intervals along each grid line.
const GRID_SAMPLES: usize = 32;
/// Where the interior grid lines sit within their cells: the golden
/// section, far from the halves and quarters symmetric curves cross at.
const GRID_OFFSET: f64 = 0.381_966_011_250_105;
/// How many times a box within a grid cell is halved, at most, while it
/// could hold a critical point of `g`: how hard the search for small
/// closed loops tries.
const BOX_DEPTH: usize = 1;
/// Newton steps bringing a point onto the silhouette.
const CORRECTOR_ITERATIONS: usize = 24;
/// Bisection steps locating a seed between two samples.
const BISECTIONS: usize = 60;
/// Trace steps per face, at most: how hard a trace tries.
const MAX_TRACE_STEPS: usize = 4000;
/// Times a trace step is halved, at most, before the trace gives up.
const STEP_HALVINGS: usize = 6;
/// Steps a trace takes along the face's size.
const STEPS_PER_FACE: f64 = 24.0;
/// Seed for the containment tests.
const FACE_CONTAINS_SEED: u64 = 0x5114_u64;

/// A piece of a face's silhouette.
#[derive(Clone, Debug)]
pub struct Silhouette {
    pub face: FaceId,
    pub curve: NurbCurve3D,
    /// Exact silhouette points along it, as `(u, v)` on the face's surface:
    /// where [`Silhouette::exact_point`] starts from.
    samples: Vec<Vector2>,
}

/// The silhouette function `g = (S_u x S_v) . d` and its gradient.
fn contour(
    surface: &NurbSurface3D,
    d: &Vector3,
    u: Scalar,
    v: Scalar,
) -> GeopResult<(Scalar, Scalar, Scalar)> {
    let (su, sv) = surface.derivatives(u, v)?;
    let [suu, suv, svv] = surface.second_derivatives(u, v)?;
    let g = su.prod_cross(&sv).prod_dot(d);
    let gu = suu.prod_cross(&sv).add(&su.prod_cross(&suv)).prod_dot(d);
    let gv = suv.prod_cross(&sv).add(&su.prod_cross(&svv)).prod_dot(d);
    Ok((g, gu, gv))
}

/// `g`'s sign at `(u, v)`, from the unit normal — which, unlike `S_u x S_v`,
/// does not vanish at a pole — or `None` where it has no definite one.
fn sign(surface: &NurbSurface3D, d: &Vector3, u: Scalar, v: Scalar) -> Option<bool> {
    surface.normal(u, v).ok()?.prod_dot(d).definite_sign()
}

/// `(u, v)` moved onto the silhouette by Newton steps along `g`'s gradient
/// (the smallest step that zeroes `g` to first order), or why it does not
/// get there. Every iterate is sharpened: each is only the seed of the
/// next. What is returned is the step from the last seed, unsharpened: an
/// enclosure of the silhouette there as wide as `g`'s leaves it — or, where
/// the seeds only cycle round it, the last step as wide as it came out, if
/// `g` could vanish over it. A start outside the domain — a sharpened foot
/// point on its boundary can land a hair past it — is clamped into it
/// first, as every iterate is. The answer is cut to the domain: the
/// silhouette sought is on the patch, and an answer as wide as a noisy `g`
/// leaves it reaches past the patch's ends, where the surface is no
/// surface and its rows' weights are no weights.
fn correct(
    surface: &NurbSurface3D,
    d: &Vector3,
    start: (Scalar, Scalar),
) -> GeopResult<(Scalar, Scalar)> {
    let (u_lo, u_hi) = surface.domain_u();
    let (v_lo, v_hi) = surface.domain_v();
    let seed = (start.0.clamped(u_lo, u_hi), start.1.clamped(v_lo, v_hi));
    let on_domain =
        |(u, v): (Scalar, Scalar)| (u.intersect(u_lo.union(u_hi)), v.intersect(v_lo.union(v_hi)));
    let mut last = None;
    let corrected = iterate(
        seed,
        seed,
        CORRECTOR_ITERATIONS,
        Repeat::Never,
        |&(u, v)| {
            let (u, v) = (u.clamped(u_lo, u_hi), v.clamped(v_lo, v_hi));
            ((u.sharpen(), v.sharpen()), (u, v))
        },
        |&(u, v)| -> GeopResult<Step<(Scalar, Scalar)>> {
            let (g, gu, gv) = contour(surface, d, u, v)?;
            if g.could_be_equal(Scalar::ZERO) {
                return Ok(Step::Done);
            }
            let scale = g.div(gu.mul(gu).add(gv.mul(gv)))?;
            last = Some((u, v, g, gu, gv));
            Ok(Step::Next((u.sub(gu.mul(scale)), v.sub(gv.mul(scale)))))
        },
    )?;
    if corrected.end == End::Done {
        // The sharp seed is where `g` could vanish, one point of where the
        // silhouette could be: the answer is the step from it, unsharpened,
        // which holds the seed and is as wide as `g` there leaves where
        // the silhouette is. Sharp, a silhouette of a padded surface ran
        // `3.6e-12` off the face's edge it lies along, and nothing could
        // see the two were one.
        let (u, v) = corrected.seed;
        let (g, gu, gv) = contour(surface, d, u, v)?;
        let Ok(scale) = g.div(gu.mul(gu).add(gv.mul(gv))) else {
            return Ok((u, v));
        };
        return Ok(on_domain((u.sub(gu.mul(scale)), v.sub(gv.mul(scale)))));
    }
    // No sharp seed lands where `g`, rounded, could vanish — the steps
    // cycle in the last bits around it. The last step's answer, as wide
    // as the step left it, may still hold the silhouette: then it is the
    // point, as precisely as the arithmetic pins it down.
    let (u, v) = on_domain(corrected.answer);
    if contour(surface, d, u, v)?.0.could_be_equal(Scalar::ZERO) {
        return Ok((u, v));
    }
    Err(GeopError::new(format!(
        "the point does not converge onto the silhouette: from {start:?}, the last \
         (u, v, g, g_u, g_v) {last:?}, in the domain {:?} x {:?}",
        surface.domain_u(),
        surface.domain_v()
    )))
}

/// `(u, v)` on the domain's boundary line `fixed` moved onto the silhouette
/// along that line, or `None` if it does not get there.
fn correct_on_boundary(
    surface: &NurbSurface3D,
    d: &Vector3,
    (u, v): (Scalar, Scalar),
    u_fixed: bool,
) -> Option<(Scalar, Scalar)> {
    // The parameter that moves, the bounds it moves in, and the point at it.
    let (moving, (lo, hi)) = if u_fixed {
        (v, surface.domain_v())
    } else {
        (u, surface.domain_u())
    };
    let at = |t: Scalar| if u_fixed { (u, t) } else { (t, v) };
    let corrected = newton_1d(
        moving,
        (lo, hi),
        CORRECTOR_ITERATIONS,
        Arrive::AtZero,
        |t| -> GeopResult<_> {
            let (u, v) = at(t);
            let (g, gu, gv) = contour(surface, d, u, v)?;
            Ok((g, if u_fixed { gv } else { gu }))
        },
    );
    let (t, arrived) = corrected.ok()?;
    arrived.then(|| at(t))
}

/// A piece of an iso line of the domain the seeds are found along: a side
/// of a box (see the module docs).
struct GridLine {
    /// Whether `u` is fixed along it (else `v` is).
    u_fixed: bool,
    value: Scalar,
    /// The other parameter at each sample.
    params: Vec<Scalar>,
    /// Sign changes between samples: the sample indices on either side.
    brackets: Vec<(usize, usize)>,
    /// Which brackets a trace has crossed.
    traced: Vec<bool>,
    /// How far apart its ends are in space: a trace seeded on it steps no
    /// further than half that, so a loop small enough to have been found
    /// only on small boxes is traced in steps its own size.
    scale: f64,
}

impl GridLine {
    fn point(&self, t: Scalar) -> (Scalar, Scalar) {
        if self.u_fixed {
            (self.value, t)
        } else {
            (t, self.value)
        }
    }

    /// Marks traced the bracket a trace crossed at `t` along this line: the
    /// one whose ends, widened by its own width and a 32nd of the line each
    /// way, hold `t`, and the nearest of them if several do. None, if no
    /// bracket is that near — the crossing is off this piece of the line.
    fn mark(&mut self, t: f64) {
        let p = |k: usize| self.params[k].to_f64();
        let line = (p(self.params.len() - 1) - p(0)) / 32.0;
        let near = self
            .brackets
            .iter()
            .enumerate()
            .filter(|(_, (lo, hi))| {
                let reach = p(*hi) - p(*lo) + line;
                p(*lo) - reach <= t && t <= p(*hi) + reach
            })
            .min_by(|(_, a), (_, b)| {
                let distance = |&(lo, hi): &(usize, usize)| (0.5 * (p(lo) + p(hi)) - t).abs();
                distance(a).total_cmp(&distance(b))
            })
            .map(|(i, _)| i);
        if let Some(i) = near {
            self.traced[i] = true;
        }
    }
}

/// The distinct knots of `knots` within `(lo, hi)`, and its ends, in
/// order.
fn breaks(knots: &[Scalar], (lo, hi): (Scalar, Scalar)) -> Vec<Scalar> {
    let mut out = vec![lo];
    for &k in knots {
        if k.definitely_greater(*out.last().expect("a first break")) && k.definitely_less(hi) {
            out.push(k);
        }
    }
    out.push(hi);
    out
}

/// The middle of `(a, b)`, sharp: where to split is a free choice.
fn middle(a: Scalar, b: Scalar) -> Scalar {
    Scalar::fraction(a, b, Scalar::from_f64(0.5)).sharpen()
}

/// The iso lines seeds are found on (see the module docs): the grid's,
/// and through every box in a cell of it that could hold a critical point
/// of `g` where `g` could vanish, a piece of a line across that cell.
fn grid(surface: &NurbSurface3D, d: &Vector3) -> Vec<GridLine> {
    let mut fractions = vec![0.0];
    fractions.extend((0..GRID_LINES).map(|k| (k as f64 + GRID_OFFSET) / GRID_LINES as f64));
    fractions.push(1.0);
    let at = |(lo, hi): (Scalar, Scalar), f: f64| match f {
        0.0 => lo,
        1.0 => hi,
        _ => Scalar::fraction(lo, hi, Scalar::from_f64(f)),
    };
    let (du, dv) = (surface.domain_u(), surface.domain_v());
    let us: Vec<Scalar> = fractions.iter().map(|&f| at(du, f)).collect();
    let vs: Vec<Scalar> = fractions.iter().map(|&f| at(dv, f)).collect();
    let mut lines = Vec::new();
    for &u in &us {
        lines.push(sampled_line(surface, d, true, u, dv));
    }
    for &v in &vs {
        lines.push(sampled_line(surface, d, false, v, du));
    }
    for u in us.windows(2) {
        for v in vs.windows(2) {
            for (at, size) in loop_seats(surface, d, (u[0], u[1]), (v[0], v[1])) {
                let mut line = sampled_line(surface, d, true, at, (v[0], v[1]));
                line.scale = line.scale.min(size);
                lines.push(line);
            }
        }
    }
    lines
}

/// Where in the cell `cu` x `cv` of the grid a loop of the silhouette that
/// crosses none of its lines could lie: the `u` of the middle of every box
/// of it, halved to [`BOX_DEPTH`], that could hold a critical point of
/// `g`, and the box's size in space — the size of a loop found there. A loop of `g = 0` encloses one (an extremum of `g`), so a line
/// through each crosses every loop in the cell — none, if `g` cannot
/// vanish in the cell at all. Each box is within a knot span, where the
/// enclosures are the span's own.
fn loop_seats(
    surface: &NurbSurface3D,
    d: &Vector3,
    cu: (Scalar, Scalar),
    cv: (Scalar, Scalar),
) -> Vec<(Scalar, f64)> {
    let zero = Scalar::ZERO;
    let us = breaks(&surface.knot_vector_u, cu);
    let vs = breaks(&surface.knot_vector_v, cv);
    let mut stack = Vec::new();
    let mut vanishes = false;
    for u in us.windows(2) {
        for v in vs.windows(2) {
            match contour(surface, d, u[0].union(u[1]), v[0].union(v[1])) {
                Ok((g, _, _)) if !g.could_be_equal(zero) => {}
                _ => vanishes = true,
            }
            stack.push(((u[0], u[1]), (v[0], v[1]), 0));
        }
    }
    let mut seats: Vec<(Scalar, f64)> = Vec::new();
    if !vanishes {
        return seats;
    }
    while let Some(((u0, u1), (v0, v1), depth)) = stack.pop() {
        let critical = match contour(surface, d, u0.union(u1), v0.union(v1)) {
            Ok((_, gu, gv)) => gu.could_be_equal(zero) && gv.could_be_equal(zero),
            Err(_) => true,
        };
        if !critical {
            continue;
        }
        let (um, vm) = (middle(u0, u1), middle(v0, v1));
        if depth >= BOX_DEPTH {
            if !seats.iter().any(|s| s.0.could_be_equal(um)) {
                let corner =
                    |u: Scalar, v: Scalar| surface.evaluate(u, v).map(|p| p.to_f64_array());
                let size = match (corner(u0, v0), corner(u1, v1)) {
                    (Ok(a), Ok(b)) => plain::distance(a, b),
                    _ => f64::INFINITY,
                };
                seats.push((um, size));
            }
            continue;
        }
        for u in [(u0, um), (um, u1)] {
            for v in [(v0, vm), (vm, v1)] {
                stack.push((u, v, depth + 1));
            }
        }
    }
    seats
}

/// The piece `range` of the iso line where `u` (if `u_fixed`, else `v`) is
/// `value`, sampled for `g`'s sign at [`GRID_SAMPLES`] intervals.
fn sampled_line(
    surface: &NurbSurface3D,
    d: &Vector3,
    u_fixed: bool,
    value: Scalar,
    (lo, hi): (Scalar, Scalar),
) -> GridLine {
    let params: Vec<Scalar> = (0..=GRID_SAMPLES)
        .map(|k| match k {
            0 => lo,
            k if k == GRID_SAMPLES => hi,
            k => Scalar::fraction(lo, hi, Scalar::from_f64(k as f64 / GRID_SAMPLES as f64)),
        })
        .collect();
    let mut line = GridLine {
        u_fixed,
        value,
        params,
        brackets: Vec::new(),
        traced: Vec::new(),
        scale: f64::INFINITY,
    };
    let at = |t: Scalar| {
        let (u, v) = line.point(t);
        surface.evaluate(u, v).map(|p| p.to_f64_array())
    };
    if let (Ok(a), Ok(b)) = (at(lo), at(hi)) {
        line.scale = plain::distance(a, b);
    }
    let brackets: Vec<(usize, usize)> = sign_changes(0..line.params.len(), |k| {
        let (u, v) = line.point(line.params[k]);
        Ok::<_, Infallible>(sign(surface, d, u, v))
    })
    .flatten()
    .map(|(j, k, _)| (j, k))
    .collect();
    line.traced = vec![false; brackets.len()];
    line.brackets = brackets;
    line
}

/// The seed in `line`'s bracket `(lo, hi)`: bisected on `g`'s sign, then
/// moved onto the silhouette.
fn seed(
    surface: &NurbSurface3D,
    d: &Vector3,
    line: &GridLine,
    (lo, hi): (usize, usize),
) -> Option<(Scalar, Scalar)> {
    let sign_at = |t| -> GeopResult<_> {
        let (u, v) = line.point(t);
        Ok(sign(surface, d, u, v))
    };
    let bracket = sign_changes([line.params[lo], line.params[hi]], sign_at)
        .bisected(BISECTIONS, false)
        .next()?
        .ok()?;
    // Where the sign cannot be told, the silhouette is there.
    let (a, b) = match bracket.undecided {
        Some(mid) => (mid, mid),
        None => (bracket.lo, bracket.hi),
    };
    let t = a.add(b).div(Scalar::TWO).ok()?.sharpen();
    correct(surface, d, line.point(t)).ok()
}

/// Marks traced every bracket the step from `a` to `b` along the
/// silhouette crosses or ends on.
///
/// Which brackets a trace has passed is only bookkeeping — it decides which
/// seeds are left to trace — but a bracket missed is a silhouette traced a
/// second time, a copy no search on the paper can tell from the first. Two
/// ways a step passes a grid line without its ends lying on either side:
///
/// - It ends on it. A trace that leaves the domain lands on its boundary,
///   inside the enclosure of the bound but not necessarily on the side of
///   its midpoint the step came from: asked of the intervals.
/// - It dips across and back. A silhouette whose extremum lies a hair past a
///   grid line crosses the line twice between the ends of one step. So the
///   step is taken as the cubic leaving `a` and arriving at `b` along the
///   silhouette's tangents there (Hermite), which follows the silhouette to
///   fourth order in the step, and its crossings are the line's.
fn mark_crossings(
    surface: &NurbSurface3D,
    d: &Vector3,
    lines: &mut [GridLine],
    a: (Scalar, Scalar),
    b: (Scalar, Scalar),
) {
    let (au, av, bu, bv) = (a.0.to_f64(), a.1.to_f64(), b.0.to_f64(), b.1.to_f64());
    let chord = [bu - au, bv - av];
    let length = chord[0].hypot(chord[1]);
    // The tangent at `p`, as long as the chord and along it; the chord
    // itself where there is none.
    let tangent = |(u, v): (Scalar, Scalar)| -> [f64; 2] {
        let Ok((_, gu, gv)) = contour(surface, d, u, v) else {
            return chord;
        };
        let t = [-gv.to_f64(), gu.to_f64()];
        let norm = t[0].hypot(t[1]);
        let along = t[0] * chord[0] + t[1] * chord[1];
        if norm == 0.0 || along == 0.0 || !norm.is_finite() {
            return chord;
        }
        let scale = length / norm * along.signum();
        [t[0] * scale, t[1] * scale]
    };
    let (ta, tb) = (tangent(a), tangent(b));
    // The step's cubic in coordinate `k` (0 for `u`, 1 for `v`), as power
    // coefficients in `s` over `[0, 1]`.
    let cubic = |k: usize| -> [f64; 4] {
        let (x0, x1) = if k == 0 { (au, bu) } else { (av, bv) };
        [
            x0,
            ta[k],
            3.0 * (x1 - x0) - 2.0 * ta[k] - tb[k],
            2.0 * (x0 - x1) + ta[k] + tb[k],
        ]
    };
    let at = |c: &[f64; 4], s: f64| ((c[3] * s + c[2]) * s + c[1]) * s + c[0];
    for line in lines.iter_mut() {
        let (k, other) = if line.u_fixed { (0, 1) } else { (1, 0) };
        for (end, p) in [(a, (au, av)), (b, (bu, bv))] {
            let (x, y) = if line.u_fixed {
                (end.0, p.1)
            } else {
                (end.1, p.0)
            };
            if x.could_be_equal(line.value) {
                line.mark(y);
            }
        }
        let mut here = cubic(k);
        here[0] -= line.value.to_f64();
        let across = cubic(other);
        // The cubic's turning points split `[0, 1]` into pieces it is
        // monotone along; a piece whose ends lie on either side of the line
        // crosses it once, bisected.
        let (qa, qb, qc) = (3.0 * here[3], 2.0 * here[2], here[1]);
        let mut ends = vec![0.0, 1.0];
        if qa == 0.0 {
            if qb != 0.0 {
                ends.push(-qc / qb);
            }
        } else {
            let disc = qb * qb - 4.0 * qa * qc;
            if disc >= 0.0 {
                let root = disc.sqrt();
                ends.extend([(-qb - root) / (2.0 * qa), (-qb + root) / (2.0 * qa)]);
            }
        }
        ends.retain(|s| (0.0..=1.0).contains(s));
        ends.sort_by(f64::total_cmp);
        // Where the cubic vanishes, it has no side.
        let side = |s| {
            let x = at(&here, s);
            Ok::<_, GeopError>((x != 0.0 && !x.is_nan()).then_some(x > 0.0))
        };
        for bracket in sign_changes(ends, side)
            .bisected(BISECTIONS, false)
            .flatten()
        {
            let s = bracket.undecided.unwrap_or(0.5 * (bracket.lo + bracket.hi));
            line.mark(at(&across, s));
        }
    }
}

/// How a march along the silhouette ended.
enum Ending {
    /// It came back to where it started.
    Closed,
    /// It left the domain, or could not go on.
    Open,
}

/// Marches from `start` along the silhouette, the way `forward` says,
/// pushing each point onto `points` (not `start` itself).
#[allow(clippy::too_many_arguments)]
fn march(
    surface: &NurbSurface3D,
    d: &Vector3,
    start: (Scalar, Scalar),
    forward: bool,
    step: Scalar,
    lines: &mut [GridLine],
    points: &mut Vec<(Scalar, Scalar)>,
) -> GeopResult<Ending> {
    let (u_lo, u_hi) = surface.domain_u();
    let (v_lo, v_hi) = surface.domain_v();
    let start_point = surface.evaluate(start.0, start.1)?;
    let mut current = start;
    let mut travelled = Scalar::ZERO;
    repeat(MAX_TRACE_STEPS, || -> GeopResult<_> {
        let (u, v) = current;
        let Ok((_, gu, gv)) = contour(surface, d, u, v) else {
            return Ok(repeat::Step::Done(Ending::Open));
        };
        let (mut tu, mut tv) = (gv.neg(), gu);
        if !forward {
            (tu, tv) = (tu.neg(), tv.neg());
        }
        let (su, sv) = surface.derivatives(u, v)?;
        let speed = su.prod_scalar(tu).add(&sv.prod_scalar(tv)).norm();
        // A step whose corrector fails, or lands further off than twice
        // the step, is retried at half the length: how far to step is a
        // free choice, and a tight turn of the silhouette needs short ones.
        let from = surface.evaluate(u, v)?;
        let mut tried = step;
        let mut halvings = 0;
        let next = loop {
            let Ok(scale) = tried.div(speed) else {
                return Ok(repeat::Step::Done(Ending::Open));
            };
            let (du, dv) = (tu.mul(scale), tv.mul(scale));
            let (pu, pv) = (u.add(du).sharpen(), v.add(dv).sharpen());
            let leaves_u = pu.definitely_less(u_lo) || pu.definitely_greater(u_hi);
            let leaves_v = pv.definitely_less(v_lo) || pv.definitely_greater(v_hi);
            if leaves_u || leaves_v {
                // The boundary the step reaches first, and the fraction of the
                // step to it: the point there, then onto the silhouette along
                // that boundary.
                let to = |x: Scalar, dx: Scalar, lo: Scalar, hi: Scalar| -> (f64, Scalar) {
                    let bound = if dx.to_f64() < 0.0 { lo } else { hi };
                    ((bound.sub(x).to_f64() / dx.to_f64()).clamp(0.0, 1.0), bound)
                };
                let (fu, bound_u) = to(u, du, u_lo, u_hi);
                let (fv, bound_v) = to(v, dv, v_lo, v_hi);
                let u_first = leaves_u && (!leaves_v || fu <= fv);
                // On the boundary, exactly: the step's own arithmetic lands a
                // rounding off it — a hair inside, on no grid line — and a
                // trace ending there does not mark the bracket the boundary
                // line has where it ends, which then seeds the same silhouette
                // again.
                let landed = if u_first {
                    let f = Scalar::from_f64(fu);
                    (bound_u, v.add(dv.mul(f)).sharpen().clamped(v_lo, v_hi))
                } else {
                    let f = Scalar::from_f64(fv);
                    (u.add(du.mul(f)).sharpen().clamped(u_lo, u_hi), bound_v)
                };
                if let Some(end) = correct_on_boundary(surface, d, landed, u_first) {
                    mark_crossings(surface, d, lines, current, end);
                    points.push(end);
                }
                return Ok(repeat::Step::Done(Ending::Open));
            }
            let landed = correct(surface, d, (pu, pv)).ok().filter(|next| {
                surface
                    .evaluate(next.0, next.1)
                    .is_ok_and(|p| p.sub(&from).norm().definitely_less(tried.mul(Scalar::TWO)))
            });
            match landed {
                Some(next) => break next,
                None if halvings < STEP_HALVINGS => {
                    tried = tried.div(Scalar::TWO)?;
                    halvings += 1;
                }
                None => return Ok(repeat::Step::Done(Ending::Open)),
            }
        };
        mark_crossings(surface, d, lines, current, next);
        points.push(next);
        travelled = travelled.add(tried);
        let here = surface.evaluate(next.0, next.1)?;
        if travelled.definitely_greater(step.mul(Scalar::TWO))
            && here.sub(&start_point).norm().definitely_less(step)
        {
            mark_crossings(surface, d, lines, next, start);
            points.push(start);
            return Ok(repeat::Step::Done(Ending::Closed));
        }
        current = next;
        Ok(repeat::Step::Again)
    })?
    .done_or(|steps| {
        GeopError::new(format!(
            "the silhouette trace from {start:?} did not end within {steps} steps of {step:?}"
        ))
    })
}

/// The silhouette through `seed`, traced both ways, as `(u, v)` points.
fn trace(
    surface: &NurbSurface3D,
    d: &Vector3,
    seed: (Scalar, Scalar),
    step: Scalar,
    lines: &mut [GridLine],
) -> GeopResult<Vec<(Scalar, Scalar)>> {
    let mut forward = Vec::new();
    let points = if let Ending::Closed = march(surface, d, seed, true, step, lines, &mut forward)? {
        let mut points = vec![seed];
        points.extend(forward);
        points
    } else {
        let mut backward = Vec::new();
        march(surface, d, seed, false, step, lines, &mut backward)?;
        backward.reverse();
        backward.push(seed);
        backward.extend(forward);
        backward
    };
    thin(surface, points, step)
}

/// `points` without those crowding the one before: where a curve is
/// sampled is a free choice, and samples a sliver apart — a step that
/// reached the domain's boundary all but, then the boundary itself — make
/// the fit through them ill posed. The ends stay.
fn thin(
    surface: &NurbSurface3D,
    points: Vec<(Scalar, Scalar)>,
    step: Scalar,
) -> GeopResult<Vec<(Scalar, Scalar)>> {
    let close = step.div(Scalar::from_f64(8.0))?;
    let at = |p: &(Scalar, Scalar)| surface.evaluate(p.0, p.1);
    let count = points.len();
    let mut kept: Vec<(Scalar, Scalar)> = Vec::with_capacity(count);
    for (k, p) in points.into_iter().enumerate() {
        if let Some(last) = kept.last()
            && at(last)?.sub(&at(&p)?).norm().definitely_less(close)
        {
            if k + 1 == count && kept.len() > 1 {
                kept.pop();
            } else {
                continue;
            }
        }
        kept.push(p);
    }
    Ok(kept)
}

/// The size of `surface`: the diagonal of its control points' box.
fn extent(surface: &NurbSurface3D) -> Scalar {
    let b = surface.aabb();
    Scalar::from_f64(plain::diagonal(
        b.map(|x| x.lower().to_f64()),
        b.map(|x| x.upper().to_f64()),
    ))
}

/// Whether `surface` is swept straight along `d` in one of its parameter
/// directions — an extruded or turned cylinder seen end on — so that its
/// normal is perpendicular to `d` everywhere and the whole face is seen
/// edge on: then its outline is its edges, and it has no silhouette curve.
///
/// Swept along `d` to the sheet's resolution ([`COARSE_SUBDIVISION`]) is
/// swept along it: every control point of a row, of one weight, lies on
/// the paper within it of the row's first, so the face lies within it of
/// the curve its first points span, which its edges draw — and so do its
/// silhouettes. A cylinder in a part placed by a pose rounded in its last
/// bits is seen `4e-14` off its axis: its silhouette is a generator that
/// projects to a point, and `g`, a rounding off zero all over the face,
/// put the corrector's answer for it anywhere round the arc, `150`
/// wide in a domain of `1`.
fn swept_along(surface: &NurbSurface3D, d: &Vector3) -> bool {
    let (nu, nv) = (surface.num_u, surface.num_v);
    let point = |i: usize, j: usize| surface.control_points[i * nv + j];
    let cartesian = |q: geop_core_math::vector::Vector4| -> Option<Vector3> {
        let w = q[3];
        Some(Vector3::from_array([
            q[0].div(w).ok()?,
            q[1].div(w).ok()?,
            q[2].div(w).ok()?,
        ]))
    };
    let along = |a: geop_core_math::vector::Vector4, b: geop_core_math::vector::Vector4| {
        let (Some(pa), Some(pb)) = (cartesian(a), cartesian(b)) else {
            return false;
        };
        a[3].could_be_equal(b[3])
            && !(pb.sub(&pa).prod_cross(d).norm())
                .definitely_greater(Scalar::from_f64(COARSE_SUBDIVISION))
    };
    let along_v = (0..nu).all(|i| (1..nv).all(|j| along(point(i, 0), point(i, j))));
    let along_u = (0..nv).all(|j| (1..nu).all(|i| along(point(0, j), point(i, j))));
    along_u || along_v
}

impl Silhouette {
    /// The point of this silhouette near `near`, exactly on it: `near`
    /// projected onto the surface from the nearest exact sample, then moved
    /// onto the silhouette.
    pub fn exact_point(&self, model: &Model, d: &Vector3, near: &Vector3) -> GeopResult<Vector3> {
        let surface = &model.get_face(self.face)?.surface;
        let ctx = |e: GeopError| {
            e.with_context(format!(
                "Silhouette::exact_point(face={}, near={near:?})",
                self.face
            ))
        };
        let distance = |uv: &Vector2| -> f64 {
            surface
                .evaluate(uv[0], uv[1])
                .map(|p| p.sub(near).norm_sq().to_f64())
                .unwrap_or(f64::INFINITY)
        };
        let seed = self
            .samples
            .iter()
            .min_by(|a, b| distance(a).total_cmp(&distance(b)))
            .ok_or_else(|| GeopError::new("a silhouette without samples"))
            .with_context(&ctx)?;
        let (u, v) = surface
            .project(*near, seed[0], seed[1], CORRECTOR_ITERATIONS)
            .with_context(&ctx)?;
        let (u, v) = correct(surface, d, (u.sharpen(), v.sharpen())).with_context(&ctx)?;
        surface.evaluate(u, v).with_context(&ctx)
    }
}

/// Samples along an edge looking for where it crosses a silhouette.
const EDGE_SAMPLES: usize = 32;

/// Where the edge `curve`, bounding the face `face_id` along the coedge
/// `coedge_id`, crosses the face's silhouette seen along `d`: the
/// parameters where the face's normal there turns perpendicular to `d`.
///
/// Projected, an edge touches the silhouette of a face it bounds
/// tangentially where it crosses it — a contact no search on the paper can
/// isolate, so it is found here, in space: `g`'s sign sampled along the
/// edge, and each change bisected. Only isolated crossings: an edge along
/// which the face is seen edge on throughout has none.
pub fn edge_crossings(
    model: &Model,
    face_id: FaceId,
    coedge_id: CoedgeId,
    curve: &NurbCurve3D,
    d: &Vector3,
) -> GeopResult<Vec<Scalar>> {
    let surface = &model.get_face(face_id)?.surface;
    if surface.as_plane()?.is_some() || swept_along(surface, d) {
        return Ok(Vec::new());
    }
    let coedge = model.get_coedge(coedge_id)?;
    let (t0, t1) = curve.domain();
    let (s0, s1) = coedge.pcurve.domain();
    let edge_domain = model.get_edge(coedge.edge()?)?.curve.domain();
    // Where on the surface the edge's point at `t` is: projected from the
    // pcurve's point as far along it — a seed, a free choice.
    let foot = |t: Scalar| -> Option<(Scalar, Scalar)> {
        let f = t
            .sub(edge_domain.0)
            .div(edge_domain.1.sub(edge_domain.0))
            .ok()?;
        let f = match coedge.sense {
            geop_core_topology::Sense::Forward => f,
            geop_core_topology::Sense::Reversed => Scalar::ONE.sub(f),
        };
        let seed = coedge.pcurve.evaluate(Scalar::fraction(s0, s1, f)).ok()?;
        let point = curve.evaluate(t).ok()?;
        surface
            .project(point, seed[0], seed[1], CORRECTOR_ITERATIONS)
            .ok()
    };
    let sign_at = |t: Scalar| -> Option<bool> {
        let (u, v) = foot(t)?;
        sign(surface, d, u, v)
    };
    let params = (0..=EDGE_SAMPLES).map(|k| match k {
        0 => t0,
        k if k == EDGE_SAMPLES => t1,
        k => Scalar::fraction(t0, t1, Scalar::from_f64(k as f64 / EDGE_SAMPLES as f64)),
    });
    sign_changes(params, |t| -> GeopResult<_> { Ok(sign_at(t)) })
        .bisected(BISECTIONS, false)
        .map(|bracket| {
            let bracket = bracket?;
            Ok(bracket
                .undecided
                .unwrap_or_else(|| bracket.lo.union(bracket.hi)))
        })
        .collect()
}

/// The silhouette of `face_id` seen along the unit vector `d`, as curves on
/// the face: none for a plane, or for a face seen edge on throughout.
pub fn face_silhouettes(
    model: &Model,
    face_id: FaceId,
    d: &Vector3,
) -> GeopResult<Vec<Silhouette>> {
    let ctx = |e: GeopError| e.with_context(format!("face_silhouettes(face={face_id}, d={d:?})"));
    let surface = &model.get_face(face_id).with_context(&ctx)?.surface;
    if surface.as_plane().with_context(&ctx)?.is_some() || swept_along(surface, d) {
        return Ok(Vec::new());
    }
    let step = extent(surface)
        .div(Scalar::from_f64(STEPS_PER_FACE))
        .with_context(&ctx)?;
    let mut lines = grid(surface, d);
    let mut traces = Vec::new();
    for l in 0..lines.len() {
        for b in 0..lines[l].brackets.len() {
            if lines[l].traced[b] {
                continue;
            }
            lines[l].traced[b] = true;
            let Some(start) = seed(surface, d, &lines[l], lines[l].brackets[b]) else {
                continue;
            };
            let half = Scalar::from_f64(lines[l].scale / 2.0);
            let step = if half.definitely_less(step) {
                half
            } else {
                step
            };
            let traced = trace(surface, d, start, step, &mut lines).map_err(|e| {
                e.with_context(format!(
                    "tracing the silhouette of {surface:?} from {start:?} in steps of {step:?}"
                ))
            });
            traces.push(traced.with_context(&ctx)?);
        }
    }
    let mut silhouettes = Vec::new();
    for points in traces {
        let trimmed = trim(model, face_id, d, &points).map_err(|e| {
            e.with_context(format!("trimming the silhouette traced through {points:?}"))
        });
        silhouettes.extend(trimmed.with_context(&ctx)?);
    }
    Ok(silhouettes)
}

/// The pieces of the traced silhouette `points` inside the face, as 3-D
/// curves through exact silhouette points.
fn trim(
    model: &Model,
    face_id: FaceId,
    d: &Vector3,
    points: &[(Scalar, Scalar)],
) -> GeopResult<Vec<Silhouette>> {
    let surface = &model.get_face(face_id)?.surface;
    let mut uv: Vec<Vector2> = Vec::with_capacity(points.len());
    for &(u, v) in points {
        let p = Vector2::from_array([u, v]);
        if uv.last().is_none_or(|q: &Vector2| !q.could_be_equal(&p)) {
            uv.push(p);
        }
    }
    if uv.len() < 2 {
        return Ok(Vec::new());
    }
    let fit = NurbCurve2D::interpolate(&uv, (uv.len() - 1).min(3))
        .map_err(|e| e.with_context(format!("fitting the traced silhouette through {uv:?}")))?;
    let (t0, t1) = fit.domain();

    // Where the trace crosses or runs along the face's trim. A trace that
    // left the domain ends exactly on its boundary, so the trim meets it
    // there, at an end it already has: only a meeting definitely inside the
    // trace is a cut. Cutting at the midpoint of the box found at an end
    // left a sliver beyond it, which at a cone's apex lay a hair inside the
    // face and drew a silhouette piece no longer than the rounding.
    //
    // A face whose boundary is its silhouette by design — a half cylinder
    // seen from above — is seen a hair off that in a placed part turned by
    // a pose rounded in its last bits: the trace runs `1e-14` beside the
    // trim all along it, and the exact search could only cover that with
    // pieces as small. Should it fail, the two are compared again to the
    // sheet's resolution, as the lines of a view are (see
    // `hidden_lines::meetings`): a stretch running that close is one.
    let coarse = Scalar::from_f64(COARSE_SUBDIVISION);
    let mut meetings = Vec::new();
    for coedge_id in model.iterate_face_coedges(face_id) {
        let pcurve = &model.get_coedge(coedge_id)?.pcurve;
        let (overlaps, crossings) = curve_curve_overlaps_and_crossings(
            &fit,
            pcurve,
            Scalar::ZERO,
            MAX_NODES,
            min_subdivision_size(),
        )
        .or_else(|error| {
            curve_curve_overlaps_and_crossings(&fit, pcurve, coarse, MAX_NODES, coarse).map_err(
                |e| {
                    e.with_context(format!(
                        "searched again to the sheet's resolution, after: {error}"
                    ))
                },
            )
        })
        .map_err(|e| {
            e.with_context(format!(
                "where the traced silhouette {fit:?} meets the pcurve {pcurve:?} of coedge \
                 {coedge_id}"
            ))
        })?;
        meetings.extend(crossings.iter().map(|(s, _)| *s));
        meetings.extend(overlaps.iter().flat_map(|o| [o.start.t, o.end.t]));
    }
    let mut cuts = vec![t0];
    cuts.extend(cut_points(meetings, (t0, t1)));
    cuts.push(t1);

    let mut silhouettes = Vec::new();
    for pair in cuts.windows(2) {
        let (a, b) = (pair[0], pair[1]);
        let mid = a.add(b).div(Scalar::TWO)?.sharpen();
        let at = fit.evaluate(mid)?;
        let inside = face_contains(
            model,
            face_id,
            at[0],
            at[1],
            MAX_NODES,
            min_subdivision_size(),
            FACE_CONTAINS_SEED,
        )?;
        if inside != PointClassification::Inside {
            continue;
        }
        // Exact points along the piece, as many as the trace had there and
        // at least enough for a cubic.
        let share = b.sub(a).div(t1.sub(t0))?.to_f64();
        let count = ((share * uv.len() as f64).ceil() as usize).max(4);
        let mut samples = Vec::with_capacity(count + 1);
        for k in 0..=count {
            let t = match k {
                0 => a,
                k if k == count => b,
                k => Scalar::fraction(a, b, Scalar::from_f64(k as f64 / count as f64)),
            };
            let q = fit.evaluate(t)?;
            let (u, v) =
                correct(surface, d, (q[0].sharpen(), q[1].sharpen())).unwrap_or((q[0], q[1]));
            samples.push(Vector2::from_array([u, v]));
        }
        let mut stations: Vec<(Vector3, Vector2)> = Vec::with_capacity(samples.len());
        for s in &samples {
            let p = surface.evaluate(s[0], s[1])?;
            if stations
                .last()
                .is_none_or(|(q, _): &(Vector3, Vector2)| !q.could_be_equal(&p))
            {
                stations.push((p, *s));
            }
        }
        if stations.len() < 2 {
            continue;
        }
        // A cubic through the stations, widened to enclose the silhouette
        // between them too: two faces whose silhouettes are one line on
        // the paper — mirror images seen along the mirror's normal — are
        // then seen to run along each other, where two interpolants of the
        // one line, each drifting from it its own way, cross and touch
        // along all of it and no search could cut them.
        let locate = |(_, from): &(Vector3, Vector2),
                      (_, to): &(Vector3, Vector2),
                      frac: Scalar|
         -> GeopResult<(Vector3, Vector2)> {
            let blend = |k: usize| Scalar::fraction(from[k], to[k], frac).sharpen();
            let (u, v) = correct(surface, d, (blend(0), blend(1)))?;
            Ok((surface.evaluate(u, v)?, Vector2::from_array([u, v])))
        };
        let ctx = |e: GeopError| {
            e.with_context(format!(
                "fitting the silhouette of face {face_id} through {stations:?}"
            ))
        };
        let curve = NurbCurve3D::fit_enclosing(&stations, locate, &mut None, &ctx)?;
        silhouettes.push(Silhouette {
            face: face_id,
            curve,
            samples,
        });
    }
    Ok(silhouettes)
}

#[cfg(test)]
mod tests {
    use super::*;
    use geop_core_math::vector::Vector4;
    use geop_core_topology::build::BodySpec;

    /// A bicubic sheet over `[0, 1]²` of 60 knot spans each way, rising
    /// gently along `x` (`z = x / 20`), with one control point near
    /// `(0.48, 0.48)` raised into a bump: a sheet whose only silhouette,
    /// seen along `x`, is a small loop round the bump's flank, where its
    /// slope cancels the rise — a loop that fits between two iso lines a
    /// twelfth of the domain apart.
    fn bumpy_sheet() -> (Model, FaceId) {
        let spans = 60;
        let n = spans + 3;
        let s = Scalar::from_f64;
        let mut knots = vec![Scalar::ZERO; 4];
        knots.extend((1..spans).map(|k| s(k as f64 / spans as f64)));
        knots.extend([Scalar::ONE; 4]);
        let mut control_points = Vec::with_capacity(n * n);
        for i in 0..n {
            for j in 0..n {
                let (x, y) = (i as f64 / (n - 1) as f64, j as f64 / (n - 1) as f64);
                let bump = if (i, j) == (30, 30) { 0.08 } else { 0.0 };
                control_points.push(Vector4::from_array([
                    s(x),
                    s(y),
                    s(x / 20.0 + bump),
                    s(1.0),
                ]));
            }
        }
        let surface = NurbSurface3D::try_new(3, 3, control_points, knots.clone(), knots).unwrap();
        let mut model = Model::new();
        let built = model
            .build_body(BodySpec::sheet_on(surface).unwrap())
            .unwrap();
        (model, built.faces[0])
    }

    /// The bump's loop is found, though it crosses no iso line a sampling
    /// grid would look along, and traced once.
    #[test]
    fn a_small_loop_between_grid_lines_is_found() {
        let (model, face) = bumpy_sheet();
        let d = Vector3::from_array([Scalar::ONE, Scalar::ZERO, Scalar::ZERO]);
        let pieces = face_silhouettes(&model, face, &d).unwrap();
        assert_eq!(pieces.len(), 1, "{pieces:?}");
        let curve = &pieces[0].curve;
        let (start, end) = (curve.start().unwrap(), curve.end().unwrap());
        assert!(start.could_be_equal(&end), "{start:?} {end:?}");
    }
}