From 3692bc9914d5a1f99f43ff8488a396e2008a5972 Mon Sep 17 00:00:00 2001 From: Alexander Schrab Date: Sun, 25 Oct 2026 23:21:20 +0100 Subject: [PATCH] perf: B-spline curve evaluation into a caller buffer (outputs bit-identical) BSplineCurve::evaluate returned a fresh Vec from crv_derivs on every call. On the #948 loft about 180 profile samples are that allocation, almost all from FindCurveDistanceExtrema::evaluate (two evaluations of 3 entries per step of every curve-cache distance solve). - nurbs_crv::crv_derivs_into writes C^(k)(u) into the caller's slice; crv_derivs allocates and calls it. - BSplineCurve::evaluate_into holds the old evaluate body (clamping, zero-tangent repair) over a slice; evaluate allocates and calls it, and keeps NLib's order error for num_derivs >= GW_MAX_DERIV without allocating an output of that size. - Curve::evaluate_into dispatches the B-spline kinds (B-spline, hyperbola, ellipse, parabola) to it; other kinds fall back to evaluate and fill missing entries with Vec3::UNDEF, as at() reads them. - FindCurveDistanceExtrema::evaluate uses two [Vec3; 3] buffers. A new test compares evaluate_into with evaluate bit for bit, errors included, for 0 to GW_MAX_DERIV + 1 derivatives, with and without the zero-tangent repair (one curve has a zero first derivative). Loft driver, 3 interleaved runs, P-cores, medians: user instructions -6.1 % and cycles -4.8 % / -4.6 % (20 / 40 sections); all non-timing output identical. Co-Authored-By: Claude Opus 5.5 --- crates/solidcore/src/curve/curve/p2.rs | 45 +++++++++ crates/solidcore/src/curve/curve_cache.rs | 6 +- .../solidcore/src/nurbs/b_spline_curve/p1.rs | 97 ++++++++++++++++--- crates/solidcore/src/nurbs/nurbs_crv.rs | 42 +++++--- 4 files changed, 165 insertions(+), 25 deletions(-) diff --git a/crates/solidcore/src/curve/curve/p2.rs b/crates/solidcore/src/curve/curve/p2.rs index a3c0efbbe..05db3bb79 100644 --- a/crates/solidcore/src/curve/curve/p2.rs +++ b/crates/solidcore/src/curve/curve/p2.rs @@ -1274,6 +1274,51 @@ impl Curve { } } + /// [`Curve::evaluate`] into the caller's slice (`num_derivs = out.len() - 1`): + /// the entries the evaluation returns are written, the others are set to + /// [`Vec3::UNDEF`] (what [`at`] reads for a missing entry). The B-spline kinds + /// evaluate without a heap allocation ([`BSplineCurve::evaluate_into`]). + /// + /// # Errors + /// + /// As [`Curve::evaluate`]; `SM_ERR` for an empty `out`. + pub fn evaluate_into( + &self, + k: &Kernel, + t: f64, + from_left: bool, + non_zero_tangents: bool, + out: &mut [Vec3], + ) -> Result<()> { + let base = match &self.kind { + CurveKind::BSpline(c) => Some(c), + CurveKind::Hyperbola(c) => Some(c.base()), + CurveKind::Ellipse(c) => Some(c.base()), + CurveKind::Parabola(c) => Some(c.base()), + CurveKind::Composite(_) + | CurveKind::Hermite(_) + | CurveKind::Conic(_) + | CurveKind::Projected(_) + | CurveKind::Iso(_) + | CurveKind::Offset(_) + | CurveKind::CrvInVolume(_) + | CurveKind::CrvOnSurf(_) + | CurveKind::Line(_) => None, + }; + if let Some(b) = base { + return b.evaluate_into(t, from_left, non_zero_tangents, out); + } + let num_derivs = out + .len() + .checked_sub(1) + .ok_or_else(|| SolidError::sm_err(OP_EVAL))?; + let pv = self.evaluate(k, t, num_derivs, from_left, non_zero_tangents)?; + for (i, slot) in out.iter_mut().enumerate() { + *slot = at(&pv, i); + } + Ok(()) + } + /// The point at `t` (`EvaluatePoint`). /// /// # Errors diff --git a/crates/solidcore/src/curve/curve_cache.rs b/crates/solidcore/src/curve/curve_cache.rs index 45eb084ab..58f085abf 100644 --- a/crates/solidcore/src/curve/curve_cache.rs +++ b/crates/solidcore/src/curve/curve_cache.rs @@ -1563,8 +1563,10 @@ impl EvalFunction for FindCurveDistanceExtrema<'_> { out.f_prime = 0.0; out.found_answer = false; - let pv1 = self.curve1.evaluate(k, t, 2, true, true)?; - let pv2 = self.curve2.evaluate(k, t, 2, true, true)?; + let mut pv1 = [Vec3::ZERO; 3]; + let mut pv2 = [Vec3::ZERO; 3]; + self.curve1.evaluate_into(k, t, true, true, &mut pv1)?; + self.curve2.evaluate_into(k, t, true, true, &mut pv2)?; let (a0, a1, a2) = (at(&pv1, 0), at(&pv1, 1), at(&pv1, 2)); let (b0, b1, b2) = (at(&pv2, 0), at(&pv2, 1), at(&pv2, 2)); diff --git a/crates/solidcore/src/nurbs/b_spline_curve/p1.rs b/crates/solidcore/src/nurbs/b_spline_curve/p1.rs index d987c0ebe..8b2019efc 100644 --- a/crates/solidcore/src/nurbs/b_spline_curve/p1.rs +++ b/crates/solidcore/src/nurbs/b_spline_curve/p1.rs @@ -138,8 +138,8 @@ use crate::nurbs::iges::{iges_126_curve, Entity126}; use crate::nurbs::knots::Side; use crate::nurbs::math::{CPoint, MatrixType, RealMatrix, NOW, NOZ}; use crate::nurbs::nurbs_crv::{ - allocate_nurb_curve, copy_nurb_curve, crv_derivs, crv_eval, crv_remove_knots, extract_frame, - fixup_knot_vector, get_knots, to_euclid, + allocate_nurb_curve, copy_nurb_curve, crv_derivs, crv_derivs_into, crv_eval, crv_remove_knots, + extract_frame, fixup_knot_vector, get_knots, to_euclid, MAX_DERIV, }; const OP_DATA: &str = "nurbs::b_spline_curve::p1::BSplineCurve::from_data"; @@ -965,6 +965,43 @@ impl BSplineCurve { if num_derivs == 0 { return Ok(vec![self.evaluate_point(t)?]); } + if num_derivs >= MAX_DERIV { + // (NLib's order error, after the NURB check, without allocating an output + // of that size) + let cur = self.nurb_ref()?; + let side = if from_left { Side::Left } else { Side::Right }; + return crv_derivs(cur, t, side, num_derivs).map_err(SolidError::nl_ser(OP_EVAL)); + } + let mut pv = vec![Vec3::ZERO; num_derivs + 1]; + self.evaluate_into(t, from_left, non_zero_tangents, &mut pv)?; + Ok(pv) + } + + /// [`BSplineCurve::evaluate`] into the caller's slice: `out[k]` is derivative `k`, + /// `num_derivs = out.len() - 1` (below `GW_MAX_DERIV`), without a heap allocation + /// (except for the rare zero-tangent re-evaluation). + /// + /// # Errors + /// + /// As [`BSplineCurve::evaluate`]; `SM_ERR` for an empty `out`. + pub fn evaluate_into( + &self, + t: f64, + from_left: bool, + non_zero_tangents: bool, + out: &mut [Vec3], + ) -> Result<()> { + let num_derivs = out + .len() + .checked_sub(1) + .ok_or_else(|| SolidError::sm_err(OP_EVAL))?; + // pass simple calls along to EvaluatePoint + if num_derivs == 0 { + if let Some(slot) = out.first_mut() { + *slot = self.evaluate_point(t)?; + } + return Ok(()); + } // natural extent let cur = self.nurb_ref()?; @@ -978,19 +1015,16 @@ impl BSplineCurve { } let side = if from_left { Side::Left } else { Side::Right }; - let mut pv = crv_derivs(cur, t, side, num_derivs).map_err(SolidError::nl_ser(OP_EVAL))?; - if pv.len() != num_derivs + 1 { - return Err(SolidError::sm_err(OP_EVAL)); - } + crv_derivs_into(cur, t, side, out).map_err(SolidError::nl_ser(OP_EVAL))?; // When asked, find a direction for zero-length first derivatives from the // 2nd derivative value. if non_zero_tangents { - let p0 = pv.first().copied().unwrap_or(Vec3::UNDEF); + let p0 = out.first().copied().unwrap_or(Vec3::UNDEF); let scaled_zero = EFF_ZERO * (1.0 + p0.max_dimension()); let new_len = 1.1 * scaled_zero; - let d1 = pv.get(1).copied().unwrap_or(Vec3::UNDEF); + let d1 = out.get(1).copied().unwrap_or(Vec3::UNDEF); if num_derivs >= 1 && d1.length_squared() < new_len * new_len { if num_derivs < 2 { // recurse so the next section sets the 1st derivative; the call @@ -999,11 +1033,11 @@ impl BSplineCurve { let vals = self .evaluate(t, 2, from_left, false) .unwrap_or_else(|_| vec![Vec3::UNDEF; 3]); - if let (Some(slot), Some(&v)) = (pv.get_mut(1), vals.get(1)) { + if let (Some(slot), Some(&v)) = (out.get_mut(1), vals.get(1)) { *slot = v; } } else { - let d2 = pv.get(2).copied().unwrap_or(Vec3::UNDEF); + let d2 = out.get(2).copied().unwrap_or(Vec3::UNDEF); let len = d2.length(); if len > new_len { let mut d1 = (new_len / len) * d2; @@ -1012,14 +1046,14 @@ impl BSplineCurve { if t > self.natural_interval().mid() { d1 *= -1.0; } - if let Some(slot) = pv.get_mut(1) { + if let Some(slot) = out.get_mut(1) { *slot = d1; } } } } } - Ok(pv) + Ok(()) } /// Tessellates the curve over `interval` (`Tessellate`): bisection with the chord @@ -1824,6 +1858,45 @@ mod tests { BSplineCurve::from_nurb(Curve::new(pw, 3, u).unwrap(), 3, false) } + #[test] + fn evaluate_into_matches_evaluate_bit_for_bit() { + let bits = |v: &[Vec3]| -> Vec<[u64; 3]> { + v.iter() + .map(|p| [p.x.to_bits(), p.y.to_bits(), p.z.to_bits()]) + .collect() + }; + // the second curve has a zero first derivative at t = 0 (repeated point) + let pw = vec![ + CPoint::new(0.0, 0.0, 0.0, NOW), + CPoint::new(0.0, 0.0, 0.0, NOW), + CPoint::new(3.0, 2.0, 1.0, NOW), + CPoint::new(4.0, 0.0, 0.0, NOW), + ]; + let u = vec![0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0]; + let cusp = BSplineCurve::from_nurb(Curve::new(pw, 3, u).unwrap(), 3, false); + for c in [cubic(), cusp] { + for t in [-0.25, 0.0, 0.3, 0.5, 0.75, 1.0, 1.5] { + for from_left in [true, false] { + for nzt in [true, false] { + for nd in 0..MAX_DERIV + 2 { + let full = c.evaluate(t, nd, from_left, nzt); + let mut out = vec![Vec3::new(7.0, 7.0, 7.0); nd + 1]; + let r = c.evaluate_into(t, from_left, nzt, &mut out); + match full { + Ok(full) => { + assert!(r.is_ok()); + assert_eq!(bits(&full), bits(&out)); + } + Err(e) => assert_eq!(r.err(), Some(e)), + } + } + } + } + } + } + assert!(cubic().evaluate_into(0.5, true, false, &mut []).is_err()); + } + #[test] fn empty_curve_has_no_nurb() { let mut c = BSplineCurve::new(2); diff --git a/crates/solidcore/src/nurbs/nurbs_crv.rs b/crates/solidcore/src/nurbs/nurbs_crv.rs index 963f984dd..437b99d3e 100644 --- a/crates/solidcore/src/nurbs/nurbs_crv.rs +++ b/crates/solidcore/src/nurbs/nurbs_crv.rs @@ -20,7 +20,7 @@ //! | [`crv_get_epolygon`] | `sm_CrvGetEPolygon` | //! | [`basis_derivs`], [`BasisDerivs`] | `sm_BasisDerivs` | //! | [`crv_rational_basis_derivs`] | `sm_CrvRationalBasisDerivs` | -//! | [`crv_derivs`] | `sm_CrvDerivs` | +//! | [`crv_derivs`], [`crv_derivs_into`] | `sm_CrvDerivs` | //! | [`basis_eval`] | `sm_BasisEval` | //! | [`crv_eval`] | `sm_CrvEval` | //! | `to_euclid` (crate) | `TO_EUCLID` (`SmNurbsCrv.h`) | @@ -979,19 +979,40 @@ pub fn crv_rational_basis_derivs( /// [`basis_derivs`] and with [`NlibCode::Index`] when the span's control points are /// missing. pub fn crv_derivs(cur: &Curve, u: f64, side: Side, der: usize) -> Result> { - const OP: &str = OP_CDERIVS; if der >= MAX_DERIV { - return Err(nlib(OP, NlibCode::MaxDerivative)); + return Err(nlib(OP_CDERIVS, NlibCode::MaxDerivative)); } + let mut cd = vec![Vec3::ZERO; der + 1]; + crv_derivs_into(cur, u, side, &mut cd)?; + Ok(cd) +} + +/// [`crv_derivs`] into the caller's slice: `out[k] = C^(k)(u)` for `k < out.len()` +/// (`der = out.len() - 1`), without a heap allocation. +/// +/// Errors as [`crv_derivs`]; an empty `out` is [`NlibCode::MaxDerivative`] as well. +pub fn crv_derivs_into(cur: &Curve, u: f64, side: Side, out: &mut [Vec3]) -> Result<()> { + const OP: &str = OP_CDERIVS; + let der = match out.len().checked_sub(1) { + Some(der) if der < MAX_DERIV => der, + _ => return Err(nlib(OP, NlibCode::MaxDerivative)), + }; if cur.degree() < SMALL_DEGREE { - crv_derivs_n::(cur, u, side, der) + crv_derivs_n::(cur, u, side, der, out) } else { - crv_derivs_n::(cur, u, side, der) + crv_derivs_n::(cur, u, side, der, out) } } -/// [`crv_derivs`] with scratch width `N` (see [`basis_derivs_into`]); `der < MAX_DERIV`. -fn crv_derivs_n(cur: &Curve, u: f64, side: Side, der: usize) -> Result> { +/// [`crv_derivs_into`] with scratch width `N` (see [`basis_derivs_into`]); +/// `der < MAX_DERIV`, `out.len() == der + 1`. +fn crv_derivs_n( + cur: &Curve, + u: f64, + side: Side, + der: usize, + out: &mut [Vec3], +) -> Result<()> { const OP: &str = OP_CDERIVS; let pu = cur.degree(); let mut bd = [[0.0_f64; N]; MAX_DERIV]; @@ -1017,8 +1038,7 @@ fn crv_derivs_n(cur: &Curve, u: f64, side: Side, der: usize) -> let pts = pts.get(..npts).ok_or_else(|| oob(OP))?; // C^(k) = Sum_i(BasisVal[k][i] * P[i]). - let mut cd = Vec::with_capacity(der + 1); - for row in bd.get(..=der).ok_or_else(|| oob(OP))? { + for (row, slot) in bd.get(..=der).ok_or_else(|| oob(OP))?.iter().zip(out) { let mut c = Vec3::new(0.0, 0.0, 0.0); for (i, q) in pts.iter().enumerate().take(pu + 1) { let b = row.get(i).copied().ok_or_else(|| oob(OP))?; @@ -1026,9 +1046,9 @@ fn crv_derivs_n(cur: &Curve, u: f64, side: Side, der: usize) -> c.y += b * q.y; c.z += b * q.z; } - cd.push(c); + *slot = c; } - Ok(cd) + Ok(()) } /// All non-vanishing basis functions of degree `p` at `u` (`sm_BasisEval`):