mirror of
https://github.com/bherbruck/solidcraft.git
synced 2026-11-05 11:35:55 +00:00
kernel: spline edges on part of their curve load as that piece
SMLib bounds an intersection edge on a piece of its spline (seed 32 intersect:
a block and a cylinder, the top disc's edges), and the truck decode refused
such an edge ("not supported yet: a spline or intersection edge on part of its
curve"), since truck bounds an edge by its curve's whole range.
bspline::segment cuts a B-spline to a sub-interval of its domain exactly:
Boehm insertion of both ends until each has multiplicity p, then the knots and
control points between them, clamped on the interval. The same points at the
same parameters, so the edge keeps its interval and the record its numbers.
load_edge_curve loads a spline or rational spline edge on part of its curve as
that piece; an intersection edge on part of its curve is still refused. Our
own records (interval = the curve's range) load as before.
Tests: segments of a rational cubic with a double knot, inside a span, across
knots, from or to a domain end, match the whole curve within 1e-12; both
spline kinds load on a piece with that range and the same points, and beyond
the domain are refused; SMLib's block ∩ cylinder decodes (it was refused) with
π r² h within 1e-3.
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01QWhMwdqJYQKQq8RCLSL1LT
This commit is contained in:
@@ -343,6 +343,46 @@ pub fn insert_knot(p: usize, knots: &[f64], ctrl: &[HPoint], u: f64) -> SplineRe
|
||||
Ok((nk, q))
|
||||
}
|
||||
|
||||
/// The piece of a B-spline (degree `p`, knots, homogeneous control points) on `[a, b]`, a
|
||||
/// sub-interval of its domain, as a clamped B-spline of its own whose domain is exactly
|
||||
/// `[a, b]`: the same points at the same parameters. Exact up to rounding: Boehm insertion of
|
||||
/// `a` and `b` until each has multiplicity `p` (NURBS Book §5.3, curve splitting), then the
|
||||
/// knots and control points between them.
|
||||
pub fn segment(p: usize, knots: &[f64], ctrl: &[HPoint], a: f64, b: f64) -> SplineResult<(Vec<f64>, Vec<HPoint>)> {
|
||||
let n = ctrl.len();
|
||||
if n <= p || knots.len() != n + p + 1 {
|
||||
return invalid("knot and control point counts do not match the degree");
|
||||
}
|
||||
let (lo, hi) = (knots.get(p).copied().unwrap_or(f64::NAN), knots.get(n).copied().unwrap_or(f64::NAN));
|
||||
if !(lo <= a && a < b && b <= hi) {
|
||||
return invalid(format!("[{a}, {b}] is not a piece of the domain [{lo}, {hi}]"));
|
||||
}
|
||||
let (mut k, mut c) = (knots.to_vec(), ctrl.to_vec());
|
||||
for t in [a, b] {
|
||||
// At most p insertions each.
|
||||
for _ in 0..p {
|
||||
if k.iter().filter(|&&x| x == t).count() >= p {
|
||||
break;
|
||||
}
|
||||
(k, c) = insert_knot(p, &k, &c, t)?;
|
||||
}
|
||||
}
|
||||
// The curve at a (from the right) is control point j0 − p, j0 the last knot equal to a;
|
||||
// at b (from the left) control point j1 − 1, j1 the first knot equal to b.
|
||||
let j0 = k.iter().rposition(|&x| x == a).filter(|&j| j >= p).ok_or_else(|| SplineError::Invalid("segment start".into()))?;
|
||||
let j1 = k.iter().position(|&x| x == b).filter(|&j| j > j0).ok_or_else(|| SplineError::Invalid("segment end".into()))?;
|
||||
let points = c.get(j0 - p..j1).ok_or_else(|| SplineError::Invalid("segment control points".into()))?.to_vec();
|
||||
let inner = k.get(j0 + 1..j1).ok_or_else(|| SplineError::Invalid("segment knots".into()))?;
|
||||
let mut out = Vec::with_capacity(inner.len() + 2 * p + 2);
|
||||
out.extend(std::iter::repeat_n(a, p + 1));
|
||||
out.extend_from_slice(inner);
|
||||
out.extend(std::iter::repeat_n(b, p + 1));
|
||||
if out.len() != points.len() + p + 1 {
|
||||
return invalid("segment counts");
|
||||
}
|
||||
Ok((out, points))
|
||||
}
|
||||
|
||||
/// The binomial coefficient `(n k)` (zero for `k > n`).
|
||||
pub(crate) fn binomial(n: usize, k: usize) -> f64 {
|
||||
if k > n {
|
||||
@@ -714,6 +754,35 @@ pub(crate) mod tests {
|
||||
assert_eq!(find_span(&knots, p, n, f64::NAN), None);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn a_segment_is_the_same_curve_on_its_interval() {
|
||||
use crate::brep::{Curve, eval::CurveEval};
|
||||
let p = 3;
|
||||
let mut rng = Lcg(11);
|
||||
let ctrl: Vec<HPoint> = (0..KNOTS.len() - p - 1)
|
||||
.map(|_| HPoint::weighted(Vec3::new(rng.unit() * 10.0, rng.unit() * 10.0, rng.unit() * 10.0), 0.5 + rng.unit()))
|
||||
.collect();
|
||||
let record =
|
||||
|k: &[f64], c: &[HPoint]| Curve::RationalSpline { knots: k.to_vec(), homogeneous: c.iter().map(|h| [h.x, h.y, h.z, h.w]).collect() };
|
||||
let whole = CurveEval::new(&record(&KNOTS, &ctrl)).unwrap();
|
||||
// Inside one span, across knots (a double one included), from or to a domain end, and
|
||||
// from one existing knot to another.
|
||||
for (a, b) in [(0.1, 0.2), (0.12, 0.95), (0.0, 0.42), (0.3, 1.0), (0.3, 0.5), (0.0, 1.0)] {
|
||||
let (k, c) = segment(p, &KNOTS, &ctrl, a, b).unwrap();
|
||||
assert_eq!((k.first(), k.last()), (Some(&a), Some(&b)), "[{a}, {b}]: clamped on the interval");
|
||||
assert_eq!(k.len(), c.len() + p + 1);
|
||||
let piece = CurveEval::new(&record(&k, &c)).unwrap();
|
||||
for x in samples() {
|
||||
let t = a + (b - a) * x;
|
||||
let (q, r) = (whole.point(t).unwrap(), piece.point(t).unwrap());
|
||||
assert!(q.dist(r) < 1e-12, "[{a}, {b}] at {t}: {q:?} vs {r:?}");
|
||||
}
|
||||
}
|
||||
assert!(segment(p, &KNOTS, &ctrl, 0.5, 0.5).is_err());
|
||||
assert!(segment(p, &KNOTS, &ctrl, -0.1, 0.5).is_err());
|
||||
assert!(segment(p, &KNOTS, &ctrl, 0.2, f64::NAN).is_err());
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn single_function_knot_insertion_is_exact() {
|
||||
let mut rng = Lcg(7);
|
||||
|
||||
@@ -473,8 +473,9 @@ fn save_surface(s: &mt::Surface, d: usize) -> Result<b::Surface> {
|
||||
}
|
||||
|
||||
/// An edge's curve for truck, which bounds an edge by its curve's range: the record itself
|
||||
/// where the interval is that range (all this adapter writes), else the piece on the interval
|
||||
/// for the kinds that have one in closed form.
|
||||
/// where the interval is that range (all this adapter writes), else the piece on the interval:
|
||||
/// in closed form for lines and conics, cut out of a spline by knot insertion (the same points
|
||||
/// at the same parameters).
|
||||
fn load_edge_curve(c: b::Curve, interval: [f64; 2]) -> Result<mt::Curve> {
|
||||
match c {
|
||||
b::Curve::Line { start, end } if interval != [0., 1.] => {
|
||||
@@ -486,10 +487,35 @@ fn load_edge_curve(c: b::Curve, interval: [f64; 2]) -> Result<mt::Curve> {
|
||||
b::Curve::Line { .. } => load_curve(c),
|
||||
c => match c.natural() {
|
||||
b::Natural::Range(a, z) if [a, z] == interval => load_curve(c),
|
||||
_ => Err(unsupported("a spline or intersection edge on part of its curve")),
|
||||
b::Natural::Range(..) if matches!(c, b::Curve::Spline { .. } | b::Curve::RationalSpline { .. }) => {
|
||||
load_curve(spline_piece(&c, interval)?)
|
||||
}
|
||||
_ => Err(unsupported("an intersection edge on part of its curve")),
|
||||
},
|
||||
}
|
||||
}
|
||||
|
||||
/// A spline record's piece on `interval` (inside its domain), as a spline record of the same
|
||||
/// kind whose domain is the interval ([`solidcraft_geom::bspline::segment`]).
|
||||
fn spline_piece(c: &b::Curve, [a, z]: [f64; 2]) -> Result<b::Curve> {
|
||||
use solidcraft_geom::bspline::{HPoint, segment};
|
||||
let (knots, h): (&[f64], Vec<HPoint>) = match c {
|
||||
b::Curve::Spline { knots, points } => (knots, points.iter().map(|q| HPoint { x: q[0], y: q[1], z: q[2], w: 1.0 }).collect()),
|
||||
b::Curve::RationalSpline { knots, homogeneous } => {
|
||||
(knots, homogeneous.iter().map(|q| HPoint { x: q[0], y: q[1], z: q[2], w: q[3] }).collect())
|
||||
}
|
||||
_ => return Err(invalid("not a spline")),
|
||||
};
|
||||
let degree = knots.len().checked_sub(h.len() + 1).ok_or_else(|| invalid("spline degree"))?;
|
||||
let (knots, h) = segment(degree, knots, &h, a, z).map_err(invalid)?;
|
||||
Ok(match c {
|
||||
b::Curve::Spline { .. } => b::Curve::Spline {
|
||||
knots,
|
||||
points: h.iter().map(|q| q.project().map(|v| [v.x, v.y, v.z]).ok_or_else(|| invalid("spline piece"))).collect::<Result<_>>()?,
|
||||
},
|
||||
_ => b::Curve::RationalSpline { knots, homogeneous: h.iter().map(|q| [q.x, q.y, q.z, q.w]).collect() },
|
||||
})
|
||||
}
|
||||
fn load_curve(c: b::Curve) -> Result<mt::Curve> {
|
||||
Ok(match c {
|
||||
b::Curve::Line { start, end } => mt::Curve::Line(mt::Line(p(start), p(end))),
|
||||
|
||||
@@ -9,12 +9,38 @@ fn refuses_domain_or_axis_changes_instead_of_approximating() {
|
||||
let plane = b::Surface::Plane { origin: [1e12; 3], u_axis: [1e-6, 0., 0.], v_axis: [0., 1., 0.] };
|
||||
assert!(compatible_surface(&plane).is_err());
|
||||
assert!(clamped(&[0., 0., 1., 1.], 2).is_ok());
|
||||
// A spline edge on part of its curve would need truck to cut it: refused, not refitted.
|
||||
let spline = b::Curve::Spline { knots: vec![0., 0., 1., 1.], points: vec![[0.; 3], [1.; 3]] };
|
||||
assert!(load_edge_curve(spline.clone(), [0., 0.5]).is_err());
|
||||
assert!(load_edge_curve(spline, [0., 1.]).is_ok());
|
||||
}
|
||||
|
||||
/// A spline edge on part of its curve (SMLib bounds an intersection edge on a piece of its
|
||||
/// spline) loads as the piece cut out by knot insertion, not a refit: the same points at the
|
||||
/// same parameters, over exactly the interval.
|
||||
#[test]
|
||||
fn spline_edges_on_part_of_their_curve_load_as_that_piece() {
|
||||
use mt::ParametricCurve;
|
||||
let knots = vec![0., 0., 0., 0., 0.4, 0.4, 0.7, 1., 1., 1., 1.];
|
||||
let points: Vec<[f64; 3]> = (0..7).map(|i| [f64::from(i), (f64::from(i) * 1.3).sin() * 4., f64::from(i * i) * 0.2]).collect();
|
||||
let weights = [1., 0.8, 1.3, 0.9, 1.1, 0.7, 1.];
|
||||
let rational = b::Curve::RationalSpline {
|
||||
knots: knots.clone(),
|
||||
homogeneous: points.iter().zip(weights).map(|(q, w)| [q[0] * w, q[1] * w, q[2] * w, w]).collect(),
|
||||
};
|
||||
for record in [b::Curve::Spline { knots, points }, rational] {
|
||||
let exact = solidcraft_geom::brep::eval::CurveEval::new(&record).unwrap();
|
||||
for interval in [[0.1, 0.35], [0.25, 0.9], [0., 0.4], [0.4, 1.]] {
|
||||
let piece = load_edge_curve(record.clone(), interval).unwrap();
|
||||
assert_eq!(piece.range_tuple(), (interval[0], interval[1]), "{interval:?}");
|
||||
for k in 0..=20 {
|
||||
let t = interval[0] + (interval[1] - interval[0]) * f64::from(k) / 20.;
|
||||
let (q, r) = (piece.subs(t), exact.point(t).unwrap());
|
||||
assert!((q - p([r.x, r.y, r.z])).magnitude() < 1e-12, "{interval:?} at {t}");
|
||||
}
|
||||
}
|
||||
assert!(load_edge_curve(record, [0.5, 1.2]).is_err(), "beyond the curve's domain");
|
||||
}
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn analytic_records_decode_as_their_mathematical_loci_and_normals() {
|
||||
use mt::{MetricSpace, ParametricSurface, ParametricSurface3D};
|
||||
@@ -371,3 +397,26 @@ fn smlib_torus_decodes_without_meshing_and_plans_with_a_sphere_beside_it() {
|
||||
assert!((ours - exact).abs() < 1e-3 * exact, "{name}: decoded volume {ours}, exact {exact}");
|
||||
}
|
||||
}
|
||||
|
||||
/// SMLib's intersection of a block and a cylinder bounds the top disc with spline edges on
|
||||
/// part of their curves: the body decodes (it was refused) with the cylinder's volume inside
|
||||
/// the block, π r² h.
|
||||
#[cfg(feature = "smlib")]
|
||||
#[test]
|
||||
fn smlib_spline_edges_on_part_of_their_curve_decode() {
|
||||
use solidcraft_smlib::BooleanOp;
|
||||
let block = crate::box_solid(Vec3::new(7.5, 15., 0.), Vec3::new(25., 32.5, 27.5)).unwrap();
|
||||
let post = crate::cylinder(Vec3::new(12.5, 20., 10.), Vec3::Z, 5., 30.).unwrap();
|
||||
let (a, c) = (crate::native::smlib_body(&block).unwrap(), crate::native::smlib_body(&post).unwrap());
|
||||
let common = a.boolean(&c, BooleanOp::Intersect).unwrap().unwrap();
|
||||
let record = common.to_native().unwrap();
|
||||
assert!(
|
||||
record.edges.iter().any(|e| matches!(e.curve.natural(), b::Natural::Range(lo, hi) if [lo, hi] != e.interval)
|
||||
&& matches!(e.curve, b::Curve::Spline { .. } | b::Curve::RationalSpline { .. })),
|
||||
"a spline edge on part of its curve"
|
||||
);
|
||||
let body = crate::native::body_from_smlib(&common).unwrap();
|
||||
let exact = std::f64::consts::PI * 25. * 17.5;
|
||||
let ours = crate::measure(&body).unwrap().volume;
|
||||
assert!((ours - exact).abs() < 1e-3 * exact, "decoded volume {ours}, exact {exact}");
|
||||
}
|
||||
|
||||
Reference in New Issue
Block a user