geop_core_geometry/nurb_curve/mod.rs
1mod area;
2pub use area::loop_signed_area;
3mod closest;
4mod compatible;
5mod construct;
6pub use construct::bezier_knots;
7mod evaluate;
8mod extend;
9mod extrema;
10mod fit;
11pub use fit::WidestBetween;
12mod helix;
13pub use helix::{Handedness, Helix, HelixRow, arc_span, cos_sin, helix_rows, sqrt2_over_2};
14mod interpolate;
15mod join;
16mod length;
17pub use length::Length;
18mod offset;
19mod refine;
20pub(crate) use interpolate::basis_funs;
21pub(crate) use interpolate::solve_interpolation_system;
22pub use interpolate::true_point_fractions;
23pub use refine::ParameterRefinable;
24mod reverse;
25mod split;
26mod sweep;
27mod tangent;
28mod translate;
29mod within;
30pub(crate) use translate::transform_control_point;
31pub use translate::{embed_point, flatten_point};
32pub use within::Within;
33
34use std::fmt::Display;
35
36use crate::{aabb::compute_aabb, knot_insertion::pinned_clamped_end};
37use geop_core_math::{
38 geop_error::{GeopError, GeopResult},
39 scalars::Scalar,
40 vector::Vector,
41};
42
43/// Dehomogenize `control_points` (`D`-dimensional homogeneous, last component
44/// the weight) into `C`-dimensional Cartesian points (`C` = `D - 1`, passed
45/// explicitly since Rust's stable const generics can't express `D - 1` in a
46/// single generic parameter's bound).
47///
48/// Infallible: every weight of a `NurbCurve`/`NurbSurface` is definitely
49/// positive by construction (see `NurbCurve::try_new`), so the division
50/// cannot fail. `pub(crate)`: surfaces and the intersection searches need
51/// the exact same points.
52pub(crate) fn dehomogenize<const D: usize, const C: usize>(
53 control_points: &[Vector<D>],
54) -> Vec<Vector<C>> {
55 control_points
56 .iter()
57 .map(|p| {
58 let inv_w = Scalar::ONE
59 .div(p[D - 1])
60 .expect("weights are definitely positive by construction");
61 let mut pt = Vector::<C>::zero();
62 for c in 0..C {
63 pt[c] = p[c].mul(inv_w);
64 }
65 pt
66 })
67 .collect()
68}
69
70/// A NURBS curve whose control points live in `D`-dimensional homogeneous space.
71///
72/// - `D = 4`: 3-D curve — control points are `(wx, wy, wz, w)`.
73/// - `D = 3`: 2-D curve (pcurve) — control points are `(wu, wv, w)`.
74#[derive(Clone, Debug, PartialEq)]
75pub struct NurbCurve<const D: usize> {
76 pub degree: usize,
77 pub control_points: Vec<Vector<D>>,
78 pub knot_vector: Vec<Scalar>,
79 /// Cached axis-aligned bounding box of the (dehomogenized) control
80 /// points — see [`crate::aabb::compute_aabb`]. Every constructor here
81 /// fills this in once, so the intersection search's per-node prefilter
82 /// (`aabb_could_overlap`, in `intersection::curve_curve`/
83 /// `curve_surface`) never has to recompute it.
84 pub(crate) aabb: [Scalar; 3],
85}
86
87/// 3-D NURBS curve (homogeneous control points in ℝ⁴).
88pub type NurbCurve3D = NurbCurve<4>;
89
90/// 2-D parameter-space NURBS curve / pcurve (homogeneous control points in ℝ³).
91pub type NurbCurve2D = NurbCurve<3>;
92
93impl<const D: usize> NurbCurve<D> {
94 pub fn try_new(
95 degree: usize,
96 control_points: Vec<Vector<D>>,
97 knot_vector: Vec<Scalar>,
98 ) -> GeopResult<Self> {
99 let n = control_points.len();
100 if knot_vector.len() != n + degree + 1 {
101 return Err(GeopError::new(format!(
102 "Invalid knot vector length: expected {}, got {}",
103 n + degree + 1,
104 knot_vector.len()
105 )));
106 }
107 // Every weight must be definitely positive. That is what makes
108 // `W(t) = Σ w_i N_i(t) > 0` on the whole domain, which the convex hull
109 // property and `contains::curve`'s division-free per-axis functions
110 // `X_k(t) - p_k W(t)` rely on. Zero or negative weights have no use
111 // case here, so they are rejected at construction rather than
112 // surfacing as a division by zero far downstream. `split` forms
113 // convex combinations of existing control points and so preserves
114 // this.
115 for p in &control_points {
116 if !p[D - 1].definitely_greater(Scalar::ZERO) {
117 return Err(GeopError::new(format!(
118 "NurbCurve::try_new: control point {p:?} has a weight that is not definitely positive"
119 )));
120 }
121 }
122 let aabb = compute_aabb(&control_points);
123 Ok(Self {
124 degree,
125 control_points,
126 knot_vector,
127 aabb,
128 })
129 }
130
131 /// Refresh the cached [`Self::aabb`] from the current `control_points`.
132 ///
133 /// Every constructor in this module keeps `aabb` in sync automatically,
134 /// but `control_points` is a `pub` field and at least one caller outside
135 /// this crate (`Model::reverse_face`, mirroring a pcurve's control
136 /// points in place to flip a face) legitimately mutates it directly
137 /// rather than building a new curve — call this afterwards or the
138 /// cached box silently goes stale and the intersection search's
139 /// `aabb_could_overlap` prefilter starts pruning real overlaps.
140 pub fn recompute_aabb(&mut self) {
141 self.aabb = compute_aabb(&self.control_points);
142 }
143
144 /// The box around the control points, and so around the curve (the
145 /// convex hull property): one enclosure per axis, holding every control
146 /// point's coordinate there. Only the first `D - 1` axes are the
147 /// curve's; a 2-D curve's third is [`Scalar::ENTIRE`].
148 pub fn aabb(&self) -> &[Scalar; 3] {
149 &self.aabb
150 }
151
152 /// Whether the boxes around two curves' control points could overlap:
153 /// `false` proves the curves apart.
154 pub fn could_overlap(&self, other: &Self) -> bool {
155 crate::aabb::aabb_could_overlap(&self.aabb, &other.aabb, D - 1)
156 }
157
158 /// Valid parameter range `(start_t, end_t)` of this curve.
159 pub fn domain(&self) -> (Scalar, Scalar) {
160 let p = self.degree;
161 let n = self.control_points.len() - 1;
162 (self.knot_vector[p], self.knot_vector[n + 1])
163 }
164
165 /// The diagonal of the box around the control points: by the convex
166 /// hull property, a length the whole curve fits within (see
167 /// `NurbSurface::size`).
168 pub fn size(&self) -> GeopResult<Scalar> {
169 self.aabb
170 .iter()
171 .take(D - 1)
172 .fold(Scalar::ZERO, |sum, axis| {
173 sum.add(axis.width().mul(axis.width()))
174 })
175 .sqrt()
176 }
177
178 /// The middle of [`Self::domain`], sharp: a parameter inside the curve
179 /// for a sample, a seed or a split point. Any value in the enclosure of
180 /// the middle serves as well, so it is a free choice, and sharpening it
181 /// loses nothing (see AGENTS.md).
182 pub fn mid_parameter(&self) -> GeopResult<Scalar> {
183 let (t0, t1) = self.domain();
184 Ok(t0.add(t1).div(Scalar::TWO)?.sharpen())
185 }
186
187 /// Where the curve stops being one polynomial piece: the ends of its
188 /// domain and every distinct knot between them (see
189 /// [`crate::spline::breakpoints`]).
190 pub fn breakpoints(&self) -> Vec<Scalar> {
191 crate::spline::breakpoints(&self.knot_vector, self.degree, self.control_points.len())
192 }
193
194 /// If `hat` pins the parameter sharply onto an end of the domain that
195 /// is clamped (see [`crate::knot_insertion::pinned_clamped_end`]), the
196 /// curve's point there: its end control point, in `C = D - 1`
197 /// Cartesian coordinates.
198 pub(crate) fn pinned_end<const C: usize>(&self, hat: Scalar) -> Option<Vector<C>> {
199 let n = self.control_points.len();
200 let first = pinned_clamped_end(hat, &self.knot_vector, n, self.degree)?;
201 let end = self.control_points[if first { 0 } else { n - 1 }];
202 Some(dehomogenize::<D, C>(&[end])[0])
203 }
204
205 pub fn domain_as_scalar(&self) -> Scalar {
206 let (s, e) = self.domain();
207 s.union(e)
208 }
209
210 /// A degenerate, maximally-unsharp curve: a single control point whose
211 /// every coordinate is [`Scalar::ENTIRE`], over a domain that accepts
212 /// any parameter. `evaluate()` at any `t` returns `ENTIRE` in every
213 /// coordinate, so it `could_be_equal`s any point — a placeholder for
214 /// geometry that is not yet known.
215 pub fn everything() -> Self {
216 let mut cp = Vector::<D>::everything();
217 cp[D - 1] = Scalar::ONE;
218 NurbCurve {
219 degree: 0,
220 control_points: vec![cp],
221 knot_vector: vec![Scalar::ENTIRE, Scalar::ENTIRE],
222 // Matches-anything, same as every other coordinate of this
223 // placeholder — no bounding box has been established yet.
224 aabb: [Scalar::ENTIRE; 3],
225 }
226 }
227}
228
229impl<const D: usize> Display for NurbCurve<D> {
230 fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
231 let at = |t: Scalar| {
232 let mut p = [Scalar::ZERO; D];
233 match self.evaluate_into(t, &mut p[..D - 1]) {
234 Ok(()) => {
235 let coordinates: Vec<String> =
236 p[..D - 1].iter().map(Scalar::to_string).collect();
237 format!("Vector{}({})", D - 1, coordinates.join(", "))
238 }
239 Err(_) => "N/A".to_string(),
240 }
241 };
242 let (start, end) = self.domain();
243 write!(f, "NurbCurve({} -> {})", at(start), at(end))
244 }
245}
246
247#[cfg(test)]
248mod tests {
249 use super::NurbCurve;
250 use geop_core_math::{scalars::Scalar, vector::Vector4};
251
252 /// Weights must be definitely positive: zero and negative are rejected.
253 #[test]
254 fn non_positive_weight_is_rejected() {
255 let f = Scalar::from_f64;
256 for w in [0., -1.] {
257 let result = NurbCurve::<4>::try_new(
258 1,
259 vec![
260 Vector4::from_array([f(0.), f(0.), f(0.), f(w)]),
261 Vector4::from_array([f(1.), f(0.), f(0.), f(1.)]),
262 ],
263 vec![f(0.), f(0.), f(1.), f(1.)],
264 );
265 assert!(result.is_err(), "weight {w} accepted");
266 }
267 }
268}