geop_core_geometry/nurb_surface/
split.rs

1use geop_core_math::{geop_error::GeopResult, scalars::Scalar, vector::Vector};
2
3use super::NurbSurface;
4use crate::{aabb::compute_aabb, knot_insertion};
5
6/// A parameter direction of a tensor-product surface.
7#[derive(Clone, Copy)]
8enum Dir {
9    U,
10    V,
11}
12
13impl<S: Scalar, const D: usize> NurbSurface<S, D> {
14    fn degree_in(&self, dir: Dir) -> usize {
15        match dir {
16            Dir::U => self.degree_u,
17            Dir::V => self.degree_v,
18        }
19    }
20
21    fn knots_in(&self, dir: Dir) -> &Vec<S> {
22        match dir {
23            Dir::U => &self.knot_vector_u,
24            Dir::V => &self.knot_vector_v,
25        }
26    }
27
28    /// The control net as rows running along `dir` (see `knot_insertion`):
29    /// the `num_v` columns `P[·][j]` for `U`, the `num_u` rows `P[i][·]` for
30    /// `V`.
31    fn rows_along(&self, dir: Dir) -> Vec<Vec<Vector<S, D>>> {
32        let (nu, nv) = (self.num_u, self.num_v);
33        match dir {
34            Dir::U => (0..nv)
35                .map(|j| (0..nu).map(|i| self.control_points[i * nv + j]).collect())
36                .collect(),
37            Dir::V => self.control_points.chunks(nv).map(<[_]>::to_vec).collect(),
38        }
39    }
40
41    /// This surface with the `dir` knots and rows replaced — the inverse of
42    /// [`Self::rows_along`].
43    fn with_rows(&self, dir: Dir, knots: Vec<S>, rows: Vec<Vec<Vector<S, D>>>) -> Self {
44        let along = rows[0].len();
45        let (num_u, num_v, control_points): (usize, usize, Vec<Vector<S, D>>) = match dir {
46            Dir::U => (
47                along,
48                rows.len(),
49                (0..along)
50                    .flat_map(|i| rows.iter().map(move |col| col[i]))
51                    .collect(),
52            ),
53            Dir::V => (rows.len(), along, rows.into_iter().flatten().collect()),
54        };
55        let (knot_vector_u, knot_vector_v) = match dir {
56            Dir::U => (knots, self.knot_vector_v.clone()),
57            Dir::V => (self.knot_vector_u.clone(), knots),
58        };
59        NurbSurface {
60            degree_u: self.degree_u,
61            degree_v: self.degree_v,
62            num_u,
63            num_v,
64            aabb: compute_aabb(&control_points),
65            control_points,
66            knot_vector_u,
67            knot_vector_v,
68        }
69    }
70
71    fn split_in(&self, dir: Dir, t: S) -> GeopResult<(Self, Self)> {
72        let mut knots = self.knots_in(dir).clone();
73        let mut rows = self.rows_along(dir);
74        let (right_knots, right_rows) =
75            knot_insertion::split(&mut knots, &mut rows, self.degree_in(dir), t)?;
76        Ok((
77            self.with_rows(dir, knots, rows),
78            self.with_rows(dir, right_knots, right_rows),
79        ))
80    }
81
82    fn restrict_in(&self, dir: Dir, t0: S, t1: S) -> GeopResult<Self> {
83        let mut knots = self.knots_in(dir).clone();
84        let mut rows = self.rows_along(dir);
85        knot_insertion::restrict(&mut knots, &mut rows, self.degree_in(dir), t0, t1)?;
86        Ok(self.with_rows(dir, knots, rows))
87    }
88
89    /// Split at parameter `t` in the **u** direction.
90    ///
91    /// `t` must lie strictly inside the u domain.  Returns `(left, right)`.
92    ///
93    /// `t` is used exactly as given — **not** sharpened here, see `NurbCurve::split`'s own
94    /// doc comment for why (this is the same Boehm-insertion construction,
95    /// one dimension up: an unsharpened `t` carried into the new knot
96    /// vector lets `alpha = (t - e) / (s - e)` blow up over repeated splits
97    /// as `s - e` shrinks while `t`'s own width doesn't). Every current caller
98    /// therefore sharpens its own midpoint before calling; a caller splitting
99    /// at a *located* parameter must validate the sharpened value it is about
100    /// to use, not the wide one it started from.
101    pub fn split_u(&self, t: S) -> GeopResult<(NurbSurface<S, D>, NurbSurface<S, D>)> {
102        self.split_in(Dir::U, t)
103    }
104
105    /// Split at parameter `t` in the **v** direction.
106    ///
107    /// `t` must lie strictly inside the v domain.  Returns `(left, right)`.
108    /// `t` is used exactly as given — **not** sharpened here, same as [`Self::split_u`].
109    pub fn split_v(&self, t: S) -> GeopResult<(NurbSurface<S, D>, NurbSurface<S, D>)> {
110        self.split_in(Dir::V, t)
111    }
112
113    fn split_in_mid(&self, dir: Dir) -> GeopResult<(Self, Self)> {
114        let knots = self.knots_in(dir);
115        let n = match dir {
116            Dir::U => self.num_u,
117            Dir::V => self.num_v,
118        };
119        let (t0, t1) = (knots[self.degree_in(dir)], knots[n]);
120        // A self-chosen subdivision point: any value in the interval cuts
121        // it equally well, so sharpening loses no accuracy and keeps
122        // repeated splits from compounding width (see AGENTS.md).
123        let mid = t0.add(t1).div(S::TWO)?.sharpen();
124        self.split_in(dir, mid)
125    }
126
127    /// Split at the midpoint of the **u** domain.
128    pub fn split_u_mid(&self) -> GeopResult<(NurbSurface<S, D>, NurbSurface<S, D>)> {
129        self.split_in_mid(Dir::U)
130    }
131
132    /// Split at the midpoint of the **v** domain.
133    pub fn split_v_mid(&self) -> GeopResult<(NurbSurface<S, D>, NurbSurface<S, D>)> {
134        self.split_in_mid(Dir::V)
135    }
136
137    /// This surface restricted to `u ∈ [u0, u1]`, `v ∈ [v0, v1]`, with the
138    /// same cut rule as `NurbCurve::sub_curve`: only bounds strictly inside
139    /// the domain are cut at, so the result covers at least the requested
140    /// box ∩ domain. Bounds are used exactly as given and must be narrow.
141    pub fn sub_surface(&self, (u0, u1): (S, S), (v0, v1): (S, S)) -> GeopResult<Self> {
142        self.restrict_in(Dir::U, u0, u1)?
143            .restrict_in(Dir::V, v0, v1)
144    }
145}
146
147#[cfg(test)]
148mod tests {
149    use super::super::NurbSurface;
150    use geop_core_math::for_all_scalars;
151    use geop_core_math::{scalars::Scalar, vector::Vector4};
152
153    fn pt<S: Scalar>(x: f64, y: f64, z: f64, w: f64) -> Vector4<S> {
154        Vector4::from_array([
155            S::from_f64(x),
156            S::from_f64(y),
157            S::from_f64(z),
158            S::from_f64(w),
159        ])
160    }
161
162    fn bilinear<S: Scalar>() -> NurbSurface<S, 4> {
163        let f = S::from_f64;
164        NurbSurface::try_new(
165            1,
166            1,
167            vec![
168                pt(0., 0., 0., 1.),
169                pt(0., 1., 0., 1.),
170                pt(1., 0., 0., 1.),
171                pt(1., 1., 0., 1.),
172            ],
173            vec![f(0.), f(0.), f(1.), f(1.)],
174            vec![f(0.), f(0.), f(1.), f(1.)],
175        )
176        .unwrap()
177    }
178
179    fn check_split_u_junction_matches<S: Scalar>() {
180        let orig = bilinear::<S>();
181        let (left, right) = orig.split_u(S::from_f64(0.5)).unwrap();
182        for v_val in [0.0_f64, 0.5, 1.0] {
183            let v = S::from_f64(v_val);
184            let o = orig.evaluate(S::from_f64(0.5), v).unwrap();
185            let l = left.evaluate(S::from_f64(0.5), v).unwrap();
186            let r = right.evaluate(S::from_f64(0.5), v).unwrap();
187            for c in 0..3 {
188                assert!(
189                    o[c].could_be_equal(l[c]),
190                    "split_u left mismatch v={v_val} c={c}"
191                );
192                assert!(
193                    o[c].could_be_equal(r[c]),
194                    "split_u right mismatch v={v_val} c={c}"
195                );
196            }
197        }
198    }
199    #[test]
200    fn split_u_junction_matches() {
201        for_all_scalars!(check_split_u_junction_matches);
202    }
203
204    fn check_split_u_left_start_matches<S: Scalar>() {
205        let orig = bilinear::<S>();
206        let (left, _) = orig.split_u(S::from_f64(0.5)).unwrap();
207        for v_val in [0.0_f64, 0.5, 1.0] {
208            let v = S::from_f64(v_val);
209            let o = orig.evaluate(S::ZERO, v).unwrap();
210            let l = left.evaluate(S::ZERO, v).unwrap();
211            for c in 0..3 {
212                assert!(o[c].could_be_equal(l[c]), "c={c} v={v_val}");
213            }
214        }
215    }
216    #[test]
217    fn split_u_left_start_matches() {
218        for_all_scalars!(check_split_u_left_start_matches);
219    }
220
221    fn check_split_u_right_end_matches<S: Scalar>() {
222        let orig = bilinear::<S>();
223        let (_, right) = orig.split_u(S::from_f64(0.5)).unwrap();
224        for v_val in [0.0_f64, 0.5, 1.0] {
225            let v = S::from_f64(v_val);
226            let o = orig.evaluate(S::ONE, v).unwrap();
227            let r = right.evaluate(S::ONE, v).unwrap();
228            for c in 0..3 {
229                assert!(o[c].could_be_equal(r[c]), "c={c} v={v_val}");
230            }
231        }
232    }
233    #[test]
234    fn split_u_right_end_matches() {
235        for_all_scalars!(check_split_u_right_end_matches);
236    }
237
238    fn check_split_u_at_boundary_returns_err<S: Scalar>() {
239        let s = bilinear::<S>();
240        assert!(s.split_u(S::ZERO).is_err());
241        assert!(s.split_u(S::ONE).is_err());
242    }
243    #[test]
244    fn split_u_at_boundary_returns_err() {
245        for_all_scalars!(check_split_u_at_boundary_returns_err);
246    }
247
248    fn check_split_v_junction_matches<S: Scalar>() {
249        let orig = bilinear::<S>();
250        let (left, right) = orig.split_v(S::from_f64(0.5)).unwrap();
251        for u_val in [0.0_f64, 0.5, 1.0] {
252            let u = S::from_f64(u_val);
253            let o = orig.evaluate(u, S::from_f64(0.5)).unwrap();
254            let l = left.evaluate(u, S::from_f64(0.5)).unwrap();
255            let r = right.evaluate(u, S::from_f64(0.5)).unwrap();
256            for c in 0..3 {
257                assert!(
258                    o[c].could_be_equal(l[c]),
259                    "split_v left mismatch u={u_val} c={c}"
260                );
261                assert!(
262                    o[c].could_be_equal(r[c]),
263                    "split_v right mismatch u={u_val} c={c}"
264                );
265            }
266        }
267    }
268    #[test]
269    fn split_v_junction_matches() {
270        for_all_scalars!(check_split_v_junction_matches);
271    }
272
273    fn check_split_v_left_start_matches<S: Scalar>() {
274        let orig = bilinear::<S>();
275        let (left, _) = orig.split_v(S::from_f64(0.5)).unwrap();
276        for u_val in [0.0_f64, 0.5, 1.0] {
277            let u = S::from_f64(u_val);
278            let o = orig.evaluate(u, S::ZERO).unwrap();
279            let l = left.evaluate(u, S::ZERO).unwrap();
280            for c in 0..3 {
281                assert!(o[c].could_be_equal(l[c]), "c={c} u={u_val}");
282            }
283        }
284    }
285    #[test]
286    fn split_v_left_start_matches() {
287        for_all_scalars!(check_split_v_left_start_matches);
288    }
289
290    fn check_split_v_right_end_matches<S: Scalar>() {
291        let orig = bilinear::<S>();
292        let (_, right) = orig.split_v(S::from_f64(0.5)).unwrap();
293        for u_val in [0.0_f64, 0.5, 1.0] {
294            let u = S::from_f64(u_val);
295            let o = orig.evaluate(u, S::ONE).unwrap();
296            let r = right.evaluate(u, S::ONE).unwrap();
297            for c in 0..3 {
298                assert!(o[c].could_be_equal(r[c]), "c={c} u={u_val}");
299            }
300        }
301    }
302    #[test]
303    fn split_v_right_end_matches() {
304        for_all_scalars!(check_split_v_right_end_matches);
305    }
306
307    fn check_split_v_at_boundary_returns_err<S: Scalar>() {
308        let s = bilinear::<S>();
309        assert!(s.split_v(S::ZERO).is_err());
310        assert!(s.split_v(S::ONE).is_err());
311    }
312    #[test]
313    fn split_v_at_boundary_returns_err() {
314        for_all_scalars!(check_split_v_at_boundary_returns_err);
315    }
316
317    fn check_split_u_domains<S: Scalar>() {
318        let orig = bilinear::<S>();
319        let (left, right) = orig.split_u(S::from_f64(0.5)).unwrap();
320        let (lu_min, lu_max) = left.domain_u();
321        let (ru_min, ru_max) = right.domain_u();
322        assert!(lu_min.could_be_equal(S::ZERO));
323        assert!(lu_max.could_be_equal(S::from_f64(0.5)));
324        assert!(ru_min.could_be_equal(S::from_f64(0.5)));
325        assert!(ru_max.could_be_equal(S::ONE));
326    }
327    #[test]
328    fn split_u_domains() {
329        for_all_scalars!(check_split_u_domains);
330    }
331
332    fn check_split_v_domains<S: Scalar>() {
333        let orig = bilinear::<S>();
334        let (left, right) = orig.split_v(S::from_f64(0.5)).unwrap();
335        let (lv_min, lv_max) = left.domain_v();
336        let (rv_min, rv_max) = right.domain_v();
337        assert!(lv_min.could_be_equal(S::ZERO));
338        assert!(lv_max.could_be_equal(S::from_f64(0.5)));
339        assert!(rv_min.could_be_equal(S::from_f64(0.5)));
340        assert!(rv_max.could_be_equal(S::ONE));
341    }
342    #[test]
343    fn split_v_domains() {
344        for_all_scalars!(check_split_v_domains);
345    }
346}
geop_core_geometry/nurb_surface/split.rs

geop_core_geometry/nurb_surface/
split.rs
