ogeom_mesh/

discretize.rs

//! Turning a curve into a polyline within a stated deflection.
//!
//! The first half of tessellation, and the half that decides whether the second
//! half can succeed: a face's triangulation is built from the polylines of its
//! boundary edges, so two faces meeting along an edge produce a watertight join
//! only if they discretize that edge to the *same* points. That is why
//! discretization is a property of the edge rather than of the face, and why it
//! is a separate step from triangulating the face itself.
//!
//! # What deflection means
//!
//! The chord deflection is the greatest distance between the polyline and the
//! curve it approximates. It is a length, so it scales with the model, and it
//! is the number a caller actually has an opinion about: "no visible error at
//! this zoom", "within machining tolerance".
//!
//! The angular deflection bounds how far the tangent may turn across one
//! segment. Without it a nearly straight curve gets two points and a
//! near-circular one gets far too few near its flattest part. Chord error
//! alone does not notice a long, gently curving span.

use ogeom_core::{OgeomResult, Tolerances, ogeom_bail};
use ogeom_geom::{Curve, Curve3d};
use ogeom_math::Point;

/// How finely to approximate a curve.
///
/// Only the chord and the angle subdivide; there is no cap on a segment's or
/// a triangle's length. A straight edge is drawn as its two ends and a planar
/// face as the triangulation of its outline, however long either is. A
/// curved edge or face is subdivided until it meets both tolerances. Inside
/// a curved face the grid the tolerances chose may also be spread so that
/// no cell is many times longer than it is wide: a ratio of that grid's own
/// spacing, which keeps the triangles between its rows from bridging
/// columns, not a length.
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct Deflection {
    /// The greatest allowed distance from the polyline to the curve.
    pub chord: f64,
    /// The greatest allowed turn in the tangent across one segment, in radians.
    pub angular: f64,
    /// A floor on the number of segments for a curve that bends.
    ///
    /// A closed curve needs at least three to enclose anything, and a nearly
    /// straight arc of one would otherwise collapse to a chord that misses the
    /// bulge entirely: the midpoint test is a sample, and one sample can be
    /// placed exactly where the curve happens to cross its own chord.
    ///
    /// It does *not* apply to a straight curve, which is exactly represented by
    /// its endpoints. Splitting a line adds points that no tolerance asked for,
    /// and every one of them becomes a vertex in every face the edge bounds.
    pub min_segments: usize,
    /// A ceiling, so a pathological curve cannot exhaust memory.
    ///
    /// Reaching it is reported rather than passed off as success; see
    /// [`Polyline::deflection_met`].
    pub max_segments: usize,
}

impl Default for Deflection {
    fn default() -> Self {
        Self {
            // A tenth of a millimetre at unit scale: invisible on screen. On
            // a small curved part the angular limit below rules instead, and
            // a mesh at this default can be a part in a hundred off in
            // volume. Mass properties integrate a face on its exact surface
            // wherever its pcurves allow, and mesh only what they cannot.
            chord: 1e-1,
            // Half a radian (twenty-eight degrees, a circle in thirteen
            // segments): coarse enough for a viewer that lets the user
            // tighten it. The interior of a face is held to it as its edges
            // are, so a tighter default is paid for on every curved face.
            angular: 0.5,
            min_segments: 2,
            max_segments: 4096,
        }
    }
}

impl Deflection {
    /// A deflection with a given chord tolerance and the default angular one.
    ///
    /// # Errors
    ///
    /// [`OgeomError::Construction`](ogeom_core::OgeomError::Construction) if `chord` is
    /// not finite and positive.
    pub fn with_chord(chord: f64) -> OgeomResult<Self> {
        if !chord.is_finite() || chord <= 0.0 {
            ogeom_bail!(
                Construction,
                "chord deflection {chord} must be finite and positive"
            );
        }
        Ok(Self {
            chord,
            ..Self::default()
        })
    }

    /// A deflection held to `chord` and `angular` alone, with the default
    /// segment bounds.
    ///
    /// The tessellation a measurement is taken against: flat faces stay as
    /// few triangles as their outlines allow, and curved ones are drawn to
    /// the chord. Every deflection subdivides this way (see [`Deflection`]);
    /// this names both tolerances and checks them.
    ///
    /// # Errors
    ///
    /// [`OgeomError::Construction`](ogeom_core::OgeomError::Construction) if `chord` or
    /// `angular` is not finite and positive.
    pub fn curvature_only(chord: f64, angular: f64) -> OgeomResult<Self> {
        let deflection = Self {
            chord,
            angular,
            ..Self::default()
        };
        deflection.validate()?;
        Ok(deflection)
    }

    /// A deflection scaled to a model of the given size.
    ///
    /// Relative tolerances are what a caller usually means: "a thousandth of the
    /// part" is a statement that survives the part being modelled in metres
    /// rather than millimetres, and an absolute default does not.
    ///
    /// # Errors
    ///
    /// [`OgeomError::Construction`](ogeom_core::OgeomError::Construction) if `size` or
    /// `fraction` is not finite and positive.
    pub fn relative(size: f64, fraction: f64) -> OgeomResult<Self> {
        if !size.is_finite() || size <= 0.0 || !fraction.is_finite() || fraction <= 0.0 {
            ogeom_bail!(
                Construction,
                "relative deflection needs a positive size and fraction, got {size} and {fraction}"
            );
        }
        Self::with_chord(size * fraction)
    }

    /// Check that the settings are usable.
    ///
    /// # Errors
    ///
    /// [`OgeomError::Construction`](ogeom_core::OgeomError::Construction) if a tolerance
    /// is not positive, or the segment bounds are inconsistent.
    pub fn validate(&self) -> OgeomResult<()> {
        if !self.chord.is_finite() || self.chord <= 0.0 {
            ogeom_bail!(
                Construction,
                "chord deflection {} must be positive",
                self.chord
            );
        }
        if !self.angular.is_finite() || self.angular <= 0.0 {
            ogeom_bail!(
                Construction,
                "angular deflection {} must be positive",
                self.angular
            );
        }
        if self.min_segments == 0 {
            ogeom_bail!(Construction, "a polyline needs at least one segment");
        }
        if self.max_segments < self.min_segments {
            ogeom_bail!(
                Construction,
                "segment ceiling {} is below the floor {}",
                self.max_segments,
                self.min_segments
            );
        }
        Ok(())
    }
}

/// A curve approximated by points, with the parameters they came from.
///
/// The parameters are kept, not discarded, because the pcurve of the same edge
/// on an adjacent face has to be sampled at exactly these values for the two
/// faces to meet.
#[derive(Debug, Clone, PartialEq)]
pub struct Polyline {
    /// The points, in order along the curve.
    pub points: Vec<Point>,
    /// The parameter each point came from.
    pub parameters: Vec<f64>,
    /// Whether the requested deflection was actually achieved.
    ///
    /// `false` when the segment ceiling was reached first. Reporting it beats
    /// silently returning a coarser polyline than asked for, which would make
    /// a downstream tolerance claim untrue with nothing to show why.
    pub deflection_met: bool,
}

impl Polyline {
    /// Number of segments.
    #[must_use]
    pub fn segment_count(&self) -> usize {
        self.points.len().saturating_sub(1)
    }

    /// Total length of the polyline.
    ///
    /// An underestimate of the curve's own length, since a chord is shorter
    /// than the arc it spans, and one that improves as the deflection tightens.
    #[must_use]
    pub fn length(&self) -> f64 {
        self.points.windows(2).map(|w| w[0].distance(w[1])).sum()
    }

    /// Whether the polyline returns to where it started.
    #[must_use]
    pub fn is_closed(&self, tol: Tolerances) -> bool {
        match (self.points.first(), self.points.last()) {
            (Some(a), Some(b)) => self.points.len() > 2 && a.is_equal(*b, tol),
            _ => false,
        }
    }
}

/// Whether a curve is a straight line, seeing through any trimming.
///
/// A line is its own polyline, so subdividing it is pure cost: the extra points
/// land on the curve, satisfy every tolerance, and then propagate into the
/// triangulation of every face the edge bounds.
#[must_use]
pub fn is_straight(curve: &Curve) -> bool {
    match curve {
        Curve::Line(_) => true,
        Curve::Trimmed(t) => is_straight(t.basis()),
        _ => false,
    }
}

/// Whether a planar curve is a straight line in parameter space.
#[must_use]
pub fn is_straight_planar(curve: &ogeom_geom::PlanarCurve) -> bool {
    match curve {
        ogeom_geom::PlanarCurve::Line(_) => true,
        ogeom_geom::PlanarCurve::Trimmed(t) => is_straight_planar(t.basis()),
        _ => false,
    }
}

/// Approximate `curve` over `range` within `deflection`.
///
/// Adaptive bisection: split a segment whenever its midpoint is further from the
/// chord than allowed, or the tangent turns too far across it. Uniform sampling
/// is the obvious alternative and wastes points on the straight parts of a curve
/// while still missing the tight ones; the whole difficulty of tessellation is
/// that curvature is not uniform.
///
/// # Errors
///
/// [`OgeomError::Construction`](ogeom_core::OgeomError::Construction) if the deflection
/// settings are unusable or the range is empty;
/// [`OgeomError::Domain`](ogeom_core::OgeomError::Domain) if the range leaves the curve.
pub fn discretize(
    curve: &Curve,
    range: (f64, f64),
    deflection: Deflection,
    tol: Tolerances,
) -> OgeomResult<Polyline> {
    deflection.validate()?;
    let (lo, hi) = range;
    if !lo.is_finite() || !hi.is_finite() || hi <= lo + tol.parametric() {
        ogeom_bail!(Construction, "range [{lo}, {hi}] is empty");
    }

    // Start from the floor, so a closed curve is never approximated by a single
    // chord from a point back to itself. A straight curve is exempt: its
    // endpoints already represent it exactly.
    let start = if is_straight(curve) {
        1
    } else {
        deflection.min_segments
    };
    let parameters: Vec<f64> = (0..=start)
        .map(|i| {
            #[allow(clippy::cast_precision_loss)]
            let t = i as f64 / start as f64;
            lo + (hi - lo) * t
        })
        .collect();
    if let Some(uniform) = circle_halvings(curve, (lo, hi), start, deflection) {
        return halved(curve, parameters, uniform, tol);
    }
    bisect(curve, (lo, hi), parameters, deflection, tol)
}

/// The seed parameters halved `times` over, by the midpoints the bisection
/// would take, and the curve's points there.
fn halved(curve: &Curve, seeds: Vec<f64>, times: u32, tol: Tolerances) -> OgeomResult<Polyline> {
    let mut parameters = seeds;
    for _ in 0..times {
        let mut finer = Vec::with_capacity(parameters.len() * 2);
        for pair in parameters.windows(2) {
            finer.push(pair[0]);
            finer.push(f64::midpoint(pair[0], pair[1]));
        }
        finer.extend(parameters.last().copied());
        parameters = finer;
    }
    let points = parameters
        .iter()
        .map(|u| curve.point_at(*u, tol))
        .collect::<OgeomResult<_>>()?;
    Ok(Polyline {
        points,
        parameters,
        deflection_met: true,
    })
}

/// The leftmost-first bisection from the seed parameters: each segment
/// split while [`needs_split`] asks it to, up to the segment cap. A curve
/// the cap stops part-way is split again by rounds instead, so the points
/// the cap allows spread over the whole curve rather than crowding its
/// start and leaving the rest at its seeds.
fn bisect(
    curve: &Curve,
    whole: (f64, f64),
    parameters: Vec<f64>,
    deflection: Deflection,
    tol: Tolerances,
) -> OgeomResult<Polyline> {
    let seeds = parameters.clone();
    let line = bisect_leftmost(curve, whole, parameters, deflection, tol)?;
    if line.deflection_met || line.parameters.len() < deflection.max_segments {
        return Ok(line);
    }
    bisect_by_rounds(curve, whole, seeds, deflection, tol)
}

/// Every segment that [`needs_split`] asks to split, split at once, round
/// after round, until none asks or the next round would pass the segment
/// cap: the polyline is as fine everywhere as the cap lets it be, never
/// fine at one end and bare at the other.
fn bisect_by_rounds(
    curve: &Curve,
    whole: (f64, f64),
    parameters: Vec<f64>,
    deflection: Deflection,
    tol: Tolerances,
) -> OgeomResult<Polyline> {
    let mut points: Vec<Point> = parameters
        .iter()
        .map(|u| curve.point_at(*u, tol))
        .collect::<OgeomResult<_>>()?;
    let mut parameters = parameters;
    loop {
        let mut splits: Vec<(usize, f64, Point)> = Vec::new();
        for i in 0..parameters.len().saturating_sub(1) {
            let (t0, t1) = (parameters[i], parameters[i + 1]);
            if let Some((mid, on_curve)) = needs_split(
                curve,
                (t0, t1),
                (points[i], points[i + 1]),
                whole,
                deflection,
                tol,
            )? && mid > t0
                && mid < t1
            {
                splits.push((i, mid, on_curve));
            }
        }
        if splits.is_empty() {
            return Ok(Polyline {
                points,
                parameters,
                deflection_met: true,
            });
        }
        if parameters.len() + splits.len() > deflection.max_segments {
            return Ok(Polyline {
                points,
                parameters,
                deflection_met: false,
            });
        }
        let mut next_parameters = Vec::with_capacity(parameters.len() + splits.len());
        let mut next_points = Vec::with_capacity(parameters.len() + splits.len());
        let mut pending = splits.into_iter().peekable();
        for i in 0..parameters.len() {
            next_parameters.push(parameters[i]);
            next_points.push(points[i]);
            if let Some(&(at, mid, on_curve)) = pending.peek()
                && at == i
            {
                next_parameters.push(mid);
                next_points.push(on_curve);
                pending.next();
            }
        }
        parameters = next_parameters;
        points = next_points;
    }
}

/// [`bisect`] leftmost first, stopping where the segment cap is reached.
fn bisect_leftmost(
    curve: &Curve,
    (lo, hi): (f64, f64),
    mut parameters: Vec<f64>,
    deflection: Deflection,
    tol: Tolerances,
) -> OgeomResult<Polyline> {
    let mut points: Vec<Point> = parameters
        .iter()
        .map(|u| curve.point_at(*u, tol))
        .collect::<OgeomResult<_>>()?;

    // Leftmost-first subdivision with a worklist: each segment is settled
    // before the walk moves right, and a split pushes its halves back for
    // re-examination. The splits come in the same order a rescan from zero
    // would make them, so the same floats come out, without re-asking every
    // settled segment on every pass or shifting the vectors per split.
    let mut met = true;
    let seeds = core::mem::take(&mut parameters);
    let seed_points = core::mem::take(&mut points);
    parameters.push(seeds[0]);
    points.push(seed_points[0]);
    // The unsettled boundary, rightmost at the bottom; the top is always the
    // segment just right of the settled prefix.
    let mut pending: Vec<(f64, Point)> = seeds[1..]
        .iter()
        .copied()
        .zip(seed_points[1..].iter().copied())
        .rev()
        .collect();
    let mut splitting = true;
    while let Some((t1, p1)) = pending.pop() {
        let t0 = *parameters.last().unwrap_or(&t1);
        let p0 = points.last().copied().unwrap_or(p1);
        if splitting
            && let Some((mid, on_curve)) =
                needs_split(curve, (t0, t1), (p0, p1), (lo, hi), deflection, tol)?
        {
            // The count the cap sees is every point currently alive, settled
            // and pending together.
            if parameters.len() + pending.len() + 1 > deflection.max_segments {
                met = false;
                splitting = false;
            } else if mid <= t0 || mid >= t1 {
                // A split that does not actually divide the interval means
                // the parameters have reached the resolution of f64;
                // refining further would loop without improving anything.
                met = false;
                splitting = false;
            } else {
                pending.push((t1, p1));
                pending.push((mid, on_curve));
                continue;
            }
        }
        parameters.push(t1);
        points.push(p1);
    }

    Ok(Polyline {
        points,
        parameters,
        deflection_met: met,
    })
}

/// How many times a circular arc's seed segments are halved, all of them
/// alike: on a circle every segment of one depth spans the same angle, so
/// [`needs_split`] answers them all the same, and the depth it settles at
/// is the first whose span passes. `None` for any other curve, where the
/// segment cap would stop the halving part-way, or where a span stands
/// within rounding of a test's limit.
fn circle_halvings(
    curve: &Curve,
    (lo, hi): (f64, f64),
    seeds: usize,
    deflection: Deflection,
) -> Option<u32> {
    let Curve::Circle(circle) = curve else {
        return None;
    };
    let radius = circle.circle().radius();
    let whole = hi - lo;
    let mut segments = seeds;
    for depth in 0..48 {
        #[allow(clippy::cast_precision_loss, reason = "a segment count")]
        let span = whole / segments as f64;
        let sag = radius * (1.0 - (span / 2.0).cos());
        let reach = 2.0 * radius * (span / 2.0).sin();
        // Where a test that decides stands within rounding of its limit, the
        // bisection's own segments, each rounded its own way, fall either
        // side of it, and only the bisection says which.
        let near = |x: f64, limit: f64| (x - limit).abs() <= limit.abs() * 1e-9;
        if near(sag, deflection.chord) {
            return None;
        }
        let split = if sag > deflection.chord {
            true
        } else {
            if near(reach, deflection.chord) {
                return None;
            }
            let short = reach < deflection.chord && {
                if near(span, whole / 16.0) {
                    return None;
                }
                span < whole / 16.0
            };
            if short {
                false
            } else {
                if near(span, deflection.angular) {
                    return None;
                }
                span > deflection.angular
            }
        };
        if !split {
            return Some(depth);
        }
        segments *= 2;
        if segments + 1 > deflection.max_segments {
            return None;
        }
    }
    None
}

/// Where one segment that violates either tolerance splits: its middle
/// parameter and the curve's point there. `None` for a segment within both.
fn needs_split(
    curve: &Curve,
    parameters: (f64, f64),
    ends: (Point, Point),
    whole: (f64, f64),
    deflection: Deflection,
    tol: Tolerances,
) -> OgeomResult<Option<(f64, Point)>> {
    let (a, b) = parameters;
    let mid = f64::midpoint(a, b);
    let on_curve = curve.point_at(mid, tol)?;
    let split = Some((mid, on_curve));

    // Chord error, measured at the midpoint. Not a bound on the true maximum
    // deviation, which would need the curve's second derivative over the span;
    // it is the standard estimate, and bisection drives it down regardless.
    let chord = ogeom_math::Axis::through(ends.0, ends.1, tol).map_or_else(
        |_| ends.0.distance(on_curve),
        |axis| axis.distance_to(on_curve),
    );
    if chord > deflection.chord {
        return Ok(split);
    }

    // Angular error: how far the tangent turns across the segment. Chord error
    // alone misses a long, gently curving span, which is exactly where a
    // silhouette goes visibly polygonal.
    //
    // Not over a segment that is both shorter than the chord tolerance and
    // a small fraction of the whole edge. A fitted edge often ends in a
    // hook a few microns long (the fit overshooting its vertex and turning
    // back), and the tangent turns through a right angle across it at
    // every scale. Asked of it, the angular test bisects the hook down to
    // the resolution of the parameter and hands the face a fan of hairs at
    // one corner, each a fin off the surface. A turn across a span under
    // the deflection the caller accepted is below what they can see; the
    // fraction keeps the test on a hole a fifth of a millimetre across,
    // every segment of which is under the chord and a thirteenth of the
    // whole.
    if ends.0.distance(ends.1) <= deflection.chord && (b - a) <= (whole.1 - whole.0) / 16.0 {
        return Ok(None);
    }
    let (Ok(start), Ok(end)) = (curve.tangent_at(a, tol), curve.tangent_at(b, tol)) else {
        // A cusp has no tangent to compare; the chord test still governs.
        return Ok(None);
    };
    Ok(split.filter(|_| start.angle(end) > deflection.angular))
}

/// Discretize a pcurve with its chord tolerance measured *in space*,
/// through the surface that gives its chart a metric.
///
/// [`discretize_planar`] measures in parameter units because it has no
/// surface to convert through; this is the version that does. Each candidate
/// segment is lifted to the surface and the sagitta measured between world
/// points, so one chord tolerance means one thing whatever the chart's
/// scale: a quarter-turn on a large cylinder refines further than the same
/// quarter-turn on a small one.
///
/// # Errors
///
/// As [`discretize_planar`].
pub fn discretize_on_surface(
    curve: &ogeom_geom::PlanarCurve,
    range: (f64, f64),
    surface: &ogeom_geom::SurfaceGeometry,
    deflection: Deflection,
    tol: Tolerances,
) -> OgeomResult<(Vec<ogeom_math::Point2>, Vec<f64>)> {
    use ogeom_geom::Curve2d;
    use ogeom_geom::Surface as _;

    deflection.validate()?;
    let (lo, hi) = range;
    if !lo.is_finite() || !hi.is_finite() || hi <= lo + tol.parametric() {
        ogeom_bail!(Construction, "range [{lo}, {hi}] is empty");
    }
    let lift = |uv: ogeom_math::Point2| -> OgeomResult<ogeom_math::Point> {
        surface.point_at(uv.x, uv.y, tol)
    };

    let start = if is_straight_planar(curve) {
        1
    } else {
        deflection.min_segments
    };
    let mut parameters: Vec<f64> = (0..=start)
        .map(|i| {
            #[allow(clippy::cast_precision_loss)]
            let t = i as f64 / start as f64;
            lo + (hi - lo) * t
        })
        .collect();
    let mut points: Vec<ogeom_math::Point2> = parameters
        .iter()
        .map(|u| curve.point_at(*u, tol))
        .collect::<OgeomResult<_>>()?;
    let mut lifted: Vec<ogeom_math::Point> = points
        .iter()
        .map(|uv| lift(*uv))
        .collect::<OgeomResult<_>>()?;

    while points.len() <= deflection.max_segments {
        let mut split_at = None;
        for i in 0..points.len() - 1 {
            let mid = f64::midpoint(parameters[i], parameters[i + 1]);
            let on_curve = curve.point_at(mid, tol)?;
            let in_space = lift(on_curve)?;
            // Sagitta in world units: the lifted midpoint against the lifted
            // chord.
            let (a, b) = (lifted[i], lifted[i + 1]);
            let chord_vector = b - a;
            let length = chord_vector.magnitude();
            let sagitta = if length <= tol.confusion() {
                in_space.distance(a)
            } else {
                (in_space - a).cross(chord_vector).magnitude() / length
            };
            if sagitta > deflection.chord {
                split_at = Some((i, mid, on_curve, in_space));
                break;
            }
        }
        let Some((i, mid, on_curve, in_space)) = split_at else {
            break;
        };
        parameters.insert(i + 1, mid);
        points.insert(i + 1, on_curve);
        lifted.insert(i + 1, in_space);
    }
    Ok((points, parameters))
}

/// Approximate a planar curve in a surface's parameter space.
///
/// The deflection is measured in parameter units here, not in space, so a caller
/// wanting a spatial tolerance has to convert through the surface's own scale:
/// the two differ by orders of magnitude near a pole. [`discretize_on_surface`]
/// is the version that measures through the surface. This exists for boundary
/// work in parameter space; for a face's actual boundary, discretize the edge's
/// 3D curve and evaluate the pcurve at those parameters instead, which is what
/// keeps adjacent faces watertight.
///
/// # Errors
///
/// As [`discretize`].
pub fn discretize_planar(
    curve: &ogeom_geom::PlanarCurve,
    range: (f64, f64),
    deflection: Deflection,
    tol: Tolerances,
) -> OgeomResult<(Vec<ogeom_math::Point2>, Vec<f64>)> {
    use ogeom_geom::Curve2d;

    deflection.validate()?;
    let (lo, hi) = range;
    if !lo.is_finite() || !hi.is_finite() || hi <= lo + tol.parametric() {
        ogeom_bail!(Construction, "range [{lo}, {hi}] is empty");
    }

    let start = if is_straight_planar(curve) {
        1
    } else {
        deflection.min_segments
    };
    let mut parameters: Vec<f64> = (0..=start)
        .map(|i| {
            #[allow(clippy::cast_precision_loss)]
            let t = i as f64 / start as f64;
            lo + (hi - lo) * t
        })
        .collect();
    let mut points: Vec<ogeom_math::Point2> = parameters
        .iter()
        .map(|u| curve.point_at(*u, tol))
        .collect::<OgeomResult<_>>()?;

    while points.len() <= deflection.max_segments {
        let mut split_at = None;
        for i in 0..points.len() - 1 {
            let mid = f64::midpoint(parameters[i], parameters[i + 1]);
            let on_curve = curve.point_at(mid, tol)?;
            let chord = ogeom_math::Axis2::through(points[i], points[i + 1], tol).map_or_else(
                |_| points[i].distance(on_curve),
                |axis| axis.distance_to(on_curve),
            );
            if chord > deflection.chord {
                split_at = Some(i);
                break;
            }
        }
        let Some(i) = split_at else { break };
        let mid = f64::midpoint(parameters[i], parameters[i + 1]);
        if mid <= parameters[i] || mid >= parameters[i + 1] {
            break;
        }
        parameters.insert(i + 1, mid);
        points.insert(i + 1, curve.point_at(mid, tol)?);
    }

    Ok((points, parameters))
}

#[cfg(test)]
#[allow(clippy::unwrap_used)]
mod tests {
    use super::*;
    use approx::assert_relative_eq;
    use ogeom_geom::{BSplineCurve, CircleCurve, LineCurve};
    use ogeom_math::{Circle, Frame, KnotVector};

    const T: Tolerances = Tolerances::millimetres();

    fn circle(radius: f64) -> Curve {
        CircleCurve::new(Circle::new(Frame::WORLD, radius, T).unwrap()).into()
    }

    /// A circle is halved uniformly in closed form, and comes out exactly
    /// as the bisection would make it: the same parameters, bit for bit,
    /// over radii, arcs and deflections either side of each test's limit.
    #[test]
    fn a_circle_halves_as_the_bisection_would() {
        for radius in [1e-3, 0.37, 1.0, 2.5, 40.0, 900.0] {
            for (lo, hi) in [(0.0, core::f64::consts::TAU), (0.3, 1.9), (-2.0, 0.05)] {
                for chord in [1e-4, 1e-3, 0.01, 0.1] {
                    for angular in [0.05, 0.2, 0.5] {
                        let deflection = Deflection {
                            chord,
                            angular,
                            ..Deflection::default()
                        };
                        let curve = circle(radius);
                        let fast = discretize(&curve, (lo, hi), deflection, T).unwrap();
                        let seeds: Vec<f64> = (0..=deflection.min_segments)
                            .map(|i| {
                                #[allow(clippy::cast_precision_loss)]
                                let t = i as f64 / deflection.min_segments as f64;
                                lo + (hi - lo) * t
                            })
                            .collect();
                        let slow = bisect(&curve, (lo, hi), seeds, deflection, T).unwrap();
                        assert_eq!(
                            fast.parameters, slow.parameters,
                            "r {radius} ({lo}, {hi}) chord {chord} angular {angular}"
                        );
                        assert_eq!(fast.deflection_met, slow.deflection_met);
                    }
                }
            }
        }
    }

    /// The greatest distance from the curve to the polyline, sampled densely.
    fn worst_error(curve: &Curve, line: &Polyline) -> f64 {
        let mut worst: f64 = 0.0;
        for window in line.parameters.windows(2) {
            for k in 1..16 {
                let t = f64::from(k) / 16.0;
                let u = window[0] + (window[1] - window[0]) * t;
                let on_curve = curve.point_at(u, T).unwrap();
                let a = curve.point_at(window[0], T).unwrap();
                let b = curve.point_at(window[1], T).unwrap();
                let chord = ogeom_math::Axis::through(a, b, T)
                    .map_or(0.0, |axis| axis.distance_to(on_curve));
                worst = worst.max(chord);
            }
        }
        worst
    }

    #[test]
    fn straightness_sees_through_a_trim() {
        // The exemption has to survive trimming, or an edge built from a
        // trimmed line (which is what a solid's edges usually are) pays the
        // floor anyway and the saving evaporates.
        use ogeom_geom::TrimmedCurve;
        let line: Curve = LineCurve::segment(Point::ORIGIN, Point::new(10.0, 0.0, 0.0), T)
            .unwrap()
            .into();
        assert!(is_straight(&line));
        assert!(is_straight(
            &TrimmedCurve::new(line, 2.0, 8.0, T).unwrap().into()
        ));
        assert!(!is_straight(&circle(1.0)));
    }

    #[test]
    fn a_line_needs_no_more_than_the_minimum_segments() {
        // A straight curve has no chord error and no turn, so refinement must
        // stop immediately rather than subdividing to the ceiling.
        let curve: Curve = LineCurve::segment(Point::ORIGIN, Point::new(100.0, 0.0, 0.0), T)
            .unwrap()
            .into();
        let line = discretize(&curve, curve.domain(), Deflection::default(), T).unwrap();
        assert_eq!(line.segment_count(), 1, "a line is its own polyline");
        assert!(line.deflection_met);
        assert_relative_eq!(line.length(), 100.0, epsilon = 1e-9);
    }

    #[test]
    fn a_circle_is_refined_until_the_chord_tolerance_is_met() {
        let curve = circle(10.0);
        for chord in [1.0_f64, 0.1, 0.01, 0.001] {
            let deflection = Deflection {
                chord,
                ..Deflection::default()
            };
            let line = discretize(&curve, curve.domain(), deflection, T).unwrap();
            assert!(line.deflection_met, "gave up at chord {chord}");
            assert!(
                worst_error(&curve, &line) <= chord * 1.5,
                "chord {chord}: worst error {}",
                worst_error(&curve, &line)
            );
        }
    }

    #[test]
    fn a_tighter_tolerance_always_gives_at_least_as_many_segments() {
        let curve = circle(10.0);
        let mut previous = 0;
        for chord in [2.0_f64, 1.0, 0.5, 0.1, 0.01] {
            let line = discretize(
                &curve,
                curve.domain(),
                Deflection {
                    chord,
                    ..Deflection::default()
                },
                T,
            )
            .unwrap();
            assert!(
                line.segment_count() >= previous,
                "chord {chord} gave fewer segments than a looser one"
            );
            previous = line.segment_count();
        }
    }

    #[test]
    fn a_polylines_length_underestimates_the_curve_and_converges_to_it() {
        // Each chord is shorter than the arc it spans, so the polyline is always
        // short, and refining closes the gap.
        let radius = 10.0;
        let curve = circle(radius);
        let exact = core::f64::consts::TAU * radius;

        let coarse = discretize(
            &curve,
            curve.domain(),
            Deflection {
                chord: 1.0,
                ..Deflection::default()
            },
            T,
        )
        .unwrap();
        let fine = discretize(
            &curve,
            curve.domain(),
            Deflection {
                chord: 1e-4,
                ..Deflection::default()
            },
            T,
        )
        .unwrap();

        assert!(coarse.length() < exact);
        assert!(fine.length() < exact);
        assert!(fine.length() > coarse.length());
        assert_relative_eq!(fine.length(), exact, max_relative = 1e-3);
    }

    #[test]
    fn the_angular_tolerance_catches_what_the_chord_one_misses() {
        // A large circle with a loose chord tolerance: the chord error over a
        // whole quadrant is huge in absolute terms but the *ratio* to the radius
        // is what a viewer sees, and the tangent turn is what bounds it.
        let curve = circle(1000.0);
        let chord_only = Deflection {
            chord: 50.0,
            angular: 10.0,
            ..Deflection::default()
        };
        let with_angle = Deflection {
            chord: 50.0,
            angular: 0.1,
            ..Deflection::default()
        };

        let loose = discretize(&curve, curve.domain(), chord_only, T).unwrap();
        let tight = discretize(&curve, curve.domain(), with_angle, T).unwrap();
        assert!(
            tight.segment_count() > loose.segment_count(),
            "the angular limit did nothing: {} vs {}",
            tight.segment_count(),
            loose.segment_count()
        );

        // Every segment really does turn by less than the limit.
        for window in tight.parameters.windows(2) {
            let a = curve.tangent_at(window[0], T).unwrap();
            let b = curve.tangent_at(window[1], T).unwrap();
            assert!(a.angle(b) <= 0.1 + 1e-9);
        }
    }

    #[test]
    fn a_closed_curve_gets_enough_segments_to_enclose_something() {
        // With a minimum of one, a full circle would come out as a single chord
        // from a point back to itself: zero length, zero area, and no error
        // reported anywhere.
        let curve = circle(5.0);
        let line = discretize(
            &curve,
            curve.domain(),
            Deflection {
                chord: 1e6,
                angular: 1e6,
                ..Deflection::default()
            },
            T,
        )
        .unwrap();
        assert!(line.segment_count() >= 2);
        assert!(line.length() > 0.0);
        assert!(line.is_closed(T));
    }

    #[test]
    fn reaching_the_ceiling_is_reported_rather_than_passed_off_as_success() {
        // A downstream tolerance claim built on a polyline that never met its
        // own would be untrue, with nothing to show why.
        let curve = circle(10.0);
        let line = discretize(
            &curve,
            curve.domain(),
            Deflection {
                chord: 1e-12,
                angular: 1e-12,
                min_segments: 2,
                max_segments: 16,
            },
            T,
        )
        .unwrap();
        assert!(!line.deflection_met);
        assert!(line.segment_count() <= 20);
    }

    #[test]
    fn a_curve_the_ceiling_stops_is_sampled_all_along() {
        // A wavy spline drawn to a chord far finer than the ceiling allows:
        // the points the ceiling does allow spread over the whole curve,
        // or the end the bisection never reached is drawn as a few long
        // chords across its waves, and a face bounded by it meshes past
        // its rim.
        let control: Vec<Point> = (0..40)
            .map(|i| Point::new(f64::from(i), if i % 2 == 0 { 0.0 } else { 1.0 }, 0.0))
            .collect();
        let knots = KnotVector::clamped_uniform(3, control.len()).unwrap();
        let curve: Curve = BSplineCurve::new(knots, control, T).unwrap().into();
        let (lo, hi) = curve.domain();
        let line = discretize(
            &curve,
            (lo, hi),
            Deflection {
                chord: 1e-9,
                angular: 0.5,
                min_segments: 2,
                max_segments: 256,
            },
            T,
        )
        .unwrap();
        assert!(!line.deflection_met);
        let widest = line
            .parameters
            .windows(2)
            .map(|w| w[1] - w[0])
            .fold(0.0_f64, f64::max);
        assert!(
            widest <= (hi - lo) / 64.0,
            "a chord spans {widest} of {}",
            hi - lo
        );
    }

    #[test]
    fn parameters_are_kept_alongside_the_points() {
        // Two faces meeting along an edge must sample its pcurves at exactly
        // these values, or the join is not watertight. Discarding them would
        // make that impossible to arrange.
        let curve = circle(3.0);
        let line = discretize(&curve, curve.domain(), Deflection::default(), T).unwrap();
        assert_eq!(line.points.len(), line.parameters.len());
        for (u, p) in line.parameters.iter().zip(&line.points) {
            assert!(curve.point_at(*u, T).unwrap().is_equal(*p, T));
        }
        // And they increase along the curve.
        assert!(line.parameters.windows(2).all(|w| w[1] > w[0]));
    }

    #[test]
    fn a_spline_is_refined_where_it_curves_and_not_where_it_does_not() {
        // The reason for adaptive rather than uniform sampling: a curve that is
        // straight for half its length and tight for the rest should not pay
        // for the tight part everywhere.
        let control = vec![
            Point::new(0.0, 0.0, 0.0),
            Point::new(10.0, 0.0, 0.0),
            Point::new(20.0, 0.0, 0.0),
            Point::new(21.0, 8.0, 0.0),
            Point::new(22.0, 0.0, 0.0),
        ];
        let curve: Curve = BSplineCurve::new(
            KnotVector::clamped_uniform(3, control.len()).unwrap(),
            control,
            T,
        )
        .unwrap()
        .into();

        let line = discretize(
            &curve,
            curve.domain(),
            Deflection {
                chord: 0.05,
                ..Deflection::default()
            },
            T,
        )
        .unwrap();
        assert!(line.deflection_met);

        // Segments in the straight first half are longer than in the curved
        // second half.
        let mid = line.points.len() / 2;
        let mean = |points: &[Point]| {
            let gaps: Vec<f64> = points.windows(2).map(|w| w[0].distance(w[1])).collect();
            #[allow(clippy::cast_precision_loss)]
            let count = gaps.len() as f64;
            gaps.iter().sum::<f64>() / count
        };
        let (early, late) = (mean(&line.points[..mid]), mean(&line.points[mid..]));
        assert!(early > late, "uniform spacing: {early} vs {late}");
    }

    #[test]
    fn discretizing_part_of_a_curve_covers_only_that_part() {
        let curve = circle(4.0);
        let line = discretize(&curve, (1.0, 2.0), Deflection::default(), T).unwrap();
        assert_relative_eq!(line.parameters[0], 1.0);
        assert_relative_eq!(line.parameters[line.parameters.len() - 1], 2.0);
        assert!(!line.is_closed(T));
        assert_relative_eq!(line.length(), 4.0, max_relative = 1e-2);
    }

    #[test]
    fn unusable_settings_are_refused() {
        let curve = circle(1.0);
        let bad = [
            Deflection {
                chord: 0.0,
                ..Deflection::default()
            },
            Deflection {
                chord: f64::NAN,
                ..Deflection::default()
            },
            Deflection {
                angular: -1.0,
                ..Deflection::default()
            },
            Deflection {
                min_segments: 0,
                ..Deflection::default()
            },
            Deflection {
                min_segments: 10,
                max_segments: 5,
                ..Deflection::default()
            },
        ];
        for deflection in bad {
            assert!(
                discretize(&curve, curve.domain(), deflection, T).is_err(),
                "accepted {deflection:?}"
            );
        }
        assert!(discretize(&curve, (1.0, 1.0), Deflection::default(), T).is_err());
        assert!(Deflection::with_chord(-1.0).is_err());
        assert!(Deflection::relative(0.0, 0.001).is_err());
        assert!(Deflection::relative(100.0, 0.001).is_ok());
    }

    #[test]
    fn a_relative_deflection_scales_with_the_model() {
        // "A thousandth of the part" survives the part being modelled in metres
        // rather than millimetres; an absolute default does not.
        let small = Deflection::relative(1.0, 0.001).unwrap();
        let large = Deflection::relative(1000.0, 0.001).unwrap();
        assert_relative_eq!(large.chord, small.chord * 1000.0);
    }

    #[test]
    fn a_planar_curve_discretizes_in_parameter_space() {
        let curve: ogeom_geom::PlanarCurve = ogeom_geom::Circle2d::new(
            ogeom_math::Circle2::centred(ogeom_math::Point2::ORIGIN, 5.0, T).unwrap(),
        )
        .into();
        let (points, parameters) = discretize_planar(
            &curve,
            (0.0, core::f64::consts::TAU),
            Deflection {
                chord: 0.05,
                ..Deflection::default()
            },
            T,
        )
        .unwrap();
        assert_eq!(points.len(), parameters.len());
        assert!(points.len() > 20, "only {} points", points.len());
        for p in &points {
            assert_relative_eq!(p.to_vector().magnitude(), 5.0, epsilon = 1e-12);
        }
    }
}
#[cfg(test)]
#[allow(clippy::unwrap_used)]
mod on_surface_tests {
    use super::*;
    use ogeom_core::Tolerances;
    use ogeom_geom::{Circle2d, PlaneSurface, Surface as _};
    use ogeom_math::{Circle2, Frame, Plane, Point2};

    const T: Tolerances = Tolerances::millimetres();

    #[test]
    fn the_spatial_chord_scales_with_the_surface_not_the_chart() {
        // Two circles in a plane's chart, radii 5 and 100: one chord
        // tolerance, measured in space, refines the big one further and
        // holds both to the same sagitta.
        let plane: ogeom_geom::SurfaceGeometry =
            PlaneSurface::over(Plane::new(Frame::WORLD), (-200.0, 200.0), (-200.0, 200.0))
                .unwrap()
                .into();
        let deflection = Deflection {
            chord: 0.05,
            ..Deflection::default()
        };
        let counts: Vec<usize> = [5.0, 100.0]
            .iter()
            .map(|&radius| {
                let circle: ogeom_geom::PlanarCurve =
                    Circle2d::new(Circle2::centred(Point2::new(0.0, 0.0), radius, T).unwrap())
                        .into();
                let (points, parameters) = discretize_on_surface(
                    &circle,
                    (0.0, core::f64::consts::TAU),
                    &plane,
                    deflection,
                    T,
                )
                .unwrap();
                // Every lifted segment's sagitta is within the chord,
                // measured at the segment's own parameter midpoint.
                for (pair, params) in points.windows(2).zip(parameters.windows(2)) {
                    use ogeom_geom::Curve2d as _;
                    let a = plane.point_at(pair[0].x, pair[0].y, T).unwrap();
                    let b = plane.point_at(pair[1].x, pair[1].y, T).unwrap();
                    let mid = circle
                        .point_at(f64::midpoint(params[0], params[1]), T)
                        .unwrap();
                    let on = plane.point_at(mid.x, mid.y, T).unwrap();
                    let chord = b - a;
                    let sagitta = (on - a).cross(chord).magnitude() / chord.magnitude();
                    assert!(
                        sagitta <= deflection.chord * 1.5,
                        "sagitta {sagitta} at radius {radius}"
                    );
                }
                points.len()
            })
            .collect();
        assert!(
            counts[1] > counts[0] * 2,
            "the larger circle refines further: {counts:?}"
        );
    }
}