geop_core_geometry/nurb_surface/
mod.rs

1mod carry;
2pub use carry::UvMap;
3mod closest;
4mod curvature;
5mod evaluate;
6mod fill;
7pub use fill::Fill;
8mod fit_pcurve;
9mod flat;
10mod join;
11pub use flat::FlatPatch;
12mod gordon;
13mod net;
14mod normal;
15mod offset;
16mod offset_fit;
17mod patch;
18pub use patch::{Neighbour, TILT};
19mod project;
20pub use project::MovingFoot;
21mod reverse;
22mod revolve;
23mod sphere;
24mod split;
25mod translate;
26
27use std::fmt::Display;
28
29use crate::aabb::compute_aabb;
30use geop_core_math::{
31    geop_error::{GeopError, GeopResult},
32    scalars::Scalar,
33    vector::Vector,
34};
35
36/// A NURBS surface patch whose control points live in `D`-dimensional homogeneous space.
37///
38/// In practice `D = 4` for all 3-D surfaces (`(wx, wy, wz, w)` control points).
39///
40/// Control points are stored row-major: `control_points[i * num_v + j]` is the
41/// point at u-index `i` (0 ≤ i < `num_u`) and v-index `j` (0 ≤ j < `num_v`).
42#[derive(Clone, Debug, PartialEq)]
43pub struct NurbSurface<const D: usize> {
44    pub degree_u: usize,
45    pub degree_v: usize,
46    pub num_u: usize,
47    pub num_v: usize,
48    pub control_points: Vec<Vector<D>>,
49    pub knot_vector_u: Vec<Scalar>,
50    pub knot_vector_v: Vec<Scalar>,
51    /// Cached axis-aligned bounding box of the (dehomogenized) control
52    /// points — see [`crate::aabb::compute_aabb`] and `NurbCurve`'s own
53    /// `aabb` field doc comment.
54    pub(crate) aabb: [Scalar; 3],
55}
56
57/// 3-D NURBS surface (homogeneous control points in ℝ⁴).
58pub type NurbSurface3D = NurbSurface<4>;
59
60impl<const D: usize> NurbSurface<D> {
61    pub fn try_new(
62        degree_u: usize,
63        degree_v: usize,
64        control_points: Vec<Vector<D>>,
65        knot_vector_u: Vec<Scalar>,
66        knot_vector_v: Vec<Scalar>,
67    ) -> GeopResult<Self> {
68        let len_u = knot_vector_u.len();
69        let len_v = knot_vector_v.len();
70
71        if len_u < degree_u + 2 {
72            return Err(GeopError::new(
73                "NurbSurface: knot_vector_u too short for the given degree",
74            ));
75        }
76        if len_v < degree_v + 2 {
77            return Err(GeopError::new(
78                "NurbSurface: knot_vector_v too short for the given degree",
79            ));
80        }
81
82        let num_u = len_u - degree_u - 1;
83        let num_v = len_v - degree_v - 1;
84
85        if control_points.len() != num_u * num_v {
86            return Err(GeopError::new(format!(
87                "NurbSurface: expected {} control points ({}×{}), got {}",
88                num_u * num_v,
89                num_u,
90                num_v,
91                control_points.len()
92            )));
93        }
94
95        // See `NurbCurve::try_new`'s identical check: every weight must be
96        // definitely positive, so `W(u, v) > 0` on the whole domain.
97        for p in &control_points {
98            if !p[D - 1].definitely_greater(Scalar::ZERO) {
99                return Err(GeopError::new(format!(
100                    "NurbSurface::try_new: control point {p:?} has a weight that is not definitely positive"
101                )));
102            }
103        }
104
105        let aabb = compute_aabb(&control_points);
106        Ok(Self {
107            degree_u,
108            degree_v,
109            num_u,
110            num_v,
111            control_points,
112            knot_vector_u,
113            knot_vector_v,
114            aabb,
115        })
116    }
117
118    /// Refresh the cached [`Self::aabb`] from the current `control_points`
119    /// — see [`crate::nurb_curve::NurbCurve::recompute_aabb`]'s identical
120    /// doc comment for why this is needed at all (`control_points` is
121    /// `pub`, and code outside this crate does mutate it in place).
122    pub fn recompute_aabb(&mut self) {
123        self.aabb = compute_aabb(&self.control_points);
124    }
125
126    /// This surface widened by `±pad` in every coordinate, everywhere: every
127    /// control point by that much (times its weight, in homogeneous form),
128    /// as [`crate::nurb_curve::NurbCurve::widen`] widens a curve — the
129    /// basis is non-negative and sums to one, so the widened patch holds
130    /// every point within `pad` of the patch along each axis.
131    pub fn widen(&mut self, pad: Scalar) {
132        for cp in &mut self.control_points {
133            let w = cp[D - 1];
134            for k in 0..D - 1 {
135                let d = pad.mul(w);
136                cp[k] = cp[k].sub(d).union(cp[k].add(d));
137            }
138        }
139        self.recompute_aabb();
140    }
141
142    /// The box around the control points, and so around the patch (see
143    /// [`crate::nurb_curve::NurbCurve::aabb`]).
144    pub fn aabb(&self) -> &[Scalar; 3] {
145        &self.aabb
146    }
147
148    /// Whether the two patches could share a point: `false` proves they do
149    /// not, as the boxes around their control points — which hold each
150    /// whole patch, by the convex hull property — are apart.
151    pub fn could_overlap(&self, other: &Self) -> bool {
152        crate::aabb::aabb_could_overlap(&self.aabb, &other.aabb, D - 1)
153    }
154
155    /// The diagonal of the box around the control points: by the convex
156    /// hull property, a length the whole patch fits within — the size of
157    /// the feature it is part of, for choices that should scale with it.
158    pub fn size(&self) -> GeopResult<Scalar> {
159        self.aabb
160            .iter()
161            .take(D - 1)
162            .fold(Scalar::ZERO, |sum, axis| {
163                sum.add(axis.width().mul(axis.width()))
164            })
165            .sqrt()
166    }
167
168    /// Valid parameter range in the u direction: `(u_min, u_max)`.
169    pub fn domain_u(&self) -> (Scalar, Scalar) {
170        (
171            self.knot_vector_u[self.degree_u],
172            self.knot_vector_u[self.num_u],
173        )
174    }
175
176    /// Valid parameter range in the v direction: `(v_min, v_max)`.
177    pub fn domain_v(&self) -> (Scalar, Scalar) {
178        (
179            self.knot_vector_v[self.degree_v],
180            self.knot_vector_v[self.num_v],
181        )
182    }
183
184    /// Where the surface stops being one polynomial piece along `u`: the
185    /// ends of its domain and every distinct knot between them.
186    pub fn breakpoints_u(&self) -> Vec<Scalar> {
187        crate::spline::breakpoints(&self.knot_vector_u, self.degree_u, self.num_u)
188    }
189
190    /// Like [`Self::breakpoints_u`], along `v`.
191    pub fn breakpoints_v(&self) -> Vec<Scalar> {
192        crate::spline::breakpoints(&self.knot_vector_v, self.degree_v, self.num_v)
193    }
194
195    /// Whether `other` could be the very same patch: equal degrees, and
196    /// knots and homogeneous control points that `could_be_equal` pairwise.
197    ///
198    /// A test of the representation, so it is exact where it answers `true`
199    /// — two patches built the same way, e.g. two copies of one solid — and
200    /// says nothing about the same surface parametrized differently.
201    pub fn could_be_equal(&self, other: &Self) -> bool {
202        self.degree_u == other.degree_u
203            && self.degree_v == other.degree_v
204            && self.num_u == other.num_u
205            && self.num_v == other.num_v
206            && self
207                .knot_vector_u
208                .iter()
209                .zip(&other.knot_vector_u)
210                .all(|(a, b)| a.could_be_equal(*b))
211            && self
212                .knot_vector_v
213                .iter()
214                .zip(&other.knot_vector_v)
215                .all(|(a, b)| a.could_be_equal(*b))
216            && self
217                .control_points
218                .iter()
219                .zip(&other.control_points)
220                .all(|(a, b)| a.could_be_equal(b))
221    }
222
223    /// The control row at index `k` along `u` (`u_fixed`: the points
224    /// `P[k][·]`, running along `v`), or else the column `P[·][k]` along
225    /// `u`.
226    pub(crate) fn control_row(&self, u_fixed: bool, k: usize) -> Vec<Vector<D>> {
227        let nv = self.num_v;
228        if u_fixed {
229            self.control_points[k * nv..(k + 1) * nv].to_vec()
230        } else {
231            (0..self.num_u)
232                .map(|i| self.control_points[i * nv + k])
233                .collect()
234        }
235    }
236
237    /// Whether this is the [`NurbSurface::everything`] placeholder — the
238    /// unsharp stand-in a face carries before it is given real geometry.
239    ///
240    /// A finished solid must have none: a placeholder face has no position,
241    /// so nothing can be classified against it, and it will silently swallow
242    /// any containment or intersection query it is handed (every comparison
243    /// against `ENTIRE` succeeds). Construction code that splits faces off a
244    /// starting placeholder has to consume the last one rather than leave it
245    /// behind, and this is how a test says so.
246    pub fn is_everything(&self) -> bool {
247        self.num_u == 1
248            && self.num_v == 1
249            && !self.control_points[0][0].is_sharp()
250            && self.knot_vector_u.iter().all(|k| !k.is_sharp())
251    }
252
253    /// A degenerate, maximally-unsharp surface: a single 1×1 control point
254    /// whose every coordinate is [`Scalar::ENTIRE`], over a domain that
255    /// accepts any `(u, v)`. `evaluate()` anywhere returns `ENTIRE` in every
256    /// coordinate, so it `could_be_equal`s any point — a placeholder for
257    /// geometry that is not yet known.
258    pub fn everything() -> Self {
259        let mut cp = Vector::<D>::everything();
260        cp[D - 1] = Scalar::ONE;
261        NurbSurface {
262            degree_u: 0,
263            degree_v: 0,
264            num_u: 1,
265            num_v: 1,
266            control_points: vec![cp],
267            knot_vector_u: vec![Scalar::ENTIRE, Scalar::ENTIRE],
268            knot_vector_v: vec![Scalar::ENTIRE, Scalar::ENTIRE],
269            // Matches-anything, same as every other coordinate of this
270            // placeholder — no bounding box has been established yet.
271            aabb: [Scalar::ENTIRE; 3],
272        }
273    }
274}
275
276impl Display for NurbSurface3D {
277    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
278        let p00 = self
279            .evaluate(self.knot_vector_u[0], self.knot_vector_v[0])
280            .map(|p| p.to_string())
281            .unwrap_or_else(|_| "N/A".to_string());
282        let p01 = self
283            .evaluate(self.knot_vector_u[0], self.knot_vector_v[self.num_v])
284            .map(|p| p.to_string())
285            .unwrap_or_else(|_| "N/A".to_string());
286        let p10 = self
287            .evaluate(self.knot_vector_u[self.num_u], self.knot_vector_v[0])
288            .map(|p| p.to_string())
289            .unwrap_or_else(|_| "N/A".to_string());
290        let p11 = self
291            .evaluate(
292                self.knot_vector_u[self.num_u],
293                self.knot_vector_v[self.num_v],
294            )
295            .map(|p| p.to_string())
296            .unwrap_or_else(|_| "N/A".to_string());
297
298        write!(
299            f,
300            "NurbSurface((0, 0) -> {}, (1, 0) -> {}, (1, 1) -> {}, (0, 1) -> {})",
301            p00, p10, p11, p01
302        )
303    }
304}
305
306#[cfg(test)]
307mod tests {
308    use super::NurbSurface;
309    use geop_core_math::{scalars::Scalar, vector::Vector4};
310
311    /// Weights must be definitely positive: zero and negative are rejected.
312    #[test]
313    fn non_positive_weight_is_rejected() {
314        let f = Scalar::from_f64;
315        let pt = |x, y, w| Vector4::from_array([f(x), f(y), f(0.), f(w)]);
316        for w in [0., -1.] {
317            let result = NurbSurface::<4>::try_new(
318                1,
319                1,
320                vec![
321                    pt(0., 0., w),
322                    pt(0., 1., 1.),
323                    pt(1., 0., 1.),
324                    pt(1., 1., 1.),
325                ],
326                vec![f(0.), f(0.), f(1.), f(1.)],
327                vec![f(0.), f(0.), f(1.), f(1.)],
328            );
329            assert!(result.is_err(), "weight {w} accepted");
330        }
331    }
332}
geop_core_geometry/nurb_surface/mod.rs

geop_core_geometry/nurb_surface/
mod.rs
