geop_core_geometry/nurb_surface/
fit_pcurve.rs

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

geop_core_geometry/nurb_surface/
fit_pcurve.rs
