masi/solvecraft
1
0
Fork 0
mirror of https://github.com/bherbruck/solidcraft.git synced 2026-11-05 19:45:54 +00:00
Code Issues Packages Projects Releases 1 Wiki Activity

kernel: SMLib sweeps along drawn 3D sketch splines, exactly

Browse Source 
A drawn 3D spline reached the kernel only as its sampled polyline (Wire.pts), so SMLib refused
to sweep along it ("not supported yet: SMLib sweep along a smooth path of points"). The doc
holds the spline's fit points; the curve is the centripetal Catmull-Rom spline that
fit_polyline samples. Each span of it is a cubic (Barry-Goldman), and its derivative at a fit
point depends only on that point, its neighbours and the two knot intervals, so neighbouring
spans share it: the spline is exactly one cubic Bézier per span, C1 across the spans.

- geom: catmull_rom holds the curve: samples() (fit_polyline's arithmetic, moved verbatim;
  sketch::fit_polyline now calls it) and beziers(), the exact spans; ArcLength on CurveN
  (Gauss-Legendre per knot span, Newton for the parameter at a length).
- kernel: PathSeg::Cubic, a cubic Bézier in space (the 3D counterpart of PathSeg::Curve, which
  is planar), with its arms wherever PathSeg is matched (start/end, tangents, length, point at
  a length, turn, SMLib curve as a degree-3 NURBS, halving). PointPath / PathPiece carry a path
  through 3D sketch curves: polylines as points, drawn splines as their fit points (reversed
  when the path runs them backwards); PointPath::segs gives SMLib the spans and the lines;
  sweep_point_path_with sweeps it (SMLib through sweep_with, so corners, twist, scale, rails and
  orientation take the existing routes). A smooth polyline piece (a projected or linked curve,
  which the kernel only has as points) is still refused on SMLib.
- doc: path_points builds the PointPath (a wire with three or more fit points is a spline
  piece); solid sweep and pipe call sweep_point_path_with; rails, surface sweeps and curve
  networks take its points as before.

Checked in scratch (not committed): the spans pass through every fit point and every
fit_polyline sample at its parameter (1e-12), C1 at the joints and round a closed spline;
samples() is bit-identical to the old fit_polyline over 1180 random open and closed inputs;
SMLib sweeps of a radius-2 circle along an open spline and round a closed one match
pi r^2 L (L by independent quadrature of the curve's definition) to 1e-7 and 3e-8.
sketch/wire.rs: the existing closed-spline test reads SPAN_SAMPLES from its new home.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_017itGh5fpMEousnZ91cYcVt
This commit is contained in:
jyviko
2026-10-22 00:18:13 -04:00
parent 1bc0df6d03
commit d1eb8b20a0
12 changed files with 524 additions and 106 deletions
Download Patch File Download Diff File Expand all files Collapse all files
crates/doc/src/eval.rs
+58 -48
View File
@@ -5700,12 +5700,13 @@ fn feature_tools_with_constructor(
Some(_) => Vec::new(),
None => path_segments(ps, path, start)?,
};
let path_start = points.as_ref().and_then(|p| p.first().copied()).or_else(|| segs.first().map(kernel::PathSeg::start)).unwrap_or(start);
let path_start =
points.as_ref().and_then(|p| p.points().first().copied()).or_else(|| segs.first().map(kernel::PathSeg::start)).unwrap_or(start);
let rails = rail_polylines(doc, st, rails, f.component, path_start)?;
let ctl = sweep_control(vals, orientation, twist, scale, rails, *rail_fit)?;
// A path through 3D sketch curves follows them as points.
if let Some(pts) = points {
return regions.iter().map(|r| Ok(tool(r, kernel::sweep_path_with(&ss.plane, r, &pts, &ctl)?))).collect();
// A path through 3D sketch curves: its drawn splines exactly, the rest as points.
if let Some(path) = points {
return regions.iter().map(|r| Ok(tool(r, kernel::sweep_point_path_with(&ss.plane, r, &path, &ctl)?))).collect();
}
regions.iter().map(|r| Ok(tool(r, kernel::sweep_with(&ss.plane, r, &segs, &ctl)?))).collect()
}
@@ -5718,11 +5719,17 @@ fn feature_tools_with_constructor(
if !(r > 1e-6) {
return Err(DocError::Invalid("pipe diameter must be positive".into()));
}
if let Some(pts) = path_points(ps, path, Vec3::new(f64::NAN, f64::NAN, f64::NAN))? {
if let Some(along) = path_points(ps, path, Vec3::new(f64::NAN, f64::NAN, f64::NAN))? {
let pts = along.points();
let (Some(a), Some(b)) = (pts.first().copied(), pts.get(1).copied()) else { return Err(DocError::Invalid("empty path".into())) };
let t = (b - a).normalized().ok_or_else(|| DocError::Invalid("path direction".into()))?;
let plane = Plane::from_normal(a, t).ok_or_else(|| DocError::Invalid("pipe profile plane".into()))?;
return Ok(unnamed(vec![kernel::sweep_path(&plane, &pipe_region(vals, r, wall)?, &pts)?]));
return Ok(unnamed(vec![kernel::sweep_point_path_with(
&plane,
&pipe_region(vals, r, wall)?,
&along,
&kernel::SweepControl::default(),
)?]));
}
let segs = path_segments(ps, path, Vec3::new(f64::NAN, f64::NAN, f64::NAN))?;
let first = segs.first().ok_or_else(|| DocError::Invalid("empty path".into()))?;
@@ -5814,7 +5821,7 @@ pub(crate) fn rail_polylines(doc: &Document, st: &ModelState, rails: &[crate::Sw
let mapped = path_in_frame(doc, rs, consumer)?;
let rs = mapped.as_ref();
match path_points(rs, &rail.curves, start)? {
Some(pts) => Ok(pts),
Some(path) => Ok(path.points()),
None => Ok(kernel::path_polyline(&path_segments(rs, &rail.curves, start)?)?),
}
})
@@ -5900,20 +5907,24 @@ fn path_in_frame<'a>(doc: &Document, sketch: &'a SolvedSketch, consumer: u64) ->
Ok(Cow::Owned(mapped))
}
/// A path that uses 3D sketch curves (wires), as one chain of world points from the end nearest
/// `start` (no start: the first curve runs toward the second). `None` when the path is all
/// planar curves.
fn path_points(ps: &SolvedSketch, ids: &[String], start: Vec3) -> Result<Option<Vec<Vec3>>> {
/// A path that uses 3D sketch curves (wires), as one chain of pieces from the end nearest
/// `start` (no start: the first curve runs toward the second): each drawn 3D spline as its fit
/// points (the kernel has its exact curve), every other curve as world points. `None` when the
/// path is all planar curves.
fn path_points(ps: &SolvedSketch, ids: &[String], start: Vec3) -> Result<Option<kernel::PointPath>> {
if !ids.iter().any(|id| ps.sketch.wire_index(id).is_some()) {
return Ok(None);
}
if ids.len() > 1000 {
return Err(DocError::Invalid("the path needs 1…1000 curves".into()));
}
let mut pieces: Vec<Vec<Vec3>> = Vec::with_capacity(ids.len());
let mut pieces: Vec<kernel::PathPiece> = Vec::with_capacity(ids.len());
for id in ids {
if let Some(wi) = ps.sketch.wire_index(id) {
pieces.push(ps.sketch.wires.get(wi).map(|w| w.pts.clone()).unwrap_or_default());
pieces.push(match ps.sketch.wires.get(wi) {
Some(w) if w.fit.len() >= 3 => kernel::PathPiece::Spline { fit: w.fit.clone(), reversed: false },
w => kernel::PathPiece::Polyline(w.map(|w| w.pts.clone()).unwrap_or_default()),
});
} else {
let ci = ps.sketch.curve_index(id).ok_or_else(|| DocError::Unknown(format!("path curve `{id}`")))?;
let mut pts: Vec<Vec3> = Vec::new();
@@ -5925,42 +5936,40 @@ fn path_points(ps: &SolvedSketch, ids: &[String], start: Vec3) -> Result<Option<
}
}
}
pieces.push(pts);
pieces.push(kernel::PathPiece::Polyline(pts));
}
}
// Pieces meet within the sketch's tolerance (projected curves are only that close).
let tol = 1e-3;
let ends = |p: &Vec<Vec3>| (p.first().copied().unwrap_or_default(), p.last().copied().unwrap_or_default());
if let Some(f) = pieces.first().cloned() {
let (a, b) = ends(&f);
let flip = if start.is_finite() {
b.dist(start) < a.dist(start)
} else {
pieces.get(1).is_some_and(|n| {
let (na, nb) = ends(n);
na.dist(a).min(nb.dist(a)) < na.dist(b).min(nb.dist(b))
})
};
if flip && let Some(p0) = pieces.first_mut() {
p0.reverse();
let ends = |p: &kernel::PathPiece| {
let pts = p.points();
(pts.first().copied().unwrap_or_default(), pts.last().copied().unwrap_or_default())
};
let mut pieces = pieces.into_iter();
let Some(first) = pieces.next() else { return Ok(None) };
let rest: Vec<kernel::PathPiece> = pieces.collect();
let (a, b) = ends(&first);
let flip = if start.is_finite() {
b.dist(start) < a.dist(start)
} else {
rest.first().is_some_and(|n| {
let (na, nb) = ends(n);
na.dist(a).min(nb.dist(a)) < na.dist(b).min(nb.dist(b))
})
};
let first = if flip { first.reversed() } else { first };
let mut last = ends(&first).1;
let mut chain = vec![first];
for p in rest {
let (a, b) = ends(&p);
let (p, a, b) = if b.dist(last) < a.dist(last) { (p.reversed(), b, a) } else { (p, a, b) };
if a.dist(last) > tol {
return Err(DocError::Invalid("the path curves are not connected end to end".into()));
}
chain.push(p);
last = b;
}
let mut chain: Vec<Vec3> = Vec::new();
for mut p in pieces {
if let Some(last) = chain.last().copied() {
let (a, b) = ends(&p);
if b.dist(last) < a.dist(last) {
p.reverse();
}
if p.first().is_some_and(|a| a.dist(last) > tol) {
return Err(DocError::Invalid("the path curves are not connected end to end".into()));
}
chain.extend(p.into_iter().skip(1));
} else {
chain = p;
}
}
Ok(Some(chain))
Ok(Some(kernel::PointPath { pieces: chain }))
}
fn path_segments(ps: &SolvedSketch, ids: &[String], start: Vec3) -> Result<Vec<kernel::PathSeg>> {
@@ -5968,8 +5977,8 @@ fn path_segments(ps: &SolvedSketch, ids: &[String], start: Vec3) -> Result<Vec<k
return Err(DocError::Invalid("the path needs 1…1000 curves".into()));
}
// A path through 3D sketch curves is followed as a 3D polyline.
if let Some(pts) = path_points(ps, ids, start)? {
return Ok(pts.windows(2).map(|w| kernel::PathSeg::Line { a: w[0], b: w[1] }).collect());
if let Some(path) = path_points(ps, ids, start)? {
return Ok(path.points().windows(2).map(|w| kernel::PathSeg::Line { a: w[0], b: w[1] }).collect());
}
let mut segs: Vec<solidcraft_geom::Seg2> = Vec::new();
for id in ids {
@@ -6278,7 +6287,7 @@ fn seg_len(s: &kernel::PathSeg) -> f64 {
kernel::PathSeg::Line { a, b } => (b - a).len(),
kernel::PathSeg::Arc { a, center, angle, .. } => (a - center).len() * angle.abs(),
// Patterns flatten their paths; a curve left is measured along its chord.
kernel::PathSeg::Curve { plane, seg } => plane.to_world(seg.start()).dist(plane.to_world(seg.end())),
kernel::PathSeg::Curve { .. } | kernel::PathSeg::Cubic { .. } => s.start().dist(s.end()),
}
}
@@ -6302,8 +6311,8 @@ fn seg_at(s: &kernel::PathSeg, t: f64) -> (Vec3, Vec3) {
let v = rotate_vec(a - center, k, s * t / r);
(center + v, k.cross(v).normalized().unwrap_or(Vec3::X) * s)
}
kernel::PathSeg::Curve { plane, seg } => {
let (a, b) = (plane.to_world(seg.start()), plane.to_world(seg.end()));
kernel::PathSeg::Curve { .. } | kernel::PathSeg::Cubic { .. } => {
let (a, b) = (s.start(), s.end());
let d = (b - a).normalized().unwrap_or(Vec3::X);
(a + d * t, d)
}
@@ -6319,6 +6328,7 @@ fn reverse_seg(s: &kernel::PathSeg) -> kernel::PathSeg {
kernel::PathSeg::Arc { a: end, center, axis: k * -1.0, angle }
}
kernel::PathSeg::Curve { plane, seg } => kernel::PathSeg::Curve { plane, seg: seg.reversed() },
kernel::PathSeg::Cubic { p: [a, b, c, d] } => kernel::PathSeg::Cubic { p: [d, c, b, a] },
}
}
crates/doc/src/eval_surface.rs
+3 -2
View File
@@ -371,7 +371,8 @@ pub(super) fn eval(doc: &Document, vals: &BTreeMap<String, Value>, f: &Feature,
let pts: Vec<Vec2> = chains.iter().flatten().flat_map(|s| s.polyline(1e-2)).collect();
let start = ss.plane.to_world(pts.iter().fold(Vec2::default(), |a, p| a + *p) * (1.0 / pts.len().max(1) as f64));
let body = match path_points(ps, path, start)? {
Some(pts) => {
Some(path) => {
let pts = path.points();
let rails = crate::eval::rail_polylines(doc, st, rails, f.component, pts.first().copied().unwrap_or(start))?;
let ctl = crate::eval::sweep_control(vals, orientation, twist, scale, rails, *rail_fit)?;
kernel::sweep_sheet_points_with(&ss.plane, &chains, &pts, &ctl)?
@@ -608,7 +609,7 @@ fn network_sources(
// 3D sketch curves are followed through their points; planar ones exactly.
let mapped = path_in_frame(doc, raw, f.component)?;
match path_points(&mapped, &c.curves, Vec3::new(f64::NAN, f64::NAN, f64::NAN))? {
Some(pts) => NetSource::Points(pts),
Some(path) => NetSource::Points(path.points()),
None => {
let ss = sketch_in_feature_frame(doc, raw, f.component)?;
let mut chains = curve_chains(&ss, &c.curves)?.into_iter();
crates/geom/src/bspline.rs
+97
View File
@@ -1021,6 +1021,103 @@ impl<const N: usize> CurveN<N> {
}
Some((x, dx))
}
/// The speed in space (the length of the first three coordinates' derivative) at `t`.
fn speed(&self, t: f64) -> Option<f64> {
let (_, d) = self.eval(t)?;
Some(d.iter().take(3).map(|x| x * x).sum::<f64>().sqrt())
}
/// The length in space (the first three coordinates) from `a` to `b`: 8-point Gauss–Legendre.
fn length_between(&self, a: f64, b: f64) -> Option<f64> {
let (mid, half) = (0.5 * (a + b), 0.5 * (b - a));
GAUSS_8.iter().try_fold(0.0, |sum, (x, w)| Some(sum + w * half * self.speed(mid + half * x)?))
}
}
/// Nodes and weights of the 8-point Gauss–Legendre rule on [-1, 1].
const GAUSS_8: [(f64, f64); 8] = [
(-0.960_289_856_497_536_2, 0.101_228_536_290_376_26),
(-0.796_666_477_413_626_7, 0.222_381_034_453_374_47),
(-0.525_532_409_916_329, 0.313_706_645_877_887_3),
(-0.183_434_642_495_649_8, 0.362_683_783_378_362),
(0.183_434_642_495_649_8, 0.362_683_783_378_362),
(0.525_532_409_916_329, 0.313_706_645_877_887_3),
(0.796_666_477_413_626_7, 0.222_381_034_453_374_47),
(0.960_289_856_497_536_2, 0.101_228_536_290_376_26),
];
/// Pieces of each knot span that [`ArcLength`] integrates over.
const ARC_PIECES_PER_SPAN: usize = 8;
/// Length in space along a [`CurveN`] (its first three coordinates), for points evenly spaced
/// by length: each knot span split into [`ARC_PIECES_PER_SPAN`] pieces, each integrated with
/// 8-point Gauss–Legendre (a polynomial curve's speed is smooth within a span, so this is exact
/// to rounding for the cubics it is used on).
#[derive(Clone, Debug, PartialEq)]
pub struct ArcLength {
/// Each piece's start and end parameter and the length before it.
pieces: Vec<[f64; 3]>,
total: f64,
}
impl ArcLength {
/// The table of `c` (`None` for a curve that does not evaluate or has no finite length).
pub fn of<const N: usize>(c: &CurveN<N>) -> Option<ArcLength> {
let (p, n) = (c.degree, c.ctrl.len());
if n <= p || c.knots.len() != n + p + 1 || n > MAX_KNOTS {
return None;
}
let mut pieces = Vec::new();
let mut total = 0.0;
for w in c.knots.get(p..=n)?.windows(2) {
let (a, b) = (w[0], w[1]);
if !(b > a) {
continue;
}
for i in 0..ARC_PIECES_PER_SPAN {
let (s, e) = (a + (b - a) * i as f64 / ARC_PIECES_PER_SPAN as f64, a + (b - a) * (i + 1) as f64 / ARC_PIECES_PER_SPAN as f64);
pieces.push([s, e, total]);
total += c.length_between(s, e)?;
}
}
(total.is_finite() && !pieces.is_empty()).then_some(ArcLength { pieces, total })
}
/// The whole length.
pub fn total(&self) -> f64 {
self.total
}
/// The parameter of `c` (the curve the table is of) at length `s` along it (clamped to the
/// curve): Newton's method within the piece holding it, kept inside that piece by bisection.
pub fn param_at<const N: usize>(&self, c: &CurveN<N>, s: f64) -> Option<f64> {
if !s.is_finite() {
return None;
}
let s = s.clamp(0.0, self.total);
let k = self.pieces.partition_point(|q| q[2] <= s).saturating_sub(1);
let [a, b, before] = *self.pieces.get(k)?;
let tol = 1e-13 * (1.0 + self.total);
let (mut lo, mut hi) = (a, b);
let piece = self.pieces.get(k + 1).map_or(self.total, |q| q[2]) - before;
let mut u = if piece > 0.0 { a + (b - a) * ((s - before) / piece).clamp(0.0, 1.0) } else { a };
for _ in 0..60 {
let f = before + c.length_between(a, u)? - s;
if f.abs() <= tol {
break;
}
if f > 0.0 {
hi = u;
} else {
lo = u;
}
let v = c.speed(u)?;
let next = if v > 0.0 { u - f / v } else { f64::NAN };
u = if next.is_finite() && next > lo && next < hi { next } else { 0.5 * (lo + hi) };
}
Some(u)
}
}
#[cfg(test)]
crates/geom/src/catmull_rom.rs
+113
View File
@@ -0,0 +1,113 @@
//! The centripetal Catmull–Rom spline through fit points: the smooth curve a 3D sketch draws
//! through the points picked for it (no overshoot loops, no cusps).
//!
//! Between fit points `P1` and `P2` the curve is the cubic of the Barry–Goldman pyramid over the
//! four points `P0…P3` around them, at the knots `t0 = 0`, `t(i+1) = t(i) + |P(i+1) - P(i)|^½`,
//! run from `t1` to `t2` (Barry and Goldman, "A recursive evaluation algorithm for a class of
//! Catmull–Rom splines", SIGGRAPH 1988; the centripetal exponent ½ after Yuksel, Schaefer and
//! Keyser, "Parameterization and applications of Catmull–Rom curves", CAD 43(7), 2011). Its
//! derivative at `P1`,
//!
//! `m1 = (P1 - P0)/(t1 - t0) - (P2 - P0)/(t2 - t0) + (P2 - P1)/(t2 - t1)`,
//!
//! depends only on `P1`, its two neighbours and the two knot intervals between them, so the spans
//! on either side of a fit point leave it with the same derivative: in the parameter that runs on
//! from span to span by each span's knot interval, the spline is C1, exactly one cubic Bézier per
//! span ([`beziers`]). [`samples`] is the polyline through [`SPAN_SAMPLES`] points per span of the
//! same curve.
//!
//! Ends: an open spline mirrors its end points (`P(-1) = 2 P(0) - P(1)`, and so at the far end),
//! which leaves it leaving each end straight towards the next fit point. A spline of more than
//! three points whose last point is its first (within 1e-9) is closed: its neighbours wrap round,
//! so it is as smooth where it closes as anywhere else.
use crate::Vec3;
/// Samples per span of [`samples`].
pub const SPAN_SAMPLES: usize = 24;
/// One span: the four points around it and their knots (`p[1]` to `p[2]` is the span).
struct Span {
p: [Vec3; 4],
t: [f64; 4],
}
/// Whether the spline through `fit` closes on itself (see the module notes).
pub fn is_closed(fit: &[Vec3]) -> bool {
let n = fit.len();
n > 3 && fit.first().zip(fit.last()).is_some_and(|(a, b)| a.dist(*b) < 1e-9)
}
/// The spans of the spline through `fit` (three or more points), in order.
fn spans(fit: &[Vec3]) -> Vec<Span> {
let n = fit.len();
let (Some(&first), Some(&second), Some(&last), Some(&before_last)) = (fit.first(), fit.get(1), fit.last(), fit.get(n.wrapping_sub(2))) else {
return Vec::new();
};
let closed = is_closed(fit);
let at = |i: isize| -> Vec3 {
match usize::try_from(i) {
Ok(i) if i < n => fit.get(i).copied().unwrap_or(last),
Ok(_) if closed => second,
Err(_) if closed => before_last,
// Mirrored end tangents.
Ok(_) => last * 2.0 - before_last,
Err(_) => first * 2.0 - second,
}
};
let tj = |a: Vec3, b: Vec3| a.dist(b).sqrt().max(1e-9);
(0..n.saturating_sub(1))
.map(|k| {
let k = k as isize;
let (p0, p1, p2, p3) = (at(k - 1), at(k), at(k + 1), at(k + 2));
let t0 = 0.0;
let t1 = t0 + tj(p0, p1);
let t2 = t1 + tj(p1, p2);
let t3 = t2 + tj(p2, p3);
Span { p: [p0, p1, p2, p3], t: [t0, t1, t2, t3] }
})
.collect()
}
/// The polyline of the spline through `fit`: the point itself (one point), the line (two), or
/// [`SPAN_SAMPLES`] points per span after the first point, the last at the last fit point.
pub fn samples(fit: &[Vec3]) -> Vec<Vec3> {
if fit.len() <= 2 {
return fit.to_vec();
}
let mut out = fit.first().copied().into_iter().collect::<Vec<_>>();
for Span { p: [p0, p1, p2, p3], t: [t0, t1, t2, t3] } in spans(fit) {
for s in 1..=SPAN_SAMPLES {
let t = t1 + (t2 - t1) * s as f64 / SPAN_SAMPLES as f64;
let a1 = p0 * ((t1 - t) / (t1 - t0)) + p1 * ((t - t0) / (t1 - t0));
let a2 = p1 * ((t2 - t) / (t2 - t1)) + p2 * ((t - t1) / (t2 - t1));
let a3 = p2 * ((t3 - t) / (t3 - t2)) + p3 * ((t - t2) / (t3 - t2));
let b1 = a1 * ((t2 - t) / (t2 - t0)) + a2 * ((t - t0) / (t2 - t0));
let b2 = a2 * ((t3 - t) / (t3 - t1)) + a3 * ((t - t1) / (t3 - t1));
out.push(b1 * ((t2 - t) / (t2 - t1)) + b2 * ((t - t1) / (t2 - t1)));
}
}
out
}
/// The spline through `fit` (three or more points) exactly, as one cubic Bézier per span (from
/// fit point to fit point): `P1`, `P1 + m1 h/3`, `P2 - m2 h/3`, `P2` (`h = t2 - t1`, `m` the
/// derivatives in the module notes), its parameter `(t - t1) / h`. Neighbouring spans meet with
/// the same derivative in `t`, so the spans make a C1 cubic B-spline with double interior knots.
/// Empty for fewer than three points (a point or a line); `None` for a spline that does not come
/// out finite.
pub fn beziers(fit: &[Vec3]) -> Option<Vec<[Vec3; 4]>> {
if fit.len() < 3 {
return Some(Vec::new());
}
let out: Vec<[Vec3; 4]> = spans(fit)
.into_iter()
.map(|Span { p: [p0, p1, p2, p3], t: [t0, t1, t2, t3] }| {
let h = t2 - t1;
let m1 = (p1 - p0) * (1.0 / (t1 - t0)) - (p2 - p0) * (1.0 / (t2 - t0)) + (p2 - p1) * (1.0 / h);
let m2 = (p2 - p1) * (1.0 / h) - (p3 - p1) * (1.0 / (t3 - t1)) + (p3 - p2) * (1.0 / (t3 - t2));
[p1, p1 + m1 * (h / 3.0), p2 - m2 * (h / 3.0), p2]
})
.collect();
out.iter().flatten().all(|p| p.is_finite()).then_some(out)
}
crates/geom/src/lib.rs
+1
View File
@@ -15,6 +15,7 @@ pub mod blend;
pub mod brep;
pub mod bspline;
mod bvh;
pub mod catmull_rom;
mod chain;
pub mod compare;
pub mod env;
crates/kernel/src/build.rs
+27 -1
View File
@@ -2337,6 +2337,11 @@ pub enum PathSeg {
plane: Plane,
seg: Seg2,
},
/// A cubic Bézier in space with control points `p`: a span of a drawn 3D sketch spline
/// ([`crate::PathPiece::Spline`]).
Cubic {
p: [Vec3; 4],
},
}
impl PathSeg {
@@ -2345,6 +2350,17 @@ impl PathSeg {
match *self {
PathSeg::Line { a, .. } | PathSeg::Arc { a, .. } => a,
PathSeg::Curve { plane, seg } => plane.to_world(seg.start()),
PathSeg::Cubic { p } => p[0],
}
}
/// Where the segment ends.
pub fn end(&self) -> Vec3 {
match *self {
PathSeg::Line { b, .. } => b,
PathSeg::Arc { a, center, axis, angle } => axis.normalized().map_or(a, |k| rotate_about(a, center, k, angle)),
PathSeg::Curve { plane, seg } => plane.to_world(seg.end()),
PathSeg::Cubic { p } => p[3],
}
}
@@ -2363,6 +2379,9 @@ impl PathSeg {
PathSeg::Curve { plane, seg } => {
Some((plane.dir_to_world(seg.start_tangent()).normalized()?, plane.dir_to_world(seg.end_tangent()).normalized()?))
}
PathSeg::Cubic { p } => {
Some((crate::point_path::cubic_derivative(&p, 0.0).normalized()?, crate::point_path::cubic_derivative(&p, 1.0).normalized()?))
}
}
}
@@ -2374,6 +2393,12 @@ impl PathSeg {
PathSeg::Curve { plane, seg } => {
seg.polyline(1e-2).windows(2).map(|q| PathSeg::Line { a: plane.to_world(q[0]), b: plane.to_world(q[1]) }).collect()
}
PathSeg::Cubic { p } => std::iter::once(p[0])
.chain(crate::point_path::cubic_polyline(&p, 1e-2))
.collect::<Vec<_>>()
.windows(2)
.map(|q| PathSeg::Line { a: q[0], b: q[1] })
.collect(),
_ => vec![*s],
})
.collect()
@@ -2384,6 +2409,7 @@ impl PathSeg {
match *self {
PathSeg::Line { b, .. } => vec![b],
PathSeg::Curve { plane, seg } => seg.polyline(tol).into_iter().skip(1).map(|q| plane.to_world(q)).collect(),
PathSeg::Cubic { p } => crate::point_path::cubic_polyline(&p, tol),
PathSeg::Arc { a, center, axis, angle } => {
let Some(k) = axis.normalized() else { return vec![a] };
let r = (a - center).len();
@@ -2521,7 +2547,7 @@ pub fn sweep(plane: &Plane, region: &Region2, path: &[PathSeg]) -> Result<Body>
builder::rsweep(&cur, p3(center), v3(ax), mt::Rad(angle))
}
// The path was flattened: no curves are left.
PathSeg::Curve { .. } => return Err(KernelError::Invalid("a free-form path curve".into())),
PathSeg::Curve { .. } | PathSeg::Cubic { .. } => return Err(KernelError::Invalid("a free-form path curve".into())),
};
let next = part
.extract_boundaries()
crates/kernel/src/lib.rs
+2
View File
@@ -111,6 +111,7 @@ mod planar_faces;
#[cfg(test)]
mod planar_faces_tests;
mod plane_fit;
mod point_path;
mod polybool;
mod polyhedron;
mod prism;
@@ -207,6 +208,7 @@ pub use ops::{
BoolOp, boolean, can_join_all_touching_faces, can_join_touching_faces, compound_disjoint_extrusions, half_space, is_planar, join_apart,
place_matrix, split_by_plane, transform, transform_matrix, union_apart, with_attempt_budget,
};
pub use point_path::{PathPiece, PointPath, sweep_point_path_with};
pub use polybool::{faceted_boolean, planar_boolean};
pub use polyhedron::{HalfSpace, convex_polyhedron, draft, offset_faces, shell};
pub use prism::{PLANAR_FIT_TOL, extrude_planar_regions, planar_regions_boolean};
crates/kernel/src/point_path.rs
+176
View File
@@ -0,0 +1,176 @@
//! Sweep paths through 3D sketch curves ([`PointPath`]): polylines (linked, projected and
//! included curves, which reach the kernel as their points) and drawn 3D splines, which it gets
//! as their fit points and so exactly ([`solidcraft_geom::catmull_rom`]). Truck sweeps the
//! path's polyline ([`PointPath::points`]: the splines' samples); SMLib follows the path's
//! segments ([`PointPath::segs`]), each spline span exactly as a [`PathSeg::Cubic`].
use solidcraft_geom::bspline::{ArcLength, CurveN};
use solidcraft_geom::catmull_rom;
use solidcraft_geom::{Plane, Region2, Vec3};
use crate::body::Body;
use crate::{KernelError, PathSeg, Result, SweepControl};
/// One curve of a [`PointPath`], in the path's direction of travel.
#[derive(Clone, Debug, PartialEq)]
pub enum PathPiece {
/// Points joined by straight lines.
Polyline(Vec<Vec3>),
/// A drawn 3D spline: the centripetal Catmull–Rom spline through `fit`, from its last fit
/// point back to its first when `reversed`.
Spline { fit: Vec<Vec3>, reversed: bool },
}
impl PathPiece {
/// The piece's points in its direction of travel: the polyline, or the spline's samples
/// ([`catmull_rom::samples`]).
pub fn points(&self) -> Vec<Vec3> {
match self {
PathPiece::Polyline(pts) => pts.clone(),
PathPiece::Spline { fit, reversed } => {
let mut pts = catmull_rom::samples(fit);
if *reversed {
pts.reverse();
}
pts
}
}
}
/// The piece the other way round.
pub fn reversed(self) -> PathPiece {
match self {
PathPiece::Polyline(mut pts) => {
pts.reverse();
PathPiece::Polyline(pts)
}
PathPiece::Spline { fit, reversed } => PathPiece::Spline { fit, reversed: !reversed },
}
}
/// A spline of three or more fit points as its exact spans ([`PathSeg::Cubic`]), in the
/// direction of travel; `None` for a polyline (and a spline of one or two points, a point or
/// a line).
pub(crate) fn spans(&self) -> Result<Option<Vec<PathSeg>>> {
let PathPiece::Spline { fit, reversed } = self else { return Ok(None) };
let spans = catmull_rom::beziers(fit).ok_or_else(|| KernelError::Invalid("the 3D sketch spline is not finite".into()))?;
if spans.is_empty() {
return Ok(None);
}
let seg = |p: [Vec3; 4]| PathSeg::Cubic { p };
Ok(Some(if *reversed {
spans.into_iter().rev().map(|[a, b, c, d]| seg([d, c, b, a])).collect()
} else {
spans.into_iter().map(seg).collect()
}))
}
}
/// A sweep path through 3D sketch curves: its pieces end to end, each starting where the one
/// before it ended.
#[derive(Clone, Debug, Default, PartialEq)]
pub struct PointPath {
pub pieces: Vec<PathPiece>,
}
impl PointPath {
/// A path of one polyline.
pub fn polyline(pts: Vec<Vec3>) -> PointPath {
PointPath { pieces: vec![PathPiece::Polyline(pts)] }
}
/// The path as one polyline: the first piece's points, then each next piece's after its
/// first (where the one before ended).
pub fn points(&self) -> Vec<Vec3> {
let mut chain: Vec<Vec3> = Vec::new();
for p in &self.pieces {
let pts = p.points();
if chain.is_empty() {
chain = pts;
} else {
chain.extend(pts.into_iter().skip(1));
}
}
chain
}
/// The path as segments: each spline span exactly ([`PathSeg::Cubic`]), each polyline as the
/// lines between its points, starting where the piece before it ended (as [`Self::points`]
/// chains them).
pub(crate) fn segs(&self) -> Result<Vec<PathSeg>> {
let mut out: Vec<PathSeg> = Vec::new();
for piece in &self.pieces {
if let Some(spans) = piece.spans()? {
out.extend(spans);
continue;
}
let mut pts = piece.points();
if let (Some(end), Some(first)) = (out.last().map(PathSeg::end), pts.first_mut()) {
*first = end;
}
out.extend(pts.windows(2).filter(|w| w[0].dist(w[1]) > 1e-9).map(|w| PathSeg::Line { a: w[0], b: w[1] }));
}
Ok(out)
}
}
/// [`crate::sweep_path_with`] along a path through 3D sketch curves. Truck sweeps the path's
/// polyline ([`PointPath::points`]); SMLib follows its splines exactly and its polylines as lines.
pub fn sweep_point_path_with(plane: &Plane, region: &Region2, path: &PointPath, ctl: &SweepControl) -> Result<Body> {
ctl.check()?;
if crate::backend()? == crate::Backend::Smlib {
return crate::sweep_path::smlib_sweep_points(plane, region, path, ctl);
}
crate::sweep_path_with(plane, region, &path.points(), ctl)
}
/// A cubic Bézier in space (a [`PathSeg::Cubic`]): its point at `x` (0…1).
pub(crate) fn cubic_at(p: &[Vec3; 4], x: f64) -> Vec3 {
let y = 1.0 - x;
p[0] * (y * y * y) + p[1] * (3.0 * x * y * y) + p[2] * (3.0 * x * x * y) + p[3] * (x * x * x)
}
/// Its derivative at `x`.
pub(crate) fn cubic_derivative(p: &[Vec3; 4], x: f64) -> Vec3 {
let y = 1.0 - x;
(p[1] - p[0]) * (3.0 * y * y) + (p[2] - p[1]) * (6.0 * x * y) + (p[3] - p[2]) * (3.0 * x * x)
}
/// The cubic split at `x` (de Casteljau): the part before and the part after.
#[cfg_attr(not(feature = "smlib"), allow(dead_code))]
pub(crate) fn cubic_split(p: &[Vec3; 4], x: f64) -> ([Vec3; 4], [Vec3; 4]) {
let [a, b, c, d] = *p;
let (ab, bc, cd) = (a.lerp(b, x), b.lerp(c, x), c.lerp(d, x));
let (abc, bcd) = (ab.lerp(bc, x), bc.lerp(cd, x));
let m = abc.lerp(bcd, x);
([a, ab, abc, m], [m, bcd, cd, d])
}
#[cfg_attr(not(feature = "smlib"), allow(dead_code))]
fn cubic_curve(p: &[Vec3; 4]) -> CurveN<3> {
CurveN { degree: 3, knots: vec![0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0], ctrl: p.iter().map(|q| [q.x, q.y, q.z]).collect() }
}
/// Its length (`None` for a cubic that does not evaluate).
#[cfg_attr(not(feature = "smlib"), allow(dead_code))]
pub(crate) fn cubic_length(p: &[Vec3; 4]) -> Option<f64> {
ArcLength::of(&cubic_curve(p)).map(|a| a.total())
}
/// The point `f` (0…1) of the way along the cubic by length, and the unit direction there.
#[cfg_attr(not(feature = "smlib"), allow(dead_code))]
pub(crate) fn cubic_along(p: &[Vec3; 4], f: f64) -> Option<(Vec3, Vec3)> {
let c = cubic_curve(p);
let arc = ArcLength::of(&c)?;
let x = arc.param_at(&c, f * arc.total())?;
Some((cubic_at(p, x), cubic_derivative(p, x).normalized()?))
}
/// Points along the cubic, without its start, the chords within `tol` of it (by the bound on
/// its second derivative, 6 max |P(i) - 2 P(i+1) + P(i+2)|: a chord of a piece `h` long in the
/// parameter is within h^2/8 of that of the curve).
pub(crate) fn cubic_polyline(p: &[Vec3; 4], tol: f64) -> Vec<Vec3> {
let bend = (p[0] - p[1] * 2.0 + p[2]).len().max((p[1] - p[2] * 2.0 + p[3]).len()) * 6.0;
let n = ((bend / (8.0 * tol.max(1e-9))).sqrt().ceil() as usize).clamp(1, 512);
(1..=n).map(|i| cubic_at(p, i as f64 / n as f64)).collect()
}
crates/kernel/src/smlib_sweeps.rs
+20
View File
@@ -458,6 +458,11 @@ mod imp {
Curve::Arc { center: p(center), axis: p(k), reference: p(reference), radius, start: 0.0, end: turn }
}
crate::PathSeg::Curve { plane, seg } => curve(&plane, &seg)?,
crate::PathSeg::Cubic { p: ctrl } => Curve::Nurbs {
degree: 3,
points: ctrl.iter().map(|q| [q.x, q.y, q.z, 1.0]).collect(),
knots: vec![0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0],
},
})
}
@@ -534,6 +539,9 @@ mod imp {
let n = plane.normal();
if translate { moved } else { turned_plane(&moved, b, n, t0.cross(t1).dot(n).atan2(t0.dot(t1))) }
}
// How SMLib's sweep turns the profile along a curve in space is its own frame's.
crate::PathSeg::Cubic { p: ctrl } if translate => moved_plane(&pl, ctrl[3] - ctrl[0]),
crate::PathSeg::Cubic { .. } => return Err(refuse("sweep round a corner next to a 3D sketch spline, the profile turning with it")),
};
}
Ok(pl)
@@ -548,6 +556,7 @@ mod imp {
(a, crate::build::rotate_about(a, center, k, angle))
}
crate::PathSeg::Curve { plane, seg } => (plane.to_world(seg.start()), plane.to_world(seg.end())),
crate::PathSeg::Cubic { p: ctrl } => (ctrl[0], ctrl[3]),
})
}
@@ -591,6 +600,7 @@ mod imp {
length += a.dist(b).max(match *s {
crate::PathSeg::Arc { a, center, angle, .. } => (a - center).len() * angle.abs(),
crate::PathSeg::Curve { seg, .. } => seg.length(),
crate::PathSeg::Cubic { .. } => seg_len(s),
crate::PathSeg::Line { .. } => 0.0,
});
}
@@ -639,6 +649,7 @@ mod imp {
crate::PathSeg::Line { a, b } => crate::PathSeg::Line { a, b: a + (b - a) * 0.5 },
crate::PathSeg::Arc { a, center, axis, angle } => crate::PathSeg::Arc { a, center, axis, angle: angle * 0.5 },
crate::PathSeg::Curve { plane, seg } => crate::PathSeg::Curve { plane, seg: seg.split_half().0 },
crate::PathSeg::Cubic { p: ctrl } => crate::PathSeg::Cubic { p: crate::point_path::cubic_split(&ctrl, 0.5).0 },
};
let there = carried(&at_mid, &[half], translate)?;
inside.extend(seeds(&there, pieces, |q| q));
@@ -700,6 +711,7 @@ mod imp {
crate::PathSeg::Line { a, b } => a.dist(b),
crate::PathSeg::Arc { a, center, angle, .. } => angle.abs() * a.dist(center),
crate::PathSeg::Curve { seg, .. } => seg.length(),
crate::PathSeg::Cubic { p: ctrl } => crate::point_path::cubic_length(&ctrl).unwrap_or(f64::NAN),
}
}
@@ -727,6 +739,7 @@ mod imp {
let t = 0.5 * (lo + hi);
Some((plane.to_world(seg.point_at(t)), plane.dir_to_world(seg.derivative(t)).normalized()?))
}
crate::PathSeg::Cubic { p: ctrl } => crate::point_path::cubic_along(&ctrl, f),
}
}
@@ -768,6 +781,13 @@ mod imp {
.map(|(u, v)| u.cross(v).atan2(u.dot(v)).abs())
.sum()
}
crate::PathSeg::Cubic { p: ctrl } => {
let pts: Vec<Vec3> = std::iter::once(ctrl[0]).chain(crate::point_path::cubic_polyline(&ctrl, 1e-3)).collect();
pts.windows(3)
.filter_map(|w| Some(((w[1] - w[0]).normalized()?, (w[2] - w[1]).normalized()?)))
.map(|(u, v)| u.cross(v).len().atan2(u.dot(v)))
.sum()
}
}
}
crates/kernel/src/surfaces.rs
+1 -1
View File
@@ -237,7 +237,7 @@ pub fn sweep_sheet(plane: &Plane, chains: &[Vec<Seg2>], path: &[crate::PathSeg])
(part, Box::new(move |p| crate::build::rotate_about(p, center, k, angle)))
}
// The path was flattened: no curves are left.
PathSeg::Curve { .. } => return Err(KernelError::Invalid("a free-form path curve".into())),
PathSeg::Curve { .. } | PathSeg::Cubic { .. } => return Err(KernelError::Invalid("a free-form path curve".into())),
};
let next = far_wire(&part, &cur, map.as_ref())?;
sides.extend(swept_sides(&part, &seeds, &[]));
crates/kernel/src/sweep_path.rs
+16 -12
View File
@@ -530,7 +530,7 @@ pub fn sweep_path(plane: &Plane, region: &Region2, path: &[Vec3]) -> Result<Body
pub fn sweep_path_with(plane: &Plane, region: &Region2, path: &[Vec3], ctl: &SweepControl) -> Result<Body> {
ctl.check()?;
if crate::backend()? == crate::Backend::Smlib {
return smlib_sweep_points(plane, region, path, ctl);
return smlib_sweep_points(plane, region, &crate::PointPath::polyline(path.to_vec()), ctl);
}
if region.area().abs() < 1e-9 {
return Err(KernelError::Invalid("profile has no area".into()));
@@ -566,18 +566,22 @@ pub fn sweep_path_with(plane: &Plane, region: &Region2, path: &[Vec3], ctl: &Swe
ruled(plane, &profile, &pts, ctl)
}
/// [`sweep_path_with`] through SMLib: a path of points is the polyline through them, its
/// straight runs the lines between them. A line, or a polyline with corners, is swept as those
/// lines ([`crate::sweep_with`], which mitres the corners); a smooth polyline of more points stands for a curve (a drawn 3D spline) that SMLib has no
/// path for yet.
fn smlib_sweep_points(plane: &Plane, region: &Region2, path: &[Vec3], ctl: &SweepControl) -> Result<Body> {
let pts = clean(path)?;
let smooth = pts.len() > 2 && !pts.windows(3).any(|w| turn(w[1] - w[0], w[2] - w[1]) > CORNER);
if smooth {
return Err(crate::smlib_sweeps::refuse("sweep along a smooth path of points (a 3D sketch spline)"));
/// [`crate::sweep_point_path_with`] through SMLib: each drawn 3D spline exactly (its spans), each
/// polyline as the lines between its points ([`crate::PointPath::segs`]), swept by
/// [`crate::sweep_with`], which mitres the corners. A smooth polyline of more points stands for a
/// curve that reaches the kernel only as points (a projected, linked or included 3D curve), which
/// SMLib has no path for yet.
pub(crate) fn smlib_sweep_points(plane: &Plane, region: &Region2, path: &crate::PointPath, ctl: &SweepControl) -> Result<Body> {
for piece in &path.pieces {
if let crate::PathPiece::Polyline(pts) = piece {
let pts = clean(pts)?;
if pts.len() > 2 && !pts.windows(3).any(|w| turn(w[1] - w[0], w[2] - w[1]) > CORNER) {
return Err(crate::smlib_sweeps::refuse("sweep along a smooth path of points (a projected or linked 3D curve)"));
}
}
}
let segs: Vec<crate::PathSeg> = pts.windows(2).map(|w| crate::PathSeg::Line { a: w[0], b: w[1] }).collect();
if segs.len() > 1000 {
let segs = path.segs()?;
if segs.is_empty() || segs.len() > 1000 {
return Err(KernelError::Invalid("the path needs 1…1000 segments".into()));
}
crate::sweep_with(plane, region, &segs, ctl)
crates/sketch/src/wire.rs
+10 -42
View File
@@ -67,48 +67,11 @@ pub struct Wire {
/// Most fit points of a drawn 3D curve.
pub const MAX_FIT_POINTS: usize = 500;
/// Samples per span of a drawn 3D spline.
const SPAN_SAMPLES: usize = 24;
/// The polyline of a drawn 3D curve through `fit`: the point itself, the line, or a centripetal
/// Catmull–Rom spline through the points (smooth, no overshoot loops).
/// The polyline of a drawn 3D curve through `fit`: the point itself, the line, or the centripetal
/// Catmull–Rom spline through the points (smooth, no overshoot loops), sampled
/// ([`solidcraft_geom::catmull_rom`] has the curve itself, exactly).
pub fn fit_polyline(fit: &[Vec3]) -> Vec<Vec3> {
if fit.len() <= 2 {
return fit.to_vec();
}
let n = fit.len();
// A curve that ends where it starts is smooth through that point (its neighbours wrap).
let closed = n > 3 && fit[0].dist(fit[n - 1]) < 1e-9;
let at = |i: isize| -> Vec3 {
match usize::try_from(i) {
Ok(i) if i < n => fit[i],
Ok(_) if closed => fit[1],
Err(_) if closed => fit[n - 2],
// Mirrored end tangents.
Ok(_) => fit[n - 1] * 2.0 - fit[n - 2],
Err(_) => fit[0] * 2.0 - fit[1],
}
};
let mut out = vec![fit[0]];
for k in 0..n - 1 {
let k = k as isize;
let (p0, p1, p2, p3) = (at(k - 1), at(k), at(k + 1), at(k + 2));
let tj = |a: Vec3, b: Vec3| a.dist(b).sqrt().max(1e-9);
let t0 = 0.0;
let t1 = t0 + tj(p0, p1);
let t2 = t1 + tj(p1, p2);
let t3 = t2 + tj(p2, p3);
for s in 1..=SPAN_SAMPLES {
let t = t1 + (t2 - t1) * s as f64 / SPAN_SAMPLES as f64;
let a1 = p0 * ((t1 - t) / (t1 - t0)) + p1 * ((t - t0) / (t1 - t0));
let a2 = p1 * ((t2 - t) / (t2 - t1)) + p2 * ((t - t1) / (t2 - t1));
let a3 = p2 * ((t3 - t) / (t3 - t2)) + p3 * ((t - t2) / (t3 - t2));
let b1 = a1 * ((t2 - t) / (t2 - t0)) + a2 * ((t - t0) / (t2 - t0));
let b2 = a2 * ((t3 - t) / (t3 - t1)) + a3 * ((t - t1) / (t3 - t1));
out.push(b1 * ((t2 - t) / (t2 - t1)) + b2 * ((t - t1) / (t2 - t1)));
}
}
out
solidcraft_geom::catmull_rom::samples(fit)
}
/// Chord-height tolerance (mm) when sampling a drawn 3D arc into `pts`.
@@ -463,7 +426,12 @@ mod tests {
let n = pts.len();
let (start, end) = (dir(pts[0], pts[1]), dir(pts[n - 2], pts[n - 1]));
let mid = |i: usize| dir(pts[i - 1], pts[i]).dot(dir(pts[i], pts[i + 1]));
assert!(start.dot(end) >= mid(SPAN_SAMPLES) - 1e-9, "turn at the close {}, at a fit point {}", start.dot(end), mid(SPAN_SAMPLES));
assert!(
start.dot(end) >= mid(solidcraft_geom::catmull_rom::SPAN_SAMPLES) - 1e-9,
"turn at the close {}, at a fit point {}",
start.dot(end),
mid(solidcraft_geom::catmull_rom::SPAN_SAMPLES)
);
let open = fit_polyline(&ring[..12]);
assert_eq!(open.first(), ring.first());
assert_eq!(open.last(), ring.get(11));