geop_core_geometry/nurb_curve/
tangent.rs

1use geop_core_math::{
2    geop_error::{GeopError, GeopResult},
3    scalars::Scalar,
4    vector::{Vector, Vector2, Vector3},
5};
6
7use super::NurbCurve;
8use crate::spline::{
9    centered, find_spans, homogeneous_derivatives, nearer_end, rational_derivatives,
10};
11
12impl<const D: usize> NurbCurve<D> {
13    /// The Cartesian point and its derivatives up to order `n` at `t`:
14    /// `[C(t), C'(t), …, C⁽ⁿ⁾(t)]`, with `C = D − 1`. Evaluated directly,
15    /// without building any derivative curve ([`homogeneous_derivatives`],
16    /// then the quotient rule [`rational_derivatives`]).
17    ///
18    /// The derivatives are taken twice and intersected, as
19    /// [`NurbCurve::evaluate`] is: from the control points as they are, and
20    /// relative to one of them ([`centered`]). They do not depend on where
21    /// the curve lies; the first's width does, far from the origin, and the
22    /// second costs a rounding that near it is the wider one.
23    ///
24    /// Over every knot span `t` reaches, united (see [`find_spans`]).
25    pub(super) fn cartesian_derivatives<const C: usize>(
26        &self,
27        t: Scalar,
28        n: usize,
29    ) -> GeopResult<Vec<Vector<C>>> {
30        let p = self.degree;
31        let last = self.control_points.len() - 1;
32        let mut out: Option<Vec<Vector<C>>> = None;
33        for (span, piece) in find_spans(p, &self.knot_vector, last, t)? {
34            let local = &self.control_points[span - p..=span];
35            let at = |points: &[Vector<D>]| {
36                rational_derivatives::<D, C>(&homogeneous_derivatives(
37                    p,
38                    &self.knot_vector,
39                    points,
40                    span,
41                    piece,
42                    n,
43                ))
44            };
45            let mut derivatives = at(local)?;
46            let near = if nearer_end(&self.knot_vector, span, piece) {
47                p
48            } else {
49                0
50            };
51            let relative = at(&centered(local, near).0)?;
52            for (d, r) in derivatives.iter_mut().zip(&relative).skip(1) {
53                for c in 0..C {
54                    d[c] = d[c].intersect(r[c]);
55                }
56            }
57            out = Some(match out {
58                Some(other) => other
59                    .iter()
60                    .zip(&derivatives)
61                    .map(|(a, b)| a.union(b))
62                    .collect(),
63                None => derivatives,
64            });
65        }
66        out.ok_or_else(|| GeopError::new(format!("NurbCurve::tangent(t={t:?}): no knot span")))
67    }
68}
69
70impl NurbCurve<4> {
71    /// Cartesian tangent vector `C'(t)` (not normalized).
72    pub fn tangent(&self, t: Scalar) -> GeopResult<Vector3> {
73        Ok(self.cartesian_derivatives::<3>(t, 1)?[1])
74    }
75
76    /// Cartesian second derivative `C''(t)` (not normalized). Nonzero for a
77    /// rational degree-1 curve too: a line with unequal weights is not
78    /// traversed at constant speed.
79    pub fn second_derivative(&self, t: Scalar) -> GeopResult<Vector3> {
80        Ok(self.cartesian_derivatives::<3>(t, 2)?[2])
81    }
82}
83
84impl NurbCurve<3> {
85    /// Cartesian tangent vector `C'(t)` (not normalized), for a 2-D pcurve.
86    pub fn tangent(&self, t: Scalar) -> GeopResult<Vector2> {
87        Ok(self.cartesian_derivatives::<2>(t, 1)?[1])
88    }
89
90    /// Cartesian second derivative `C''(t)` (not normalized), for a 2-D
91    /// curve.
92    pub fn second_derivative(&self, t: Scalar) -> GeopResult<Vector2> {
93        Ok(self.cartesian_derivatives::<2>(t, 2)?[2])
94    }
95}
96
97#[cfg(test)]
98mod tests {
99    use crate::nurb_curve::NurbCurve;
100    use crate::test_fixtures::hom;
101    use geop_core_math::scalars::Scalar;
102
103    #[test]
104    fn line_tangent_is_constant_direction() {
105        let f = Scalar::from_f64;
106        let c = NurbCurve::try_new(
107            1,
108            vec![hom(0., 0., 0., 1.), hom(2., 0., 0., 1.)],
109            vec![f(0.), f(0.), f(1.), f(1.)],
110        )
111        .unwrap();
112
113        for &t in &[0.0, 0.25, 0.5, 0.75, 1.0] {
114            let tan = c.tangent(f(t)).unwrap();
115            assert!(tan[0].could_be_equal(Scalar::TWO));
116            assert!(tan[1].could_be_equal(Scalar::ZERO));
117            assert!(tan[2].could_be_equal(Scalar::ZERO));
118        }
119    }
120
121    /// Quadratic Bézier `(0,0,0) -> (1,1,0) -> (2,0,0)`: tangent at t=0.5
122    /// (the apex) should be purely horizontal (dy/dt = 0), pointing +x.
123    #[test]
124    fn quadratic_apex_tangent_is_horizontal() {
125        let f = Scalar::from_f64;
126        let c = NurbCurve::try_new(
127            2,
128            vec![
129                hom(0., 0., 0., 1.),
130                hom(1., 1., 0., 1.),
131                hom(2., 0., 0., 1.),
132            ],
133            vec![f(0.), f(0.), f(0.), f(1.), f(1.), f(1.)],
134        )
135        .unwrap();
136
137        let tan = c.tangent(f(0.5)).unwrap();
138        assert!(tan[0].definitely_greater(Scalar::ZERO));
139        assert!(tan[1].could_be_equal(Scalar::ZERO));
140        assert!(tan[2].could_be_equal(Scalar::ZERO));
141    }
142
143    /// At the start of that same Bézier, the tangent should point up and to
144    /// the right (+x, +y), matching the initial control-polygon direction.
145    #[test]
146    fn quadratic_start_tangent_direction() {
147        let f = Scalar::from_f64;
148        let c = NurbCurve::try_new(
149            2,
150            vec![
151                hom(0., 0., 0., 1.),
152                hom(1., 1., 0., 1.),
153                hom(2., 0., 0., 1.),
154            ],
155            vec![f(0.), f(0.), f(0.), f(1.), f(1.), f(1.)],
156        )
157        .unwrap();
158
159        let tan = c.tangent(f(0.0)).unwrap();
160        assert!(tan[0].definitely_greater(Scalar::ZERO));
161        assert!(tan[1].definitely_greater(Scalar::ZERO));
162    }
163
164    /// Degree-1 line with weights 1 and 2: `x(t) = 2t / (1 + t)`, so
165    /// `x' = 2 / (1+t)²` and `x'' = −4 / (1+t)³` — nonzero although the
166    /// curve is straight.
167    #[test]
168    fn rational_line_derivatives() {
169        let f = Scalar::from_f64;
170        let c = NurbCurve::try_new(
171            1,
172            vec![hom(0., 0., 0., 1.), hom(2., 0., 0., 2.)],
173            vec![f(0.), f(0.), f(1.), f(1.)],
174        )
175        .unwrap();
176        for t in [0., 0.3, 1.] {
177            let d1 = c.tangent(f(t)).unwrap();
178            let d2 = c.second_derivative(f(t)).unwrap();
179            assert!(d1[0].could_be_equal(f(2. / (1. + t).powi(2))), "{d1:?}");
180            assert!(d2[0].could_be_equal(f(-4. / (1. + t).powi(3))), "{d2:?}");
181            assert!(d1[1].could_be_equal(Scalar::ZERO) && d2[1].could_be_equal(Scalar::ZERO));
182        }
183    }
184
185    /// On the exact rational unit quarter circle `|C| = 1`, so
186    /// `C·C' = 0` and `C·C'' + |C'|² = 0` for every `t`.
187    #[test]
188    fn quarter_circle_derivatives() {
189        let f = Scalar::from_f64;
190        let w = std::f64::consts::FRAC_1_SQRT_2;
191        let c = NurbCurve::try_new(
192            2,
193            vec![hom(1., 0., 0., 1.), hom(w, w, 0., w), hom(0., 1., 0., 1.)],
194            vec![f(0.), f(0.), f(0.), f(1.), f(1.), f(1.)],
195        )
196        .unwrap();
197        for t in [0., 0.2, 0.5, 0.9, 1.] {
198            let p = c.evaluate(f(t)).unwrap();
199            let d1 = c.tangent(f(t)).unwrap();
200            let d2 = c.second_derivative(f(t)).unwrap();
201            assert!(p.prod_dot(&d1).could_be_equal(Scalar::ZERO), "t={t}");
202            assert!(
203                p.prod_dot(&d2)
204                    .add(d1.prod_dot(&d1))
205                    .could_be_equal(Scalar::ZERO),
206                "t={t}: {d1:?} {d2:?}"
207            );
208            assert!(d1.norm().definitely_greater(Scalar::ZERO));
209        }
210    }
211
212    // ── 2-D (pcurve) tangent ──────────────────────────────────────────────────
213
214    fn pt2(x: f64, y: f64, w: f64) -> geop_core_math::vector::Vector3 {
215        geop_core_math::vector::Vector3::from_array([
216            Scalar::from_f64(x),
217            Scalar::from_f64(y),
218            Scalar::from_f64(w),
219        ])
220    }
221
222    #[test]
223    fn line_2d_tangent_is_constant_direction() {
224        let f = Scalar::from_f64;
225        let c: crate::nurb_curve::NurbCurve2D = NurbCurve::try_new(
226            1,
227            vec![pt2(0., 0., 1.), pt2(2., 0., 1.)],
228            vec![f(0.), f(0.), f(1.), f(1.)],
229        )
230        .unwrap();
231
232        for &t in &[0.0, 0.25, 0.5, 0.75, 1.0] {
233            let tan = c.tangent(f(t)).unwrap();
234            assert!(tan[0].could_be_equal(Scalar::TWO));
235            assert!(tan[1].could_be_equal(Scalar::ZERO));
236        }
237    }
238
239    /// Quadratic Bézier `(0,0) -> (1,1) -> (2,0)` in 2-D: tangent at t=0.5
240    /// (the apex) should be purely horizontal (dy/dt = 0), pointing +x.
241    #[test]
242    fn quadratic_2d_apex_tangent_is_horizontal() {
243        let f = Scalar::from_f64;
244        let c: crate::nurb_curve::NurbCurve2D = NurbCurve::try_new(
245            2,
246            vec![pt2(0., 0., 1.), pt2(1., 1., 1.), pt2(2., 0., 1.)],
247            vec![f(0.), f(0.), f(0.), f(1.), f(1.), f(1.)],
248        )
249        .unwrap();
250
251        let tan = c.tangent(f(0.5)).unwrap();
252        assert!(tan[0].definitely_greater(Scalar::ZERO));
253        assert!(tan[1].could_be_equal(Scalar::ZERO));
254    }
255}
geop_core_geometry/nurb_curve/tangent.rs

geop_core_geometry/nurb_curve/
tangent.rs
