geop_core_geometry/nurb_curve/
evaluate.rs

1use geop_core_math::{
2    geop_error::{GeopError, GeopResult, WithContext},
3    scalars::Scalar,
4    vector::{Vector2, Vector3},
5};
6
7use super::NurbCurve;
8use crate::{
9    knot_insertion::pinned_clamped_end,
10    spline::{centered, de_boor, find_spans, nearer_end},
11};
12
13impl<const D: usize> NurbCurve<D> {
14    /// The Cartesian point at `t`, its `D − 1` coordinates written to `out`:
15    /// the union of its value on every knot span `t` reaches (see
16    /// [`find_spans`]).
17    ///
18    /// On each span it is evaluated twice and intersected — both are
19    /// enclosures of the one point: `A / W` as it is, and relative to a
20    /// control point of the span ([`centered`]), `origin + A / W`. The first
21    /// divides two enclosures interval arithmetic cannot correlate, so its
22    /// width grows with `|A|`: an interval `t` gives a point as wide as the
23    /// curve is far from the origin, rather than as wide as the stretch of
24    /// curve it covers. The second is as wide as that stretch, but moving
25    /// away and back costs a rounding of the origin, which near the origin
26    /// is wider than the first — and wider than the control points, which
27    /// containment searches clip against, say the curve is.
28    ///
29    /// At a clamped end, pinned sharply, the curve is its end control point.
30    pub(super) fn evaluate_into(&self, t: Scalar, out: &mut [Scalar]) -> GeopResult<()> {
31        let p = self.degree;
32        let n = self.control_points.len();
33        let inverse = |w: Scalar, span: usize| -> GeopResult<Scalar> {
34            let ctx = |e: GeopError| {
35                e.with_context(format!(
36                    "NurbCurve::evaluate(t={t:?}, span={span}): degree={}, knot_vector={:?}, control_points={:?}",
37                    self.degree, self.knot_vector, self.control_points
38                ))
39            };
40            if w.could_be_equal(Scalar::ZERO) {
41                return Err(ctx(GeopError::new(format!(
42                    "weight is zero at evaluation point (homogeneous de_boor result w={w:?})"
43                ))));
44            }
45            Scalar::ONE.div(w).with_context(&ctx)
46        };
47        if let Some(first) = pinned_clamped_end(t, &self.knot_vector, n, p) {
48            let end = self.control_points[if first { 0 } else { n - 1 }];
49            let inv_w = inverse(end[D - 1], if first { p } else { n - 1 })?;
50            for (c, o) in out.iter_mut().enumerate() {
51                *o = end[c].mul(inv_w);
52            }
53            return Ok(());
54        }
55        let mut reached = false;
56        for (span, piece) in find_spans(p, &self.knot_vector, n - 1, t)? {
57            let local = &self.control_points[span - p..=span];
58            let knots = &self.knot_vector[span - p..];
59            let absolute = de_boor(p, knots, local, piece, p);
60            let near = if nearer_end(&self.knot_vector, span, piece) {
61                p
62            } else {
63                0
64            };
65            let (moved, origin) = centered(local, near);
66            let relative = de_boor(p, knots, &moved, piece, p);
67            // The weights are not moved: one `W` for both.
68            let inv_w = inverse(absolute[D - 1], span)?;
69            for (c, o) in out.iter_mut().enumerate() {
70                let near = origin[c].add(relative[c].mul(inv_w));
71                let value = absolute[c].mul(inv_w).intersect(near);
72                *o = if reached { o.union(value) } else { value };
73            }
74            reached = true;
75        }
76        if !reached {
77            return Err(GeopError::new(format!(
78                "NurbCurve::evaluate(t={t:?}): no knot span"
79            )));
80        }
81        Ok(())
82    }
83}
84
85impl<const D: usize> NurbCurve<D> {
86    /// The weight at `t`: the homogeneous coordinate the curve's point there
87    /// is divided by — the union of its value on every knot span `t`
88    /// reaches. A polynomial curve's is one everywhere; a rational one's says
89    /// how its parameter runs: its points, multiplied by it, are a
90    /// polynomial in `t` on each span, which is what makes a rational
91    /// quadratic an exact circle.
92    pub fn weight(&self, t: Scalar) -> GeopResult<Scalar> {
93        let p = self.degree;
94        let n = self.control_points.len();
95        let mut weight: Option<Scalar> = None;
96        for (span, piece) in find_spans(p, &self.knot_vector, n - 1, t)? {
97            let local = &self.control_points[span - p..=span];
98            let knots = &self.knot_vector[span - p..];
99            let w = de_boor(p, knots, local, piece, p)[D - 1];
100            weight = Some(weight.map_or(w, |u| u.union(w)));
101        }
102        weight.ok_or_else(|| GeopError::new(format!("NurbCurve::weight(t={t:?}): no knot span")))
103    }
104}
105
106// ── 3-D curve: evaluate returns Vector3 ─────────────────────────────────────
107
108impl NurbCurve<4> {
109    /// Evaluate the 3-D NURBS curve at `t`, returning a Cartesian `Vector3`.
110    pub fn evaluate(&self, t: Scalar) -> GeopResult<Vector3> {
111        let mut out = [Scalar::ZERO; 3];
112        self.evaluate_into(t, &mut out)?;
113        Ok(Vector3::from_array(out))
114    }
115
116    /// Where it starts and where it ends.
117    pub fn ends(&self) -> GeopResult<(Vector3, Vector3)> {
118        let (t0, t1) = self.domain();
119        Ok((self.evaluate(t0)?, self.evaluate(t1)?))
120    }
121}
122
123// ── 2-D curve (pcurve): evaluate returns Vector2 ─────────────────────────────
124
125impl NurbCurve<3> {
126    /// Evaluate the 2-D pcurve at `t`, returning a Cartesian `Vector2`.
127    pub fn evaluate(&self, t: Scalar) -> GeopResult<Vector2> {
128        let mut out = [Scalar::ZERO; 2];
129        self.evaluate_into(t, &mut out)?;
130        Ok(Vector2::from_array(out))
131    }
132
133    /// Where it starts and where it ends.
134    pub fn ends(&self) -> GeopResult<(Vector2, Vector2)> {
135        let (t0, t1) = self.domain();
136        Ok((self.evaluate(t0)?, self.evaluate(t1)?))
137    }
138}
139
140#[cfg(test)]
141mod tests {
142    use crate::nurb_curve::NurbCurve;
143    use crate::test_fixtures::hom;
144    use geop_core_math::{scalars::Scalar, vector::Vector3};
145
146    fn line_curve() -> NurbCurve<4> {
147        NurbCurve::try_new(
148            1,
149            vec![hom(0.0, 0.0, 0.0, 1.0), hom(1.0, 0.0, 0.0, 1.0)],
150            vec![
151                Scalar::from_f64(0.0),
152                Scalar::from_f64(0.0),
153                Scalar::from_f64(1.0),
154                Scalar::from_f64(1.0),
155            ],
156        )
157        .unwrap()
158    }
159
160    #[test]
161    fn line_at_start() {
162        let c = line_curve();
163        let p = c.evaluate(Scalar::ZERO).unwrap();
164        assert!(p[0].could_be_equal(Scalar::ZERO));
165        assert!(p[1].could_be_equal(Scalar::ZERO));
166    }
167
168    #[test]
169    fn line_at_end() {
170        let c = line_curve();
171        let p = c.evaluate(Scalar::ONE).unwrap();
172        assert!(p[0].could_be_equal(Scalar::ONE));
173        assert!(p[1].could_be_equal(Scalar::ZERO));
174    }
175
176    #[test]
177    fn line_at_midpoint() {
178        let c = line_curve();
179        let p = c.evaluate(Scalar::from_f64(0.5)).unwrap();
180        assert!(p[0].could_be_equal(Scalar::from_f64(0.5)));
181    }
182
183    #[test]
184    fn out_of_domain_returns_err() {
185        let c = line_curve();
186        assert!(c.evaluate(Scalar::from_f64(-0.1)).is_err());
187        assert!(c.evaluate(Scalar::from_f64(1.1)).is_err());
188    }
189
190    #[test]
191    fn quadratic_midpoint() {
192        let curve = NurbCurve::try_new(
193            2,
194            vec![
195                hom(0.0, 0.0, 0.0, 1.0),
196                hom(0.5, 0.0, 0.0, 1.0),
197                hom(1.0, 0.0, 0.0, 1.0),
198            ],
199            vec![
200                Scalar::from_f64(0.0),
201                Scalar::from_f64(0.0),
202                Scalar::from_f64(0.0),
203                Scalar::from_f64(1.0),
204                Scalar::from_f64(1.0),
205                Scalar::from_f64(1.0),
206            ],
207        )
208        .unwrap();
209        let p = curve.evaluate(Scalar::from_f64(0.5)).unwrap();
210        assert!(p[0].could_be_equal(Scalar::from_f64(0.5)));
211        assert!(p[1].could_be_equal(Scalar::ZERO));
212    }
213
214    /// Regression check for a `weight is zero at evaluation point` failure
215    /// observed from `remesh` on a degree-1, weight-1-constant, [0,0,1,1]
216    /// curve — data that on paper cannot produce a near-zero interpolated
217    /// weight (both endpoint weights are exactly 1). Reproduces the exact
218    /// control points/knots/`t` from that failure's error context to check
219    /// whether `evaluate` itself is at fault, independent of `remesh`.
220    #[test]
221    fn weight_one_constant_line_does_not_report_zero_weight() {
222        let curve = NurbCurve::try_new(
223            1,
224            vec![hom(0.5, 0.5, 0.5, 1.0), hom(0.5, 0.5, -0.5, 1.0)],
225            vec![
226                Scalar::from_f64(0.0),
227                Scalar::from_f64(0.0),
228                Scalar::from_f64(1.0),
229                Scalar::from_f64(1.0),
230            ],
231        )
232        .unwrap();
233        let p = curve.evaluate(Scalar::from_f64(0.972)).unwrap();
234        assert!(
235            p[2].could_be_equal(Scalar::from_f64(0.5 - 1.0 * 0.972)),
236            "p={p:?}"
237        );
238    }
239
240    /// A quarter circle of radius 3 in the plane `z = 50`, centred at
241    /// `(1000, -500)`: the kernel's rational arc, far from the origin.
242    fn far_arc() -> NurbCurve<4> {
243        let h = std::f64::consts::FRAC_1_SQRT_2;
244        let (cx, cy, z) = (1000.0, -500.0, 50.0);
245        NurbCurve::try_new(
246            2,
247            [(3.0, 0.0, 1.0), (3.0, 3.0, h), (0.0, 3.0, 1.0)]
248                .iter()
249                .map(|&(x, y, w)| hom((cx + x) * w, (cy + y) * w, z * w, w))
250                .collect(),
251            [0., 0., 0., 1., 1., 1.].map(Scalar::from_f64).to_vec(),
252        )
253        .unwrap()
254    }
255
256    /// Evaluated over a short parameter box, the arc far from the origin is
257    /// as wide as the stretch of it the box covers (about 0.05 long), and as
258    /// flat as the plane it lies in. Evaluated as `A / W` in absolute
259    /// coordinates its width grew with its distance from the origin instead.
260    #[test]
261    fn an_arc_far_from_the_origin_is_as_wide_as_its_stretch() {
262        let f = Scalar::from_f64;
263        let c = far_arc();
264        let t = f(0.33).union(f(0.34));
265        let p = c.evaluate(t).unwrap();
266        for k in 0..2 {
267            assert!(p[k].width().to_f64() < 0.2, "{p:?}");
268        }
269        assert!(p[2].could_be_equal(f(50.0)), "{p:?}");
270        assert!(p[2].width().to_f64() < 1e-6, "{p:?}");
271        // Its tangent too lies in the plane, and is as wide as it turns.
272        let d = c.tangent(t).unwrap();
273        assert!(d[2].could_be_equal(Scalar::ZERO), "{d:?}");
274        assert!(d[2].width().to_f64() < 1e-6, "{d:?}");
275        assert!(d[0].width().to_f64() < 1.0, "{d:?}");
276    }
277
278    /// The same for a pcurve: a quarter circle in `(u, v)` around
279    /// `(1000, -500)`, as on a surface parameterized in millimetres.
280    #[test]
281    fn a_pcurve_far_from_the_origin_is_as_wide_as_its_stretch() {
282        let f = Scalar::from_f64;
283        let h = std::f64::consts::FRAC_1_SQRT_2;
284        let c = NurbCurve::<3>::try_new(
285            2,
286            [(3.0, 0.0, 1.0), (3.0, 3.0, h), (0.0, 3.0, 1.0)]
287                .iter()
288                .map(|&(x, y, w)| {
289                    geop_core_math::vector::Vector3::from_array([
290                        f((1000.0 + x) * w),
291                        f((-500.0 + y) * w),
292                        f(w),
293                    ])
294                })
295                .collect(),
296            [0., 0., 0., 1., 1., 1.].map(f).to_vec(),
297        )
298        .unwrap();
299        let t = f(0.33).union(f(0.34));
300        let p = c.evaluate(t).unwrap();
301        for k in 0..2 {
302            assert!(p[k].width().to_f64() < 0.2, "{p:?}");
303        }
304        let d = c.tangent(t).unwrap();
305        assert!(d[0].width().to_f64() < 1.0, "{d:?}");
306    }
307
308    #[test]
309    fn everything_matches_any_point() {
310        let c = NurbCurve::<4>::everything();
311        let p = c.evaluate(Scalar::from_f64(7.0)).unwrap();
312        assert!(p[0].could_be_equal(Scalar::from_f64(123.456)));
313        assert!(p[1].could_be_equal(Scalar::from_f64(-9.0)));
314        assert!(p[2].could_be_equal(Scalar::ZERO));
315    }
316
317    /// Over an interval reaching across many knot spans, a curve is the union
318    /// of its spans' parts, not one span's cubic extrapolated over all of
319    /// them: a quarter circle of radius 2.25 interpolated through 49 points
320    /// came out 57 wide over its whole domain, and over an eighth of it
321    /// eight times as wide as that eighth.
322    #[test]
323    fn a_wide_parameter_encloses_the_curve_tightly() {
324        let r = 2.25;
325        let points: Vec<_> = (0..=48)
326            .map(|k| {
327                let a = std::f64::consts::FRAC_PI_2 * k as f64 / 48.0;
328                Vector3::from_array([r * a.cos(), r * a.sin(), 0.0].map(Scalar::from_f64))
329            })
330            .collect();
331        let curve = NurbCurve::<4>::interpolate(&points, 3).unwrap();
332        let (t0, t1) = curve.domain();
333        for (lo, hi) in [
334            (t0, t1),
335            (t0, Scalar::from_f64(0.125)),
336            (Scalar::from_f64(0.3), Scalar::from_f64(0.55)),
337        ] {
338            let box_ = curve.evaluate(lo.union(hi)).unwrap();
339            for k in 0..=64 {
340                let t = match k {
341                    0 => lo,
342                    64 => hi,
343                    _ => Scalar::interpolate(lo, hi, Scalar::from_ratio(k, 64).unwrap()).sharpen(),
344                };
345                let p = curve.evaluate(t).unwrap();
346                assert!(box_.could_be_equal(&p), "{t:?}: {p:?} outside {box_:?}");
347            }
348            // The arc between `lo` and `hi` spans at most its length on
349            // either axis; each span's part, evaluated over all of it in
350            // interval arithmetic, comes out somewhat wider than that.
351            let length = r * std::f64::consts::FRAC_PI_2 * (hi.to_f64() - lo.to_f64());
352            for c in 0..2 {
353                let width = box_[c].width().to_f64();
354                assert!(width < 1.5 * length, "{c}: {width} over {lo:?}..{hi:?}");
355            }
356        }
357    }
358}
geop_core_geometry/nurb_curve/evaluate.rs

geop_core_geometry/nurb_curve/
evaluate.rs
