geop_core_geometry/nurb_surface/

fit_pcurve.rs

use crate::{
    nurb_curve::{NurbCurve, NurbCurve2D, true_point_fractions},
    spline::interior_knots,
};
use geop_core_math::{
    geop_error::{GeopError, GeopResult, WithContext},
    scalars::Scalar,
    solvers::refine::{Refine, Refined, refine},
    vector::Vector2,
};

use super::NurbSurface;

/// How many intervals a curve is split into for projection to start with,
/// shared out among its smooth pieces by their share of its domain (but at
/// least [`MIN_PIECE_SAMPLES`] each). Each piece's pcurve passes exactly
/// through its samples and is widened to enclose the trace between them
/// (see `NurbCurve2D::interpolate_enclosing`).
///
/// A cubic interpolant's drift falls as `h^4` on a smooth trace, and since
/// the pcurve carries that drift as width, this decides how *wide* the
/// pcurve is, not whether it is right. At 9 samples the drift reached ~2e-5
/// on strongly curved patches; 48 keeps it well inside the accuracy
/// `validation::numerical_accuracy` holds every entity to. Where it does
/// not, a piece is sampled more densely (see [`DRIFT_PER_RESOLUTION`]).
///
/// This bounds effort, not correctness: each sample is one Newton foot-point
/// projection (plus one per true point between samples, see
/// `true_point_fractions`), and a pcurve fitted through more of them is
/// strictly narrower.
const SAMPLES: usize = 48;

/// The fewest intervals a smooth piece of a curve is sampled in.
const MIN_PIECE_SAMPLES: usize = 8;

/// The most intervals a smooth piece of a curve is sampled in. A pcurve
/// still drifting then has a trace that is not smooth where the curve's
/// knots say nothing (a seam of the surface): it is returned as it is, wide
/// and honest.
const MAX_PIECE_SAMPLES: usize = 1536;

/// A piece is sampled twice as densely while its pcurve is wider than
/// `min_subdivision_size` over this (and than twice its widest sample, a
/// width no sampling removes): a width that small is lost in any
/// search resolving to that size, which is what every question later asked
/// of the pcurve is. Like `min_subdivision_size` itself, it decides effort,
/// not correctness: the pcurve encloses its trace however wide it is.
const DRIFT_PER_RESOLUTION: i64 = 16;

impl NurbSurface<4> {
    /// The `(u, v)` trace of `curve` across this surface: sample the curve,
    /// Newton-project each sample onto the surface (each projection seeded
    /// from the previous one's result, so the walk stays continuous), and
    /// fit a pcurve through the results.
    ///
    /// The pcurve is fitted piece by piece, the curve broken at each knot
    /// where it is less than twice continuously differentiable — the joints
    /// of the kernel's rational arcs and helices, which are only `C1` — and
    /// the pieces joined. The fit is a cubic spline, `C2` at its own knots:
    /// across a jump in the trace's curvature it cannot follow, and its
    /// drift there shrinks only as `h^2` with the spacing of the samples. A
    /// helix of pitch 1 fitted across its joints drifted past the accuracy
    /// limit. Within a piece, the samples are doubled while the pcurve is
    /// wider than the target (see [`DRIFT_PER_RESOLUTION`]).
    ///
    /// The walk's first seed has to be found globally, because Newton only
    /// polishes a foot point it is already near: seeded from anywhere else
    /// on a curved patch it settles on whichever local foot point is closest
    /// — including one clamped against a domain bound, where the residual is
    /// merely orthogonal to the boundary — and every later sample, seeded
    /// from its predecessor, follows it there. So the start is `pin_start`
    /// when there is one (the face's own authoritative `(u, v)`), and is
    /// otherwise located by [`Self::locate`] (`max_nodes` /
    /// `min_subdivision_size` bound its search). A curve that does not start on this surface is
    /// an error: there is no trace of it to fit.
    ///
    /// `pin_start` / `pin_end` override the projected `(u, v)` of the first
    /// and last sample. Pass them whenever the curve's endpoint is a place
    /// this surface's face *already* has a coedge for: that coedge's own
    /// pcurve endpoint is the authoritative `(u, v)` there, and an
    /// independently re-projected one lands a hair away from it — enough to
    /// break the exact `could_be_equal` continuity a face's boundary loop
    /// requires between one coedge's pcurve end and the next one's start.
    /// The endpoints are free to pin without disturbing the rest of the
    /// curve because `interpolate` produces a clamped B-spline, which
    /// passes exactly through each sample.
    pub fn fit_pcurve(
        &self,
        curve: &NurbCurve<4>,
        pin_start: Option<Vector2>,
        pin_end: Option<Vector2>,
        max_nodes: usize,
        min_subdivision_size: Scalar,
    ) -> GeopResult<NurbCurve2D> {
        let ctx = |e: GeopError| {
            e.with_context(format!(
                "NurbSurface::fit_pcurve: curve={curve:?}, domain_u={:?}, domain_v={:?}",
                self.domain_u(),
                self.domain_v(),
            ))
        };

        let (t0, t1) = curve.domain();
        // Seed the first projection where the curve actually starts on this
        // surface (see above); every later one seeds from its predecessor.
        // Sharp, since a seed is a free choice.
        let (seed_u, seed_v) = match pin_start {
            Some(pin) => (pin[0], pin[1]),
            None => {
                let start = curve.evaluate(t0).with_context(&ctx)?;
                self.locate(&start, max_nodes, min_subdivision_size)
                    .with_context(&ctx)?
                    .ok_or_else(|| {
                        ctx(GeopError::new(format!(
                            "the curve starts at {start:?}, which is not on this surface"
                        )))
                    })?
            }
        };
        let mut seed = (seed_u.sharpen(), seed_v.sharpen());

        // The foot point of the curve at `frac` of the way from `a` to `b`,
        // projected from `seed`, which then moves on to it. Seeded from a
        // sharp value (any point inside the previous iterate is an equally
        // valid starting guess), but the projection's own enclosure is
        // what's returned: these `(u, v)` end up in the fitted pcurve, which
        // is later compared against the edge's 3D points, so narrowing them
        // here would claim precision the projection did not have.
        let project_at = |seed: &mut (Scalar, Scalar),
                          (a, b): (Scalar, Scalar),
                          frac: Scalar|
         -> GeopResult<Vector2> {
            let t = a.add(b.sub(a).mul(frac));
            let p = curve.evaluate(t)?;
            let (u, v) = self.project(p, seed.0, seed.1, NEWTON_ITERATIONS)?;
            *seed = (u.sharpen(), v.sharpen());
            Ok(Vector2::from_array([u, v]))
        };

        // Where the curve is less than `C2`: the ends of its domain, and
        // every interior knot of multiplicity at least `degree − 1`.
        let p = curve.degree;
        let mut breaks = vec![t0];
        breaks.extend(
            interior_knots(&curve.knot_vector, p, curve.control_points.len())
                .into_iter()
                .filter(|&(_, m)| m + 1 >= p)
                .map(|(k, _)| k),
        );
        breaks.push(t1);

        let target = min_subdivision_size
            .div(Scalar::from_i64(DRIFT_PER_RESOLUTION))
            .with_context(&ctx)?
            .upper();
        let length = t1.sub(t0).to_f64();
        let mut pieces = Vec::with_capacity(breaks.len() - 1);
        // Where the next piece starts: the pin, then where the last one
        // ended — the same point projected again would land a hair away.
        let mut start = pin_start;
        for (k, piece) in breaks.windows(2).enumerate() {
            let (a, b) = (piece[0], piece[1]);
            let end = if k + 2 == breaks.len() { pin_end } else { None };
            let piece_seed = seed;
            // How many samples is a free choice: the piece's share of 48.
            let share = b.sub(a).to_f64() / length;
            let first =
                ((SAMPLES as f64 * share).ceil() as usize).clamp(MIN_PIECE_SAMPLES, SAMPLES);
            // The cap ends it: `intervals` doubles each step, so it passes
            // any cap within as many steps as it has bits.
            let refined = refine(first, usize::BITS as usize, |intervals| -> GeopResult<_> {
                // Each attempt walks the piece afresh, from its start.
                seed = piece_seed;
                // The samples the pcurve passes through, and — walking the
                // same path, so every projection seeds from its neighbour —
                // the true trace between each consecutive pair, which the
                // pcurve is widened to enclose
                // (`NurbCurve2D::interpolate_enclosing`): an interpolant
                // drifts from its trace between samples, and that drift is
                // part of what the pcurve honestly knows about where the
                // trace is.
                let count = Scalar::from_i64(intervals as i64);
                let mut uvs = Vec::with_capacity(intervals + 1);
                let mut between = Vec::with_capacity(intervals);
                for i in 0..=intervals {
                    let frac = Scalar::from_i64(i as i64).div(count).with_context(&ctx)?;
                    uvs.push(project_at(&mut seed, (a, b), frac).with_context(&ctx)?);
                    if i < intervals {
                        let fractions = true_point_fractions(i, intervals);
                        let mut inside = Vec::with_capacity(fractions.len());
                        for &(num, den) in fractions {
                            let frac =
                                Scalar::from_ratio(i as i64 * den + num, intervals as i64 * den)
                                    .with_context(&ctx)?;
                            inside.push(project_at(&mut seed, (a, b), frac).with_context(&ctx)?);
                        }
                        between.push(inside);
                    }
                }
                // Still projected, as the walk's seed.
                if let Some(at) = start {
                    uvs[0] = at;
                }
                if let Some(pin) = end {
                    *uvs.last_mut().expect("uvs is never empty") = pin;
                }
                let fitted =
                    NurbCurve2D::interpolate_enclosing(&uvs, &between, 3).with_context(&ctx)?;
                // Denser samples narrow the pcurve only by its drift: as
                // wide as its samples are, it stays however many there are.
                // So it is refined while it is wider than the target and
                // more than twice as wide as its widest sample.
                let sampled = uvs
                    .iter()
                    .chain(between.iter().flatten())
                    .fold(Scalar::ZERO, |w, p| w.max(p[0].width()).max(p[1].width()));
                let fitted_width = fitted.control_points.iter().fold(Scalar::ZERO, |w, cp| {
                    w.max(cp[0].width()).max(cp[1].width())
                });
                let drifting = fitted_width.definitely_greater(target)
                    && fitted_width.definitely_greater(sampled.add(sampled));
                if !drifting || intervals * 2 > MAX_PIECE_SAMPLES {
                    return Ok(Refine::Done((fitted, uvs.last().copied())));
                }
                Ok(Refine::Finer(intervals * 2))
            })?;
            let Refined::Done((fitted, last)) = refined else {
                return Err(ctx(GeopError::new(
                    "the samples of a piece doubled past every cap",
                )));
            };
            start = last;
            pieces.push(fitted);
        }
        NurbCurve2D::join(&pieces).with_context(&ctx)
    }

    /// The pcurve of `curve`, which lies on this surface, where the
    /// surface is affine in `(u, v)` — degree one both ways, its four
    /// control points a parallelogram of one weight: then `(u, v)` are an
    /// affine function of the point, and the pcurve is `curve` itself with
    /// each control point mapped to them, of the same weights and knots.
    /// Exact, where a fit only encloses. None for any other surface.
    pub fn affine_pcurve(&self, curve: &NurbCurve<4>) -> GeopResult<Option<NurbCurve2D>> {
        let cps = &self.control_points;
        if self.degree_u != 1
            || self.degree_v != 1
            || self.num_u != 2
            || self.num_v != 2
            || cps.iter().any(|p| !p[3].is_exactly(cps[0][3]))
        {
            return Ok(None);
        }
        let point = |k: usize| cps[k].dehomogenized();
        // `cps[i * 2 + j]` is the corner at the `i`-th end of `u` and the
        // `j`-th of `v`.
        let origin = point(0)?;
        let along_u = point(2)?.sub(&origin);
        let along_v = point(1)?.sub(&origin);
        let skew = point(3)?.sub(&origin).sub(&along_u).sub(&along_v);
        if (0..3).any(|k| !skew[k].could_be_equal(Scalar::ZERO)) {
            return Ok(None);
        }
        let (uu, uv, vv) = (
            along_u.prod_dot(&along_u),
            along_u.prod_dot(&along_v),
            along_v.prod_dot(&along_v),
        );
        let det = uu.mul(vv).sub(uv.mul(uv));
        if !det.definitely_greater(Scalar::ZERO) {
            return Ok(None);
        }
        let ((u0, u1), (v0, v1)) = (self.domain_u(), self.domain_v());
        let control_points = curve
            .control_points
            .iter()
            .map(|q| {
                let weight = q[3];
                let d = q.dehomogenized()?.sub(&origin);
                let (du, dv) = (along_u.prod_dot(&d), along_v.prod_dot(&d));
                let s = vv.mul(du).sub(uv.mul(dv)).div(det)?;
                let t = uu.mul(dv).sub(uv.mul(du)).div(det)?;
                let u = u0.add(s.mul(u1.sub(u0)));
                let v = v0.add(t.mul(v1.sub(v0)));
                Ok(geop_core_math::vector::Vector3::from_array([
                    u.mul(weight),
                    v.mul(weight),
                    weight,
                ]))
            })
            .collect::<GeopResult<Vec<_>>>()?;
        NurbCurve::try_new(curve.degree, control_points, curve.knot_vector.clone()).map(Some)
    }
}

/// Newton iteration count for each sample's foot-point projection. Unlike
/// `max_nodes`/`min_subdivision_size` this doesn't decide whether a search
/// converges to a correct-or-error answer, only how tightly a
/// fixed-iteration projection tracks its target, so it's a constant rather
/// than a threaded parameter.
const NEWTON_ITERATIONS: usize = 20;

#[cfg(test)]
mod tests {
    use geop_core_math::{
        primitives::CoordinateSystem,
        scalars::Scalar,
        vector::{Vector3, Vector4},
    };

    use crate::{
        nurb_curve::{Handedness, Helix, NurbCurve3D},
        nurb_surface::NurbSurface,
        shape::Axis,
    };

    /// The pcurve of the kernel's own helix — radius 1, pitch 1, almost a
    /// full turn, so its quarter-turn spans meet at three joints where it is
    /// only `C1` — on the cylinder it lies on: as narrow as a smooth trace's,
    /// and on the helix at both ends. Fitted as one cubic across the joints
    /// it drifted by about 1e-4, past what validation allows any entity to
    /// carry (`a_strip_turning_more_than_once_is_cut_along_meridians` in
    /// `geop-ops-step` had to use a pitch of a quarter).
    #[test]
    fn a_helix_is_fitted_between_its_joints() {
        let f = Scalar::from_f64;
        let v = |x: f64, y: f64, z: f64| Vector3::from_array([f(x), f(y), f(z)]);
        let helix = NurbCurve3D::helix(
            &CoordinateSystem::world_at(v(0.0, 0.0, 0.0)),
            Scalar::ONE,
            &Helix::new(Scalar::ONE, 0.9, Handedness::Right),
        )
        .unwrap();
        let line = NurbCurve3D::try_new(
            1,
            vec![
                Vector4::from_array([f(1.0), f(0.0), f(-0.5), Scalar::ONE]),
                Vector4::from_array([f(1.0), f(0.0), f(1.5), Scalar::ONE]),
            ],
            vec![f(0.0), f(0.0), f(1.0), f(1.0)],
        )
        .unwrap();
        // Its seam a little behind where the helix starts.
        let axis = Axis::try_new(v(0.0, 0.0, 0.0), v(0.0, 0.0, 1.0)).unwrap();
        let cylinder =
            NurbSurface::revolve(&line, &axis, -0.5, std::f64::consts::TAU - 0.5).unwrap();
        let pcurve = cylinder
            .fit_pcurve(&helix, None, None, 5000, f(1e-4))
            .unwrap();
        let widest = pcurve
            .control_points
            .iter()
            .flat_map(|cp| (0..2).map(|c| cp[c].width().to_f64()))
            .fold(0.0, f64::max);
        assert!(widest < 1e-5, "widest control point {widest:e}");
        let ((s0, s1), (t0, t1)) = (pcurve.domain(), helix.domain());
        for (s, t) in [(s0, t0), (s1, t1)] {
            let uv = pcurve.evaluate(s).unwrap();
            let on = cylinder.evaluate(uv[0], uv[1]).unwrap();
            assert!(on.could_be_equal(&helix.evaluate(t).unwrap()), "{on:?}");
        }
    }
}