geop_core_geometry/nurb_curve/
tangent.rs

1use geop_core_math::{
2    geop_error::GeopResult,
3    scalars::Scalar,
4    vector::{Vector2, Vector3},
5};
6
7use super::NurbCurve;
8
9impl<S: Scalar> NurbCurve<S, 4> {
10    /// Cartesian tangent vector `C'(t)` (not normalized).
11    ///
12    /// Evaluates the homogeneous derivative curve from [`NurbCurve::derivative`]
13    /// and applies the quotient rule `C'(t) = (A'(t) − w'(t)·C(t)) / w(t)`,
14    /// where `(A(t), w(t))` is `self`'s own homogeneous point at `t`.
15    pub fn tangent(&self, t: S) -> GeopResult<Vector3<S>> {
16        let pos = self.evaluate(t)?;
17
18        let span = self.find_knot_span(t)?;
19        let hw = self.de_boor(t, span);
20        let w = hw[3];
21
22        let deriv = self.derivative()?;
23        let dspan = deriv.find_knot_span(t)?;
24        let dhw = deriv.de_boor(t, dspan);
25
26        let mut result = Vector3::zero();
27        for c in 0..3 {
28            result[c] = dhw[c].sub(dhw[3].mul(pos[c])).div(w)?;
29        }
30        Ok(result)
31    }
32
33    /// Cartesian second derivative `C''(t)` (not normalized).
34    ///
35    /// Applies the quotient rule twice: `C''(t) = (A''(t) − 2·w'(t)·C'(t) −
36    /// w''(t)·C(t)) / w(t)`, where `(A(t), w(t))` is `self`'s own homogeneous
37    /// point at `t` and `A''(t), w''(t)` come from differentiating the
38    /// homogeneous curve twice via [`NurbCurve::derivative`].
39    ///
40    /// A curve of degree `< 2` has no well-defined second derivative of its
41    /// *homogeneous* representation (differentiating twice underflows the
42    /// degree) — but geometrically a degree-`0`/`1` curve is a straight
43    /// line, whose Cartesian second derivative is genuinely zero, so that's
44    /// what's returned instead of an error.
45    pub fn second_derivative(&self, t: S) -> GeopResult<Vector3<S>> {
46        if self.degree < 2 {
47            return Ok(Vector3::zero());
48        }
49
50        let pos = self.evaluate(t)?;
51        let tan = self.tangent(t)?;
52
53        let span = self.find_knot_span(t)?;
54        let w = self.de_boor(t, span)[3];
55
56        let d1 = self.derivative()?;
57        let d1span = d1.find_knot_span(t)?;
58        let wp = d1.de_boor(t, d1span)[3];
59
60        let d2 = d1.derivative()?;
61        let d2span = d2.find_knot_span(t)?;
62        let d2hw = d2.de_boor(t, d2span);
63        let wpp = d2hw[3];
64
65        let mut result = Vector3::zero();
66        for c in 0..3 {
67            result[c] = d2hw[c]
68                .sub(wp.mul(S::TWO).mul(tan[c]))
69                .sub(wpp.mul(pos[c]))
70                .div(w)?;
71        }
72        Ok(result)
73    }
74}
75
76impl<S: Scalar> NurbCurve<S, 3> {
77    /// Cartesian tangent vector `C'(t)` (not normalized), for a 2-D pcurve.
78    ///
79    /// Same quotient-rule derivation as the 3-D [`NurbCurve::<S, 4>::tangent`].
80    pub fn tangent(&self, t: S) -> GeopResult<Vector2<S>> {
81        let pos = self.evaluate(t)?;
82
83        let span = self.find_knot_span(t)?;
84        let hw = self.de_boor(t, span);
85        let w = hw[2];
86
87        let deriv = self.derivative()?;
88        let dspan = deriv.find_knot_span(t)?;
89        let dhw = deriv.de_boor(t, dspan);
90
91        let mut result = Vector2::zero();
92        for c in 0..2 {
93            result[c] = dhw[c].sub(dhw[2].mul(pos[c])).div(w)?;
94        }
95        Ok(result)
96    }
97}
98
99#[cfg(test)]
100mod tests {
101    use crate::nurb_curve::NurbCurve;
102    use geop_core_math::for_all_scalars;
103    use geop_core_math::{scalars::Scalar, vector::Vector4};
104
105    fn pt<S: Scalar>(x: f64, y: f64, z: f64, w: f64) -> Vector4<S> {
106        Vector4::from_array([
107            S::from_f64(x),
108            S::from_f64(y),
109            S::from_f64(z),
110            S::from_f64(w),
111        ])
112    }
113
114    fn check_line_tangent_is_constant_direction<S: Scalar>() {
115        let f = S::from_f64;
116        let c = NurbCurve::try_new(
117            1,
118            vec![pt(0., 0., 0., 1.), pt(2., 0., 0., 1.)],
119            vec![f(0.), f(0.), f(1.), f(1.)],
120        )
121        .unwrap();
122
123        for &t in &[0.0, 0.25, 0.5, 0.75, 1.0] {
124            let tan = c.tangent(f(t)).unwrap();
125            assert!(tan[0].could_be_equal(S::TWO));
126            assert!(tan[1].could_be_equal(S::ZERO));
127            assert!(tan[2].could_be_equal(S::ZERO));
128        }
129    }
130    #[test]
131    fn line_tangent_is_constant_direction() {
132        for_all_scalars!(check_line_tangent_is_constant_direction);
133    }
134
135    /// Quadratic Bézier `(0,0,0) -> (1,1,0) -> (2,0,0)`: tangent at t=0.5
136    /// (the apex) should be purely horizontal (dy/dt = 0), pointing +x.
137    fn check_quadratic_apex_tangent_is_horizontal<S: Scalar>() {
138        let f = S::from_f64;
139        let c = NurbCurve::try_new(
140            2,
141            vec![pt(0., 0., 0., 1.), pt(1., 1., 0., 1.), pt(2., 0., 0., 1.)],
142            vec![f(0.), f(0.), f(0.), f(1.), f(1.), f(1.)],
143        )
144        .unwrap();
145
146        let tan = c.tangent(f(0.5)).unwrap();
147        assert!(tan[0].definitely_greater(S::ZERO));
148        assert!(tan[1].could_be_equal(S::ZERO));
149        assert!(tan[2].could_be_equal(S::ZERO));
150    }
151    #[test]
152    fn quadratic_apex_tangent_is_horizontal() {
153        for_all_scalars!(check_quadratic_apex_tangent_is_horizontal);
154    }
155
156    /// At the start of that same Bézier, the tangent should point up and to
157    /// the right (+x, +y), matching the initial control-polygon direction.
158    fn check_quadratic_start_tangent_direction<S: Scalar>() {
159        let f = S::from_f64;
160        let c = NurbCurve::try_new(
161            2,
162            vec![pt(0., 0., 0., 1.), pt(1., 1., 0., 1.), pt(2., 0., 0., 1.)],
163            vec![f(0.), f(0.), f(0.), f(1.), f(1.), f(1.)],
164        )
165        .unwrap();
166
167        let tan = c.tangent(f(0.0)).unwrap();
168        assert!(tan[0].definitely_greater(S::ZERO));
169        assert!(tan[1].definitely_greater(S::ZERO));
170    }
171    #[test]
172    fn quadratic_start_tangent_direction() {
173        for_all_scalars!(check_quadratic_start_tangent_direction);
174    }
175
176    // ── 2-D (pcurve) tangent ──────────────────────────────────────────────────
177
178    fn pt2<S: Scalar>(x: f64, y: f64, w: f64) -> geop_core_math::vector::Vector3<S> {
179        geop_core_math::vector::Vector3::from_array([
180            S::from_f64(x),
181            S::from_f64(y),
182            S::from_f64(w),
183        ])
184    }
185
186    fn check_line_2d_tangent_is_constant_direction<S: Scalar>() {
187        let f = S::from_f64;
188        let c: crate::nurb_curve::NurbCurve2D<S> = NurbCurve::try_new(
189            1,
190            vec![pt2(0., 0., 1.), pt2(2., 0., 1.)],
191            vec![f(0.), f(0.), f(1.), f(1.)],
192        )
193        .unwrap();
194
195        for &t in &[0.0, 0.25, 0.5, 0.75, 1.0] {
196            let tan = c.tangent(f(t)).unwrap();
197            assert!(tan[0].could_be_equal(S::TWO));
198            assert!(tan[1].could_be_equal(S::ZERO));
199        }
200    }
201    #[test]
202    fn line_2d_tangent_is_constant_direction() {
203        for_all_scalars!(check_line_2d_tangent_is_constant_direction);
204    }
205
206    /// Quadratic Bézier `(0,0) -> (1,1) -> (2,0)` in 2-D: tangent at t=0.5
207    /// (the apex) should be purely horizontal (dy/dt = 0), pointing +x.
208    fn check_quadratic_2d_apex_tangent_is_horizontal<S: Scalar>() {
209        let f = S::from_f64;
210        let c: crate::nurb_curve::NurbCurve2D<S> = NurbCurve::try_new(
211            2,
212            vec![pt2(0., 0., 1.), pt2(1., 1., 1.), pt2(2., 0., 1.)],
213            vec![f(0.), f(0.), f(0.), f(1.), f(1.), f(1.)],
214        )
215        .unwrap();
216
217        let tan = c.tangent(f(0.5)).unwrap();
218        assert!(tan[0].definitely_greater(S::ZERO));
219        assert!(tan[1].could_be_equal(S::ZERO));
220    }
221    #[test]
222    fn quadratic_2d_apex_tangent_is_horizontal() {
223        for_all_scalars!(check_quadratic_2d_apex_tangent_is_horizontal);
224    }
225}
geop_core_geometry/nurb_curve/tangent.rs

geop_core_geometry/nurb_curve/
tangent.rs
