geop_core_geometry/nurb_surface/fit_pcurve.rs
1use crate::{
2 contains::surface::surface_could_contain,
3 nurb_curve::{NurbCurve, NurbCurve2D, true_point_fractions},
4};
5use geop_core_math::{
6 geop_error::{GeopError, GeopResult, WithContext},
7 scalars::Scalar,
8 vector::Vector2,
9};
10
11use super::NurbSurface;
12
13/// How many intervals the curve is split into for projection; the fitted
14/// pcurve passes exactly through all `SAMPLES + 1` points and is widened to
15/// enclose the trace between them (see `NurbCurve2D::interpolate_enclosing`).
16///
17/// A cubic interpolant's drift falls as `h^4`, and since the pcurve carries
18/// that drift as width, this decides how *wide* the pcurve is, not whether it
19/// is right. At 9 samples the drift reached ~2e-5 on strongly curved patches
20/// and exceeded 1e-4 on the worst of them — too wide for the accuracy
21/// `validation::numerical_accuracy` holds every entity to. 48 keeps it well
22/// inside that.
23///
24/// This bounds effort, not correctness: each sample is one Newton foot-point
25/// projection (plus one per true point between samples, see
26/// `true_point_fractions`),
27/// and a pcurve fitted through more of them is strictly narrower.
28const SAMPLES: usize = 48;
29
30impl<S: Scalar> NurbSurface<S, 4> {
31 /// The `(u, v)` trace of `curve` across this surface: sample the curve,
32 /// Newton-project each sample onto the surface (each projection seeded
33 /// from the previous one's result, so the walk stays continuous), and
34 /// fit a pcurve through the results.
35 ///
36 /// The walk's first seed has to be found globally, because Newton only
37 /// polishes a foot point it is already near: seeded from anywhere else
38 /// on a curved patch it settles on whichever local foot point is closest
39 /// — including one clamped against a domain bound, where the residual is
40 /// merely orthogonal to the boundary — and every later sample, seeded
41 /// from its predecessor, follows it there. So the start is `pin_start`
42 /// when there is one (the face's own authoritative `(u, v)`), and is
43 /// otherwise isolated by [`surface_could_contain`] (`max_nodes` /
44 /// `min_subdivision_size` bound that search) — subdivide to isolate,
45 /// then Newton to refine. A curve that does not start on this surface is
46 /// an error: there is no trace of it to fit.
47 ///
48 /// `pin_start` / `pin_end` override the projected `(u, v)` of the first
49 /// and last sample. Pass them whenever the curve's endpoint is a place
50 /// this surface's face *already* has a coedge for: that coedge's own
51 /// pcurve endpoint is the authoritative `(u, v)` there, and an
52 /// independently re-projected one lands a hair away from it — enough to
53 /// break the exact `could_be_equal` continuity a face's boundary loop
54 /// requires between one coedge's pcurve end and the next one's start.
55 /// The endpoints are free to pin without disturbing the rest of the
56 /// curve because `interpolate` produces a clamped B-spline, which
57 /// passes exactly through each sample.
58 pub fn fit_pcurve(
59 &self,
60 curve: &NurbCurve<S, 4>,
61 pin_start: Option<Vector2<S>>,
62 pin_end: Option<Vector2<S>>,
63 max_nodes: usize,
64 min_subdivision_size: S,
65 ) -> GeopResult<NurbCurve2D<S>> {
66 let ctx = |e: GeopError| {
67 e.with_context(format!(
68 "NurbSurface::fit_pcurve: curve={curve:?}, domain_u={:?}, domain_v={:?}",
69 self.domain_u(),
70 self.domain_v(),
71 ))
72 };
73
74 let (t0, t1) = curve.domain();
75 // Seed the first projection where the curve actually starts on this
76 // surface (see above); every later one seeds from its predecessor.
77 // Sharp, since a seed is a free choice.
78 let (seed_u, seed_v) = match pin_start {
79 Some(pin) => (pin[0], pin[1]),
80 None => {
81 let start = curve.evaluate(t0).with_context(&ctx)?;
82 surface_could_contain(self, &start, max_nodes, min_subdivision_size)
83 .with_context(&ctx)?
84 .ok_or_else(|| {
85 ctx(GeopError::new(format!(
86 "the curve starts at {start:?}, which is not on this surface"
87 )))
88 })?
89 }
90 };
91 let (mut seed_u, mut seed_v) = (seed_u.sharpen(), seed_v.sharpen());
92
93 // The foot point of the curve at `frac` of its domain. Seeded from a
94 // sharp value (any point inside the previous iterate is an equally
95 // valid starting guess), but the projection's own enclosure is what's
96 // returned: these `(u, v)` end up in the fitted pcurve, which is later
97 // compared against the edge's 3D points, so narrowing them here would
98 // claim precision the projection did not have.
99 let mut project_at = |frac: S| -> GeopResult<Vector2<S>> {
100 let t = t0.add(t1.sub(t0).mul(frac));
101 let p = curve.evaluate(t)?;
102 let (u, v) = self.project(p, seed_u, seed_v, NEWTON_ITERATIONS)?;
103 seed_u = u.sharpen();
104 seed_v = v.sharpen();
105 Ok(Vector2::from_array([u, v]))
106 };
107
108 // The samples the pcurve passes through, and — walking the same path,
109 // so every projection seeds from its neighbour — the true trace
110 // between each consecutive pair, which the pcurve is widened to
111 // enclose (`NurbCurve2D::interpolate_enclosing`): an interpolant
112 // drifts from its trace between samples, and that drift is part of
113 // what the pcurve honestly knows about where the trace is.
114 let samples = S::from_i64(SAMPLES as i64);
115 let mut uvs = Vec::with_capacity(SAMPLES + 1);
116 let mut between = Vec::with_capacity(SAMPLES);
117 for i in 0..=SAMPLES {
118 uvs.push(
119 project_at(S::from_i64(i as i64).div(samples).with_context(&ctx)?)
120 .with_context(&ctx)?,
121 );
122 if i < SAMPLES {
123 let fractions = true_point_fractions(i, SAMPLES);
124 let mut inside = Vec::with_capacity(fractions.len());
125 for &(a, b) in fractions {
126 let frac =
127 S::from_ratio(i as i64 * b + a, SAMPLES as i64 * b).with_context(&ctx)?;
128 inside.push(project_at(frac).with_context(&ctx)?);
129 }
130 between.push(inside);
131 }
132 }
133
134 if let Some(pin) = pin_start {
135 uvs[0] = pin;
136 }
137 if let Some(pin) = pin_end {
138 *uvs.last_mut().expect("uvs is never empty") = pin;
139 }
140
141 NurbCurve2D::interpolate_enclosing(&uvs, &between, 3).with_context(&ctx)
142 }
143}
144
145/// Newton iteration count for each sample's foot-point projection. Unlike
146/// `max_nodes`/`min_subdivision_size` this doesn't decide whether a search
147/// converges to a correct-or-error answer, only how tightly a
148/// fixed-iteration projection tracks its target, so it's a constant rather
149/// than a threaded parameter.
150const NEWTON_ITERATIONS: usize = 20;