geop_core_geometry/nurb_curve/refine.rs
1use geop_core_math::{
2 geop_error::{GeopError, GeopResult},
3 scalars::Scalar,
4 vector::Vector,
5};
6
7use super::NurbCurve;
8use geop_core_math::solvers::newton;
9
10impl<const D: usize> NurbCurve<D> {
11 /// Refine `t` — a parameter enclosure produced by a subdivision search —
12 /// into the tightest enclosure of the parameter at which this curve
13 /// passes through `point`: Gauss-Newton on `|C(t) - P|²`, by
14 /// [`newton::refine`], which says why and how honestly.
15 ///
16 /// A search's parameter is as wide as the search's tolerance; the
17 /// refined one is as narrow as the curve and the point pin it down, and
18 /// still an honest enclosure. [`Self::split_parameter_at_point`] cuts
19 /// within it.
20 ///
21 /// A point nowhere near the curve, or a vanishing tangent, leaves `t`
22 /// as it is.
23 pub fn refine_parameter_at_point<const C: usize>(&self, t: Scalar, point: &Vector<C>) -> Scalar
24 where
25 NurbCurve<D>: ParameterRefinable<C>,
26 {
27 let [t] = newton::refine([t], [self.domain()], |[current]| {
28 let position = self.evaluate_cartesian(current).ok()?;
29 let tangent = self.tangent_cartesian(current).ok()?;
30 // The step that zeroes the directional residual along the
31 // tangent.
32 let residual = position.sub(point);
33 let numerator = residual.prod_dot(&tangent);
34 let denominator = tangent.prod_dot(&tangent);
35 let step = numerator.div(denominator).ok()?;
36 Some([current.sub(step)])
37 });
38 t
39 }
40
41 /// The parameter to [`split`](Self::split) this curve at so a piece
42 /// ends at `point`, from a parameter `t` that isolates it: `t` refined
43 /// ([`Self::refine_parameter_at_point`]), then the sharp middle of it.
44 ///
45 /// Where in the refined enclosure the curve is cut is a free choice.
46 /// `point` is an enclosure, and every parameter in the refined one
47 /// gives a piece that lies on the curve and ends within `point`, so
48 /// any of them serves; the check below asks it of the one chosen.
49 /// Cutting at the whole enclosure instead keeps *every* such piece in
50 /// one: the junction's uncertainty along the curve becomes width in the
51 /// pieces' control points (through `α = (t - a) / (b - a)` in Boehm
52 /// insertion) and in their end knots. The next cut's parameter is then
53 /// found on the wider piece, comes out wider, and widens the next
54 /// pieces more: a quarter arc of radius 4 cut forty times went from
55 /// `3e-15` wide to `0.08` (three times wider per cut), and a disc's rim
56 /// cut by thirty-two slots ended `4e-10` off its cap. A sharp cut adds
57 /// the rounding of one subdivision per cut, nothing more.
58 ///
59 /// A knot in the refined enclosure is the good choice there: the curve
60 /// is cut at it as at the knot it already has, where cutting beside it
61 /// would leave a piece a sliver of a span, a hair long, its control
62 /// points one point and its end without a tangent (a chamfer's contact
63 /// edge cut where its run-out span begins, as a fillet then found).
64 ///
65 /// Refused if `t` could be an end of the domain — then `point` is an
66 /// end of the curve, nothing to cut — or if the curve at the chosen
67 /// parameter does not reach `within`, the enclosure `point` was chosen
68 /// in: the vertex. Not `point` itself: a pcurve running from `u = 1`
69 /// has parameters a hair `1e-16` apart near its end, and a `point` of
70 /// `u = 1e-4`, sharper than that, no parameter reaches.
71 pub fn split_parameter_at_point<const C: usize>(
72 &self,
73 t: Scalar,
74 point: &Vector<C>,
75 within: &Vector<C>,
76 ) -> GeopResult<Scalar>
77 where
78 NurbCurve<D>: ParameterRefinable<C>,
79 {
80 let (lo, hi) = self.domain();
81 let refined = self.refine_parameter_at_point(t, point);
82 if !refined.definitely_greater(lo) || !refined.definitely_less(hi) {
83 return Err(GeopError::new(format!(
84 "the parameter {refined:?} of {point:?} could be an end of the domain ({lo:?}, {hi:?}): the point is an end of the curve, not inside it"
85 )));
86 }
87 let cut = self
88 .knot_vector
89 .iter()
90 .find(|k| {
91 k.could_be_equal(refined) && k.definitely_greater(lo) && k.definitely_less(hi)
92 })
93 .map_or_else(|| refined.sharpen(), |k| k.sharpen());
94 let reached = self.evaluate_cartesian(cut)?;
95 if !reached.could_be_equal(within) {
96 return Err(GeopError::new(format!(
97 "the curve at {cut:?}, the middle of the refined parameter {refined:?}, is at {reached:?}, not within {within:?}"
98 )));
99 }
100 Ok(cut)
101 }
102}
103
104/// Bridges a `NurbCurve<D>`'s Cartesian evaluation (`D = 4` -> 3-D points,
105/// `D = 3` -> 2-D pcurve points) so [`NurbCurve::refine_parameter_at_point`]
106/// can be written once for both — stable Rust's const generics cannot express
107/// `C = D - 1` directly.
108pub trait ParameterRefinable<const C: usize> {
109 fn evaluate_cartesian(&self, t: Scalar) -> GeopResult<Vector<C>>;
110 /// The Cartesian tangent, via each dimension's own `tangent` (the
111 /// quotient rule over the homogeneous `derivative()`).
112 fn tangent_cartesian(&self, t: Scalar) -> GeopResult<Vector<C>>;
113}
114
115impl ParameterRefinable<3> for NurbCurve<4> {
116 fn evaluate_cartesian(&self, t: Scalar) -> GeopResult<Vector<3>> {
117 self.evaluate(t)
118 }
119 fn tangent_cartesian(&self, t: Scalar) -> GeopResult<Vector<3>> {
120 self.tangent(t)
121 }
122}
123
124impl ParameterRefinable<2> for NurbCurve<3> {
125 fn evaluate_cartesian(&self, t: Scalar) -> GeopResult<Vector<2>> {
126 self.evaluate(t)
127 }
128 fn tangent_cartesian(&self, t: Scalar) -> GeopResult<Vector<2>> {
129 self.tangent(t)
130 }
131}
132
133#[cfg(test)]
134mod tests {
135 use crate::nurb_curve::NurbCurve;
136 use crate::test_fixtures::hom;
137 use geop_core_math::scalars::Scalar;
138
139 /// Degree-2 arc, so the parametrization is genuinely nonlinear.
140 fn arc() -> NurbCurve<4> {
141 let f = Scalar::from_f64;
142 NurbCurve::try_new(
143 2,
144 vec![
145 hom(0.0, 0.0, 0.0, 1.),
146 hom(1.0, 2.0, 0.0, 1.),
147 hom(2.0, 0.0, 0.0, 1.),
148 ],
149 vec![f(0.), f(0.), f(0.), f(1.), f(1.), f(1.)],
150 )
151 .unwrap()
152 }
153
154 /// A search-width parameter refines to a far tighter one that still
155 /// encloses the true parameter, and still lands on the same point.
156 #[test]
157 fn refines_a_wide_parameter() {
158 let curve = arc();
159 let exact = Scalar::from_f64(0.375);
160 let target = curve.evaluate(exact).unwrap();
161
162 // What a 1e-4 subdivision search would hand back.
163 let wide = Scalar::from_f64(0.3745).union(Scalar::from_f64(0.3755));
164 let refined = curve.refine_parameter_at_point(wide, &target);
165
166 assert!(
167 refined.could_be_equal(exact),
168 "refined parameter must still enclose the true one: {refined:?}"
169 );
170 let landed = curve.evaluate(refined).unwrap();
171 for c in 0..3 {
172 assert!(landed[c].could_be_equal(target[c]), "coord {c} moved");
173 }
174 }
175
176 /// A point nowhere near the curve must not drag the parameter somewhere
177 /// arbitrary — the incoming enclosure comes back untouched.
178 #[test]
179 fn off_curve_point_falls_back() {
180 let curve = arc();
181 let target = geop_core_math::vector::Vector3::from_array([
182 Scalar::from_f64(50.0),
183 Scalar::from_f64(50.0),
184 Scalar::from_f64(50.0),
185 ]);
186 let wide = Scalar::from_f64(0.3745).union(Scalar::from_f64(0.3755));
187 let refined = curve.refine_parameter_at_point(wide, &target);
188 assert!(refined.could_be_equal(wide));
189 }
190}