diff --git a/crates/doc/src/eval.rs b/crates/doc/src/eval.rs index c05eaf2c..58d01bf4 100644 --- a/crates/doc/src/eval.rs +++ b/crates/doc/src/eval.rs @@ -1106,10 +1106,14 @@ fn path_segments(ps: &SolvedSketch, ids: &[String], start: Vec3) -> Result 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()) diff --git a/crates/doc/src/project.rs b/crates/doc/src/project.rs index 281f998f..9f94cfff 100644 --- a/crates/doc/src/project.rs +++ b/crates/doc/src/project.rs @@ -424,7 +424,7 @@ fn tilted_circle(plane: &Plane, pts: &[Vec3]) -> Option> { 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> { } 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; diff --git a/crates/engine/src/cmd/sketch_modify_tests.rs b/crates/engine/src/cmd/sketch_modify_tests.rs index 64525d59..eab71904 100644 --- a/crates/engine/src/cmd/sketch_modify_tests.rs +++ b/crates/engine/src/cmd/sketch_modify_tests.rs @@ -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}"); +} diff --git a/crates/geom/src/profile.rs b/crates/geom/src/profile.rs index 0d5ee390..7197c90c 100644 --- a/crates/geom/src/profile.rs +++ b/crates/geom/src/profile.rs @@ -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() + } } } } diff --git a/crates/geom/src/tests.rs b/crates/geom/src/tests.rs index bd6fd93c..a8be2b49 100644 --- a/crates/geom/src/tests.rs +++ b/crates/geom/src/tests.rs @@ -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); + } +} diff --git a/crates/kernel/src/build.rs b/crates/kernel/src/build.rs index e4527917..bc50205d 100644 --- a/crates/kernel/src/build.rs +++ b/crates/kernel/src/build.rs @@ -41,6 +41,9 @@ pub(crate) fn wire(plane: &Plane, lp: &Loop2) -> Result { 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 { 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::>>()?; 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); } diff --git a/crates/kernel/src/freeform.rs b/crates/kernel/src/freeform.rs new file mode 100644 index 00000000..53c651ce --- /dev/null +++ b/crates/kernel/src/freeform.rs @@ -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 { + 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, + } +} diff --git a/crates/kernel/src/lib.rs b/crates/kernel/src/lib.rs index f5ab4068..f3c21e50 100644 --- a/crates/kernel/src/lib.rs +++ b/crates/kernel/src/lib.rs @@ -16,6 +16,7 @@ mod blend; mod body; mod build; mod coplanar; +mod freeform; mod heal; mod helix; mod loopblend; diff --git a/crates/sketch/src/curves.rs b/crates/sketch/src/curves.rs index 130739e2..debfbadc 100644 --- a/crates/sketch/src/curves.rs +++ b/crates/sketch/src/curves.rs @@ -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> { + 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 { 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)]; diff --git a/crates/sketch/src/lib.rs b/crates/sketch/src/lib.rs index 4eb775f0..349a2362 100644 --- a/crates/sketch/src/lib.rs +++ b/crates/sketch/src/lib.rs @@ -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}; diff --git a/crates/sketch/src/model.rs b/crates/sketch/src/model.rs index 4f35a9c1..a96e8d5f 100644 --- a/crates/sketch/src/model.rs +++ b/crates/sketch/src/model.rs @@ -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> { + 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 = pts.iter().map(pt).collect::>()?; + 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 = pts.iter().map(pt).collect::>()?; + 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 { diff --git a/crates/sketch/src/profiles.rs b/crates/sketch/src/profiles.rs index e1b7525c..a54fc8cd 100644 --- a/crates/sketch/src/profiles.rs +++ b/crates/sketch/src/profiles.rs @@ -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, } +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 { match *sh { @@ -68,6 +72,41 @@ fn end_param(sh: &Shape, p: Vec2, tol: f64) -> Option { } } +/// Put exact free-form segments back where a loop runs over all the pieces of one, in order. +fn restore_exact(segs: Vec, tags: &[Option<(Piece, bool)>], exact: &std::collections::HashMap>) -> Vec { + 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 = (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 = run.iter().filter_map(|x| tags.get(*x).copied().flatten().map(|((_, _, p, _), _)| p)).collect(); + let want: Vec = 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 { let tol = MERGE_TOL * 10.0; // Analytic shapes of the profile curves. let mut shapes: Vec<(usize, Shape)> = Vec::new(); + let mut tags: Vec> = Vec::new(); + let mut exact: std::collections::HashMap> = 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 { 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 { 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 { } 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 { } let mut segs = Vec::new(); let mut ids = Vec::new(); + let mut ptags: Vec> = 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 { 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 { 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)); diff --git a/crates/sketch/src/shape.rs b/crates/sketch/src/shape.rs index 86d5b763..56446b18 100644 --- a/crates/sketch/src/shape.rs +++ b/crates/sketch/src/shape.rs @@ -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, }