geop_core_geometry/nurb_curve/refine.rs
1use geop_core_math::{geop_error::GeopResult, scalars::Scalar, vector::Vector};
2
3use super::NurbCurve;
4
5/// Gauss-Newton foot-point iterations. Quadratic convergence means a handful
6/// is plenty; unlike `max_nodes` this cannot change whether a correct answer
7/// is found, only how tightly an already-isolated one is pinned down.
8const ITERATIONS: usize = 12;
9
10impl<S: Scalar, const D: usize> NurbCurve<S, 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`.
14 ///
15 /// # Why this exists
16 ///
17 /// Subdivision is a *global* method: it reliably finds and separates
18 /// every solution, and (via the leaf-count signal the intersection
19 /// searches rely on) recognizes coincidence even for a partial overlap.
20 /// What it is bad at is the last few digits — it converges one bit per
21 /// split, so squeezing a parameter down to machine accuracy would take
22 /// ~50 levels, which is exponentially more work than the ~7 needed to
23 /// isolate the solution in the first place.
24 ///
25 /// Newton is the opposite: useless for finding solutions, unbeatable for
26 /// polishing one that is already isolated, converging quadratically. So
27 /// the two compose — subdivide to isolate, then refine here.
28 ///
29 /// This is what lets [`NurbCurve::split`] be called without sharpening.
30 /// A `min_subdivision_size`-wide `t` cannot be fed to Boehm insertion:
31 /// its width flows into `alpha = (t - e) / (s - e)`, whose denominator
32 /// shrinks with every successive split while the width does not, so the
33 /// sub-curves' control points widen without bound. The old answer was to
34 /// `sharpen` the parameter at each call site, which moved the split to
35 /// the interval's midpoint rather than the point actually located — a
36 /// silent geometric error of `|t_mid - t*| x |C'(t)|`, and the reason an
37 /// edge endpoint could land ~1e-8 from the vertex it is anchored to. A
38 /// refined parameter is narrow *and* still an honest enclosure, so it
39 /// needs no sharpening and introduces no such error.
40 ///
41 /// # Honesty of the result
42 ///
43 /// Every iterate except the last is sharpened, which is legitimate: it is
44 /// only a seed for the next step, and any value inside it is an equally
45 /// good one. The final step is left unsharpened, so the returned width is
46 /// the honest statement of how precisely `point` pins down a parameter
47 /// (see "Sharpen only where the value is a free choice" in `AGENTS.md` —
48 /// this is exactly the rule `NurbSurface::project` follows).
49 ///
50 /// The result is then intersected with the incoming `t`: both are valid
51 /// enclosures of the same parameter, so their intersection is too, and is
52 /// tighter than either. If they turn out to be disjoint, Newton has
53 /// wandered out of the box the search proved the solution lies in — the
54 /// incoming enclosure is returned unchanged rather than trusting the
55 /// refinement. The same fallback covers a vanishing tangent (the
56 /// Gauss-Newton denominator could be zero), so this never turns a usable
57 /// answer into a failure.
58 pub fn refine_parameter_at_point<const C: usize>(
59 &self,
60 t: S,
61 point: &Vector<S, C>,
62 ) -> GeopResult<S>
63 where
64 NurbCurve<S, D>: ParameterRefinable<S, C>,
65 {
66 let (lo, hi) = self.domain();
67
68 let mut current = t.sharpen();
69 for iteration in 0..ITERATIONS {
70 let position = self.evaluate_cartesian(current)?;
71 let tangent = self.tangent_cartesian(current)?;
72
73 // Gauss-Newton on |C(t) - P|^2: the step that zeroes the
74 // directional residual along the tangent.
75 let residual = position.sub(point);
76 let numerator = residual.prod_dot(&tangent);
77 let denominator = tangent.prod_dot(&tangent);
78 let Ok(step) = numerator.div(denominator) else {
79 return Ok(t);
80 };
81
82 let next = current.sub(step);
83 let next = if iteration + 1 == ITERATIONS {
84 next
85 } else {
86 next.sharpen()
87 };
88 // Clamp rather than bail. A foot point can legitimately step
89 // outside the domain on the way in — bailing there silently
90 // returns the wide search parameter, which is indistinguishable
91 // from a successful refinement to every caller and reintroduces
92 // exactly the fat sub-curves this exists to prevent.
93 current = if next.definitely_less(lo) {
94 lo
95 } else if next.definitely_greater(hi) {
96 hi
97 } else {
98 next
99 };
100 }
101
102 if !current.could_be_equal(t) {
103 return Ok(t);
104 }
105 Ok(t.intersect(current))
106 }
107}
108
109/// Bridges a `NurbCurve<S, D>`'s Cartesian evaluation (`D = 4` -> 3-D points,
110/// `D = 3` -> 2-D pcurve points) so [`NurbCurve::refine_parameter_at_point`]
111/// can be written once for both — stable Rust's const generics cannot express
112/// `C = D - 1` directly.
113pub trait ParameterRefinable<S: Scalar, const C: usize> {
114 fn evaluate_cartesian(&self, t: S) -> GeopResult<Vector<S, C>>;
115 /// The Cartesian tangent, via each dimension's own `tangent`. Not the
116 /// `derivative()` curve: that is the *homogeneous* derivative, whose
117 /// weight component is zero for a non-rational curve, so evaluating it
118 /// as a rational curve fails outright.
119 fn tangent_cartesian(&self, t: S) -> GeopResult<Vector<S, C>>;
120}
121
122impl<S: Scalar> ParameterRefinable<S, 3> for NurbCurve<S, 4> {
123 fn evaluate_cartesian(&self, t: S) -> GeopResult<Vector<S, 3>> {
124 self.evaluate(t)
125 }
126 fn tangent_cartesian(&self, t: S) -> GeopResult<Vector<S, 3>> {
127 self.tangent(t)
128 }
129}
130
131impl<S: Scalar> ParameterRefinable<S, 2> for NurbCurve<S, 3> {
132 fn evaluate_cartesian(&self, t: S) -> GeopResult<Vector<S, 2>> {
133 self.evaluate(t)
134 }
135 fn tangent_cartesian(&self, t: S) -> GeopResult<Vector<S, 2>> {
136 self.tangent(t)
137 }
138}
139
140#[cfg(test)]
141mod tests {
142 use crate::nurb_curve::NurbCurve;
143 use geop_core_math::for_all_scalars;
144 use geop_core_math::{scalars::Scalar, vector::Vector4};
145
146 fn pt<S: Scalar>(x: f64, y: f64, z: f64) -> Vector4<S> {
147 Vector4::from_array([S::from_f64(x), S::from_f64(y), S::from_f64(z), S::ONE])
148 }
149
150 /// Degree-2 arc, so the parametrization is genuinely nonlinear.
151 fn arc<S: Scalar>() -> NurbCurve<S, 4> {
152 let f = S::from_f64;
153 NurbCurve::try_new(
154 2,
155 vec![pt(0.0, 0.0, 0.0), pt(1.0, 2.0, 0.0), pt(2.0, 0.0, 0.0)],
156 vec![f(0.), f(0.), f(0.), f(1.), f(1.), f(1.)],
157 )
158 .unwrap()
159 }
160
161 /// A search-width parameter refines to a far tighter one that still
162 /// encloses the true parameter, and still lands on the same point.
163 fn check_refines_a_wide_parameter<S: Scalar>() {
164 let curve = arc::<S>();
165 let exact = S::from_f64(0.375);
166 let target = curve.evaluate(exact).unwrap();
167
168 // What a 1e-4 subdivision search would hand back.
169 let wide = S::from_f64(0.3745).union(S::from_f64(0.3755));
170 let refined = curve.refine_parameter_at_point(wide, &target).unwrap();
171
172 assert!(
173 refined.could_be_equal(exact),
174 "refined parameter must still enclose the true one: {refined:?}"
175 );
176 let landed = curve.evaluate(refined).unwrap();
177 for c in 0..3 {
178 assert!(landed[c].could_be_equal(target[c]), "coord {c} moved");
179 }
180 }
181 #[test]
182 fn refines_a_wide_parameter() {
183 for_all_scalars!(check_refines_a_wide_parameter);
184 }
185
186 /// A point nowhere near the curve must not drag the parameter somewhere
187 /// arbitrary — the incoming enclosure comes back untouched.
188 fn check_off_curve_point_falls_back<S: Scalar>() {
189 let curve = arc::<S>();
190 let target = geop_core_math::vector::Vector3::from_array([
191 S::from_f64(50.0),
192 S::from_f64(50.0),
193 S::from_f64(50.0),
194 ]);
195 let wide = S::from_f64(0.3745).union(S::from_f64(0.3755));
196 let refined = curve.refine_parameter_at_point(wide, &target).unwrap();
197 assert!(refined.could_be_equal(wide));
198 }
199 #[test]
200 fn off_curve_point_falls_back() {
201 for_all_scalars!(check_off_curve_point_falls_back);
202 }
203}