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;
8use crate::spline::{centered, de_boor, find_span, nearer_end};
9
10/// The index, among the `(p + 1) × (q + 1)` control points acting on
11/// `(u, v)` in `u` index major order, of the corner nearest `(u, v)`: the
12/// one to evaluate relative to (see [`centered`]).
13pub(super) fn corner_nearest(
14    surface: &NurbSurface<4>,
15    (span_u, u): (usize, Scalar),
16    (span_v, v): (usize, Scalar),
17) -> usize {
18    let (p, q) = (surface.degree_u, surface.degree_v);
19    let i = if nearer_end(&surface.knot_vector_u, span_u, u) {
20        p
21    } else {
22        0
23    };
24    let j = if nearer_end(&surface.knot_vector_v, span_v, v) {
25        q
26    } else {
27        0
28    };
29    i * (q + 1) + j
30}
31
32// ── 3-D surface ──────────────────────────────────────────────────────────────
33
34impl NurbSurface<4> {
35    /// The weight of the surface's homogeneous form at `(u, v)`: its
36    /// rational basis's denominator there.
37    pub(crate) fn weight(&self, u: Scalar, v: Scalar) -> GeopResult<Scalar> {
38        let (p, q) = (self.degree_u, self.degree_v);
39        let nv = self.num_v;
40        let span_u = find_span(p, &self.knot_vector_u, self.num_u - 1, u)?;
41        let span_v = find_span(q, &self.knot_vector_v, nv - 1, v)?;
42        let cols: Vec<Vector<1>> = (span_v - q..=span_v)
43            .map(|j| {
44                let col: Vec<Vector<1>> = (span_u - p..=span_u)
45                    .map(|i| Vector::from_array([self.control_points[i * nv + j][3]]))
46                    .collect();
47                de_boor(p, &self.knot_vector_u[span_u - p..], &col, u, p)
48            })
49            .collect();
50        Ok(de_boor(q, &self.knot_vector_v[span_v - q..], &cols, v, q)[0])
51    }
52
53    /// Evaluate the surface at `(u, v)`, returning a 3-D Cartesian point.
54    ///
55    /// Relative to one of the control points acting there ([`centered`]),
56    /// so an interval `(u, v)` gives a point as wide as the patch is, not as
57    /// far from the origin as it lies.
58    ///
59    /// Over a box of `(u, v)` that is also the mean value form, the point
60    /// at the box's middle plus the partials over the box times the half
61    /// widths, and the two intersected: both enclose the same points. Taken
62    /// directly, `A / w` over a box divides two enclosures that cannot be
63    /// correlated, and the quotient widens with how far the point lies from
64    /// the span's control point times how much the weight changes: a
65    /// cylinder evaluated at a foot point a micron wide came out ten
66    /// microns wide along its axis, and a shell's copy of an edge each time
67    /// eight times as wide as the edge.
68    pub fn evaluate(&self, u: Scalar, v: Scalar) -> GeopResult<Vector3> {
69        let direct = self.evaluate_at(u, v)?;
70        if u.is_sharp() && v.is_sharp() {
71            return Ok(direct);
72        }
73        // Only ever tighter: where the partials over the box cannot be had,
74        // as at a pole or past the domain, the direct enclosure stands.
75        let (um, vm) = (u.midpoint(), v.midpoint());
76        let (Ok(middle), Ok((su, sv))) = (self.evaluate_at(um, vm), self.derivatives(u, v)) else {
77            return Ok(direct);
78        };
79        let (du, dv) = (u.sub(um), v.sub(vm));
80        Ok(Vector3::from_array(std::array::from_fn(|c| {
81            let mean = middle[c].add(su[c].mul(du)).add(sv[c].mul(dv));
82            direct[c].intersect(mean)
83        })))
84    }
85
86    /// [`Self::evaluate`] directly: the homogeneous point over `(u, v)`
87    /// divided by its weight there.
88    fn evaluate_at(&self, u: Scalar, v: Scalar) -> GeopResult<Vector3> {
89        let (p, q) = (self.degree_u, self.degree_v);
90        let nv = self.num_v;
91        let span_u = find_span(p, &self.knot_vector_u, self.num_u - 1, u)?;
92        let span_v = find_span(q, &self.knot_vector_v, nv - 1, v)?;
93
94        // The `(p + 1) × (q + 1)` control points acting on the span, `u`
95        // index major; each column is evaluated at `u`, then those at `v`.
96        let local: Vec<Vector<4>> = (span_u - p..=span_u)
97            .flat_map(|i| (span_v - q..=span_v).map(move |j| i * nv + j))
98            .map(|k| self.control_points[k])
99            .collect();
100        let at = corner_nearest(self, (span_u, u), (span_v, v));
101        let (local, origin) = centered(&local, at);
102        let cols: Vec<Vector<4>> = (0..=q)
103            .map(|j| {
104                let col: Vec<Vector<4>> = (0..=p).map(|i| local[i * (q + 1) + j]).collect();
105                de_boor(p, &self.knot_vector_u[span_u - p..], &col, u, p)
106            })
107            .collect();
108        let hw = de_boor(q, &self.knot_vector_v[span_v - q..], &cols, v, q);
109
110        let w = hw[3];
111        if w.could_be_equal(Scalar::ZERO) {
112            return Err(GeopError::new(format!(
113                "weight is zero at evaluation point (u, v) = ({u:?}, {v:?}) in the domain \
114                 {:?} x {:?} (homogeneous point {hw:?})",
115                self.domain_u(),
116                self.domain_v()
117            )));
118        }
119        let inv_w = Scalar::ONE.div(w)?;
120        let mut result = Vector3::zero();
121        for c in 0..3 {
122            result[c] = origin[c].add(hw[c].mul(inv_w));
123        }
124        Ok(result)
125    }
126}
127
128#[cfg(test)]
129mod tests {
130    use super::super::NurbSurface;
131    use crate::test_fixtures::{bilinear, hom};
132    use geop_core_math::scalars::Scalar;
133
134    #[test]
135    fn bilinear_corner_00() {
136        let s = bilinear();
137        let p = s.evaluate(Scalar::ZERO, Scalar::ZERO).unwrap();
138        assert!(p[0].could_be_equal(Scalar::ZERO));
139        assert!(p[1].could_be_equal(Scalar::ZERO));
140        assert!(p[2].could_be_equal(Scalar::ZERO));
141    }
142
143    #[test]
144    fn bilinear_corner_10() {
145        let s = bilinear();
146        let p = s.evaluate(Scalar::ONE, Scalar::ZERO).unwrap();
147        assert!(p[0].could_be_equal(Scalar::ONE));
148        assert!(p[1].could_be_equal(Scalar::ZERO));
149    }
150
151    #[test]
152    fn bilinear_corner_01() {
153        let s = bilinear();
154        let p = s.evaluate(Scalar::ZERO, Scalar::ONE).unwrap();
155        assert!(p[0].could_be_equal(Scalar::ZERO));
156        assert!(p[1].could_be_equal(Scalar::ONE));
157    }
158
159    #[test]
160    fn bilinear_corner_11() {
161        let s = bilinear();
162        let p = s.evaluate(Scalar::ONE, Scalar::ONE).unwrap();
163        assert!(p[0].could_be_equal(Scalar::ONE));
164        assert!(p[1].could_be_equal(Scalar::ONE));
165    }
166
167    #[test]
168    fn bilinear_center() {
169        let s = bilinear();
170        let p = s
171            .evaluate(Scalar::from_f64(0.5), Scalar::from_f64(0.5))
172            .unwrap();
173        assert!(p[0].could_be_equal(Scalar::from_f64(0.5)));
174        assert!(p[1].could_be_equal(Scalar::from_f64(0.5)));
175        assert!(p[2].could_be_equal(Scalar::ZERO));
176    }
177
178    #[test]
179    fn bilinear_mid_u_edge() {
180        let s = bilinear();
181        let p = s.evaluate(Scalar::from_f64(0.5), Scalar::ZERO).unwrap();
182        assert!(p[0].could_be_equal(Scalar::from_f64(0.5)));
183        assert!(p[1].could_be_equal(Scalar::ZERO));
184    }
185
186    #[test]
187    fn out_of_domain_u_returns_err() {
188        let s = bilinear();
189        assert!(
190            s.evaluate(Scalar::from_f64(-0.1), Scalar::from_f64(0.5))
191                .is_err()
192        );
193        assert!(
194            s.evaluate(Scalar::from_f64(1.1), Scalar::from_f64(0.5))
195                .is_err()
196        );
197    }
198
199    #[test]
200    fn out_of_domain_v_returns_err() {
201        let s = bilinear();
202        assert!(
203            s.evaluate(Scalar::from_f64(0.5), Scalar::from_f64(-0.1))
204                .is_err()
205        );
206        assert!(
207            s.evaluate(Scalar::from_f64(0.5), Scalar::from_f64(1.1))
208                .is_err()
209        );
210    }
211
212    #[test]
213    fn quadratic_u_midpoint() {
214        let f = Scalar::from_f64;
215        let s = NurbSurface::try_new(
216            2,
217            1,
218            vec![
219                hom(0., 0., 0., 1.),
220                hom(0., 1., 0., 1.),
221                hom(0.5, 0., 1., 1.),
222                hom(0.5, 1., 1., 1.),
223                hom(1., 0., 0., 1.),
224                hom(1., 1., 0., 1.),
225            ],
226            vec![f(0.), f(0.), f(0.), f(1.), f(1.), f(1.)],
227            vec![f(0.), f(0.), f(1.), f(1.)],
228        )
229        .unwrap();
230        let p = s.evaluate(f(0.5), f(0.)).unwrap();
231        assert!(p[0].could_be_equal(f(0.5)));
232        assert!(p[1].could_be_equal(f(0.)));
233        assert!(p[2].could_be_equal(f(0.5)));
234    }
235
236    /// A quarter of a cylinder's cap (radius 3, at height 50), evaluated
237    /// over a wide parameter box near its centre: every point of it is at
238    /// height 50, and so is the enclosure — as narrow as at the origin.
239    /// Evaluated as `A / W` in absolute coordinates it spanned `[38, 94]`
240    /// in fixed point, and a cylinder 50 up along its own axis was found
241    /// to overlap its own other cap.
242    #[test]
243    fn a_cap_far_from_the_origin_stays_flat() {
244        let f = Scalar::from_f64;
245        let h = std::f64::consts::FRAC_1_SQRT_2;
246        let z = 50.0;
247        // u round the quarter circle, v from the centre out.
248        let corners = [(1.0, 0.0, 1.0), (1.0, 1.0, h), (0.0, 1.0, 1.0)];
249        let points = corners
250            .iter()
251            .flat_map(|&(x, y, w)| {
252                [
253                    hom(0.0, 0.0, z * w, w),
254                    hom(3.0 * x * w, 3.0 * y * w, z * w, w),
255                ]
256            })
257            .collect();
258        let s = NurbSurface::try_new(
259            2,
260            1,
261            points,
262            vec![f(0.), f(0.), f(0.), f(1.), f(1.), f(1.)],
263            vec![f(0.), f(0.), f(1.), f(1.)],
264        )
265        .unwrap();
266        let (u, v) = (f(0.33).union(f(0.83)), f(0.0005).union(f(0.0008)));
267        for (u, v) in [(u, v), (f(0.5), f(0.5))] {
268            let p = s.evaluate(u, v).unwrap();
269            assert!(p[2].could_be_equal(f(z)), "{p:?}");
270            assert!(p[2].width().to_f64() < 1e-6, "{p:?}");
271        }
272    }
273
274    #[test]
275    fn everything_matches_any_point() {
276        let s = NurbSurface::<4>::everything();
277        let p = s
278            .evaluate(Scalar::from_f64(0.5), Scalar::from_f64(-3.0))
279            .unwrap();
280        assert!(p[0].could_be_equal(Scalar::from_f64(123.456)));
281        assert!(p[1].could_be_equal(Scalar::from_f64(-9.0)));
282        assert!(p[2].could_be_equal(Scalar::ZERO));
283    }
284}
geop_core_geometry/nurb_surface/evaluate.rs

geop_core_geometry/nurb_surface/
evaluate.rs
