geop_core_geometry/nurb_curve/
offset.rs

1//! How far a curve fitted to a plane curve's offset is from it, proven.
2//!
3//! The offset of a plane NURBS curve `S` by `d`, `O = S + d N` with `N` its
4//! unit normal, is no NURBS: a spline fitted to it only comes near it. How
5//! near is bounded here, box by box over the fit's domain, by the mean
6//! value form of the miss `e = F − O`: its value in the middle of the box,
7//! plus its derivative over the box times the half width. That derivative
8//! is the fit's less the offset's, `O' = S' + d N'` with
9//! `N' = (S''⊥ |S'|² − S'⊥ (S'·S'')) / |S'|³`, all enclosed over the box:
10//! the bound is honest wherever it is, and tight once the boxes are small
11//! against the fit's own error.
12
13use geop_core_math::{
14    geop_error::{GeopError, GeopResult},
15    scalars::Scalar,
16    vector::Vector2,
17};
18
19use super::NurbCurve2D;
20
21/// How many times a box of the fit's domain is split in two at most: the
22/// bound converges with the square of the box's width.
23const MISS_SPLITS: usize = 12;
24
25/// `v` turned a quarter counter-clockwise.
26fn left(v: Vector2) -> Vector2 {
27    Vector2::from_array([v[1].neg(), v[0]])
28}
29
30impl NurbCurve2D {
31    /// A bound on how far `fit`, over its domain, is from this curve's
32    /// offset by `distance` to its left — `fit` at the start of its domain
33    /// standing for the offset of this curve at `t0`, at its end for that
34    /// at `t1`, affinely between — proven over boxes of `fit`'s domain,
35    /// each split in two while its bound exceeds `tolerance` (at most
36    /// [`MISS_SPLITS`] times; then the bound is only wider) and its middle
37    /// is not definitely further off than that. Fails where the curve has
38    /// no tangent.
39    pub fn offset_miss(
40        &self,
41        distance: Scalar,
42        (t0, t1): (Scalar, Scalar),
43        fit: &NurbCurve2D,
44        tolerance: Scalar,
45    ) -> GeopResult<Scalar> {
46        let ctx = |e: GeopError| {
47            e.with_context(format!(
48                "NurbCurve2D::offset_miss(distance={distance:?}, t=({t0:?}, {t1:?}))"
49            ))
50        };
51        let (s0, s1) = fit.domain();
52        let rate = t1.sub(t0).div(s1.sub(s0)).map_err(ctx)?;
53        let along = |s: Scalar| t0.add(s.sub(s0).mul(rate));
54        let offset_at = |t: Scalar| -> GeopResult<Vector2> {
55            let normal = left(self.tangent(t)?).normalize()?;
56            Ok(self.evaluate(t)?.add(&normal.prod_scalar(distance)))
57        };
58        // The bound over the box from `a` to `b`, and whether its middle
59        // alone is off by more than `tolerance`: then splitting it cannot
60        // bring it under.
61        let miss_over = |a: Scalar, b: Scalar| -> GeopResult<(Scalar, bool)> {
62            let middle = a.add(b).div(Scalar::TWO)?.sharpen();
63            let here = fit.evaluate(middle)?.sub(&offset_at(along(middle))?);
64            let s = a.union(b);
65            let t = along(s);
66            let (first, second) = (self.tangent(t)?, self.second_derivative(t)?);
67            let speed_sq = first.prod_dot(&first);
68            let speed = speed_sq.sqrt()?;
69            let turning = left(second)
70                .prod_scalar(speed_sq)
71                .sub(&left(first).prod_scalar(first.prod_dot(&second)))
72                .prod_scalar(Scalar::ONE.div(speed_sq.mul(speed))?);
73            let slope = fit
74                .tangent(s)?
75                .sub(&first.add(&turning.prod_scalar(distance)).prod_scalar(rate));
76            let half = b.sub(a).div(Scalar::TWO)?;
77            let reach = half.neg().union(half);
78            let e = here.add(&slope.prod_scalar(reach));
79            let (x, y) = (e[0].abs(), e[1].abs());
80            let off = here.prod_dot(&here).sqrt()?.definitely_greater(tolerance);
81            Ok((x.mul(x).add(y.mul(y)).sqrt()?.upper(), off))
82        };
83        let breaks = fit.breakpoints();
84        let mut stack: Vec<(Scalar, Scalar, usize)> =
85            breaks.windows(2).map(|w| (w[0], w[1], 0)).collect();
86        let mut worst = Scalar::ZERO;
87        while let Some((a, b, depth)) = stack.pop() {
88            // Over a box too wide its enclosures are too: the speed could be
89            // zero over it, and the bound is not decided there yet.
90            let miss = match miss_over(a, b) {
91                Err(_) if depth < MISS_SPLITS => None,
92                miss => Some(miss.map_err(ctx)?),
93            };
94            if miss.is_none_or(|(m, off)| m.definitely_greater(tolerance) && !off)
95                && depth < MISS_SPLITS
96            {
97                let middle = a.add(b).div(Scalar::TWO).map_err(ctx)?.sharpen();
98                stack.push((a, middle, depth + 1));
99                stack.push((middle, b, depth + 1));
100                continue;
101            }
102            worst = worst.max(miss.expect("decided at the last split").0);
103        }
104        Ok(worst)
105    }
106}
107
108#[cfg(test)]
109mod tests {
110    use super::*;
111    use crate::test_fixtures::knots;
112
113    fn h(x: f64, y: f64, w: f64) -> geop_core_math::vector::Vector3 {
114        geop_core_math::vector::Vector3::from_array([x * w, y * w, w].map(Scalar::from_f64))
115    }
116
117    /// A quarter circle of radius 2 and its offset by 0.5 inward, the
118    /// quarter circle of radius 1.5, exactly: missed by no more than the
119    /// tolerance; and the same offset by 0.6, by 0.1 at least.
120    #[test]
121    fn an_arc_offset_is_missed_by_no_more_than_it_is_off() {
122        let r2 = std::f64::consts::FRAC_1_SQRT_2;
123        let arc = |r: f64| {
124            NurbCurve2D::try_new(
125                2,
126                vec![h(r, 0., 1.), h(r, r, r2), h(0., r, 1.)],
127                knots(&[0., 0., 0., 1., 1., 1.]),
128            )
129            .unwrap()
130        };
131        let domain = (Scalar::ZERO, Scalar::ONE);
132        let tolerance = Scalar::from_f64(1e-6);
133        let exact = arc(2.).offset_miss(Scalar::from_f64(0.5), domain, &arc(1.5), tolerance);
134        let exact = exact.unwrap().to_f64();
135        assert!(exact < 1e-6, "{exact}");
136        let off = arc(2.)
137            .offset_miss(Scalar::from_f64(0.6), domain, &arc(1.5), tolerance)
138            .unwrap()
139            .to_f64();
140        assert!(off >= 0.1, "{off}");
141    }
142}
geop_core_geometry/nurb_curve/offset.rs

geop_core_geometry/nurb_curve/
offset.rs
