From b9bc75670c71fb13dae245288d51f1552595cf79 Mon Sep 17 00:00:00 2001 From: bherbruck Date: Thu, 8 Oct 2026 07:27:40 +0000 Subject: [PATCH] =?UTF-8?q?Exact=20profiles=20for=20splines=20of=20any=20d?= =?UTF-8?q?egree=20(a=20degree-4=E2=80=A67=20B=C3=A9zier=20segment=20kind,?= =?UTF-8?q?=20Seg2::Bezier,=20extruded=20as=20B-spline=20edges)=20and=20fo?= =?UTF-8?q?r=20free-form=20curves=20cut=20by=20crossings=20(exact=20sub-cu?= =?UTF-8?q?rves=20whose=20ends=20are=20moved=20to=20the=20true=20crossing,?= =?UTF-8?q?=20found=20by=20Newton;=20polylines=20only=20where=20that=20fai?= =?UTF-8?q?ls)?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01AS5Tn3xieewUSg3uxrMVdL --- ROADMAP.md | 6 +- crates/doc/src/eval.rs | 2 +- crates/engine/src/cmd/sketch_freeform.rs | 9 + crates/engine/src/cmd/sketch_modify_tests.rs | 29 +++ crates/geom/src/lib.rs | 2 +- crates/geom/src/profile.rs | 89 +++++++- crates/kernel/src/build.rs | 5 +- crates/kernel/src/freeform.rs | 6 +- crates/sketch/src/model.rs | 14 +- crates/sketch/src/profiles.rs | 211 +++++++++++++++++-- crates/sketch/src/tests.rs | 50 +++++ 11 files changed, 389 insertions(+), 34 deletions(-) diff --git a/ROADMAP.md b/ROADMAP.md index 8a58e2c2..d50cb9cb 100644 --- a/ROADMAP.md +++ b/ROADMAP.md @@ -40,7 +40,7 @@ headline counts in-scope tabs only and lists the deferred ones separately. | # | Milestone | Status | Hours left | |---|---|---|---| | M0 | Vertical slice | done | — | -| M1 | Sketch depth: done — projection/intersect/include as linked geometry (auto-project on faces, lost-reference warnings, break link), full SKETCH create/modify set (slots, polygons, tangent arcs/circles, ellipse, fit/control splines, conic, text, trim/extend/break, fillet/chamfer, offset, mirror, patterns, move/scale, blend curve, centerline), constraints incl. curvature and polygon with over-constraint refusal, driven/arc-length/linear-diameter dimensions, AutoConstrain, DXF/SVG insert, 3D sketch curves (Project To Surface, Intersection Curve, Isoparametric Curve, Spun Profile), exact free-form profiles (B-spline/NURBS edges). Left: sweeps along 3D sketch curves, degree 4+ splines as exact profiles | in progress | 4 | +| M1 | Sketch depth: done — projection/intersect/include as linked geometry (auto-project on faces, lost-reference warnings, break link), full SKETCH create/modify set (slots, polygons, tangent arcs/circles, ellipse, fit/control splines, conic, text, trim/extend/break, fillet/chamfer, offset, mirror, patterns, move/scale, blend curve, centerline), constraints incl. curvature and polygon with over-constraint refusal, driven/arc-length/linear-diameter dimensions, AutoConstrain, DXF/SVG insert, 3D sketch curves (Project To Surface, Intersection Curve, Isoparametric Curve, Spun Profile), exact free-form profiles (B-spline/NURBS edges, any spline degree, curves cut by crossings), sweeps and pipes along 3D sketch curves, canvases and decals, curvature combs | done | 4 | | M2 | Robust modelling: coplanar booleans in the kernel, holes (placed on faces and on sketch points), threads (cosmetic, then modelled), patterns, mirror, shell, draft, split, construction geometry | in progress | 16 | | M3 | General fillets and chamfers (curved edges, chains, vertex blends) | planned | 30 | | M4 | Files: STEP import (done: all 29 Fusion STEP files — [docs/step-import.md](docs/step-import.md)), 3MF export and 3MF/STL import as mesh bodies (done), DXF sketches, OBJ import | in progress | 6 | @@ -77,8 +77,8 @@ headline counts in-scope tabs only and lists the deferred ones separately. and assembly components are flattened into bodies. Closed periodic faces (torus, B-spline bands) are split in two for meshing; measure() uses a finer tolerance for small radii. - No persistent naming (fillet edges are re-found by position; projected sketch geometry too). -- Control-point splines of degree 4 and up, and free-form curves split by crossings, enter - profiles as polylines; 3D sketch curves are reference geometry (not sweep paths yet). +- Free-form curves cut at a tangent touch (no clean crossing) enter that profile as polylines. + Sweeps along 3D sketch curves follow their polyline (smooth walls when the path is smooth). - Shell and offset faces work on bodies of planes and cylinders (convex or not: filleted boxes, bosses); draft on convex planar bodies and on walls (planes, cylinders → cones) between caps; loft is ruled. diff --git a/crates/doc/src/eval.rs b/crates/doc/src/eval.rs index 2ab5aaa3..d5db0d5d 100644 --- a/crates/doc/src/eval.rs +++ b/crates/doc/src/eval.rs @@ -1223,7 +1223,7 @@ fn path_segments(ps: &SolvedSketch, ids: &[String], start: Vec3) -> Result { + solvecraft_geom::Seg2::Cubic { .. } | solvecraft_geom::Seg2::Conic { .. } | solvecraft_geom::Seg2::Bezier { .. } => { s.polyline(1e-2).windows(2).map(|q| kernel::PathSeg::Line { a: w(q[0]), b: w(q[1]) }).collect() } }) diff --git a/crates/engine/src/cmd/sketch_freeform.rs b/crates/engine/src/cmd/sketch_freeform.rs index b15c1825..0ec0e92f 100644 --- a/crates/engine/src/cmd/sketch_freeform.rs +++ b/crates/engine/src/cmd/sketch_freeform.rs @@ -153,6 +153,15 @@ fn add_pieces(sk: &mut Sketch, segs: &[Seg2], start: Option, end: Option< let x = sk.add_point(apex, None)?; CurveKind::Conic { a: prev, b: to, apex: x, rho: w / (1.0 + w) } } + Seg2::Bezier { n, p } => { + let k = (n as usize).clamp(1, solvecraft_geom::MAX_BEZIER_DEGREE); + let mut pts = vec![prev]; + for q in p.iter().take(k).skip(1) { + pts.push(sk.add_point(*q, None)?); + } + pts.push(to); + CurveKind::Spline { pts, control: true, degree: n } + } Seg2::Arc { .. } => return Err(bad("trim", "unexpected arc piece")), }; let c = sk.add_curve(kind, None)?; diff --git a/crates/engine/src/cmd/sketch_modify_tests.rs b/crates/engine/src/cmd/sketch_modify_tests.rs index 38e992e4..fad5fe75 100644 --- a/crates/engine/src/cmd/sketch_modify_tests.rs +++ b/crates/engine/src/cmd/sketch_modify_tests.rs @@ -343,6 +343,35 @@ fn free_form_profiles_extrude_exactly() { assert!(m["bodies"][0]["faces"].as_u64().unwrap() <= 5, "{m}"); } +#[test] +fn degree_five_and_cut_free_form_profiles_extrude_exactly() { + let vol = |s: &mut Session| run(s, "MeasureCommand", json!({}))["total"]["volume_mm3"].as_f64().unwrap(); + // A degree-5 control spline closed by a line: one B-spline side face. + let mut s = new_sketch(); + let sp = ids(&run(&mut s, "DrawCVMSpline5D", json!({"points": [[0, 0], [4, 10], [10, 12], [16, 12], [22, 10], [26, 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": 3})); + let v = vol(&mut s); + assert!((v - a * 3.0).abs() / (a * 3.0) < 2e-4, "{v} vs {}", a * 3.0); + let m = run(&mut s, "MeasureCommand", json!({})); + assert_eq!(m["bodies"][0]["faces"], 4, "{m}"); + // A fit spline arch cut by a line: both sides extrude with exact curved faces. + let mut s = new_sketch(); + let sp = ids(&run(&mut s, "DrawSpline", json!({"points": [[0, 0], [6, 9], [14, 9], [20, 0]]}))["curves"])[0].clone(); + run(&mut s, "DrawPolyline", json!({"points": [format!("{sp}.end"), format!("{sp}.start")]})); + run(&mut s, "DrawPolyline", json!({"points": [[10, -5], [10, 20]]})); + let areas = profiles(&s); + assert_eq!(areas.len(), 2, "{areas:?}"); + run(&mut s, "SketchStop", json!({})); + run(&mut s, "Extrude", json!({"distance": 2, "profiles": [0]})); + let v = vol(&mut s); + assert!((v - areas[0] * 2.0).abs() / (areas[0] * 2.0) < 2e-4, "{v} vs {}", areas[0] * 2.0); + let m = run(&mut s, "MeasureCommand", json!({})); + assert!(m["bodies"][0]["faces"].as_u64().unwrap() <= 7, "no polyline facets: {m}"); +} + #[test] fn trim_and_break_free_form_curves_exactly() { let pi = std::f64::consts::PI; diff --git a/crates/geom/src/lib.rs b/crates/geom/src/lib.rs index 55c05ffd..dc2cbbdd 100644 --- a/crates/geom/src/lib.rs +++ b/crates/geom/src/lib.rs @@ -13,7 +13,7 @@ mod vec; pub use mesh::{Aabb3, Mesh, MeshMeasure}; pub use plane::Plane; -pub use profile::{Loop2, Region2, Seg2}; +pub use profile::{Loop2, MAX_BEZIER_DEGREE, Region2, Seg2, bezier}; pub use vec::{Vec2, Vec3}; /// Geometric tolerance for coincidence tests (mm). diff --git a/crates/geom/src/profile.rs b/crates/geom/src/profile.rs index 7fc235e1..71b82358 100644 --- a/crates/geom/src/profile.rs +++ b/crates/geom/src/profile.rs @@ -38,6 +38,52 @@ pub enum Seg2 { b: Vec2, w: f64, }, + /// Bézier of degree `n` (4…7) over `p[0..=n]` (the rest unused): pieces of higher-degree + /// splines. + Bezier { + n: u8, + p: [Vec2; 8], + }, +} + +/// Highest degree of [`Seg2::Bezier`]. +pub const MAX_BEZIER_DEGREE: usize = 7; + +/// A Bézier segment over `pts` (2…8 points): a line, a cubic or a higher-degree piece. A +/// quadratic becomes a conic of weight 1. +pub fn bezier(pts: &[Vec2]) -> Option { + Some(match *pts { + [a, b] => Seg2::Line { a, b }, + [a, apex, b] => Seg2::Conic { a, apex, b, w: 1.0 }, + [p0, p1, p2, p3] => Seg2::Cubic { p0, p1, p2, p3 }, + _ if pts.len() <= MAX_BEZIER_DEGREE + 1 => { + let mut p = [Vec2::ZERO; 8]; + for (d, s) in p.iter_mut().zip(pts) { + *d = *s; + } + Seg2::Bezier { n: (pts.len() - 1) as u8, p } + } + _ => return None, + }) +} + +/// de Casteljau: the point at `t` and the two halves' control points. +fn casteljau(c: &[Vec2], t: f64) -> (Vec2, Vec, Vec) { + let mut w = c.to_vec(); + let mut left = Vec::with_capacity(c.len()); + let mut right = Vec::with_capacity(c.len()); + while let Some(first) = w.first().copied() { + left.push(first); + if let Some(last) = w.last().copied() { + right.push(last); + } + if w.len() == 1 { + break; + } + w = w.windows(2).map(|x| x[0].lerp(x[1], t)).collect(); + } + right.reverse(); + (left.last().copied().unwrap_or(Vec2::ZERO), left, right) } /// Gauss–Legendre nodes and weights on [0, 1] (5 points). @@ -50,6 +96,9 @@ const GL5: [(f64, f64); 5] = [ ]; impl Seg2 { + fn bez(n: u8, p: &[Vec2; 8]) -> &[Vec2] { + &p[..(n as usize).clamp(1, MAX_BEZIER_DEGREE) + 1] + } /// 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 { @@ -65,12 +114,18 @@ impl Seg2 { 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) } + Seg2::Bezier { n, ref p } => { + let c = Self::bez(n, p); + let diffs: Vec = c.windows(2).map(|w| (w[1] - w[0]) * (c.len() - 1) as f64).collect(); + casteljau(&diffs, t).0 + } } } pub fn start(&self) -> Vec2 { match *self { Seg2::Line { a, .. } | Seg2::Conic { a, .. } => a, Seg2::Cubic { p0, .. } => p0, + Seg2::Bezier { p, .. } => p[0], Seg2::Arc { center, radius, start, .. } => center + Vec2::from_angle(start) * radius, } } @@ -78,6 +133,7 @@ impl Seg2 { match *self { Seg2::Line { b, .. } | Seg2::Conic { b, .. } => b, Seg2::Cubic { p3, .. } => p3, + Seg2::Bezier { n, ref p } => Self::bez(n, p).last().copied().unwrap_or(p[0]), Seg2::Arc { center, radius, start, sweep } => center + Vec2::from_angle(start + sweep) * radius, } } @@ -94,6 +150,7 @@ impl Seg2 { 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) } + Seg2::Bezier { n, ref p } => casteljau(Self::bez(n, p), t).0, } } pub fn mid(&self) -> Vec2 { @@ -106,6 +163,10 @@ impl Seg2 { 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), + Seg2::Bezier { n, ref p } => { + let c = Self::bez(n, p); + c.iter().skip(1).find_map(|q| (*q - c[0]).normalized()).unwrap_or(Vec2::X) + } } } /// Unit tangent at the end (direction of travel). @@ -115,6 +176,11 @@ impl Seg2 { 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), + Seg2::Bezier { n, ref p } => { + let c = Self::bez(n, p); + let last = c.last().copied().unwrap_or(p[0]); + c.iter().rev().skip(1).find_map(|q| (last - *q).normalized()).unwrap_or(Vec2::X) + } } } pub fn reversed(&self) -> Seg2 { @@ -123,13 +189,19 @@ impl Seg2 { 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 }, + Seg2::Bezier { n, ref p } => { + let mut q = *p; + let k = (n as usize).clamp(1, MAX_BEZIER_DEGREE); + q[..=k].reverse(); + Seg2::Bezier { n, p: q } + } } } 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()), + Seg2::Cubic { .. } | Seg2::Conic { .. } | Seg2::Bezier { .. } => self.integrate(|s, t| s.derivative(t).len()), } } /// The segment split at parameter `t` in (0, 1) (same kind, exact; conic halves are @@ -165,6 +237,11 @@ impl Seg2 { Seg2::Conic { a: m, apex: pt(q1), b, w: q1.1 / r.1.max(1e-300).sqrt() }, ) } + Seg2::Bezier { n, ref p } => { + let (_, l, r) = casteljau(Self::bez(n, p), t); + let mk = |c: &[Vec2]| bezier(c).unwrap_or(*self); + (mk(&l), mk(&r)) + } } } /// Parameter of the point of the segment closest to `p` (sampling, then Newton). @@ -213,6 +290,7 @@ impl Seg2 { 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 }) } + Seg2::Bezier { .. } => self.split_at(0.5), } } /// ∫₀¹ f(t) dt by composite Gauss–Legendre (16 pieces of 5 points). @@ -236,7 +314,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))), + Seg2::Cubic { .. } | Seg2::Conic { .. } | Seg2::Bezier { .. } => 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. @@ -260,6 +338,13 @@ impl Seg2 { 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::Bezier { n, ref p } => { + // Second derivative ≤ n(n−1) × the largest second difference of the controls. + let c = Self::bez(n, p); + let d2 = c.windows(3).map(|w| (w[0] - w[1] * 2.0 + w[2]).len()).fold(0.0, f64::max) * (n as f64) * (n as f64 - 1.0); + let k = ((d2 / (8.0 * tol.max(1e-9))).sqrt().ceil() as usize).clamp(2, 4096); + (0..=k).map(|i| self.point_at(i as f64 / k as f64)).collect() + } } } } diff --git a/crates/kernel/src/build.rs b/crates/kernel/src/build.rs index 998e449f..2405dff9 100644 --- a/crates/kernel/src/build.rs +++ b/crates/kernel/src/build.rs @@ -41,7 +41,7 @@ pub(crate) fn wire(plane: &Plane, lp: &Loop2) -> Result { builder::circle_arc(a, b, p3(plane.to_world(s.mid()))) } } - Seg2::Cubic { .. } | Seg2::Conic { .. } => { + Seg2::Cubic { .. } | Seg2::Conic { .. } | Seg2::Bezier { .. } => { crate::freeform::edge(plane, s, a, b).ok_or_else(|| KernelError::Invalid("profile curve".into()))? } }; @@ -71,6 +71,7 @@ fn region_ok(r: &Region2) -> Result<()> { } 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, + Seg2::Bezier { n, p } => (4..=solvecraft_geom::MAX_BEZIER_DEGREE as u8).contains(&n) && p.iter().all(|q| q.is_finite()), }; if !ok { return Err(KernelError::Invalid("non-finite profile".into())); @@ -417,7 +418,7 @@ 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) { + Seg2::Cubic { .. } | Seg2::Conic { .. } | Seg2::Bezier { .. } => match crate::freeform::edge(plane, s, a, b) { Some(e) => e, None => builder::line(a, b), }, diff --git a/crates/kernel/src/freeform.rs b/crates/kernel/src/freeform.rs index 53c651ce..bb6555a1 100644 --- a/crates/kernel/src/freeform.rs +++ b/crates/kernel/src/freeform.rs @@ -1,4 +1,4 @@ -//! Free-form profile edges: cubic Béziers as B-spline curves and conics as rational quadratic +//! Free-form profile edges: Béziers (cubic and higher) as B-spline curves and conics as rational quadratic //! NURBS curves (exact ellipses, parabolas and hyperbolas from sketches). use solvecraft_geom::{Plane, Seg2}; @@ -20,6 +20,10 @@ pub(crate) fn edge(plane: &Plane, s: &Seg2, a: &mt::Vertex, b: &mt::Vertex) -> O 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::Bezier { n, ref p } => { + let inner = p.get(1..(n as usize).clamp(1, solvecraft_geom::MAX_BEZIER_DEGREE))?; + Some(mt::builder::bezier(a, b, inner.iter().map(|q| p3(plane.to_world(*q))).collect())) + } Seg2::Line { .. } | Seg2::Arc { .. } => None, } } diff --git a/crates/sketch/src/model.rs b/crates/sketch/src/model.rs index 602610bc..844bc999 100644 --- a/crates/sketch/src/model.rs +++ b/crates/sketch/src/model.rs @@ -904,20 +904,10 @@ impl Sketch { 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 => { + CurveKind::Spline { pts, control: true, degree } if (*degree as usize) <= solvecraft_geom::MAX_BEZIER_DEGREE => { 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(), - ) + Some(crate::curves::bspline_beziers(&p, k).iter().filter_map(|c| solvecraft_geom::bezier(c)).collect()) } _ => None, } diff --git a/crates/sketch/src/profiles.rs b/crates/sketch/src/profiles.rs index b19c35e5..a2ce3596 100644 --- a/crates/sketch/src/profiles.rs +++ b/crates/sketch/src/profiles.rs @@ -72,7 +72,10 @@ 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. +/// Put exact free-form segments back. A loop that runs over all the pieces of one, in order, +/// gets the whole segment; one that runs over part of it (cut by a crossing) gets the exact +/// sub-segment, its ends moved to the true crossing with the neighbouring segment. Where that +/// crossing cannot be found the loop keeps the polyline pieces. 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) { @@ -81,7 +84,7 @@ fn restore_exact(segs: Vec, tags: &[Option<(Piece, bool)>], exact: &std::c 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 runs: Vec = Vec::new(); let mut i = 0; while i < n { let k0 = key(r + i); @@ -90,21 +93,205 @@ fn restore_exact(segs: Vec, tags: &[Option<(Piece, bool)>], exact: &std::c 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 plain: Vec = run.iter().filter_map(|x| segs.get(*x).copied()).collect(); + let full = k0.and_then(|(ci, k, fwd)| exact.get(&ci).and_then(|e| e.get(k)).map(|s| if fwd { *s } else { s.reversed() })); + let whole = match (full, tags.get(run[0]).copied().flatten()) { + (Some(_), Some(((_, _, _, cnt), fwd))) 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 } + order == want } - _ => None, + _ => false, }; - match whole { - Some(s) => out.push(s), - None => out.extend(run.iter().filter_map(|x| segs.get(*x).copied())), - } + runs.push(match full { + Some(seg) if whole => Run::Exact { seg, t0: 0.0, t1: 1.0, fixed: (true, true), plain }, + Some(seg) => { + let (a, b) = (plain.first().map(Seg2::start).unwrap_or_default(), plain.last().map(Seg2::end).unwrap_or_default()); + let (t0, t1) = (seg.closest_param(a), seg.closest_param(b)); + // Ends that are the segment's own ends need no search. + let snap = |t: f64, p: Vec2, at: f64| (t - at).abs() < 1e-6 && seg.point_at(at).dist(p) < 1e-6; + let fixed = (snap(t0, a, 0.0), snap(t1, b, 1.0)); + let (t0, t1) = (if fixed.0 { 0.0 } else { t0 }, if fixed.1 { 1.0 } else { t1 }); + if t1 - t0 > 1e-9 { Run::Exact { seg, t0, t1, fixed, plain } } else { Run::Plain(plain) } + } + None => Run::Plain(plain), + }); i = j; } - out + // Junctions next to a cut exact run: find the true crossing. + let m = runs.len(); + for a in 0..m { + let b = (a + 1) % m; + let need = match (runs.get(a), runs.get(b)) { + (Some(Run::Exact { fixed, .. }), _) if !fixed.1 => true, + (_, Some(Run::Exact { fixed, .. })) if !fixed.0 => true, + _ => false, + }; + if !need { + continue; + } + let (Some(ra), Some(rb)) = (runs.get(a), runs.get(b)) else { continue }; + let p0 = ra.end(); + let crossing = refine_crossing(ra, rb, p0); + match crossing { + Some((p, ta, tb)) if p.dist(p0) < 1e-2 => { + if let Some(x) = runs.get_mut(a) { + x.set_end(p, ta); + } + if let Some(x) = runs.get_mut(b) { + x.set_start(p, tb); + } + } + _ => { + for k in [a, b] { + if let Some(x) = runs.get_mut(k) { + x.give_up(); + } + } + } + } + } + runs.into_iter().flat_map(Run::into_segs).collect() +} + +/// A stretch of a profile loop: plain segments, or (part of) an exact free-form segment. +enum Run { + Plain(Vec), + /// `seg` (in loop direction) between parameters `t0` and `t1`; `fixed` ends are exact; + /// `plain` is the polyline fallback. + Exact { + seg: Seg2, + t0: f64, + t1: f64, + fixed: (bool, bool), + plain: Vec, + }, +} + +impl Run { + fn end(&self) -> Vec2 { + match self { + Run::Plain(v) => v.last().map(Seg2::end).unwrap_or_default(), + Run::Exact { seg, t1, .. } => seg.point_at(*t1), + } + } + /// The segment at the run's end and the parameter of the end on it. + fn end_geom(&self) -> Option<(Seg2, f64)> { + match self { + Run::Plain(v) => v.last().map(|s| (*s, 1.0)), + Run::Exact { seg, t1, .. } => Some((*seg, *t1)), + } + } + /// The segment at the run's start and the parameter of the start on it. + fn start_geom(&self) -> Option<(Seg2, f64)> { + match self { + Run::Plain(v) => v.first().map(|s| (*s, 0.0)), + Run::Exact { seg, t0, .. } => Some((*seg, *t0)), + } + } + fn set_end(&mut self, p: Vec2, t: f64) { + match self { + Run::Plain(v) => { + if let Some(s) = v.last_mut() { + *s = with_end(*s, p); + } + } + Run::Exact { t1, .. } => *t1 = t, + } + } + fn set_start(&mut self, p: Vec2, t: f64) { + match self { + Run::Plain(v) => { + if let Some(s) = v.first_mut() { + *s = with_end(s.reversed(), p).reversed(); + } + } + Run::Exact { t0, .. } => *t0 = t, + } + } + fn give_up(&mut self) { + if let Run::Exact { fixed: (false, _) | (_, false), plain, .. } = self { + *self = Run::Plain(std::mem::take(plain)); + } + } + fn into_segs(self) -> Vec { + match self { + Run::Plain(v) => v, + Run::Exact { seg, t0, t1, .. } => { + if t0 <= 1e-12 && t1 >= 1.0 - 1e-12 { + return vec![seg]; + } + let s = if t1 < 1.0 - 1e-12 { seg.split_at(t1).0 } else { seg }; + if t0 <= 1e-12 { + return vec![s]; + } + // Conic halves are reparametrised: find t0's point again on the left part. + let u = local_param(&s, seg.point_at(t0), if t1 > 1e-12 { t0 / t1 } else { 0.0 }); + vec![s.split_at(u.clamp(1e-9, 1.0 - 1e-9)).1] + } + } + } +} + +/// The segment with its end moved to `p` (lines; arcs keep their circle). +fn with_end(s: Seg2, p: Vec2) -> Seg2 { + match s { + Seg2::Line { a, .. } => Seg2::Line { a, b: p }, + Seg2::Arc { center, radius, start, sweep } => { + let target = (p - center).angle(); + let end = start + sweep; + // The same end angle, wound the same way. + let delta = (target - end + std::f64::consts::PI).rem_euclid(TAU) - std::f64::consts::PI; + Seg2::Arc { center, radius, start, sweep: sweep + delta } + } + other => other, + } +} + +/// Where the end of run `a` meets the start of run `b`, near `p`: the point and the parameters +/// on the two runs' segments (Newton on A(ta) = B(tb); lines and arcs extend past their ends). +fn refine_crossing(a: &Run, b: &Run, p: Vec2) -> Option<(Vec2, f64, f64)> { + let (sa, mut ta) = a.end_geom()?; + let (sb, mut tb) = b.start_geom()?; + // Start from the parameters nearest the polyline crossing. + ta = local_param(&sa, p, ta); + tb = local_param(&sb, p, tb); + for _ in 0..50 { + let f = sa.point_at(ta) - sb.point_at(tb); + if f.len() < 1e-11 { + let pt = match (a, b) { + (Run::Exact { .. }, _) | (Run::Plain(_), Run::Plain(_)) => sa.point_at(ta), + (_, Run::Exact { .. }) => sb.point_at(tb), + }; + return (ta.is_finite() && tb.is_finite()).then_some((pt, ta, tb)); + } + let (da, db) = (sa.derivative(ta), -sb.derivative(tb)); + let det = da.cross(db); + if det.abs() < 1e-14 * da.len().max(1e-300) * db.len().max(1e-300) { + return None; + } + // Solve [da db] (dta, dtb)ᵀ = −f. + ta -= f.cross(db) / det; + tb -= da.cross(f) / det; + } + None +} + +/// Parameter of the point of `s` nearest `p`, by Newton from `t`. +fn local_param(s: &Seg2, p: Vec2, mut t: f64) -> f64 { + for _ in 0..30 { + let (q, d) = (s.point_at(t), s.derivative(t)); + let dd = d.dot(d); + if dd < 1e-300 { + break; + } + let step = (p - q).dot(d) / dd; + t += step; + if step.abs() < 1e-15 { + break; + } + } + t } /// A uniform grid over a set of boxes (about one box per cell), for neighbour queries. @@ -161,7 +348,7 @@ 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 { .. } => { + Seg2::Cubic { .. } | Seg2::Conic { .. } | Seg2::Bezier { .. } => { 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) diff --git a/crates/sketch/src/tests.rs b/crates/sketch/src/tests.rs index cd45b950..a7a2ca51 100644 --- a/crates/sketch/src/tests.rs +++ b/crates/sketch/src/tests.rs @@ -407,3 +407,53 @@ fn long_connected_chain_solves() { assert_eq!(rep.dof, 0); assert!(ms < 10_000, "{ms} ms"); } + +/// Free-form curves stay exact in profiles: a degree-5 spline closed by a line, and a fit spline +/// cut in two by a line through it (each side gets an exact sub-curve ending on the line). +#[test] +fn exact_profiles_for_high_degree_and_cut_free_form_curves() { + use solvecraft_geom::Seg2; + let exact_loop = |p: &Profile| p.region.outer.segs.iter().all(|s| !matches!(s, Seg2::Line { .. }) || s.length() > 1.0); + let closed = |p: &Profile| { + let s = &p.region.outer.segs; + s.iter().zip(s.iter().cycle().skip(1)).all(|(a, b)| a.end().dist(b.start()) < 1e-9) + }; + // Degree 5 over six points, closed by a line. + let mut sk = Sketch::default(); + let ids: Vec = [v(0.0, 0.0), v(4.0, 10.0), v(10.0, 12.0), v(16.0, 12.0), v(22.0, 10.0), v(26.0, 0.0)] + .iter() + .map(|p| sk.add_point(*p, None).unwrap()) + .collect(); + sk.add_curve(CurveKind::Spline { pts: ids.clone(), control: true, degree: 5 }, None).unwrap(); + sk.add_line_pts(ids[5], ids[0], None).unwrap(); + let pr = find_profiles(&sk); + assert_eq!(pr.len(), 1); + assert!(pr[0].region.outer.segs.iter().any(|s| matches!(s, Seg2::Bezier { n: 5, .. })), "{:?}", pr[0].region.outer.segs); + assert!(closed(&pr[0])); + // A fit spline arch over a base line, cut by a vertical line at x = 10. + let mut sk = Sketch::default(); + let ids: Vec = [v(0.0, 0.0), v(6.0, 9.0), v(14.0, 9.0), v(20.0, 0.0)].iter().map(|p| sk.add_point(*p, None).unwrap()).collect(); + sk.add_curve(CurveKind::Spline { pts: ids.clone(), control: false, degree: 3 }, None).unwrap(); + sk.add_line_pts(ids[3], ids[0], None).unwrap(); + let (a, b) = (sk.add_point(v(10.0, -5.0), None).unwrap(), sk.add_point(v(10.0, 20.0), None).unwrap()); + sk.add_line_pts(a, b, None).unwrap(); + let pr = find_profiles(&sk); + assert_eq!(pr.len(), 2, "{pr:?}"); + let whole: f64 = pr.iter().map(|p| p.area).sum(); + for p in &pr { + assert!(closed(p), "{:?}", p.region.outer.segs); + assert!(exact_loop(p), "polyline pieces left: {:?}", p.region.outer.segs); + // The crossing lies on the line x = 10. + assert!( + p.region + .outer + .segs + .iter() + .any(|s| matches!(s, Seg2::Cubic { .. }) && ((s.start().x - 10.0).abs() < 1e-9 || (s.end().x - 10.0).abs() < 1e-9)) + ); + } + // Same total area as the uncut arch. + sk.curves.pop(); + let uncut: f64 = find_profiles(&sk).iter().map(|p| p.area).sum(); + assert!((whole - uncut).abs() < 1e-6 * uncut, "{whole} vs {uncut}"); +}