Skip to content

perf: B-spline curve evaluation into a caller buffer (outputs bit-identical) #1089

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
bherbruck merged 1 commit into from
Oct 26, 2026
Merged

perf: B-spline curve evaluation into a caller buffer (outputs bit-identical) #1089

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
45 changes: 45 additions & 0 deletions crates/solidcore/src/curve/curve/p2.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
6 changes: 4 additions & 2 deletions crates/solidcore/src/curve/curve_cache.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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));

Expand Down
97 changes: 85 additions & 12 deletions crates/solidcore/src/nurbs/b_spline_curve/p1.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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";
Expand Down Expand Up @@ -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()?;
Expand All @@ -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
Expand All @@ -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;
Expand All @@ -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
Expand Down Expand Up @@ -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);
Expand Down
42 changes: 31 additions & 11 deletions crates/solidcore/src/nurbs/nurbs_crv.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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`) |
Expand Down Expand Up @@ -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<Vec<Vec3>> {
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::<SMALL_DEGREE>(cur, u, side, der)
crv_derivs_n::<SMALL_DEGREE>(cur, u, side, der, out)
} else {
crv_derivs_n::<MAX_DEGREE>(cur, u, side, der)
crv_derivs_n::<MAX_DEGREE>(cur, u, side, der, out)
}
}

/// [`crv_derivs`] with scratch width `N` (see [`basis_derivs_into`]); `der < MAX_DERIV`.
fn crv_derivs_n<const N: usize>(cur: &Curve, u: f64, side: Side, der: usize) -> Result<Vec<Vec3>> {
/// [`crv_derivs_into`] with scratch width `N` (see [`basis_derivs_into`]);
/// `der < MAX_DERIV`, `out.len() == der + 1`.
fn crv_derivs_n<const N: usize>(
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];
Expand All @@ -1017,18 +1038,17 @@ fn crv_derivs_n<const N: usize>(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))?;
c.x += b * q.x;
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`):
Expand Down
Loading