geop_core_geometry/nurb_curve/
derivative.rs

1use geop_core_math::{
2    geop_error::{GeopError, GeopResult},
3    scalars::Scalar,
4    vector::Vector,
5};
6
7use super::NurbCurve;
8use crate::aabb::compute_aabb;
9
10impl<S: Scalar, const D: usize> NurbCurve<S, D> {
11    /// Return the derivative of this curve as a new `NurbCurve` of degree `p − 1`.
12    ///
13    /// The returned curve is the derivative of the *homogeneous* B-spline.
14    /// Evaluating it at `t` yields the homogeneous tangent vector.  To obtain
15    /// the Cartesian tangent `C′(t)`, apply the quotient rule:
16    ///
17    /// ```text
18    /// C′(t) = (A′(t) − w′(t)·C(t)) / w(t)
19    /// ```
20    ///
21    /// Errors if `self.degree == 0`.
22    pub fn derivative(&self) -> GeopResult<NurbCurve<S, D>> {
23        let p = self.degree;
24        if p == 0 {
25            return Err(GeopError::new(
26                "derivative is undefined for a degree-0 curve",
27            ));
28        }
29
30        let n = self.control_points.len() - 1;
31        let u = &self.knot_vector;
32        let p_s = S::from_i64(p as i64);
33
34        let mut deriv_pts: Vec<Vector<S, D>> = Vec::with_capacity(n);
35        for i in 0..n {
36            let den = u[i + p + 1].sub(u[i + 1]);
37            let mut d = Vector::<S, D>::zero();
38            if den.definitely_not_equal(S::ZERO) {
39                let scale = p_s.div(den).unwrap_or(S::ZERO);
40                for c in 0..D {
41                    d[c] = self.control_points[i + 1][c]
42                        .sub(self.control_points[i][c])
43                        .mul(scale);
44                }
45            }
46            deriv_pts.push(d);
47        }
48
49        let deriv_knots: Vec<S> = u[1..=n + p].to_vec();
50
51        // Computed the same way every other constructor here does (see
52        // `compute_aabb`'s own doc comment) — this matches what
53        // `NurbCurve::evaluate` already does with a derivative curve's own
54        // "weight" component (it has no separate rational-vs-derivative
55        // special case; it just divides by the last coordinate like it
56        // would for any other curve), so this stays consistent with that
57        // rather than falling back to a match-anything placeholder that
58        // would silently disable any future caller relying on `.aabb`
59        // (not just today's intersection-search prefilters).
60        let aabb = compute_aabb(&deriv_pts);
61
62        Ok(NurbCurve {
63            degree: p - 1,
64            control_points: deriv_pts,
65            knot_vector: deriv_knots,
66            aabb,
67        })
68    }
69}
70
71#[cfg(test)]
72mod tests {
73    use crate::nurb_curve::NurbCurve;
74    use geop_core_math::for_all_scalars;
75    use geop_core_math::{scalars::Scalar, vector::Vector4};
76
77    fn pt<S: Scalar>(x: f64, y: f64, z: f64, w: f64) -> Vector4<S> {
78        Vector4::from_array([
79            S::from_f64(x),
80            S::from_f64(y),
81            S::from_f64(z),
82            S::from_f64(w),
83        ])
84    }
85
86    fn check_degree_zero_errors<S: Scalar>() {
87        let f = S::from_f64;
88        let c = NurbCurve::try_new(
89            1,
90            vec![pt(0., 0., 0., 1.), pt(1., 0., 0., 1.)],
91            vec![f(0.), f(0.), f(1.), f(1.)],
92        )
93        .unwrap();
94        let d = c.derivative().unwrap();
95        assert!(d.derivative().is_err());
96    }
97    #[test]
98    fn degree_zero_errors() {
99        for_all_scalars!(check_degree_zero_errors);
100    }
101
102    fn check_structural_properties<S: Scalar>() {
103        let f = S::from_f64;
104        let c = NurbCurve::try_new(
105            3,
106            vec![
107                pt(0., 0., 0., 1.),
108                pt(1., 1., 0., 1.),
109                pt(2., 1., 0., 1.),
110                pt(3., 0., 0., 1.),
111            ],
112            vec![f(0.), f(0.), f(0.), f(0.), f(1.), f(1.), f(1.), f(1.)],
113        )
114        .unwrap();
115
116        let d1 = c.derivative().unwrap();
117        assert_eq!(d1.degree, 2);
118        assert_eq!(d1.control_points.len(), 3);
119        assert_eq!(
120            d1.knot_vector.len(),
121            d1.control_points.len() + d1.degree + 1
122        );
123
124        let d2 = d1.derivative().unwrap();
125        assert_eq!(d2.degree, 1);
126        assert_eq!(d2.control_points.len(), 2);
127        assert_eq!(
128            d2.knot_vector.len(),
129            d2.control_points.len() + d2.degree + 1
130        );
131
132        let d3 = d2.derivative().unwrap();
133        assert_eq!(d3.degree, 0);
134        assert_eq!(d3.control_points.len(), 1);
135        assert_eq!(
136            d3.knot_vector.len(),
137            d3.control_points.len() + d3.degree + 1
138        );
139    }
140    #[test]
141    fn structural_properties() {
142        for_all_scalars!(check_structural_properties);
143    }
144
145    fn check_quadratic_derivative_control_points<S: Scalar>() {
146        let f = S::from_f64;
147        let c = NurbCurve::try_new(
148            2,
149            vec![pt(0., 0., 0., 1.), pt(0.5, 0.5, 0., 1.), pt(1., 0., 0., 1.)],
150            vec![f(0.), f(0.), f(0.), f(1.), f(1.), f(1.)],
151        )
152        .unwrap();
153
154        let d = c.derivative().unwrap();
155        assert_eq!(d.degree, 1);
156        assert_eq!(d.control_points.len(), 2);
157
158        assert!(d.control_points[0][0].could_be_equal(S::ONE));
159        assert!(d.control_points[0][1].could_be_equal(S::ONE));
160        assert!(d.control_points[0][2].could_be_equal(S::ZERO));
161        assert!(d.control_points[0][3].could_be_equal(S::ZERO));
162
163        assert!(d.control_points[1][0].could_be_equal(S::ONE));
164        assert!(d.control_points[1][1].could_be_equal(S::ONE.neg()));
165        assert!(d.control_points[1][2].could_be_equal(S::ZERO));
166        assert!(d.control_points[1][3].could_be_equal(S::ZERO));
167    }
168    #[test]
169    fn quadratic_derivative_control_points() {
170        for_all_scalars!(check_quadratic_derivative_control_points);
171    }
172
173    fn check_knot_vector_is_trimmed<S: Scalar>() {
174        let f = S::from_f64;
175        let c = NurbCurve::try_new(
176            2,
177            vec![pt(0., 0., 0., 1.), pt(0.5, 0.5, 0., 1.), pt(1., 0., 0., 1.)],
178            vec![f(0.), f(0.), f(0.), f(1.), f(1.), f(1.)],
179        )
180        .unwrap();
181
182        let d = c.derivative().unwrap();
183        let expected = [0., 0., 1., 1.];
184        for (got, &exp) in d.knot_vector.iter().zip(expected.iter()) {
185            assert!(
186                got.could_be_equal(f(exp)),
187                "expected {exp}, got {}",
188                got.to_f64()
189            );
190        }
191    }
192    #[test]
193    fn knot_vector_is_trimmed() {
194        for_all_scalars!(check_knot_vector_is_trimmed);
195    }
196
197    fn check_weighted_derivative_is_evaluatable<S: Scalar>() {
198        let f = S::from_f64;
199        let c = NurbCurve::try_new(
200            1,
201            vec![pt(0., 0., 0., 1.), pt(1., 0., 0., 2.)],
202            vec![f(0.), f(0.), f(1.), f(1.)],
203        )
204        .unwrap();
205
206        let d = c.derivative().unwrap();
207        assert_eq!(d.degree, 0);
208
209        assert!(d.control_points[0][0].could_be_equal(S::ONE));
210        assert!(d.control_points[0][3].could_be_equal(S::ONE));
211
212        let a_prime = d.evaluate(f(0.5)).unwrap();
213        assert!(a_prime[0].could_be_equal(S::ONE));
214        assert!(a_prime[1].could_be_equal(S::ZERO));
215    }
216    #[test]
217    fn weighted_derivative_is_evaluatable() {
218        for_all_scalars!(check_weighted_derivative_is_evaluatable);
219    }
220}
geop_core_geometry/nurb_curve/derivative.rs

geop_core_geometry/nurb_curve/
derivative.rs
