ogeom_algo/
pcurve_gap.rs

1//! How far an edge's pcurves stand off its curve, and edges' tolerances
2//! raised to cover it.
3//!
4//! An edge's tolerance is the radius about its curve within which every
5//! description of the edge lies, each pcurve lifted through its surface
6//! included. A producer that fits or projects a pcurve knows its error
7//! only where it sampled; these measure it densely, after the fact.
8
9use ogeom_core::{OgeomResult, Tolerance, Tolerances};
10use ogeom_geom::{Curve, Curve2d as _, Curve3d as _, PlanarCurve, Surface as _, SurfaceGeometry};
11use ogeom_math::{Point, Transform, Vector};
12use ogeom_topo::{EdgeRepr, Model, Shape, ShapeType};
13
14mod closed;
15pub(crate) use closed::closed_gap;
16
17/// How much wider than the gap it was measured from an edge's tolerance is
18/// recorded: a millionth.
19const MARGIN: f64 = 1e-6;
20
21/// How far a pcurve, lifted through its surface, stands off an edge's
22/// curve: at each of many samples along both ranges alike, the distance
23/// between the two points, or where that exceeds `stated`, the distance
24/// from the lifted point to the nearest point of the curve's stretch, and
25/// the widest peaks between samples searched for their tops, so a stretch
26/// of the edge sampled anywhere stays within the answer. Within `stated`,
27/// the answer is only known to be within it.
28///
29/// Unlike [`lifted_gap`](crate::pcurve_fit::lifted_gap), which charges a
30/// pcurve for a pace differing from the curve's by more than a sample's
31/// step, this measures where the pcurve lies, as `check` does: an iso line
32/// run the other way round a closed rim stands on the rim.
33///
34/// A line, a circle or an ellipse against a line, a circle or an ellipse
35/// on a plane, or an iso line of a cylinder, a cone, a sphere or a torus,
36/// is answered in closed form where that stands within `stated`: within
37/// rounding of what the sampling measures.
38///
39/// # Errors
40///
41/// If the curve or the pcurve cannot be evaluated over its range.
42pub fn pcurve_gap(
43    curve: (&Curve, (f64, f64)),
44    pcurve: (&PlanarCurve, (f64, f64)),
45    surface: &SurfaceGeometry,
46    stated: f64,
47    tol: Tolerances,
48) -> OgeomResult<f64> {
49    let Some((gap, margin)) = closed::closed_gap(curve, pcurve, surface, stated, tol) else {
50        return sampled_gap(curve, pcurve, surface, stated, tol);
51    };
52    if cfg!(debug_assertions) {
53        let sampled = sampled_gap(curve, pcurve, surface, stated, tol);
54        assert!(
55            sampled
56                .as_ref()
57                .is_ok_and(|s| (gap - s).abs() <= margin && *s <= stated),
58            "closed-form pcurve gap {gap:e} (within {margin:e}) against sampled {sampled:?}, \
59             stated {stated:e}"
60        );
61    }
62    Ok(gap)
63}
64
65/// [`pcurve_gap`] by sampling alone.
66pub(crate) fn sampled_gap(
67    (curve, range): (&Curve, (f64, f64)),
68    (pcurve, pcurve_range): (&PlanarCurve, (f64, f64)),
69    surface: &SurfaceGeometry,
70    stated: f64,
71    tol: Tolerances,
72) -> OgeomResult<f64> {
73    // The checker's sample counts divide this one, so every point it
74    // samples is one of these.
75    const SAMPLES: u32 = 256;
76    // A sample this close to the widest may stand beside a peak wider
77    // than every sample: the stretch either side of it is searched, for
78    // at most this many of the widest.
79    const NEAR_PEAK: f64 = 0.9;
80    const PEAKS: usize = 4;
81    // The gap at `f` of the way along, the nearest point of the curve
82    // standing in where the same parameter's is past `within`.
83    let gap_beyond = |f: f64, within: f64| -> OgeomResult<Option<f64>> {
84        let uv = pcurve.point_at(pcurve_range.0 + (pcurve_range.1 - pcurve_range.0) * f, tol)?;
85        let Ok(lifted) = surface.point_at(uv.x, uv.y, tol) else {
86            return Ok(None);
87        };
88        let gap = curve
89            .point_at(range.0 + (range.1 - range.0) * f, tol)?
90            .distance(lifted);
91        Ok(Some(if gap > within {
92            gap.min(nearest_on_stretch(curve, range, lifted, tol)?)
93        } else {
94            gap
95        }))
96    };
97    let step = 1.0 / f64::from(SAMPLES);
98    let mut widest: f64 = 0.0;
99    let mut gaps = Vec::with_capacity(SAMPLES as usize + 1);
100    for i in 0..=SAMPLES {
101        // A gap no wider than the widest so far changes nothing, and is
102        // kept as found.
103        let gap = gap_beyond(f64::from(i) * step, widest.max(stated))?;
104        if let Some(gap) = gap {
105            widest = widest.max(gap);
106        }
107        gaps.push(gap);
108    }
109    // A piece of the edge is sampled between these points, so the peaks
110    // between them count: the widest local peaks near the widest sample are
111    // searched for their tops.
112    let sampled = widest;
113    let gap_of = |j: Option<usize>| j.and_then(|j| gaps.get(j).copied().flatten());
114    let mut peaks: Vec<(f64, usize)> = gaps
115        .iter()
116        .enumerate()
117        .filter_map(|(i, gap)| {
118            let gap = (*gap)?;
119            let peak = gap >= sampled * NEAR_PEAK
120                && gap_of(i.checked_sub(1)).is_none_or(|g| g <= gap)
121                && gap_of(Some(i + 1)).is_none_or(|g| g <= gap);
122            peak.then_some((gap, i))
123        })
124        .collect();
125    peaks.sort_by(|a, b| b.0.total_cmp(&a.0));
126    let gap_at = |f: f64| gap_beyond(f, stated);
127    for &(_, i) in peaks.iter().take(PEAKS) {
128        #[allow(clippy::cast_precision_loss, reason = "a sample index")]
129        let at = i as f64 * step;
130        widest = widest.max(peak_about(&gap_at, at, step)?);
131    }
132    Ok(widest)
133}
134
135/// The widest of `gap_at` within a step either side of `at`, by a
136/// golden-section search, never less than at `at`.
137fn peak_about(
138    gap_at: &impl Fn(f64) -> OgeomResult<Option<f64>>,
139    at: f64,
140    step: f64,
141) -> OgeomResult<f64> {
142    const ROUNDS: u32 = 40;
143    let ratio = (5.0_f64.sqrt() - 1.0) / 2.0;
144    let mut widest = gap_at(at)?.unwrap_or(0.0);
145    let (mut a, mut b) = ((at - step).max(0.0), (at + step).min(1.0));
146    for _ in 0..ROUNDS {
147        let (c, d) = (b - (b - a) * ratio, a + (b - a) * ratio);
148        let (gc, gd) = (gap_at(c)?.unwrap_or(0.0), gap_at(d)?.unwrap_or(0.0));
149        widest = widest.max(gc).max(gd);
150        if gc > gd {
151            b = d;
152        } else {
153            a = c;
154        }
155    }
156    Ok(widest)
157}
158
159/// The widest `|gap(f)|` for `f` over `[0, 1]`, bounded between samples
160/// rather than read only at them.
161///
162/// The stretch is cut into equal spans and at `breaks` (fractions where a
163/// piecewise curve's smoothness drops), and each span is halved until the
164/// gap at its middle stands within `settle` of the chord between its ends.
165/// Where the gap bends as a parabola across a span, its offset from that
166/// chord peaks at the middle and the chord's length peaks at an end, so the
167/// span is bounded by its widest sample plus the middle's offset. A span
168/// bounded within `floor` is left whole, nothing in it being wanted; one
169/// narrower than `finest` is a point at the parameter's resolution, bounded
170/// by its samples. A sample that cannot be had (a lift off its surface's
171/// domain) leaves its span bounded by the samples it has.
172///
173/// # Errors
174///
175/// As `gap`.
176pub(crate) fn widest_refined(
177    gap: &impl Fn(f64) -> OgeomResult<Option<Vector>>,
178    breaks: &[f64],
179    floor: f64,
180    settle: f64,
181    finest: f64,
182) -> OgeomResult<f64> {
183    const SPANS: u32 = 32;
184    let mut cuts: Vec<f64> = (0..=SPANS)
185        .map(|i| f64::from(i) / f64::from(SPANS))
186        .chain(breaks.iter().copied().filter(|f| *f > 0.0 && *f < 1.0))
187        .collect();
188    cuts.sort_by(f64::total_cmp);
189    cuts.dedup();
190    let at = cuts
191        .iter()
192        .map(|f| gap(*f))
193        .collect::<OgeomResult<Vec<Option<Vector>>>>()?;
194    let mut spans: Vec<(f64, f64, Option<Vector>, Option<Vector>)> = cuts
195        .windows(2)
196        .zip(at.windows(2))
197        .map(|(f, g)| (f[0], f[1], g[0], g[1]))
198        .collect();
199    let mut widest: f64 = 0.0;
200    while let Some((a, b, ga, gb)) = spans.pop() {
201        let middle = f64::midpoint(a, b);
202        let gm = gap(middle)?;
203        let sampled = [ga, gm, gb]
204            .iter()
205            .flatten()
206            .map(|g| g.magnitude())
207            .fold(0.0_f64, f64::max);
208        let (Some(va), Some(vm), Some(vb)) = (ga, gm, gb) else {
209            widest = widest.max(sampled);
210            continue;
211        };
212        let bent = (vm - (va + vb) * 0.5).magnitude();
213        let bound = sampled + bent;
214        if bound <= floor || bent <= settle {
215            widest = widest.max(bound);
216        } else if b - a <= finest {
217            widest = widest.max(sampled);
218        } else {
219            spans.push((a, middle, ga, gm));
220            spans.push((middle, b, gm, gb));
221        }
222    }
223    Ok(widest)
224}
225
226/// Where `curve`'s smoothness may drop inside `range`, as fractions of the
227/// way along it: a B-spline's interior knots.
228pub(crate) fn curve_breaks(curve: &Curve, range: (f64, f64)) -> Vec<f64> {
229    match curve {
230        Curve::BSpline(spline) => fractions(spline.knots(), range),
231        _ => Vec::new(),
232    }
233}
234
235/// [`curve_breaks`] for a pcurve over `range`.
236pub(crate) fn pcurve_breaks(pcurve: &PlanarCurve, range: (f64, f64)) -> Vec<f64> {
237    match pcurve {
238        PlanarCurve::BSpline(spline) => fractions(spline.knots(), range),
239        _ => Vec::new(),
240    }
241}
242
243/// The distinct knots strictly inside `range`, as fractions of the way along
244/// it from its first end.
245fn fractions(knots: &ogeom_math::KnotVector, range: (f64, f64)) -> Vec<f64> {
246    let width = range.1 - range.0;
247    if width == 0.0 {
248        return Vec::new();
249    }
250    knots
251        .distinct()
252        .into_iter()
253        .map(|(at, _)| (at - range.0) / width)
254        .filter(|f| *f > 0.0 && *f < 1.0)
255        .collect()
256}
257
258/// The distance from `target` to the nearest point of `curve` over `range`:
259/// a scan, then a golden-section search about each sample nearer than the
260/// samples beside it. A closed curve's two ends are one point, and a target
261/// just short of one end is nearest a point beside that end, though the
262/// other end's sample ties with it.
263pub(crate) fn nearest_on_stretch(
264    curve: &Curve,
265    range: (f64, f64),
266    target: Point,
267    tol: Tolerances,
268) -> OgeomResult<f64> {
269    Ok(nearest_point_on_stretch(curve, range, target, tol)?.1)
270}
271
272/// [`nearest_on_stretch`], with the parameter of the nearest point found.
273pub(crate) fn nearest_point_on_stretch(
274    curve: &Curve,
275    range: (f64, f64),
276    target: Point,
277    tol: Tolerances,
278) -> OgeomResult<(f64, f64)> {
279    const SCAN: u32 = 64;
280    let at = |t: f64| -> OgeomResult<f64> { Ok(curve.point_at(t, tol)?.distance(target)) };
281    let step = (range.1 - range.0) / f64::from(SCAN);
282    let scan = (0..=SCAN)
283        .map(|k| {
284            let t = range.0 + step * f64::from(k);
285            Ok((t, at(t)?))
286        })
287        .collect::<OgeomResult<Vec<(f64, f64)>>>()?;
288    let (lo, hi) = (range.0.min(range.1), range.0.max(range.1));
289    let ratio = (5.0_f64.sqrt() - 1.0) / 2.0;
290    let mut nearest = (range.0, f64::INFINITY);
291    for (k, &(t, d)) in scan.iter().enumerate() {
292        if d < nearest.1 {
293            nearest = (t, d);
294        }
295        let low = k == 0 || scan[k - 1].1 >= d;
296        let high = k + 1 == scan.len() || scan[k + 1].1 >= d;
297        if !(low && high) {
298            continue;
299        }
300        let (mut a, mut b) = ((t - step.abs()).max(lo), (t + step.abs()).min(hi));
301        for _ in 0..80 {
302            let (c, d) = (b - (b - a) * ratio, a + (b - a) * ratio);
303            if at(c)? < at(d)? {
304                b = d;
305            } else {
306                a = c;
307            }
308        }
309        let middle = f64::midpoint(a, b);
310        let d = at(middle)?;
311        if d < nearest.1 {
312            nearest = (middle, d);
313        }
314    }
315    Ok(nearest)
316}
317
318/// The widest gap between `edge`'s curve and any of its pcurves at the
319/// edge's own placement, measured by [`pcurve_gap`] against the edge's
320/// stated tolerance; `None` for an edge with no curve.
321///
322/// # Errors
323///
324/// [`OgeomError::Dangling`](ogeom_core::OgeomError::Dangling) if the edge
325/// or its geometry is not in `model`; as [`pcurve_gap`].
326pub fn edge_pcurve_gap(model: &Model, edge: &Shape, tol: Tolerances) -> OgeomResult<Option<f64>> {
327    edge_gap_by(model, edge, tol, pcurve_gap)
328}
329
330/// [`edge_pcurve_gap`] by sampling alone, for a reader of the gap below the
331/// stated tolerance, where the closed form is only known to be within
332/// rounding of it.
333pub(crate) fn edge_pcurve_gap_sampled(
334    model: &Model,
335    edge: &Shape,
336    tol: Tolerances,
337) -> OgeomResult<Option<f64>> {
338    edge_gap_by(model, edge, tol, sampled_gap)
339}
340
341/// [`edge_pcurve_gap`], each pcurve's gap measured by `gap`.
342pub(crate) fn edge_gap_by(
343    model: &Model,
344    edge: &Shape,
345    tol: Tolerances,
346    mut gap: impl FnMut(
347        (&Curve, (f64, f64)),
348        (&PlanarCurve, (f64, f64)),
349        &SurfaceGeometry,
350        f64,
351        Tolerances,
352    ) -> OgeomResult<f64>,
353) -> OgeomResult<Option<f64>> {
354    let Some(data) = model.node(edge).and_then(|n| n.data().as_edge()) else {
355        ogeom_core::ogeom_bail!(Dangling, "an edge is not in this model");
356    };
357    let Some(EdgeRepr::Curve3d {
358        curve,
359        location,
360        range,
361    }) = data.curve3d()
362    else {
363        return Ok(None);
364    };
365    let Some(curve) = model.geometry().curve(*curve) else {
366        ogeom_core::ogeom_bail!(Dangling, "an edge names a curve not in this model");
367    };
368    let stated = data.tolerance.get();
369    let mut widest: f64 = 0.0;
370    for repr in &data.representations {
371        let (sides, pcurve_range, surface, at) = match repr {
372            EdgeRepr::PCurve {
373                curve,
374                range,
375                surface,
376                location,
377            } => ([Some(*curve), None], *range, *surface, location),
378            EdgeRepr::Seam {
379                forward,
380                reversed,
381                range,
382                surface,
383                location,
384            } => (
385                [Some(*forward), Some(*reversed)],
386                *range,
387                *surface,
388                location,
389            ),
390            _ => continue,
391        };
392        // A pcurve at another placement describes the edge where that
393        // occurrence stands, not where the curve does.
394        if at != location {
395            continue;
396        }
397        let Some(surface) = model.geometry().surface(surface) else {
398            ogeom_core::ogeom_bail!(Dangling, "an edge names a surface not in this model");
399        };
400        for id in sides.into_iter().flatten() {
401            let Some(pcurve) = model.geometry().pcurve(id) else {
402                ogeom_core::ogeom_bail!(Dangling, "an edge names a pcurve not in this model");
403            };
404            widest = widest.max(gap(
405                (curve, *range),
406                (pcurve, pcurve_range),
407                surface,
408                stated,
409                tol,
410            )?);
411        }
412    }
413    Ok(Some(widest))
414}
415
416/// How far apart two edges' curves stand: the farthest that any of many
417/// samples along either curve's stretch lies from the nearest point of the
418/// other's, each curve placed as its occurrence is. `None` where either
419/// edge has no curve.
420///
421/// # Errors
422///
423/// [`OgeomError::Dangling`](ogeom_core::OgeomError::Dangling) if an edge or
424/// its curve is not in `model`; [`OgeomError::Domain`](ogeom_core::OgeomError::Domain)
425/// if a curve cannot be evaluated over its stretch.
426pub fn edges_apart(
427    model: &Model,
428    a: &Shape,
429    b: &Shape,
430    tol: Tolerances,
431) -> OgeomResult<Option<f64>> {
432    let mut placed = Vec::with_capacity(2);
433    for edge in [a, b] {
434        let Some(data) = model.node(edge).and_then(|n| n.data().as_edge()) else {
435            ogeom_core::ogeom_bail!(Dangling, "an edge is not in this model");
436        };
437        let Some(EdgeRepr::Curve3d { curve, range, .. }) = data.curve3d() else {
438            return Ok(None);
439        };
440        let Some(curve) = model.geometry().curve(*curve) else {
441            ogeom_core::ogeom_bail!(Dangling, "an edge names a curve not in this model");
442        };
443        placed.push((curve, *range, edge.transform(model.datums())?));
444    }
445    Ok(Some(stretches_apart(placed[0], placed[1], tol)?))
446}
447
448/// How far apart two placed stretches of curves stand: the farthest that
449/// any of many samples along either lies from the nearest point of the
450/// other, as [`edges_apart`] measures two edges.
451pub(crate) fn stretches_apart(
452    a: (&Curve, (f64, f64), Transform),
453    b: (&Curve, (f64, f64), Transform),
454    tol: Tolerances,
455) -> OgeomResult<f64> {
456    const SAMPLES: u32 = 64;
457    let placed = [a, b];
458    let mut widest: f64 = 0.0;
459    for (from, to) in [(0, 1), (1, 0)] {
460        let (curve, range, placement) = placed[from];
461        let (other, other_range, other_placement) = placed[to];
462        // The sample brought into the other curve's own frame, and the
463        // distance found there scaled back as that placement scales.
464        let back = other_placement.inverse()?;
465        let scale = other_placement.scale_factor().abs();
466        for i in 0..=SAMPLES {
467            let t = range.0 + (range.1 - range.0) * f64::from(i) / f64::from(SAMPLES);
468            let at = back.apply(placement.apply(curve.point_at(t, tol)?));
469            widest = widest.max(nearest_on_stretch(other, other_range, at, tol)? * scale);
470        }
471    }
472    Ok(widest)
473}
474
475/// Raise each of `edges`' tolerance, and its vertices', to how far its
476/// pcurves stand off its curve, where that is wider than it states.
477///
478/// # Errors
479///
480/// As [`edge_pcurve_gap`].
481pub fn state_pcurve_gaps_of(
482    model: &mut Model,
483    edges: &[Shape],
484    tol: Tolerances,
485) -> OgeomResult<()> {
486    let mut wider = Vec::new();
487    for edge in edges {
488        let Some(data) = model.node(edge).and_then(|n| n.data().as_edge()) else {
489            continue;
490        };
491        let stated = data.tolerance.get();
492        if let Some(gap) = edge_pcurve_gap(model, edge, tol)?
493            && gap > stated
494        {
495            wider.push((edge.clone(), gap));
496        }
497    }
498    for (edge, gap) in wider {
499        model.widen(&edge, Tolerance::new(gap * (1.0 + MARGIN))?)?;
500    }
501    Ok(())
502}
503
504/// [`state_pcurve_gaps_of`] over every edge of `shape`.
505///
506/// # Errors
507///
508/// As [`edge_pcurve_gap`].
509pub fn state_pcurve_gaps(model: &mut Model, shape: &Shape, tol: Tolerances) -> OgeomResult<()> {
510    let edges = ogeom_topo::explore_unique(model, shape, ShapeType::Edge)?;
511    state_pcurve_gaps_of(model, &edges, tol)
512}
513
514#[cfg(test)]
515#[allow(clippy::unwrap_used, clippy::expect_used)]
516mod tests {
517    use super::*;
518    use ogeom_geom::{BSpline2d, LineCurve, PlaneSurface};
519    use ogeom_math::{KnotVector, Plane, Point2};
520
521    const T: Tolerances = Tolerances::millimetres();
522
523    /// A point of a whole circle just short of where it closes is on the
524    /// circle, though its start, the same point as its end, ties with
525    /// the end as the nearest sample.
526    #[test]
527    fn a_point_short_of_a_closed_curve_s_end_is_on_it() {
528        use ogeom_geom::CircleCurve;
529        use ogeom_math::{Circle, Frame};
530        let circle: Curve = CircleCurve::new(Circle::new(Frame::WORLD, 0.5, T).unwrap()).into();
531        let range = (0.0, core::f64::consts::TAU);
532        let target = circle.point_at(range.1 - 0.01, T).unwrap();
533        let nearest = nearest_on_stretch(&circle, range, target, T).unwrap();
534        assert!(nearest < 1e-9, "{nearest}");
535    }
536
537    /// A pcurve rising to a corner half way between two samples, along a
538    /// straight edge on the plane: the samples either side stand short of
539    /// the corner by a quarter of a percent, and the gap is the corner's
540    /// height, where a piece of the edge sampled at its corner finds it.
541    #[test]
542    fn a_peak_between_samples_is_measured_at_its_top() {
543        let curve: Curve = LineCurve::segment(Point::ORIGIN, Point::new(1.0, 0.0, 0.0), T)
544            .unwrap()
545            .into();
546        let surface: SurfaceGeometry =
547            PlaneSurface::new(Plane::new(ogeom_math::Frame::WORLD)).into();
548        let (corner, height) = (128.5 / 256.0, 1e-3);
549        let knots = KnotVector::new(vec![0.0, 0.0, corner, 1.0, 1.0], 1).unwrap();
550        let control = vec![
551            Point2::new(0.0, 0.0),
552            Point2::new(corner, height),
553            Point2::new(1.0, 0.0),
554        ];
555        let pcurve: PlanarCurve = BSpline2d::new(knots, control, T).unwrap().into();
556        let gap = pcurve_gap(
557            (&curve, (0.0, 1.0)),
558            (&pcurve, (0.0, 1.0)),
559            &surface,
560            1e-7,
561            T,
562        )
563        .unwrap();
564        assert!(
565            (gap - height).abs() < height * 1e-6,
566            "{gap} against {height}"
567        );
568    }
569}
ogeom_algo/pcurve_gap.rs

ogeom_algo/
pcurve_gap.rs
