kernel: iso_curves, a face's isoparametric curve exactly, trimmed by its loops

iso_curves gives the isoparametric curve of a body's face through a point, along
its first or second parameter, as record curves: on a B-spline surface the
surface's own row at the point's other parameter (each control point its
column evaluated there, homogeneous, so exact for rational nets), on an
analytic surface its line or circle. It is cut where it crosses the face's
loops: each crossing solved on the loop's pcurve itself (bracketed on samples,
then bisected on the pcurve), the pieces between kept where their middle lies
inside the face (even-odd over the loops), neighbouring inside pieces joined.
A revolved surface record is not supported yet.

lib.rs (coordinator): export point_path::iso_curves.
(not built: base in progress)

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

crates/kernel/src/point_path.rs

@@ -171,6 +171,217 @@ pub fn project_onto_face(
     crate::smlib_sweeps::project_onto_face(body, face_at, plane, segs, direction)
 }
 
+/// Samples per knot span (and at least this many per pcurve) where the iso curves look for the
+/// crossings of a face's loops and classify the face's inside.
+const ISO_SAMPLES: usize = 64;
+
+/// The isoparametric curve of face `face` of `body` through the point `at` on it, the pieces of
+/// it inside the face, each as its exact record curve on its interval: along the surface's first
+/// parameter (`along_u`; around a cylinder, cone, sphere or torus) or its second (along them).
+/// On a B-spline surface the curve is the surface's row at the point's other parameter (each
+/// control point the column through it evaluated there, exactly); on an analytic surface its
+/// line or circle. The pieces are where the curve crosses the face's loops (their pcurves).
+pub fn iso_curves(body: &Body, face: usize, at: Vec3, along_u: bool) -> Result<Vec<(brep::Curve, [f64; 2])>> {
+    use solidcraft_geom::brep::eval::SurfaceEval;
+    let record = body.record();
+    let f = record.faces.get(face).ok_or_else(|| KernelError::Invalid(format!("the body has no face {face}")))?;
+    let eval = SurfaceEval::new(&f.surface).map_err(|e| KernelError::Failed(format!("iso curve: {e}")))?;
+    let [[u0, u1], [v0, v1]] = f.domain;
+    let size = body.size().max(1e-9);
+    let mid = [0.5 * (u0 + u1), 0.5 * (v0 + v1)];
+    let [mut pu, mut pv] =
+        eval.project(at, Some(mid), size * 1e-3 + 1e-6).ok_or_else(|| KernelError::Invalid("iso curve: the point is not on the face".into()))?;
+    // Into the face's box, by whole periods.
+    let period = eval.domain().period;
+    for (x, (lo, hi), p) in [(&mut pu, (u0, u1), period[0]), (&mut pv, (v0, v1), period[1])] {
+        if let Some(p) = p.filter(|p| *p > 0.0) {
+            *x = lo + (*x - lo).rem_euclid(p);
+            if *x > hi + 1e-12 {
+                *x -= p;
+            }
+        }
+    }
+    // The iso line in (u, v): `fixed` held, the other running over the box.
+    let (fixed, (lo, hi)) = if along_u { (pv, (u0, u1)) } else { (pu, (v0, v1)) };
+    let at_t = |t: f64| if along_u { [t, fixed] } else { [fixed, t] };
+    let pcurves: Vec<Vec<&brep::Curve2>> = f
+        .loops
+        .iter()
+        .filter_map(|l| match l {
+            brep::Loop::Edges { edges } => Some(edges.iter().map(|e| &e.pcurve).collect()),
+            brep::Loop::Vertex { .. } => None,
+        })
+        .collect();
+    // Where the line crosses the loops' pcurves (along it), each crossing solved on the pcurve
+    // itself, with the box's ends.
+    let across = |q: [f64; 2]| if along_u { q[1] - fixed } else { q[0] - fixed };
+    let along = |q: [f64; 2]| if along_u { q[0] } else { q[1] };
+    let mut cuts = vec![lo, hi];
+    for c in pcurves.iter().flatten() {
+        let samples = pcurve_samples(c);
+        for w in samples.windows(2) {
+            let ((ta, qa), (tb, qb)) = (w[0], w[1]);
+            let (ka, kb) = (across(qa), across(qb));
+            if (ka < 0.0) == (kb < 0.0) || ka == kb {
+                continue;
+            }
+            let (mut x, mut y) = (ta, tb);
+            for _ in 0..80 {
+                let m = 0.5 * (x + y);
+                match c.eval(m).map(across) {
+                    Some(k) if (k < 0.0) == (ka < 0.0) => x = m,
+                    Some(_) => y = m,
+                    None => break,
+                }
+            }
+            if let Some(t) = c.eval(0.5 * (x + y)).map(along).filter(|t| *t > lo && *t < hi) {
+                cuts.push(t);
+            }
+        }
+    }
+    let loops: Vec<Vec<[f64; 2]>> = pcurves.iter().map(|l| l.iter().flat_map(|c| pcurve_samples(c)).map(|(_, q)| q).collect()).collect();
+    cuts.sort_by(f64::total_cmp);
+    cuts.dedup_by(|a, b| (*a - *b).abs() <= 1e-9 * (1.0 + b.abs()));
+    let inside = |q: [f64; 2]| loops.iter().filter(|l| crosses_right(q, l)).count() % 2 == 1;
+    let mut out = Vec::new();
+    // Neighbouring inside pieces are one piece.
+    let mut run: Option<(f64, f64)> = None;
+    for w in cuts.windows(2) {
+        let (a, b) = (w[0], w[1]);
+        if !(b - a > 1e-12) {
+            continue;
+        }
+        if inside(at_t(0.5 * (a + b))) {
+            run = Some(run.map_or((a, b), |(s, _)| (s, b)));
+        } else if let Some(r) = run.take() {
+            out.push(iso_record(&f.surface, along_u, fixed, r)?);
+        }
+    }
+    if let Some(r) = run {
+        out.push(iso_record(&f.surface, along_u, fixed, r)?);
+    }
+    if out.is_empty() {
+        return Err(KernelError::Invalid("the iso curve runs outside the face".into()));
+    }
+    Ok(out)
+}
+
+/// A pcurve's parameters and points at [`ISO_SAMPLES`] per knot span, in its direction.
+fn pcurve_samples(c: &brep::Curve2) -> Vec<(f64, [f64; 2])> {
+    let Some([a, b]) = c.interval() else { return Vec::new() };
+    let spans = match c {
+        brep::Curve2::Spline { knots, .. } | brep::Curve2::RationalSpline { knots, .. } => knots.windows(2).filter(|w| w[1] > w[0]).count().max(1),
+    };
+    let n = (spans * ISO_SAMPLES).min(1 << 16);
+    (0..=n)
+        .filter_map(|i| {
+            let t = a + (b - a) * i as f64 / n as f64;
+            c.eval(t).map(|q| (t, q))
+        })
+        .collect()
+}
+
+/// Whether the ray from `q` along +u crosses the closed polyline `l` an odd number of times.
+fn crosses_right(q: [f64; 2], l: &[[f64; 2]]) -> bool {
+    let mut odd = false;
+    for w in l.windows(2) {
+        let (a, b) = (w[0], w[1]);
+        if (a[1] > q[1]) != (b[1] > q[1]) && q[0] < a[0] + (b[0] - a[0]) * (q[1] - a[1]) / (b[1] - a[1]) {
+            odd = !odd;
+        }
+    }
+    odd
+}
+
+/// The iso curve of `surface` at `fixed` (the other parameter) as a record curve on `[a, b]` of
+/// its own parameter: the surface's own line or circle, or its B-spline row.
+fn iso_record(surface: &brep::Surface, along_u: bool, fixed: f64, (a, b): (f64, f64)) -> Result<(brep::Curve, [f64; 2])> {
+    use solidcraft_geom::brep::{Curve, Frame, Surface};
+    let v = |q: [f64; 3]| Vec3::new(q[0], q[1], q[2]);
+    let d = |q: Vec3| [q.x, q.y, q.z];
+    let line = |p: Vec3, q: Vec3| (Curve::Line { start: d(p), end: d(q) }, [0.0, 1.0]);
+    let circle = |o: Vec3, x: Vec3, y: Vec3, r: f64| (Curve::Circle { frame: Frame { origin: d(o), x_axis: d(x), y_axis: d(y) }, radius: r }, [a, b]);
+    let frame = |f: &Frame| (v(f.origin), v(f.x_axis), v(f.y_axis), v(f.x_axis).cross(v(f.y_axis)));
+    Ok(match surface {
+        Surface::Plane { origin, u_axis, v_axis } => {
+            let (o, x, y) = (v(*origin), v(*u_axis), v(*v_axis));
+            let at = |t: f64| if along_u { o + x * t + y * fixed } else { o + x * fixed + y * t };
+            line(at(a), at(b))
+        }
+        Surface::PlaneThrough { origin, u_point, v_point } => {
+            let (o, x, y) = (v(*origin), v(*u_point) - v(*origin), v(*v_point) - v(*origin));
+            let at = |t: f64| if along_u { o + x * t + y * fixed } else { o + x * fixed + y * t };
+            line(at(a), at(b))
+        }
+        Surface::Cylinder { frame: f, radius } => {
+            let (o, x, y, z) = frame(f);
+            if along_u {
+                circle(o + z * fixed, x, y, *radius)
+            } else {
+                let r = (x * fixed.cos() + y * fixed.sin()) * *radius;
+                line(o + r + z * a, o + r + z * b)
+            }
+        }
+        Surface::Cone { frame: f, radius_at_origin, half_angle } => {
+            let (o, x, y, z) = frame(f);
+            let r = |h: f64| radius_at_origin + h * half_angle.tan();
+            if along_u {
+                circle(o + z * fixed, x, y, r(fixed))
+            } else {
+                let dir = x * fixed.cos() + y * fixed.sin();
+                line(o + dir * r(a) + z * a, o + dir * r(b) + z * b)
+            }
+        }
+        Surface::Sphere { frame: f, radius } => {
+            let (o, x, y, z) = frame(f);
+            if along_u {
+                circle(o + z * (radius * fixed.sin()), x, y, radius * fixed.cos())
+            } else {
+                circle(o, x * fixed.cos() + y * fixed.sin(), z, *radius)
+            }
+        }
+        Surface::Torus { frame: f, major_radius, minor_radius } => {
+            let (o, x, y, z) = frame(f);
+            if along_u {
+                circle(o + z * (minor_radius * fixed.sin()), x, y, major_radius + minor_radius * fixed.cos())
+            } else {
+                let dir = x * fixed.cos() + y * fixed.sin();
+                circle(o + dir * *major_radius, dir, z, *minor_radius)
+            }
+        }
+        Surface::Spline { u_knots, v_knots, points } => {
+            let rows: Vec<Vec<[f64; 4]>> = points.iter().map(|r| r.iter().map(|q| [q[0], q[1], q[2], 1.0]).collect()).collect();
+            spline_iso(u_knots, v_knots, &rows, along_u, fixed, a, b)?
+        }
+        Surface::RationalSpline { u_knots, v_knots, homogeneous } => spline_iso(u_knots, v_knots, homogeneous, along_u, fixed, a, b)?,
+        Surface::PlacedRevolution { .. } => {
+            return Err(KernelError::Failed("not supported yet: an iso curve of a revolved surface record".into()));
+        }
+    })
+}
+
+/// The row of a B-spline surface (homogeneous `net[i][j]`, i along u) at `fixed`: along u, each
+/// control point the v-column through it evaluated at `fixed` (along v, the u-row), on `[a, b]`.
+fn spline_iso(u_knots: &[f64], v_knots: &[f64], net: &[Vec<[f64; 4]>], along_u: bool, fixed: f64, a: f64, b: f64) -> Result<(brep::Curve, [f64; 2])> {
+    use solidcraft_geom::bspline::CurveN;
+    let bad = || KernelError::Invalid("iso curve: the surface's B-spline does not evaluate".into());
+    let (nu, nv) = (net.len(), net.first().map_or(0, Vec::len));
+    if nu * nv > MAX_CURVE_POINTS || net.iter().any(|r| r.len() != nv) {
+        return Err(bad());
+    }
+    let (pu, pv) = (u_knots.len().checked_sub(nu + 1).ok_or_else(bad)?, v_knots.len().checked_sub(nv + 1).ok_or_else(bad)?);
+    let at = |degree: usize, knots: &[f64], ctrl: Vec<[f64; 4]>| CurveN::<4> { degree, knots: knots.to_vec(), ctrl }.eval(fixed).map(|(p, _)| p);
+    let (knots, homogeneous): (Vec<f64>, Vec<[f64; 4]>) = if along_u {
+        (u_knots.to_vec(), net.iter().map(|row| at(pv, v_knots, row.clone())).collect::<Option<_>>().ok_or_else(bad)?)
+    } else {
+        (
+            v_knots.to_vec(),
+            (0..nv).map(|j| at(pu, u_knots, net.iter().filter_map(|r| r.get(j).copied()).collect())).collect::<Option<_>>().ok_or_else(bad)?,
+        )
+    };
+    Ok((brep::Curve::RationalSpline { knots, homogeneous }, [a, b]))
+}
+
 /// Most control points a record curve given to [`curve_path`] may have.
 const MAX_CURVE_POINTS: usize = 100_000;