Exact free-form profiles: geom Seg2 gains cubic Béziers and conics (rational quadratics: exact length, area, split, tangents); sketch ellipses, conics, text, fit splines and control splines up to degree 3 enter profiles as exact segments wherever a loop uses a whole segment, and the kernel builds them as B-spline / NURBS edges (an extruded ellipse is three faces with the exact volume)

Co-Authored-By: Claude Opus 5.5 <[email protected]>
Claude-Session: https://claude.ai/code/session_01AS5Tn3xieewUSg3uxrMVdL
This commit is contained in:
bherbruckandClaude Opus 5.5 committed 2026-10-08 02:16:00 +00:00
1 parent 8b4ec13d3b
commit 31030ed008
13 files changed
+447 -31

No files matched your search

+7 -3
View File
@@ -1106,10 +1106,14 @@ fn path_segments(ps: &SolvedSketch, ids: &[String], start: Vec3) -> Result<Vec<k
let n = ps.plane.normal();
Ok(segs
.iter()
.map(|s| match *s {
solvecraft_geom::Seg2::Line { a, b } => kernel::PathSeg::Line { a: w(a), b: w(b) },
.flat_map(|s| match *s {
solvecraft_geom::Seg2::Line { a, b } => vec![kernel::PathSeg::Line { a: w(a), b: w(b) }],
solvecraft_geom::Seg2::Arc { center, sweep, .. } => {
kernel::PathSeg::Arc { a: w(s.start()), center: w(center), axis: n * sweep.signum(), angle: sweep.abs() }
vec![kernel::PathSeg::Arc { a: w(s.start()), center: w(center), axis: n * sweep.signum(), angle: sweep.abs() }]
}
// Free-form path curves are followed as polylines.
solvecraft_geom::Seg2::Cubic { .. } | solvecraft_geom::Seg2::Conic { .. } => {
s.polyline(1e-2).windows(2).map(|q| kernel::PathSeg::Line { a: w(q[0]), b: w(q[1]) }).collect()
}
})
.collect())
+1 -2
View File
@@ -424,7 +424,7 @@ fn tilted_circle(plane: &Plane, pts: &[Vec3]) -> Option<Vec<LinkGeom>> {
return Some(vec![LinkGeom::Ellipse { c: plane.to_local(c), major: plane.to_local(c + u * r), minor: r * cosang }]);
}
// Sweep from the first point along the polyline's direction.
let (p0, p1) = (*pts.first()?, *pts.last()?);
let p0 = *pts.first()?;
let e1 = (p0 - c).normalized()?;
let e2 = n.cross(e1);
let ang = |p: Vec3| (p - c).dot(e2).atan2((p - c).dot(e1));
@@ -439,7 +439,6 @@ fn tilted_circle(plane: &Plane, pts: &[Vec3]) -> Option<Vec<LinkGeom>> {
}
sweep += d;
}
let _ = p1;
let k = ((sweep.abs() / std::f64::consts::FRAC_PI_2).ceil() as usize).clamp(1, 8);
let h = sweep / k as f64;
let at = |t: f64| c + (e1 * t.cos() + e2 * t.sin()) * r;
@@ -314,3 +314,31 @@ fn centerline_is_left_out_of_profiles() {
assert!(s.execute("sketch.centerline", &json!({"curves": ["nope"]})).is_err());
let _ = l;
}
#[test]
fn free_form_profiles_extrude_exactly() {
let pi = std::f64::consts::PI;
let vol = |s: &mut Session| run(s, "MeasureCommand", json!({}))["total"]["volume_mm3"].as_f64().unwrap();
// Ellipse a = 10, b = 4: π a b h.
let mut s = new_sketch();
run(&mut s, "CircleElipse", json!({"center": [0, 0], "major": [10, 0], "minor_radius": 4}));
let a = profiles(&s)[0];
assert!((a - pi * 40.0).abs() < 1e-9, "{a}");
run(&mut s, "SketchStop", json!({}));
run(&mut s, "Extrude", json!({"distance": 5}));
let v = vol(&mut s);
assert!((v - pi * 40.0 * 5.0).abs() / (pi * 200.0) < 2e-4, "{v}");
let m = run(&mut s, "MeasureCommand", json!({}));
assert_eq!(m["bodies"][0]["faces"], 3, "an exact ellipse: top, bottom and one side ({m})");
// A region closed by a fit spline and a line.
let mut s = new_sketch();
let sp = ids(&run(&mut s, "DrawSpline", json!({"points": [[0, 0], [10, 8], [20, 0]]}))["curves"])[0].clone();
run(&mut s, "DrawPolyline", json!({"points": [format!("{sp}.end"), format!("{sp}.start")]}));
let a = profiles(&s)[0];
run(&mut s, "SketchStop", json!({}));
run(&mut s, "Extrude", json!({"distance": 2}));
let v = vol(&mut s);
assert!((v - a * 2.0).abs() / (a * 2.0) < 2e-4, "{v} vs {}", a * 2.0);
let m = run(&mut s, "MeasureCommand", json!({}));
assert!(m["bodies"][0]["faces"].as_u64().unwrap() <= 5, "{m}");
}
+114 -4
View File
@@ -1,5 +1,6 @@
//! Planar profiles: closed loops of lines and circular arcs, and regions (an outer loop minus
//! holes). Sketches produce them; the kernel turns them into faces.
//! Planar profiles: closed loops of lines, circular arcs, cubic Béziers and conics (rational
//! quadratic Béziers), and regions (an outer loop minus holes). Sketches produce them; the
//! kernel turns them into faces.
use std::f64::consts::TAU;
@@ -22,18 +23,61 @@ pub enum Seg2 {
start: f64,
sweep: f64,
},
/// Cubic Bézier.
Cubic {
p0: Vec2,
p1: Vec2,
p2: Vec2,
p3: Vec2,
},
/// Rational quadratic Bézier from `a` to `b` with control `apex` of weight `w` (end
/// weights 1): exact conic arcs (ellipses, parabolas, hyperbolas).
Conic {
a: Vec2,
apex: Vec2,
b: Vec2,
w: f64,
},
}
/// Gauss–Legendre nodes and weights on [0, 1] (5 points).
const GL5: [(f64, f64); 5] = [
(0.046_910_077_030_668, 0.118_463_442_528_095),
(0.230_765_344_947_158, 0.239_314_335_249_683),
(0.5, 0.284_444_444_444_444),
(0.769_234_655_052_842, 0.239_314_335_249_683),
(0.953_089_922_969_332, 0.118_463_442_528_095),
];
impl Seg2 {
/// Derivative with respect to the parameter at `t` (Béziers and conics; lines and arcs too).
pub fn derivative(&self, t: f64) -> Vec2 {
match *self {
Seg2::Line { a, b } => b - a,
Seg2::Arc { radius, start, sweep, .. } => Vec2::from_angle(start + sweep * t).perp() * (radius * sweep),
Seg2::Cubic { p0, p1, p2, p3 } => {
let u = 1.0 - t;
(p1 - p0) * (3.0 * u * u) + (p2 - p1) * (6.0 * u * t) + (p3 - p2) * (3.0 * t * t)
}
Seg2::Conic { a, apex, b, w } => {
let u = 1.0 - t;
let (n, d) = (a * (u * u) + apex * (2.0 * w * u * t) + b * (t * t), u * u + 2.0 * w * u * t + t * t);
let (dn, dd) = (a * (-2.0 * u) + apex * (2.0 * w * (1.0 - 2.0 * t)) + b * (2.0 * t), -2.0 * u + 2.0 * w * (1.0 - 2.0 * t) + 2.0 * t);
(dn * d - n * dd) / (d * d).max(1e-300)
}
}
}
pub fn start(&self) -> Vec2 {
match *self {
Seg2::Line { a, .. } => a,
Seg2::Line { a, .. } | Seg2::Conic { a, .. } => a,
Seg2::Cubic { p0, .. } => p0,
Seg2::Arc { center, radius, start, .. } => center + Vec2::from_angle(start) * radius,
}
}
pub fn end(&self) -> Vec2 {
match *self {
Seg2::Line { b, .. } => b,
Seg2::Line { b, .. } | Seg2::Conic { b, .. } => b,
Seg2::Cubic { p3, .. } => p3,
Seg2::Arc { center, radius, start, sweep } => center + Vec2::from_angle(start + sweep) * radius,
}
}
@@ -41,6 +85,15 @@ impl Seg2 {
match *self {
Seg2::Line { a, b } => a.lerp(b, t),
Seg2::Arc { center, radius, start, sweep } => center + Vec2::from_angle(start + sweep * t) * radius,
Seg2::Cubic { p0, p1, p2, p3 } => {
let u = 1.0 - t;
p0 * (u * u * u) + p1 * (3.0 * u * u * t) + p2 * (3.0 * u * t * t) + p3 * (t * t * t)
}
Seg2::Conic { a, apex, b, w } => {
let u = 1.0 - t;
let d = u * u + 2.0 * w * u * t + t * t;
(a * (u * u) + apex * (2.0 * w * u * t) + b * (t * t)) / d.max(1e-300)
}
}
}
pub fn mid(&self) -> Vec2 {
@@ -51,6 +104,8 @@ impl Seg2 {
match *self {
Seg2::Line { a, b } => (b - a).normalized().unwrap_or(Vec2::X),
Seg2::Arc { start, sweep, .. } => Vec2::from_angle(start).perp() * sweep.signum(),
Seg2::Cubic { p0, p1, p2, .. } => (p1 - p0).normalized().or_else(|| (p2 - p0).normalized()).unwrap_or(Vec2::X),
Seg2::Conic { a, apex, b, .. } => (apex - a).normalized().or_else(|| (b - a).normalized()).unwrap_or(Vec2::X),
}
}
/// Unit tangent at the end (direction of travel).
@@ -58,20 +113,63 @@ impl Seg2 {
match *self {
Seg2::Line { a, b } => (b - a).normalized().unwrap_or(Vec2::X),
Seg2::Arc { start, sweep, .. } => Vec2::from_angle(start + sweep).perp() * sweep.signum(),
Seg2::Cubic { p1, p2, p3, .. } => (p3 - p2).normalized().or_else(|| (p3 - p1).normalized()).unwrap_or(Vec2::X),
Seg2::Conic { a, apex, b, .. } => (b - apex).normalized().or_else(|| (b - a).normalized()).unwrap_or(Vec2::X),
}
}
pub fn reversed(&self) -> Seg2 {
match *self {
Seg2::Line { a, b } => Seg2::Line { a: b, b: a },
Seg2::Arc { center, radius, start, sweep } => Seg2::Arc { center, radius, start: start + sweep, sweep: -sweep },
Seg2::Cubic { p0, p1, p2, p3 } => Seg2::Cubic { p0: p3, p1: p2, p2: p1, p3: p0 },
Seg2::Conic { a, apex, b, w } => Seg2::Conic { a: b, apex, b: a, w },
}
}
pub fn length(&self) -> f64 {
match *self {
Seg2::Line { a, b } => a.dist(b),
Seg2::Arc { radius, sweep, .. } => (radius * sweep).abs(),
Seg2::Cubic { .. } | Seg2::Conic { .. } => self.integrate(|s, t| s.derivative(t).len()),
}
}
/// The two halves of the segment (same kind, exact).
pub fn split_half(&self) -> (Seg2, Seg2) {
match *self {
Seg2::Line { a, b } => {
let m = a.lerp(b, 0.5);
(Seg2::Line { a, b: m }, Seg2::Line { a: m, b })
}
Seg2::Arc { center, radius, start, sweep } => (
Seg2::Arc { center, radius, start, sweep: sweep / 2.0 },
Seg2::Arc { center, radius, start: start + sweep / 2.0, sweep: sweep / 2.0 },
),
Seg2::Cubic { p0, p1, p2, p3 } => {
let (a, b, c) = (p0.lerp(p1, 0.5), p1.lerp(p2, 0.5), p2.lerp(p3, 0.5));
let (d, e) = (a.lerp(b, 0.5), b.lerp(c, 0.5));
let m = d.lerp(e, 0.5);
(Seg2::Cubic { p0, p1: a, p2: d, p3: m }, Seg2::Cubic { p0: m, p1: e, p2: c, p3 })
}
Seg2::Conic { a, apex, b, w } => {
// Rational de Casteljau at the middle, weights normalised to 1 at the ends.
let m = self.point_at(0.5);
let w2 = ((1.0 + w) / 2.0).max(0.0).sqrt();
let (l, r) = ((a + apex * w) / (1.0 + w), (apex * w + b) / (1.0 + w));
(Seg2::Conic { a, apex: l, b: m, w: w2 }, Seg2::Conic { a: m, apex: r, b, w: w2 })
}
}
}
/// ∫₀¹ f(t) dt by composite Gauss–Legendre (16 pieces of 5 points).
fn integrate(&self, f: impl Fn(&Seg2, f64) -> f64) -> f64 {
let n = 16;
let mut sum = 0.0;
for k in 0..n {
let (t0, h) = (k as f64 / n as f64, 1.0 / n as f64);
for (x, w) in GL5 {
sum += w * h * f(self, t0 + x * h);
}
}
sum
}
/// Contribution to the loop's signed area (shoelace integral ½∮(x dy − y dx)).
pub fn area_term(&self) -> f64 {
match *self {
@@ -81,6 +179,7 @@ impl Seg2 {
// Chord term plus the circular segment between chord and arc.
0.5 * a.cross(b) + 0.5 * radius * radius * (sweep - sweep.sin())
}
Seg2::Cubic { .. } | Seg2::Conic { .. } => self.integrate(|s, t| 0.5 * s.point_at(t).cross(s.derivative(t))),
}
}
/// Polyline approximation (including both end points) with at most `tol` chord error.
@@ -93,6 +192,17 @@ impl Seg2 {
let n = ((sweep.abs() / step.max(1e-6)).ceil() as usize).clamp(2, 4096);
(0..=n).map(|i| self.point_at(i as f64 / n as f64)).collect()
}
Seg2::Cubic { p0, p1, p2, p3 } => {
// Chord error of n pieces ≈ max second difference / (8 n²).
let bow = (p0 - p1 * 2.0 + p2).len().max((p1 - p2 * 2.0 + p3).len()) * 6.0;
let n = ((bow / (8.0 * tol.max(1e-9))).sqrt().ceil() as usize).clamp(2, 4096);
(0..=n).map(|i| self.point_at(i as f64 / n as f64)).collect()
}
Seg2::Conic { a, apex, b, w } => {
let bow = (a - apex * 2.0 + b).len() * 2.0 * w.max(1.0);
let n = ((bow / (8.0 * tol.max(1e-9))).sqrt().ceil() as usize).clamp(2, 4096);
(0..=n).map(|i| self.point_at(i as f64 / n as f64)).collect()
}
}
}
}
+55
View File
@@ -83,3 +83,58 @@ fn tangent_chains_follow_smooth_edges() {
assert_eq!(m.tangent_chain(2, 5f64.to_radians()), vec![2]);
assert_eq!(m.face_edges(1), vec![1, 2]);
}
#[test]
fn cubic_and_conic_segments_measure_exactly() {
use crate::{Loop2, Seg2, Vec2};
// An ellipse (a = 5, b = 3) as four quarter conics: area π a b.
let w = std::f64::consts::FRAC_1_SQRT_2;
let (a, b) = (5.0, 3.0);
let q = |p0: Vec2, c: Vec2, p1: Vec2| Seg2::Conic { a: p0, apex: c, b: p1, w };
let e = Loop2 {
segs: vec![
q(Vec2::new(a, 0.0), Vec2::new(a, b), Vec2::new(0.0, b)),
q(Vec2::new(0.0, b), Vec2::new(-a, b), Vec2::new(-a, 0.0)),
q(Vec2::new(-a, 0.0), Vec2::new(-a, -b), Vec2::new(0.0, -b)),
q(Vec2::new(0.0, -b), Vec2::new(a, -b), Vec2::new(a, 0.0)),
],
};
assert!((e.signed_area() - std::f64::consts::PI * a * b).abs() < 1e-9, "{}", e.signed_area());
for s in &e.segs {
let m = s.point_at(0.37);
assert!(((m.x / a).powi(2) + (m.y / b).powi(2) - 1.0).abs() < 1e-12);
}
// A cubic on a straight line has the line's length; reversing flips the area sign.
let c = Seg2::Cubic { p0: Vec2::new(0.0, 0.0), p1: Vec2::new(1.0, 0.0), p2: Vec2::new(2.0, 0.0), p3: Vec2::new(3.0, 0.0) };
assert!((c.length() - 3.0).abs() < 1e-12);
// Square with one bulging cubic side: area = 1 + the cubic's area over the chord.
let bulge = Seg2::Cubic { p0: Vec2::new(1.0, 0.0), p1: Vec2::new(1.5, 1.0 / 3.0), p2: Vec2::new(1.5, 2.0 / 3.0), p3: Vec2::new(1.0, 1.0) };
let sq = Loop2 {
segs: vec![
Seg2::Line { a: Vec2::new(0.0, 0.0), b: Vec2::new(1.0, 0.0) },
bulge,
Seg2::Line { a: Vec2::new(1.0, 1.0), b: Vec2::new(0.0, 1.0) },
Seg2::Line { a: Vec2::new(0.0, 1.0), b: Vec2::new(0.0, 0.0) },
],
};
// x(t) - 1 = 1.5 t(1-t) (Bernstein: 3·0.5·t(1-t)²+3·0.5·t²(1-t)), y = t: ∫ 1.5 t(1-t) dt = 0.25.
assert!((sq.signed_area() - 1.25).abs() < 1e-12, "{}", sq.signed_area());
let rev = bulge.reversed();
assert!((rev.area_term() + bulge.area_term()).abs() < 1e-12);
assert!(bulge.polyline(1e-4).len() > 4);
// Halves cover the same curve.
let (l, r) = bulge.split_half();
assert!(l.point_at(0.5).dist(bulge.point_at(0.25)) < 1e-12);
assert!(r.point_at(0.5).dist(bulge.point_at(0.75)) < 1e-12);
// Conic halves are re-parameterised but stay on the ellipse.
let (l2, r2) = e.segs[0].split_half();
for t in [0.1, 0.5, 0.9] {
for q in [l2.point_at(t), r2.point_at(t)] {
assert!(((q.x / a).powi(2) + (q.y / b).powi(2) - 1.0).abs() < 1e-12, "{q:?}");
}
}
for (seg, (l, r)) in [(bulge, (l, r)), (e.segs[0], (l2, r2))] {
assert!(l.end().dist(r.start()) < 1e-12);
assert!((l.area_term() + r.area_term() - seg.area_term()).abs() < 1e-12);
}
}
+11 -11
View File
@@ -41,6 +41,9 @@ pub(crate) fn wire(plane: &Plane, lp: &Loop2) -> Result<mt::Wire> {
builder::circle_arc(a, b, p3(plane.to_world(s.mid())))
}
}
Seg2::Cubic { .. } | Seg2::Conic { .. } => {
crate::freeform::edge(plane, s, a, b).ok_or_else(|| KernelError::Invalid("profile curve".into()))?
}
};
edges.push(e);
}
@@ -66,6 +69,8 @@ fn region_ok(r: &Region2) -> Result<()> {
Seg2::Arc { center, radius, start, sweep } => {
center.is_finite() && radius.is_finite() && radius > 1e-9 && start.is_finite() && sweep.is_finite()
}
Seg2::Cubic { p0, p1, p2, p3 } => p0.is_finite() && p1.is_finite() && p2.is_finite() && p3.is_finite(),
Seg2::Conic { a, apex, b, w } => a.is_finite() && apex.is_finite() && b.is_finite() && w.is_finite() && w > 1e-9,
};
if !ok {
return Err(KernelError::Invalid("non-finite profile".into()));
@@ -234,7 +239,7 @@ fn offset_loop(lp: &Loop2, d: f64) -> Result<Loop2> {
Err(KernelError::Invalid("the taper closes the profile".into()))
}
}
Seg2::Line { .. } => Err(KernelError::Failed("not supported yet: tapered profiles mixing lines and arcs".into())),
_ => Err(KernelError::Failed("not supported yet: tapered profiles mixing lines and arcs".into())),
})
.collect::<Result<Vec<_>>>()?;
return Ok(Loop2 { segs });
@@ -377,6 +382,10 @@ fn revolve_touching_axis(plane: &Plane, lp: &Loop2, axis_origin: Vec2, d: Vec2)
edges.push(match *s {
Seg2::Line { .. } => builder::line(a, b),
Seg2::Arc { .. } => builder::circle_arc(a, b, p3(plane.to_world(s.mid()))),
Seg2::Cubic { .. } | Seg2::Conic { .. } => match crate::freeform::edge(plane, s, a, b) {
Some(e) => e,
None => builder::line(a, b),
},
});
}
let mut w: mt::Wire = edges.into();
@@ -459,16 +468,7 @@ fn resample(lp: &Loop2, n: usize, angles: &[f64]) -> Loop2 {
while segs.len() < n {
let Some((k, _)) = segs.iter().enumerate().max_by(|a, b| a.1.length().total_cmp(&b.1.length())) else { break };
let s = segs.remove(k);
let (p, q) = match s {
Seg2::Line { a, b } => {
let m = a.lerp(b, 0.5);
(Seg2::Line { a, b: m }, Seg2::Line { a: m, b })
}
Seg2::Arc { center, radius, start, sweep } => (
Seg2::Arc { center, radius, start, sweep: sweep / 2.0 },
Seg2::Arc { center, radius, start: start + sweep / 2.0, sweep: sweep / 2.0 },
),
};
let (p, q) = s.split_half();
segs.insert(k, q);
segs.insert(k, p);
}
+25
View File
@@ -0,0 +1,25 @@
//! Free-form profile edges: cubic Béziers as B-spline curves and conics as rational quadratic
//! NURBS curves (exact ellipses, parabolas and hyperbolas from sketches).
use solvecraft_geom::{Plane, Seg2};
use truck_modeling as mt;
use crate::body::p3;
/// The edge for a free-form segment between two existing vertices (None for lines and arcs).
pub(crate) fn edge(plane: &Plane, s: &Seg2, a: &mt::Vertex, b: &mt::Vertex) -> Option<mt::Edge> {
match *s {
Seg2::Cubic { p1, p2, .. } => Some(mt::builder::bezier(a, b, vec![p3(plane.to_world(p1)), p3(plane.to_world(p2))])),
Seg2::Conic { apex, w, .. } => {
let (pa, pb, px) = (a.point(), b.point(), p3(plane.to_world(apex)));
let ctrl = vec![
mt::Vector4::new(pa.x, pa.y, pa.z, 1.0),
mt::Vector4::new(px.x * w, px.y * w, px.z * w, w),
mt::Vector4::new(pb.x, pb.y, pb.z, 1.0),
];
let curve = mt::NurbsCurve::new(mt::BSplineCurve::new(mt::KnotVec::bezier_knot(2), ctrl));
Some(mt::Edge::new(a, b, mt::Curve::NurbsCurve(curve)))
}
Seg2::Line { .. } | Seg2::Arc { .. } => None,
}
}
+1
View File
@@ -16,6 +16,7 @@ mod blend;
mod body;
mod build;
mod coplanar;
mod freeform;
mod heal;
mod helix;
mod loopblend;
+76
View File
@@ -163,6 +163,53 @@ pub fn end_curvature(f: &dyn Fn(f64) -> Vec2, at_end: bool) -> Vec2 {
(d2 - t * t.dot(d2)) / l2
}
/// The cubic Bézier pieces of a fit spline (exact: each span is a cubic Hermite).
pub fn fit_beziers(p: &[Vec2]) -> Vec<[Vec2; 4]> {
fit_spans(p).into_iter().map(|(p0, p1, t0, t1)| [p0, p0 + t0 / 3.0, p1 - t1 / 3.0, p1]).collect()
}
/// The Bézier pieces (each `degree + 1` points) of the clamped uniform B-spline over `ctrl`
/// (Boehm knot insertion until every interior knot has multiplicity `degree`).
pub fn bspline_beziers(ctrl: &[Vec2], degree: usize) -> Vec<Vec<Vec2>> {
let n = ctrl.len();
let k = degree.min(n.saturating_sub(1)).max(1);
let spans = n - k;
let mut knots = vec![0.0; k + 1];
for i in 1..spans {
knots.push(i as f64 / spans as f64);
}
knots.extend(std::iter::repeat_n(1.0, k + 1));
let mut pts = ctrl.to_vec();
for s in 1..spans {
let u = s as f64 / spans as f64;
for _ in 1..k {
// Span holding u (the last knot ≤ u).
let Some(span) = knots.iter().rposition(|x| *x <= u + 1e-15) else { break };
let mut q = Vec::with_capacity(pts.len() + 1);
for i in 0..=pts.len() {
let p = if i + k <= span {
pts.get(i).copied()
} else if i > span {
pts.get(i - 1).copied()
} else {
let (lo, hi) = (knots.get(i).copied().unwrap_or(0.0), knots.get(i + k).copied().unwrap_or(1.0));
let a = if hi - lo > 1e-15 { (u - lo) / (hi - lo) } else { 0.0 };
match (pts.get(i.wrapping_sub(1)), pts.get(i)) {
(Some(x), Some(y)) => Some(*x * (1.0 - a) + *y * a),
_ => None,
}
};
if let Some(p) = p {
q.push(p);
}
}
knots.insert(span + 1, u);
pts = q;
}
}
pts.windows(k + 1).step_by(k).map(|w| w.to_vec()).collect()
}
/// Polyline of a spline through (`control == false`) or over (`control`) the points.
pub fn spline_polyline(p: &[Vec2], control: bool, degree: u8) -> Vec<Vec2> {
if p.len() < 2 {
@@ -244,6 +291,35 @@ mod tests {
assert!(end_curvature(&|t| spline_point(&p, false, 3, t), false).len() < 1e-3);
}
#[test]
fn spline_bezier_pieces_match_the_spline() {
let p = [Vec2::new(0.0, 0.0), Vec2::new(10.0, 5.0), Vec2::new(20.0, -3.0), Vec2::new(30.0, 8.0), Vec2::new(35.0, 0.0)];
let bez = |c: &[Vec2], t: f64| -> Vec2 {
// de Casteljau.
let mut v = c.to_vec();
while v.len() > 1 {
v = v.windows(2).map(|w| w[0] * (1.0 - t) + w[1] * t).collect();
}
v[0]
};
for degree in [2usize, 3] {
let pieces = bspline_beziers(&p, degree);
assert_eq!(pieces.len(), p.len() - degree, "degree {degree}");
let spans = pieces.len();
for (i, c) in pieces.iter().enumerate() {
for t in [0.0, 0.3, 0.7, 1.0] {
let g = (i as f64 + t) / spans as f64;
let want = spline_point(&p, true, degree as u8, g);
assert!(bez(c, t).dist(want) < 1e-9, "degree {degree} piece {i} t {t}: {:?} vs {want:?}", bez(c, t));
}
}
}
for (i, c) in fit_beziers(&p).iter().enumerate() {
let g = (i as f64 + 0.4) / 4.0;
assert!(bez(c, 0.4).dist(spline_point(&p, false, 3, g)) < 1e-9);
}
}
#[test]
fn control_spline_is_clamped_and_conic_hits_its_ends() {
let p = [Vec2::new(0.0, 0.0), Vec2::new(0.0, 10.0), Vec2::new(10.0, 10.0), Vec2::new(10.0, 0.0)];
+3 -1
View File
@@ -18,7 +18,9 @@ mod solver;
mod text;
mod wire;
pub use curves::{MAX_SPLINE_POINTS, conic_point, ellipse_point, end_curvature, spline_end_tangent, spline_point, spline_polyline};
pub use curves::{
MAX_SPLINE_POINTS, bspline_beziers, conic_point, ellipse_point, end_curvature, fit_beziers, spline_end_tangent, spline_point, spline_polyline,
};
pub use link::{Link, LinkGeom, LinkKind, LinkSource, MAX_LINK_ENTITIES};
pub use model::{Constraint, ConstraintKind, Curve, CurveKind, SPoint, Sketch, SketchError};
pub use profiles::{Profile, find_drawn_profiles, find_profiles, merge_regions};
+46
View File
@@ -720,6 +720,52 @@ impl Sketch {
}
}
/// Exact boundary segments of a free-form curve: ellipses as four conics, conics, fit
/// splines and control splines up to degree 3 as Béziers. None for other curves (and for
/// higher-degree splines, which stay polylines).
pub fn exact_segs(&self, ci: usize) -> Option<Vec<Seg2>> {
let cu = self.curves.get(ci)?;
let pt = |i: &usize| self.point(*i);
match &cu.kind {
CurveKind::Ellipse { c, m, r } => {
let (c, m) = (pt(c)?, pt(m)?);
let u = m - c;
let v = u.perp().normalized()? * *r;
let w = std::f64::consts::FRAC_1_SQRT_2;
let corners = [c + u, c + v, c - u, c - v];
Some(
(0..4)
.map(|i| {
let (a, b) = (corners[i], corners[(i + 1) % 4]);
Seg2::Conic { a, apex: a + b - c, b, w }
})
.collect(),
)
}
CurveKind::Conic { a, b, apex, rho } => Some(vec![Seg2::Conic { a: pt(a)?, apex: pt(apex)?, b: pt(b)?, w: rho / (1.0 - rho).max(1e-9) }]),
CurveKind::Spline { pts, control: false, .. } => {
let p: Vec<Vec2> = pts.iter().map(pt).collect::<Option<_>>()?;
Some(crate::curves::fit_beziers(&p).into_iter().map(|[p0, p1, p2, p3]| Seg2::Cubic { p0, p1, p2, p3 }).collect())
}
CurveKind::Spline { pts, control: true, degree } if *degree <= 3 => {
let p: Vec<Vec2> = pts.iter().map(pt).collect::<Option<_>>()?;
let k = (*degree as usize).min(p.len().saturating_sub(1)).max(1);
Some(
crate::curves::bspline_beziers(&p, k)
.into_iter()
.filter_map(|c| match c[..] {
[a, b] => Some(Seg2::Line { a, b }),
[a, x, b] => Some(Seg2::Conic { a, apex: x, b, w: 1.0 }),
[p0, p1, p2, p3] => Some(Seg2::Cubic { p0, p1, p2, p3 }),
_ => None,
})
.collect(),
)
}
_ => None,
}
}
/// Polyline of a free-form curve (ellipse, spline, conic); other curves: their segments'
/// polylines.
pub fn polyline(&self, ci: usize) -> Vec<Vec2> {
+79 -9
View File
@@ -46,8 +46,12 @@ struct Edge {
u: usize,
v: usize,
id: String,
/// A piece of an exact free-form segment: (curve, segment, piece, pieces).
tag: Option<Piece>,
}
type Piece = (usize, usize, usize, usize);
/// Parameter of `p` at an end of the shape (lines: 0 or 1), for intersections at end points.
fn end_param(sh: &Shape, p: Vec2, tol: f64) -> Option<f64> {
match *sh {
@@ -68,6 +72,41 @@ fn end_param(sh: &Shape, p: Vec2, tol: f64) -> Option<f64> {
}
}
/// Put exact free-form segments back where a loop runs over all the pieces of one, in order.
fn restore_exact(segs: Vec<Seg2>, tags: &[Option<(Piece, bool)>], exact: &std::collections::HashMap<usize, Vec<Seg2>>) -> Vec<Seg2> {
let n = segs.len();
if n == 0 || tags.len() != n || tags.iter().all(Option::is_none) {
return segs;
}
let key = |i: usize| tags.get(i % n).copied().flatten().map(|((ci, k, _, _), fwd)| (ci, k, fwd));
// Start where a run begins so no run wraps around.
let r = (0..n).find(|i| key(*i) != key(*i + n - 1)).unwrap_or(0);
let mut out = Vec::with_capacity(n);
let mut i = 0;
while i < n {
let k0 = key(r + i);
let mut j = i + 1;
while j < n && key(r + j) == k0 && k0.is_some() {
j += 1;
}
let run: Vec<usize> = (i..j).map(|x| (r + x) % n).collect();
let whole = match (k0, tags.get(run[0]).copied().flatten()) {
(Some((ci, k, fwd)), Some(((_, _, _, cnt), _))) if run.len() == cnt => {
let order: Vec<usize> = run.iter().filter_map(|x| tags.get(*x).copied().flatten().map(|((_, _, p, _), _)| p)).collect();
let want: Vec<usize> = if fwd { (0..cnt).collect() } else { (0..cnt).rev().collect() };
if order == want { exact.get(&ci).and_then(|e| e.get(k)).map(|s| if fwd { *s } else { s.reversed() }) } else { None }
}
_ => None,
};
match whole {
Some(s) => out.push(s),
None => out.extend(run.iter().filter_map(|x| segs.get(*x).copied())),
}
i = j;
}
out
}
/// A uniform grid over a set of boxes (about one box per cell), for neighbour queries.
struct Grid {
lo: Vec2,
@@ -122,6 +161,11 @@ fn curvature(s: &Seg2) -> f64 {
match *s {
Seg2::Line { .. } => 0.0,
Seg2::Arc { radius, sweep, .. } => sweep.signum() / radius.max(1e-12),
Seg2::Cubic { .. } | Seg2::Conic { .. } => {
let (d1, h) = (s.derivative(0.0), 1e-5);
let d2 = (s.derivative(h) - d1) / h;
d1.cross(d2) / d1.len().powi(3).max(1e-300)
}
}
}
@@ -145,6 +189,8 @@ pub fn find_profiles(sk: &Sketch) -> Vec<Profile> {
let tol = MERGE_TOL * 10.0;
// Analytic shapes of the profile curves.
let mut shapes: Vec<(usize, Shape)> = Vec::new();
let mut tags: Vec<Option<Piece>> = Vec::new();
let mut exact: std::collections::HashMap<usize, Vec<Seg2>> = std::collections::HashMap::new();
for (ci, c) in sk.curves.iter().enumerate() {
if c.construction || c.centerline {
continue;
@@ -167,15 +213,33 @@ pub fn find_profiles(sk: &Sketch) -> Vec<Profile> {
shapes.push((ci, Shape::Round { c: center, r: radius, start, sweep }));
}
}
// Free-form curves take part as their polylines.
_ => {
for s in sk.segs(ci) {
if let Seg2::Line { a, b } = s {
shapes.push((ci, Shape::Line { a, b }));
// Free-form curves take part as polylines of their exact segments; loops that use
// a whole segment get the exact one back.
_ => match sk.exact_segs(ci) {
Some(ex) => {
for (k, seg) in ex.iter().enumerate() {
let poly = seg.polyline(1e-3);
let n = poly.len().saturating_sub(1);
for (j, w) in poly.windows(2).enumerate() {
if w[0].dist(w[1]) > MERGE_TOL {
shapes.push((ci, Shape::Line { a: w[0], b: w[1] }));
tags.resize(shapes.len() - 1, None);
tags.push(Some((ci, k, j, n)));
}
}
}
exact.insert(ci, ex);
}
None => {
for s in sk.segs(ci) {
if let Seg2::Line { a, b } = s {
shapes.push((ci, Shape::Line { a, b }));
}
}
}
}
},
}
tags.resize(shapes.len(), None);
if shapes.len() > MAX_EDGES {
break;
}
@@ -344,11 +408,13 @@ pub fn find_profiles(sk: &Sketch) -> Vec<Profile> {
let mut ps = vec![a];
ps.extend(ts.iter().map(|t| sh.point(*t)));
ps.push(b);
// Pieces split by other curves no longer stand for their whole piece.
let tag = if ts.is_empty() { tags.get(k).copied().flatten() } else { None };
for w in ps.windows(2) {
let (p, q) = (w[0], w[1]);
if p.dist(q) > MERGE_TOL {
let (u, v) = (node(p), node(q));
edges.push(Edge { seg: Seg2::Line { a: p, b: q }, u, v, id: id.clone() });
edges.push(Edge { seg: Seg2::Line { a: p, b: q }, u, v, id: id.clone(), tag });
}
}
}
@@ -371,7 +437,7 @@ pub fn find_profiles(sk: &Sketch) -> Vec<Profile> {
}
let seg = Seg2::Arc { center: c, radius: r, start: start + t0, sweep: t1 - t0 };
let (u, v) = (node(seg.start()), node(seg.end()));
edges.push(Edge { seg, u, v, id: id.clone() });
edges.push(Edge { seg, u, v, id: id.clone(), tag: None });
}
}
}
@@ -441,6 +507,7 @@ pub fn find_profiles(sk: &Sketch) -> Vec<Profile> {
}
let mut segs = Vec::new();
let mut ids = Vec::new();
let mut ptags: Vec<Option<(Piece, bool)>> = Vec::new();
let mut h = start;
let mut ok = true;
for _ in 0..=edges.len() * 2 {
@@ -454,6 +521,9 @@ pub fn find_profiles(sk: &Sketch) -> Vec<Profile> {
segs.push(seg);
if let Some(e) = edges.get(h / 2) {
ids.push(e.id.clone());
ptags.push(e.tag.map(|t| (t, h.is_multiple_of(2))));
} else {
ptags.push(None);
}
let back = (-seg.end_tangent()).angle();
let twin = h ^ 1;
@@ -496,7 +566,7 @@ pub fn find_profiles(sk: &Sketch) -> Vec<Profile> {
if !ok || h != start {
continue;
}
let lp = Loop2 { segs };
let lp = Loop2 { segs: restore_exact(segs, &ptags, &exact) };
if lp.signed_area() > 1e-9 {
let c = edges.get(start / 2).map(|e| find(&mut comp, e.u)).unwrap_or(0);
loops.push((lp, ids, c));
+1 -1
View File
@@ -159,7 +159,7 @@ impl Sketch {
CurveKind::Circle { c, r } => Some(Shape::Round { c: self.point(c)?, r, start: 0.0, sweep: TAU }),
CurveKind::Arc { .. } => match self.segs(ci).into_iter().next()? {
Seg2::Arc { center, radius, start, sweep } => Some(Shape::Round { c: center, r: radius, start, sweep }),
Seg2::Line { .. } => None,
_ => None,
},
_ => None,
}