Sketch constraints and dimensions: Curvature (G2, ConstraintSmooth) and Polygon constraints; redundant constraints and dimensions are refused as over-constraining (Fusion behaviour), with driven (reference) dimensions that measure instead; arc length and linear diameter (about a centerline) dimensions; sketch.glyphs (constraint glyph anchors, conflicts) and sketch.snap (coincident/midpoint/on-curve/horizontal/vertical inference); solver damping floor keeps weakly seen variables from flying off

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01AS5Tn3xieewUSg3uxrMVdL

crates/sketch/src/curves.rs

@@ -94,38 +94,73 @@ fn hermite(p0: Vec2, p1: Vec2, t0: Vec2, t1: Vec2, s: f64) -> Vec2 {
     p0 * (2.0 * s3 - 3.0 * s2 + 1.0) + t0 * (s3 - 2.0 * s2 + s) + p1 * (-2.0 * s3 + 3.0 * s2) + t1 * (s3 - s2)
 }
 
-/// Clamped uniform B-spline of `degree` over control points `p`, sampled.
-fn bspline_polyline(p: &[Vec2], degree: usize) -> Vec<Vec2> {
+/// Clamped uniform B-spline of `degree` over control points `p` at `t` in [0, 1] (de Boor).
+fn bspline_at(p: &[Vec2], degree: usize, t: f64) -> Vec2 {
     let n = p.len();
     let k = degree.min(n.saturating_sub(1)).max(1);
-    // Knots: k+1 zeros, uniform interior, k+1 ones.
     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 total = SPAN_SAMPLES * spans;
-    (0..=total)
-        .map(|i| {
-            let t = i as f64 / total as f64;
-            // De Boor.
-            let mut s = k;
-            while s + 1 < n && knots.get(s + 1).is_some_and(|x| *x <= t) {
-                s += 1;
-            }
-            let mut d: Vec<Vec2> = (0..=k).map(|j| p.get(j + s - k).copied().unwrap_or_default()).collect();
-            for r in 1..=k {
-                for j in (r..=k).rev() {
-                    let i0 = j + s - k;
-                    let (lo, hi) = (knots.get(i0).copied().unwrap_or(0.0), knots.get(i0 + k + 1 - r).copied().unwrap_or(1.0));
-                    let a = if hi - lo > 1e-15 { (t - lo) / (hi - lo) } else { 0.0 };
-                    d[j] = d[j - 1] * (1.0 - a) + d[j] * a;
-                }
-            }
-            d[k]
-        })
-        .collect()
+    let mut s = k;
+    while s + 1 < n && knots.get(s + 1).is_some_and(|x| *x <= t) {
+        s += 1;
+    }
+    let mut d: Vec<Vec2> = (0..=k).map(|j| p.get(j + s - k).copied().unwrap_or_default()).collect();
+    for r in 1..=k {
+        for j in (r..=k).rev() {
+            let i0 = j + s - k;
+            let (lo, hi) = (knots.get(i0).copied().unwrap_or(0.0), knots.get(i0 + k + 1 - r).copied().unwrap_or(1.0));
+            let a = if hi - lo > 1e-15 { (t - lo) / (hi - lo) } else { 0.0 };
+            d[j] = d[j - 1] * (1.0 - a) + d[j] * a;
+        }
+    }
+    d[k]
+}
+
+/// Clamped uniform B-spline of `degree` over control points `p`, sampled.
+fn bspline_polyline(p: &[Vec2], degree: usize) -> Vec<Vec2> {
+    let n = p.len();
+    let k = degree.min(n.saturating_sub(1)).max(1);
+    let total = SPAN_SAMPLES * (n - k);
+    (0..=total).map(|i| bspline_at(p, degree, i as f64 / total as f64)).collect()
+}
+
+/// Point of a spline at global parameter `t` in [0, 1].
+pub fn spline_point(p: &[Vec2], control: bool, degree: u8, t: f64) -> Vec2 {
+    let t = t.clamp(0.0, 1.0);
+    if p.len() < 2 {
+        return p.first().copied().unwrap_or_default();
+    }
+    if control {
+        return bspline_at(p, degree.clamp(1, 7) as usize, t);
+    }
+    let spans = fit_spans(p);
+    let n = spans.len().max(1);
+    let x = t * n as f64;
+    let i = (x.floor() as usize).min(n - 1);
+    match spans.get(i) {
+        Some((p0, p1, t0, t1)) => hermite(*p0, *p1, *t0, *t1, x - i as f64),
+        None => p[0],
+    }
+}
+
+/// Curvature vector (towards the centre of curvature, magnitude 1/radius) of a parametric
+/// curve at an end (`at_end`), by one-sided second-order differences.
+pub fn end_curvature(f: &dyn Fn(f64) -> Vec2, at_end: bool) -> Vec2 {
+    let h = 1e-3;
+    let (t0, sgn) = if at_end { (1.0, -1.0) } else { (0.0, 1.0) };
+    let (f0, f1, f2, f3) = (f(t0), f(t0 + sgn * h), f(t0 + sgn * 2.0 * h), f(t0 + sgn * 3.0 * h));
+    let d1 = (f0 * -3.0 + f1 * 4.0 - f2) / (2.0 * h) * sgn;
+    let d2 = (f0 * 2.0 - f1 * 5.0 + f2 * 4.0 - f3) / (h * h);
+    let l2 = d1.len2();
+    if l2 < 1e-24 {
+        return Vec2::ZERO;
+    }
+    let t = d1 / l2.sqrt();
+    (d2 - t * t.dot(d2)) / l2
 }
 
 /// Polyline of a spline through (`control == false`) or over (`control`) the points.
@@ -161,20 +196,20 @@ pub fn spline_end_tangent(p: &[Vec2], control: bool, at_end: bool) -> Option<Vec
     if at_end { spans.last().and_then(|s| s.3.normalized()) } else { spans.first().and_then(|s| s.2.normalized()) }
 }
 
-/// Distance from `q` to a polyline.
-pub fn polyline_dist(poly: &[Vec2], q: Vec2) -> f64 {
-    let mut best = f64::INFINITY;
+/// Signed distance from `q` to a polyline (positive on the left of the nearest segment).
+pub fn polyline_signed_dist(poly: &[Vec2], q: Vec2) -> f64 {
+    let mut best = (f64::INFINITY, 1.0);
     for w in poly.windows(2) {
         let (a, b) = (w[0], w[1]);
         let d = b - a;
         let l2 = d.len2();
         let t = if l2 > 0.0 { ((q - a).dot(d) / l2).clamp(0.0, 1.0) } else { 0.0 };
-        best = best.min(q.dist(a + d * t));
+        let dist = q.dist(a + d * t);
+        if dist < best.0 {
+            best = (dist, if d.cross(q - a) >= 0.0 { 1.0 } else { -1.0 });
+        }
     }
-    if poly.len() == 1 {
-        best = poly[0].dist(q);
-    }
-    best
+    best.0 * best.1
 }
 
 #[cfg(test)]
@@ -195,6 +230,20 @@ mod tests {
         assert!((t.x - t.y).abs() < 1e-9);
     }
 
+    #[test]
+    fn end_curvature_of_a_circle_quarter() {
+        let w = std::f64::consts::FRAC_1_SQRT_2;
+        let rho = w / (1.0 + w);
+        let (a, x, b) = (Vec2::new(10.0, 0.0), Vec2::new(10.0, 10.0), Vec2::new(0.0, 10.0));
+        let k = end_curvature(&|t| conic_point(a, x, b, rho, t), false);
+        assert!(k.dist(Vec2::new(-0.1, 0.0)) < 1e-5, "{k:?}");
+        let k1 = end_curvature(&|t| conic_point(a, x, b, rho, t), true);
+        assert!(k1.dist(Vec2::new(0.0, -0.1)) < 1e-5, "{k1:?}");
+        // Fit splines have natural (straight) ends.
+        let p = [Vec2::new(0.0, 0.0), Vec2::new(10.0, 5.0), Vec2::new(20.0, 0.0)];
+        assert!(end_curvature(&|t| spline_point(&p, false, 3, t), false).len() < 1e-3);
+    }
+
     #[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)];

crates/sketch/src/lib.rs

@@ -17,7 +17,7 @@ mod shape;
 mod solver;
 mod text;
 
-pub use curves::{MAX_SPLINE_POINTS, conic_point, ellipse_point, spline_end_tangent, spline_polyline};
+pub use curves::{MAX_SPLINE_POINTS, conic_point, ellipse_point, end_curvature, 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};

crates/sketch/src/model.rs

@@ -248,6 +248,23 @@ pub enum ConstraintKind {
         #[serde(default, skip_serializing_if = "std::ops::Not::not")]
         flip: bool,
     },
+    /// Curvature continuity (G2) where curves `a` and `b` share an end point: tangent and
+    /// equal curvature.
+    Smooth {
+        a: usize,
+        b: usize,
+    },
+    /// Length along an arc.
+    ArcLength {
+        c: usize,
+        value: f64,
+    },
+    /// Twice the distance from point `p` to line `l` (a diameter about a centerline).
+    LinearDiameter {
+        p: usize,
+        l: usize,
+        value: f64,
+    },
 }
 
 impl ConstraintKind {
@@ -264,6 +281,8 @@ impl ConstraintKind {
             | Length { value, .. }
             | Radius { value, .. }
             | Diameter { value, .. }
+            | ArcLength { value, .. }
+            | LinearDiameter { value, .. }
             | Angle { value, .. } => Some(value),
             _ => None,
         }
@@ -278,6 +297,8 @@ impl ConstraintKind {
             | Length { value, .. }
             | Radius { value, .. }
             | Diameter { value, .. }
+            | ArcLength { value, .. }
+            | LinearDiameter { value, .. }
             | Angle { value, .. } => *value = v,
             _ => {}
         }
@@ -305,6 +326,9 @@ impl ConstraintKind {
             Radius { .. } => "RadialDimension",
             Diameter { .. } => "DiameterDimension",
             Angle { .. } => "AngularDimension",
+            Smooth { .. } => "Smooth",
+            ArcLength { .. } => "ArcLengthDimension",
+            LinearDiameter { .. } => "LinearDiameterDimension",
         }
     }
     fn indices(&self) -> (Vec<usize>, Vec<usize>) {
@@ -322,6 +346,9 @@ impl ConstraintKind {
             Midpoint { p, l } | PointLineDistance { p, l, .. } => (vec![p], vec![l]),
             Symmetric { p, q, l } => (vec![p, q], vec![l]),
             Fix { p } => (vec![p], vec![]),
+            Smooth { a, b } => (vec![], vec![a, b]),
+            ArcLength { c, .. } => (vec![], vec![c]),
+            LinearDiameter { p, l, .. } => (vec![p], vec![l]),
         }
     }
     fn remap(&mut self, pmap: &dyn Fn(usize) -> Option<usize>, cmap: &dyn Fn(usize) -> Option<usize>) -> bool {
@@ -369,6 +396,15 @@ impl ConstraintKind {
                 fc(l);
             }
             Fix { p } => fp(p),
+            Smooth { a, b } => {
+                fc(a);
+                fc(b);
+            }
+            ArcLength { c, .. } => fc(c),
+            LinearDiameter { p, l, .. } => {
+                fp(p);
+                fc(l);
+            }
         }
         ok && ok2
     }
@@ -382,6 +418,9 @@ pub struct Constraint {
     /// For dimensions: the model parameter that drives the value (e.g. `d1`).
     #[serde(default, skip_serializing_if = "Option::is_none")]
     pub param: Option<String>,
+    /// A driven (reference) dimension: it measures, it does not constrain.
+    #[serde(default, skip_serializing_if = "std::ops::Not::not")]
+    pub driven: bool,
 }
 
 /// A 2D sketch. Point 0 is always the fixed sketch origin (id `origin`).
@@ -553,7 +592,7 @@ impl Sketch {
         }
         self.validate(&kind)?;
         let id = self.fresh("k");
-        self.constraints.push(Constraint { id: id.clone(), kind, param });
+        self.constraints.push(Constraint { id: id.clone(), kind, param, driven: false });
         Ok(id)
     }
 
@@ -577,8 +616,27 @@ impl Sketch {
         let need_line = |i: usize| if is_line(i) { Ok(()) } else { Err(SketchError::WrongKind(name(i), "needs a line".into())) };
         let need_round = |i: usize| if is_round(i) { Ok(()) } else { Err(SketchError::WrongKind(name(i), "needs a circle or arc".into())) };
         match *k {
-            Horizontal { l } | Vertical { l } | Length { l, .. } | Midpoint { l, .. } | Symmetric { l, .. } | PointLineDistance { l, .. } => {
-                need_line(l)?
+            Horizontal { l }
+            | Vertical { l }
+            | Length { l, .. }
+            | Midpoint { l, .. }
+            | Symmetric { l, .. }
+            | PointLineDistance { l, .. }
+            | LinearDiameter { l, .. } => need_line(l)?,
+            ArcLength { c, .. } => {
+                if !matches!(self.curves.get(c).map(|c| &c.kind), Some(CurveKind::Arc { .. })) {
+                    return Err(SketchError::WrongKind(name(c), "arc length needs an arc".into()));
+                }
+            }
+            Smooth { a, b } => {
+                let ends = |i: usize| self.curves.get(i).and_then(|c| c.kind.ends());
+                let shared = match (ends(a), ends(b)) {
+                    (Some((a0, a1)), Some((b0, b1))) => a0 == b0 || a0 == b1 || a1 == b0 || a1 == b1,
+                    _ => false,
+                };
+                if !shared || a == b {
+                    return Err(SketchError::Invalid("curvature continuity needs two curves sharing an end point".into()));
+                }
             }
             Parallel { a, b } | Perpendicular { a, b } | Collinear { a, b } | Angle { a, b, .. } => {
                 need_line(a)?;

crates/sketch/src/solver.rs

@@ -157,6 +157,43 @@ impl State<'_> {
             _ => Vec::new(),
         }
     }
+    fn ends(&self, c: usize) -> Option<(usize, usize)> {
+        self.sk.curves.get(c)?.kind.ends()
+    }
+    /// Length along an arc (counter-clockwise from its start).
+    fn arc_length(&self, c: usize) -> f64 {
+        match self.sk.curves.get(c).map(|c| &c.kind) {
+            Some(CurveKind::Arc { c: ci, a, b }) => {
+                let (cc, pa, pb) = (self.p(*ci), self.p(*a), self.p(*b));
+                let mut sw = (pb - cc).angle() - (pa - cc).angle();
+                while sw <= 1e-12 {
+                    sw += std::f64::consts::TAU;
+                }
+                cc.dist(pa) * sw
+            }
+            _ => 0.0,
+        }
+    }
+    /// Curvature vector of a curve at its point `p` (towards the centre of curvature).
+    fn curvature_at(&self, c: usize, p: usize) -> Vec2 {
+        let Some(kind) = self.sk.curves.get(c).map(|c| &c.kind) else { return Vec2::ZERO };
+        match kind {
+            CurveKind::Line { .. } | CurveKind::Ellipse { .. } => Vec2::ZERO,
+            CurveKind::Circle { .. } | CurveKind::Arc { .. } => {
+                let (cc, r) = self.round(c);
+                (cc - self.p(p)) / (r * r).max(1e-24)
+            }
+            CurveKind::Spline { pts, control, degree } => {
+                let q: Vec<Vec2> = pts.iter().map(|i| self.p(*i)).collect();
+                let at_end = pts.last() == Some(&p) && pts.first() != Some(&p);
+                crate::curves::end_curvature(&|t| crate::curves::spline_point(&q, *control, *degree, t), at_end)
+            }
+            CurveKind::Conic { a, b, apex, rho } => {
+                let (pa, pb, px) = (self.p(*a), self.p(*b), self.p(*apex));
+                crate::curves::end_curvature(&|t| crate::curves::conic_point(pa, px, pb, *rho, t), p == *b)
+            }
+        }
+    }
     /// Unit tangent of any curve at its point `p` (an end point for open free-form curves).
     fn tangent_at(&self, c: usize, p: usize) -> Option<Vec2> {
         match self.sk.curves.get(c).map(|c| &c.kind)? {
@@ -214,7 +251,7 @@ fn residuals(s: &State, k: &ConstraintKind, out: &mut Vec<f64>) {
         PointOnCurve { p, c } => {
             let pp = s.p(p);
             if s.is_freeform(c) {
-                out.push(crate::curves::polyline_dist(&s.poly(c), pp));
+                out.push(crate::curves::polyline_signed_dist(&s.poly(c), pp));
             } else if s.is_line(c) {
                 let (a, b) = s.line(c);
                 out.push(line_dist(pp, a, b));
@@ -253,15 +290,17 @@ fn residuals(s: &State, k: &ConstraintKind, out: &mut Vec<f64>) {
                 // Without a shared end point: the curves touch (closest approach zero).
                 let (f, o) = if s.is_freeform(a) { (a, b) } else { (b, a) };
                 let poly = s.poly(f);
+                // Signed gap at the closest approach (smooth through zero, unlike its size).
+                let closest = |v: &mut dyn Iterator<Item = f64>| v.fold(f64::INFINITY, |m: f64, x| if x.abs() < m.abs() { x } else { m });
                 let d = if s.is_line(o) {
                     let (p0, p1) = s.line(o);
-                    poly.iter().map(|q| line_dist(*q, p0, p1).abs()).fold(f64::INFINITY, f64::min)
+                    closest(&mut poly.iter().map(|q| line_dist(*q, p0, p1)))
                 } else if s.is_freeform(o) {
                     let other = s.poly(o);
-                    poly.iter().map(|q| crate::curves::polyline_dist(&other, *q)).fold(f64::INFINITY, f64::min)
+                    closest(&mut poly.iter().map(|q| crate::curves::polyline_signed_dist(&other, *q)))
                 } else {
                     let (cc, r) = s.round(o);
-                    poly.iter().map(|q| (q.dist(cc) - r).abs()).fold(f64::INFINITY, f64::min)
+                    closest(&mut poly.iter().map(|q| q.dist(cc) - r))
                 };
                 out.push(if d.is_finite() { d } else { 0.0 });
             } else if la || lb {
@@ -327,9 +366,63 @@ fn residuals(s: &State, k: &ConstraintKind, out: &mut Vec<f64>) {
             let target = if flip { value + std::f64::consts::PI } else { value };
             out.push(wrap_angle(ang - target) * scale.max(1e-6));
         }
+        Smooth { a, b } => {
+            let shared = match (s.ends(a), s.ends(b)) {
+                (Some((a0, a1)), Some((b0, b1))) => [a0, a1].into_iter().find(|x| *x == b0 || *x == b1),
+                _ => None,
+            };
+            let Some(p) = shared else {
+                out.extend([0.0, 0.0, 0.0]);
+                return;
+            };
+            match (s.tangent_at(a, p), s.tangent_at(b, p)) {
+                (Some(u), Some(v)) => out.push(u.cross(v) * 10.0),
+                _ => out.push(0.0),
+            }
+            let k = s.curvature_at(a, p) - s.curvature_at(b, p);
+            out.extend([k.x * 100.0, k.y * 100.0]);
+        }
+        ArcLength { c, value } => out.push(s.arc_length(c) - value),
+        LinearDiameter { p, l, value } => {
+            let (a, b) = s.line(l);
+            out.push(2.0 * line_dist(s.p(p), a, b).abs() - value);
+        }
     }
 }
 
+/// Measured value of a dimension (for driven dimensions).
+fn measure(s: &State, k: &ConstraintKind) -> Option<f64> {
+    use ConstraintKind::*;
+    Some(match *k {
+        Distance { p, q, .. } => s.p(p).dist(s.p(q)),
+        DistanceX { p, q, .. } => (s.p(q).x - s.p(p).x).abs(),
+        DistanceY { p, q, .. } => (s.p(q).y - s.p(p).y).abs(),
+        PointLineDistance { p, l, .. } => {
+            let (a, b) = s.line(l);
+            line_dist(s.p(p), a, b).abs()
+        }
+        Length { l, .. } => {
+            let (a, b) = s.line(l);
+            a.dist(b)
+        }
+        Radius { c, .. } => s.round(c).1,
+        Diameter { c, .. } => 2.0 * s.round(c).1,
+        Angle { a, b, flip, .. } => {
+            let (a0, a1) = s.line(a);
+            let (b0, b1) = s.line(b);
+            let (da, db) = (a1 - a0, b1 - b0);
+            let ang = da.cross(db).atan2(da.dot(db)) - if flip { std::f64::consts::PI } else { 0.0 };
+            ang.rem_euclid(std::f64::consts::TAU)
+        }
+        ArcLength { c, .. } => s.arc_length(c),
+        LinearDiameter { p, l, .. } => {
+            let (a, b) = s.line(l);
+            2.0 * line_dist(s.p(p), a, b).abs()
+        }
+        _ => return None,
+    })
+}
+
 /// One residual block: implicit arc radius equality, a user constraint, or a tangency at a known
 /// touch point (better conditioned than the distance form when an end point is on the curve).
 enum Block<'a> {
@@ -445,6 +538,15 @@ fn block_vars(sk: &Sketch, lay: &Layout, b: &Block) -> Vec<usize> {
                 curve(&mut v, l);
             }
             Fix { .. } => {}
+            Smooth { a, b } => {
+                curve(&mut v, a);
+                curve(&mut v, b);
+            }
+            ArcLength { c, .. } => curve(&mut v, c),
+            LinearDiameter { p, l, .. } => {
+                pt(&mut v, p);
+                curve(&mut v, l);
+            }
         },
         Block::TangentAt(_, k, p) => {
             pt(&mut v, *p);
@@ -553,11 +655,13 @@ fn lm(sk: &Sketch, lay: &Layout, blocks: &[(Block, Vec<usize>)], vars: &[usize],
             }
         }
         let mut improved = false;
+        // Damping floor: variables the residuals barely see must not take huge steps.
+        let floor = (0..n).map(|i| a.get(i, i)).fold(0.0_f64, f64::max) * 1e-4;
         for _ in 0..30 {
             let mut m = a.clone();
             for i in 0..n {
                 let d = a.get(i, i);
-                m.add(i, i, lambda * (d + 1e-9) + 1e-12);
+                m.add(i, i, lambda * (d.max(floor) + 1e-9) + 1e-12);
             }
             let rhs: Vec<f64> = g.iter().map(|v| -v).collect();
             let Some(step) = lin_solve(m, rhs) else {
@@ -625,6 +729,9 @@ pub fn solve(sk: &mut Sketch) -> SolveReport {
         }
     }
     for (i, c) in sk.constraints.iter().enumerate() {
+        if c.driven {
+            continue;
+        }
         let bl = match touch_point(sk, &c.kind) {
             Some(p) => Block::TangentAt(i, &c.kind, p),
             None => Block::User(i, &c.kind),
@@ -725,6 +832,17 @@ pub fn solve(sk: &mut Sketch) -> SolveReport {
         Some(v) => vdet(Some(v)) && vdet(Some(v + 1)),
         None => true,
     };
+    // Driven dimensions measure the solved geometry.
+    let fixed_lay = Layout { pvar: vec![None; sk.points.len()], rvar: vec![None; sk.curves.len()], n: 0 };
+    let measured: Vec<Option<f64>> =
+        sk.constraints.iter().map(|c| if c.driven { measure(&State { sk, lay: &fixed_lay, x: &[] }, &c.kind) } else { None }).collect();
+    for (c, m) in sk.constraints.iter_mut().zip(measured) {
+        if let Some(v) = m
+            && v.is_finite()
+        {
+            c.kind.set_value(v);
+        }
+    }
     let point_determined: Vec<bool> = (0..sk.points.len()).map(pdet).collect();
     let curve_determined =
         sk.curves.iter().enumerate().map(|(i, c)| c.kind.point_ids().into_iter().all(pdet) && vdet(lay.rvar.get(i).copied().flatten())).collect();