geop_core_geometry/nurb_surface/
evaluate.rs

1use geop_core_math::{
2    geop_error::{GeopError, GeopResult},
3    scalars::Scalar,
4    vector::{Vector, Vector3},
5};
6
7use super::NurbSurface;
8
9/// Find the knot span: last index k in [degree, n] where knots[k] <= t < knots[k+1].
10pub(super) fn find_span<S: Scalar>(
11    degree: usize,
12    knots: &[S],
13    n: usize,
14    t: S,
15) -> GeopResult<usize> {
16    let p = degree;
17    if t.definitely_less(knots[p]) || t.definitely_greater(knots[n + 1]) {
18        return Err(GeopError::new(format!(
19            "parameter t={} out of domain [{}, {}]",
20            t,
21            knots[p],
22            knots[n + 1]
23        )));
24    }
25    if !t.definitely_less(knots[n + 1]) {
26        for k in (p..=n).rev() {
27            if knots[k].definitely_less(knots[n + 1]) {
28                return Ok(k);
29            }
30        }
31        return Ok(p);
32    }
33    for k in p..=n {
34        if !t.definitely_less(knots[k]) && t.definitely_less(knots[k + 1]) {
35            return Ok(k);
36        }
37    }
38    Err(GeopError::new("could not find knot span"))
39}
40
41/// De Boor triangular recursion in homogeneous space; generic over CP dimension D.
42pub(super) fn de_boor<S: Scalar, const D: usize>(
43    degree: usize,
44    knots: &[S],
45    points: &[Vector<S, D>],
46    t: S,
47    span: usize,
48) -> Vector<S, D> {
49    let p = degree;
50    let mut d: Vec<Vector<S, D>> = (0..=p).map(|j| points[span - p + j]).collect();
51    for r in 1..=p {
52        for j in (r..=p).rev() {
53            let i = span - p + j;
54            let denom = knots[i + p - r + 1].sub(knots[i]);
55            let alpha = if denom.could_be_equal(S::ZERO) {
56                S::ZERO
57            } else {
58                t.sub(knots[i]).div(denom).unwrap_or(S::ZERO)
59            };
60            d[j] = Vector::interpolate(&d[j - 1], &d[j], alpha);
61        }
62    }
63    d[p]
64}
65
66// ── 3-D surface ──────────────────────────────────────────────────────────────
67
68impl<S: Scalar> NurbSurface<S, 4> {
69    /// Evaluate the surface at `(u, v)`, returning a 3-D Cartesian point.
70    pub fn evaluate(&self, u: S, v: S) -> GeopResult<Vector3<S>> {
71        let p = self.degree_u;
72        let q = self.degree_v;
73        let nu = self.num_u;
74        let nv = self.num_v;
75
76        let span_u = find_span(p, &self.knot_vector_u, nu - 1, u)?;
77
78        let mut col_pts: Vec<Vector<S, 4>> = Vec::with_capacity(nv);
79        for j in 0..nv {
80            let row: Vec<Vector<S, 4>> = (0..nu).map(|i| self.control_points[i * nv + j]).collect();
81            col_pts.push(de_boor(p, &self.knot_vector_u, &row, u, span_u));
82        }
83
84        let span_v = find_span(q, &self.knot_vector_v, nv - 1, v)?;
85        let hw = de_boor(q, &self.knot_vector_v, &col_pts, v, span_v);
86
87        let w = hw[3];
88        if w.could_be_equal(S::ZERO) {
89            return Err(GeopError::new("weight is zero at evaluation point"));
90        }
91        let inv_w = S::ONE.div(w)?;
92        let mut result = Vector3::zero();
93        for c in 0..3 {
94            result[c] = hw[c].mul(inv_w);
95        }
96        Ok(result)
97    }
98}
99
100impl<S: Scalar> geop_core_math::primitives::scene::RasterizableSurface<S> for NurbSurface<S, 4> {
101    fn eval_at(&self, u: S, v: S) -> GeopResult<Vector3<S>> {
102        self.evaluate(u, v)
103    }
104}
105
106#[cfg(test)]
107mod tests {
108    use super::super::NurbSurface;
109    use geop_core_math::for_all_scalars;
110    use geop_core_math::{scalars::Scalar, vector::Vector4};
111
112    fn pt<S: Scalar>(x: f64, y: f64, z: f64, w: f64) -> Vector4<S> {
113        Vector4::from_array([
114            S::from_f64(x),
115            S::from_f64(y),
116            S::from_f64(z),
117            S::from_f64(w),
118        ])
119    }
120
121    fn bilinear<S: Scalar>() -> NurbSurface<S, 4> {
122        let f = S::from_f64;
123        NurbSurface::try_new(
124            1,
125            1,
126            vec![
127                pt(0., 0., 0., 1.),
128                pt(0., 1., 0., 1.),
129                pt(1., 0., 0., 1.),
130                pt(1., 1., 0., 1.),
131            ],
132            vec![f(0.), f(0.), f(1.), f(1.)],
133            vec![f(0.), f(0.), f(1.), f(1.)],
134        )
135        .unwrap()
136    }
137
138    fn check_bilinear_corner_00<S: Scalar>() {
139        let s = bilinear::<S>();
140        let p = s.evaluate(S::ZERO, S::ZERO).unwrap();
141        assert!(p[0].could_be_equal(S::ZERO));
142        assert!(p[1].could_be_equal(S::ZERO));
143        assert!(p[2].could_be_equal(S::ZERO));
144    }
145    #[test]
146    fn bilinear_corner_00() {
147        for_all_scalars!(check_bilinear_corner_00);
148    }
149
150    fn check_bilinear_corner_10<S: Scalar>() {
151        let s = bilinear::<S>();
152        let p = s.evaluate(S::ONE, S::ZERO).unwrap();
153        assert!(p[0].could_be_equal(S::ONE));
154        assert!(p[1].could_be_equal(S::ZERO));
155    }
156    #[test]
157    fn bilinear_corner_10() {
158        for_all_scalars!(check_bilinear_corner_10);
159    }
160
161    fn check_bilinear_corner_01<S: Scalar>() {
162        let s = bilinear::<S>();
163        let p = s.evaluate(S::ZERO, S::ONE).unwrap();
164        assert!(p[0].could_be_equal(S::ZERO));
165        assert!(p[1].could_be_equal(S::ONE));
166    }
167    #[test]
168    fn bilinear_corner_01() {
169        for_all_scalars!(check_bilinear_corner_01);
170    }
171
172    fn check_bilinear_corner_11<S: Scalar>() {
173        let s = bilinear::<S>();
174        let p = s.evaluate(S::ONE, S::ONE).unwrap();
175        assert!(p[0].could_be_equal(S::ONE));
176        assert!(p[1].could_be_equal(S::ONE));
177    }
178    #[test]
179    fn bilinear_corner_11() {
180        for_all_scalars!(check_bilinear_corner_11);
181    }
182
183    fn check_bilinear_center<S: Scalar>() {
184        let s = bilinear::<S>();
185        let p = s.evaluate(S::from_f64(0.5), S::from_f64(0.5)).unwrap();
186        assert!(p[0].could_be_equal(S::from_f64(0.5)));
187        assert!(p[1].could_be_equal(S::from_f64(0.5)));
188        assert!(p[2].could_be_equal(S::ZERO));
189    }
190    #[test]
191    fn bilinear_center() {
192        for_all_scalars!(check_bilinear_center);
193    }
194
195    fn check_bilinear_mid_u_edge<S: Scalar>() {
196        let s = bilinear::<S>();
197        let p = s.evaluate(S::from_f64(0.5), S::ZERO).unwrap();
198        assert!(p[0].could_be_equal(S::from_f64(0.5)));
199        assert!(p[1].could_be_equal(S::ZERO));
200    }
201    #[test]
202    fn bilinear_mid_u_edge() {
203        for_all_scalars!(check_bilinear_mid_u_edge);
204    }
205
206    fn check_out_of_domain_u_returns_err<S: Scalar>() {
207        let s = bilinear::<S>();
208        assert!(s.evaluate(S::from_f64(-0.1), S::from_f64(0.5)).is_err());
209        assert!(s.evaluate(S::from_f64(1.1), S::from_f64(0.5)).is_err());
210    }
211    #[test]
212    fn out_of_domain_u_returns_err() {
213        for_all_scalars!(check_out_of_domain_u_returns_err);
214    }
215
216    fn check_out_of_domain_v_returns_err<S: Scalar>() {
217        let s = bilinear::<S>();
218        assert!(s.evaluate(S::from_f64(0.5), S::from_f64(-0.1)).is_err());
219        assert!(s.evaluate(S::from_f64(0.5), S::from_f64(1.1)).is_err());
220    }
221    #[test]
222    fn out_of_domain_v_returns_err() {
223        for_all_scalars!(check_out_of_domain_v_returns_err);
224    }
225
226    fn check_quadratic_u_midpoint<S: Scalar>() {
227        let f = S::from_f64;
228        let s = NurbSurface::try_new(
229            2,
230            1,
231            vec![
232                pt(0., 0., 0., 1.),
233                pt(0., 1., 0., 1.),
234                pt(0.5, 0., 1., 1.),
235                pt(0.5, 1., 1., 1.),
236                pt(1., 0., 0., 1.),
237                pt(1., 1., 0., 1.),
238            ],
239            vec![f(0.), f(0.), f(0.), f(1.), f(1.), f(1.)],
240            vec![f(0.), f(0.), f(1.), f(1.)],
241        )
242        .unwrap();
243        let p = s.evaluate(f(0.5), f(0.)).unwrap();
244        assert!(p[0].could_be_equal(f(0.5)));
245        assert!(p[1].could_be_equal(f(0.)));
246        assert!(p[2].could_be_equal(f(0.5)));
247    }
248    #[test]
249    fn quadratic_u_midpoint() {
250        for_all_scalars!(check_quadratic_u_midpoint);
251    }
252
253    fn check_everything_matches_any_point<S: Scalar>() {
254        let s = NurbSurface::<S, 4>::everything();
255        let p = s.evaluate(S::from_f64(0.5), S::from_f64(-3.0)).unwrap();
256        assert!(p[0].could_be_equal(S::from_f64(123.456)));
257        assert!(p[1].could_be_equal(S::from_f64(-9.0)));
258        assert!(p[2].could_be_equal(S::ZERO));
259    }
260    #[test]
261    fn everything_matches_any_point() {
262        for_all_scalars!(check_everything_matches_any_point);
263    }
264}
geop_core_geometry/nurb_surface/evaluate.rs

geop_core_geometry/nurb_surface/
evaluate.rs
