geop_core_geometry/nurb_curve/fit.rs
1//! [`NurbCurve3D::fit_enclosing`]: a cubic through stations of a curve
2//! known only point by point, widened to enclose the curve between them.
3
4use geop_core_math::{
5 geop_error::{GeopError, GeopResult, WithContext},
6 scalars::Scalar,
7 solvers::refine::{Refine, Refined, refine},
8 vector::Vector3,
9};
10
11use super::{NurbCurve, NurbCurve3D, true_point_fractions};
12
13/// The fewest legs a fitted curve runs through (see
14/// [`NurbCurve3D::fit_enclosing`]): a cubic interpolant needs points enough
15/// to follow the curve, however few stations it was given. Like
16/// [`FIT_DRIFT`] this bounds effort — every extra leg is one more point
17/// located — and decides how *wide* the fit's enclosure is, not whether it
18/// holds.
19const MIN_FITTED_LEGS: usize = 8;
20
21/// How wide a fit's enclosure may come out before its legs are split finer
22/// and it is fitted again: ten times what the usual traced branch leaves.
23/// Like [`MIN_FITTED_LEGS`] it bounds effort and decides how wide the
24/// enclosure is, not whether it holds.
25const FIT_DRIFT: f64 = 1e-5;
26
27/// The most legs a fit is refined to: the budget [`FIT_DRIFT`] is pursued
28/// with, each leg a point located and a row of the interpolation's system.
29const MAX_FITTED_LEGS: usize = 1024;
30
31/// Where a fit's widest point between stations was, of its last attempt:
32/// leg, fraction, point, width.
33pub type WidestBetween = Option<(usize, usize, Vector3, f64)>;
34
35impl NurbCurve3D {
36 /// A cubic through `stations` of a curve — each a point of it and
37 /// whatever `locate` needs to find more — that encloses the curve
38 /// between them too. `locate(from, to, frac)` is the point of the curve
39 /// a fraction `frac` of the way along the leg from station `from` to
40 /// `to`, with what it needs to go on from there.
41 ///
42 /// Every leg is first split into pieces: a curve given in only a few
43 /// stations would leave too few points for a cubic to follow it, so
44 /// every fit has at least [`MIN_FITTED_LEGS`] legs. The cubic drifts
45 /// from the curve between its points, which is enclosed as width,
46 /// honestly (see [`NurbCurve3D::interpolate_enclosing`]); while that
47 /// width is more than [`FIT_DRIFT`] and refining the legs still halves
48 /// it, the legs that drift are halved: those a halving of all, which
49 /// shrinks a cubic's drift sixteen times, would leave the widest. Where
50 /// the curve bends evenly that is every leg. Where it turns sharply in
51 /// one place — a blend's arc, tangent to a plane it is cut by at three
52 /// degrees, meets it in an ellipse nineteen times as long as it is wide,
53 /// its tip of radius `5e-3` — it is those at the tip, which legs of an
54 /// even length reach only when there are hundreds of them. Once
55 /// refining no longer halves the width, the width is the points' own,
56 /// and the coarser fit is kept.
57 ///
58 /// `widest_between` is told the widest point of the curve located
59 /// between stations of the last attempt: the curve is widened to hold
60 /// it, so one that is wide makes the whole curve wide, and this says
61 /// which.
62 pub fn fit_enclosing<T: Clone>(
63 stations: &[(Vector3, T)],
64 locate: impl Fn(&(Vector3, T), &(Vector3, T), Scalar) -> GeopResult<(Vector3, T)>,
65 widest_between: &mut WidestBetween,
66 ctx: &dyn Fn(GeopError) -> GeopError,
67 ) -> GeopResult<NurbCurve<4>> {
68 let strides = stations.len() - 1;
69 let pieces = MIN_FITTED_LEGS.div_ceil(strides).max(2);
70 let mut dense = vec![stations[0].clone()];
71 for i in 0..strides {
72 for k in 1..pieces {
73 let frac = Scalar::from_ratio(k as i64, pieces as i64).with_context(&ctx)?;
74 dense.push(locate(&stations[i], &stations[i + 1], frac).with_context(&ctx)?);
75 }
76 dense.push(stations[i + 1].clone());
77 }
78 // With each step the fit before, to tell whether refining the legs
79 // still halved its width: each step halves it or is the last, so
80 // the steps' cap is never met before [`FIT_DRIFT`] is.
81 let first = (dense, None::<(NurbCurve<4>, Scalar)>);
82 let refined = refine(
83 first,
84 usize::BITS as usize,
85 |(dense, best)| -> GeopResult<_> {
86 // The curve inside each of the legs, which
87 // `interpolate_enclosing` widens the cubic to hold.
88 let legs = dense.len() - 1;
89 let mut between = Vec::with_capacity(legs);
90 for i in 0..legs {
91 let fractions = true_point_fractions(i, legs);
92 let mut inside = Vec::with_capacity(fractions.len());
93 for (k, &(a, b)) in fractions.iter().enumerate() {
94 let frac = Scalar::from_ratio(a, b).with_context(&ctx)?;
95 let (p, _) = locate(&dense[i], &dense[i + 1], frac).with_context(&ctx)?;
96 let width = (0..3).map(|c| p[c].width().to_f64()).fold(0.0, f64::max);
97 if widest_between.as_ref().is_none_or(|w| width > w.3) {
98 *widest_between = Some((i, k, p, width));
99 }
100 inside.push(p);
101 }
102 between.push(inside);
103 }
104 let points: Vec<Vector3> = dense.iter().map(|(p, _)| *p).collect();
105 let curve =
106 NurbCurve3D::interpolate_enclosing(&points, &between, 3).with_context(&ctx)?;
107 let width = curve
108 .control_points
109 .iter()
110 .flat_map(|p| (0..3).map(move |k| p[k].width()))
111 .fold(
112 Scalar::ZERO,
113 |a, w| if w.definitely_greater(a) { w } else { a },
114 );
115 let halved = best
116 .as_ref()
117 .is_none_or(|(_, before)| before.definitely_greater(width.mul(Scalar::TWO)));
118 if !halved {
119 return Ok(Refine::Done(best.expect("a fit to compare with").0));
120 }
121 if legs >= MAX_FITTED_LEGS || !width.definitely_greater(Scalar::from_f64(FIT_DRIFT))
122 {
123 return Ok(Refine::Done(curve));
124 }
125 // The width is twice the drift it holds. Halved, a leg's
126 // drift falls sixteen times: those that would then still be
127 // the widest are all that need it.
128 let drifts =
129 NurbCurve3D::drift_by_interval(&points, &between, 3).with_context(&ctx)?;
130 let enough = width.upper().to_f64() / (2.0 * 16.0);
131 let mut finer = vec![dense[0].clone()];
132 for i in 0..legs {
133 if drifts[i] > enough {
134 let half = Scalar::from_ratio(1, 2).with_context(&ctx)?;
135 finer.push(locate(&dense[i], &dense[i + 1], half).with_context(&ctx)?);
136 }
137 finer.push(dense[i + 1].clone());
138 }
139 Ok(Refine::Finer((finer, Some((curve, width)))))
140 },
141 )?;
142 let Refined::Done(curve) = refined else {
143 return Err(ctx(GeopError::new(
144 "the fitted curve's legs were split past every cap",
145 )));
146 };
147 Ok(curve)
148 }
149}
150
151#[cfg(test)]
152mod tests {
153 use super::*;
154
155 /// A quarter of an ellipse nineteen times as long as it is wide, marched
156 /// in legs of one length, ten times its tip's radius of curvature at the
157 /// tip: a blend's arc cut by a plane at three degrees to its contact.
158 /// Only the few legs at the tip drift. Halving them all would soon pass
159 /// the legs a fit may have; those few bring the width down.
160 #[test]
161 fn a_fit_refines_the_legs_where_the_curve_turns_sharply() {
162 let (a, b) = (1.9, 0.1);
163 let at = |t: f64| {
164 Vector3::from_array([
165 Scalar::from_f64(a * t.cos()),
166 Scalar::from_f64(b * t.sin()),
167 Scalar::ZERO,
168 ])
169 };
170 // Stations a length apart along it, from a table of its length.
171 let table = 20000;
172 let top = std::f64::consts::FRAC_PI_2;
173 let mut lengths = vec![0.0];
174 for k in 1..=table {
175 let (t0, t1) = (
176 top * (k - 1) as f64 / table as f64,
177 top * k as f64 / table as f64,
178 );
179 let d =
180 ((a * (t1.cos() - t0.cos())).powi(2) + (b * (t1.sin() - t0.sin())).powi(2)).sqrt();
181 lengths.push(lengths[k - 1] + d);
182 }
183 let strides = 108;
184 let stations: Vec<(Vector3, f64)> = (0..=strides)
185 .map(|k| {
186 let want = lengths[table] * k as f64 / strides as f64;
187 let i = lengths.partition_point(|&l| l < want).min(table);
188 let t = top * i as f64 / table as f64;
189 (at(t), t)
190 })
191 .collect();
192 let curve = NurbCurve3D::fit_enclosing(
193 &stations,
194 |from, to, frac| {
195 let t = from.1 + frac.to_f64() * (to.1 - from.1);
196 Ok((at(t), t))
197 },
198 &mut None,
199 &|e| e,
200 )
201 .unwrap();
202 let widest = (curve.control_points.iter())
203 .flat_map(|p| (0..3).map(move |k| p[k].width().to_f64()))
204 .fold(0.0, f64::max);
205 assert!(widest <= FIT_DRIFT, "{widest:e}");
206 // Two legs a stride, and a few more at the tip.
207 assert!(
208 curve.control_points.len() < 300,
209 "{}",
210 curve.control_points.len()
211 );
212 }
213}