Hidden lines: a traced intersection is drawn only where its spline lies on both surfaces

A leader off either surface is refused like a curve with no spline to follow. Tests: a traced rim
drawn on the circle and plane within tolerance, the same leader 1 mm off refused, and kernel bodies
with arcs through the production prepare against closed forms.

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

crates/geom/src/hlr/curve.rs

@@ -5,7 +5,7 @@
 
 use super::HlrError;
 use crate::Vec3;
-use crate::brep::eval::CurveEval;
+use crate::brep::eval::{CurveEval, SurfaceEval};
 use crate::brep::{Curve, Edge, Frame};
 use crate::section::CurveAccuracy;
 
@@ -44,6 +44,24 @@ pub(crate) fn frame_axes(f: &Frame) -> (Vec3, Vec3, Vec3) {
     (v3(f.origin), x, y)
 }
 
+/// Whether every sample of `c`, a few to each of its breaks, lies within `tolerance` of both
+/// `surfaces`.
+fn on_both_surfaces(c: &C3, surfaces: &[Box<crate::brep::Surface>; 2], tolerance: f64) -> bool {
+    /// Most samples checked; a longer curve's are spread evenly over it.
+    const MAX_SAMPLES: usize = 512;
+    let [Ok(a), Ok(b)] = [SurfaceEval::new(&surfaces[0]), SurfaceEval::new(&surfaces[1])] else { return false };
+    let breaks = c.breaks();
+    let samples: Vec<f64> = breaks.windows(2).flat_map(|w| (0..4).map(move |i| w[0] + (w[1] - w[0]) * f64::from(i) / 4.0)).chain(breaks.last().copied()).collect();
+    let stride = (samples.len() / MAX_SAMPLES).max(1);
+    samples.iter().step_by(stride).all(|t| {
+        let p = c.at(*t);
+        // A plane's `project` answers with the foot of a point however far off, so the distance
+        // is measured here.
+        let on = |s: &SurfaceEval| s.project(p, None, tolerance).and_then(|uv| s.point(uv)).is_some_and(|q| q.dist(p) <= tolerance);
+        p.is_finite() && on(&a) && on(&b)
+    })
+}
+
 impl C3 {
     /// The bounded curve of an edge. A line, circle or ellipse with a NURBS form has its bounds
     /// on the form ([`Edge::parameters`]): a line's ends are evaluated there, a conic's angles come
@@ -51,7 +69,7 @@ impl C3 {
     pub fn of_edge(e: &Edge) -> Result<(C3, CurveAccuracy), HlrError> {
         let exact = |c: C3| (c, CurveAccuracy::Analytic);
         match (&e.curve, &e.nurbs) {
-            (Curve::SurfaceIntersection { leader, tolerance, .. }, form) => {
+            (Curve::SurfaceIntersection { surfaces, leader, tolerance }, form) => {
                 // Procedural: traced on demand from its leader. A spline to follow (its form, or
                 // the leader itself) is a fitted curve, good to the record's tolerance; a polyline
                 // leader carries no such curve.
@@ -59,7 +77,13 @@ impl C3 {
                     (Some(n @ (Curve::Spline { .. } | Curve::RationalSpline { .. })), _) | (None, n @ Curve::Spline { .. }) => n,
                     _ => return Err(HlrError::Unsupported("surface intersection curve")),
                 };
-                Ok((C3::bounded(fitted, e.interval)?, CurveAccuracy::Bounded { max_deviation_mm: *tolerance }))
+                let c = C3::bounded(fitted, e.interval)?;
+                // What is drawn is that spline: it has to lie on both surfaces, or the tolerance
+                // the record claims says nothing and the drawing would show a wrong line.
+                if !on_both_surfaces(&c, surfaces, *tolerance) {
+                    return Err(HlrError::Unsupported("surface intersection curve"));
+                }
+                Ok((c, CurveAccuracy::Bounded { max_deviation_mm: *tolerance }))
             }
             (Curve::Line { .. }, Some(form)) => {
                 let eval = CurveEval::new(form).map_err(|_| HlrError::Invalid("edge form"))?;

crates/geom/src/hlr_tests.rs

@@ -2682,17 +2682,21 @@ fn an_edge_is_cut_to_its_interval_on_its_unbounded_curve() {
     let (part, _) = C3::of_edge(&bounded_edge(spline.clone(), None, [0.1, 0.2])).unwrap();
     assert_eq!(part.breaks(), vec![0.1, 0.2]);
     assert!(part.at(0.15).dist(whole.at(0.15)) < 1e-12);
-    // A traced intersection with a spline leader is that spline, good to its tolerance; with a
-    // polyline leader there is no curve to follow.
+    // A traced intersection with a spline leader on its surfaces is that spline, good to its
+    // tolerance; with a polyline leader there is no curve to follow.
     let plane = Box::new(brep::Surface::Plane { origin: [0.; 3], u_axis: [1., 0., 0.], v_axis: [0., 1., 0.] });
     let traced = |leader: brep::Curve| brep::Curve::SurfaceIntersection {
         surfaces: [plane.clone(), plane.clone()],
         leader: Box::new(leader),
         tolerance: 1e-4,
     };
-    let (c, acc) = C3::of_edge(&bounded_edge(traced(spline), None, [0.25, 0.75])).unwrap();
+    let flat = brep::Curve::Spline { knots: vec![0., 0., 0., 0., 0.5, 1., 1., 1., 1.], points: vec![[0., 0., 0.], [1., 3., 0.], [2., -3., 0.], [3., 3., 0.], [4., 0., 0.]] };
+    let (whole, _) = C3::of_edge(&bounded_edge(flat.clone(), None, [0., 1.])).unwrap();
+    let (c, acc) = C3::of_edge(&bounded_edge(traced(flat), None, [0.25, 0.75])).unwrap();
     assert_eq!(acc, CurveAccuracy::Bounded { max_deviation_mm: 1e-4 });
     assert!(c.at(0.5).dist(whole.at(0.5)) < 1e-12);
+    // That spline leaves the plane: it is not their intersection.
+    assert!(matches!(C3::of_edge(&bounded_edge(traced(spline), None, [0.25, 0.75])), Err(HlrError::Unsupported("surface intersection curve"))));
     let leader = brep::Curve::Polyline { points: vec![[0.; 3], [1., 0., 0.]] };
     assert!(matches!(C3::of_edge(&bounded_edge(traced(leader), None, [0., 1.])), Err(HlrError::Unsupported("surface intersection curve"))));
 }
@@ -2776,3 +2780,45 @@ fn an_edge_no_face_uses_is_drawn_like_any_other() {
     let d = hidden_lines(&[def], &[identity(0, [0.; 3])], &view, &Options::default(), &never).unwrap();
     assert!(edge_line(&d, 0, 12).pieces.iter().all(|p| p.visibility == Visibility::Visible), "{d:?}");
 }
+
+/// A traced intersection is drawn from its spline leader only where that spline lies on both
+/// surfaces, so the drawn line is the true intersection to the record's tolerance. The top rim of
+/// a cylinder of radius 10 and height 10, stored as the meeting of its side and the plane z = 10,
+/// its leader four cubic quarter circles (2.8e-3 from the circle at most): seen from above the
+/// drawn rim is that circle, and from the side the plane, within the tolerance. The same leader
+/// 1 mm off the plane is refused: it would draw a wrong line.
+#[test]
+fn a_traced_rim_is_drawn_on_the_true_intersection_or_refused() {
+    let tolerance = 5e-3;
+    let k = 10. * 4. / 3. * (std::f64::consts::FRAC_PI_8).tan();
+    let knots = vec![0., 0., 0., 0., 1., 1., 1., 2., 2., 2., 3., 3., 3., 4., 4., 4., 4.];
+    let quarter = [[10., 0.], [10., k], [k, 10.], [0., 10.], [-k, 10.], [-10., k], [-10., 0.], [-10., -k], [-k, -10.], [0., -10.], [k, -10.], [10., -k], [10., 0.]];
+    let traced = |z: f64| {
+        let mut body = cylinder(10., 10.);
+        let side = Box::new(body.shells[0].faces[0].surface.clone());
+        let top = Box::new(Surface::Plane { origin: [0., 0., 10.], u_axis: [1., 0., 0.], v_axis: [0., 1., 0.] });
+        let points = quarter.iter().map(|[x, y]| [*x, *y, z]).collect();
+        let leader = Box::new(Curve::Spline { knots: knots.clone(), points });
+        body.shells[0].edges[1].curve = Curve::SurfaceIntersection { surfaces: [side, top], leader, tolerance };
+        body
+    };
+    let o = Options::default();
+    let def = prepare(&traced(10.), &o).unwrap();
+    let rim_points = |view: &View| -> Vec<[f64; 2]> {
+        let d = hidden_lines(std::slice::from_ref(&def), &[identity(0, [0.; 3])], view, &o, &never).unwrap();
+        let rim: Vec<_> = d.lines.iter().filter(|l| matches!(l.source, Source::Edge { edge: 1, .. })).flat_map(|l| &l.pieces).flat_map(|p| p.points.iter().copied()).collect();
+        assert!(rim.len() > 16, "{rim:?}");
+        rim
+    };
+    let above = rim_points(&View::new(Vec3::new(0., 0., -1.), Vec3::new(0., 1., 0.)).unwrap());
+    assert!(above.iter().all(|p| (p[0].hypot(p[1]) - 10.).abs() <= tolerance), "{above:?}");
+    let spread = |at: usize| above.iter().map(|p| p[at]).fold(f64::MIN, f64::max) - above.iter().map(|p| p[at]).fold(f64::MAX, f64::min);
+    assert!(spread(0) > 19.9 && spread(1) > 19.9, "the whole circle");
+    let side = rim_points(&View::new(Vec3::new(0., 1., 0.), Vec3::new(0., 0., 1.)).unwrap());
+    let height = side.first().map_or(f64::NAN, |p| p[1]);
+    let lowest = hidden_lines(std::slice::from_ref(&def), &[identity(0, [0.; 3])], &View::new(Vec3::new(0., 1., 0.), Vec3::new(0., 0., 1.)).unwrap(), &o, &never).unwrap();
+    let bottom = lowest.lines.iter().filter(|l| matches!(l.source, Source::Edge { edge: 0, .. })).flat_map(|l| &l.pieces).flat_map(|p| p.points.iter()).next().map_or(f64::NAN, |p| p[1]);
+    assert!(side.iter().all(|p| (p[1] - height).abs() <= tolerance) && ((height - bottom) - 10.).abs() <= tolerance, "{side:?} over {bottom}");
+    let off = prepare(&traced(11.), &o).err();
+    assert!(matches!(off, Some(HlrError::Unsupported("surface intersection curve"))), "{off:?}");
+}

crates/kernel/tests/hidden_lines.rs

@@ -0,0 +1,92 @@
+//! Hidden lines of bodies the kernel makes, through the production `prepare`, against closed forms.
+
+use std::f64::consts::{FRAC_PI_3, PI};
+
+use solidcraft_geom::hlr::{Line, Options, Placement, Source, View, Visibility, hidden_lines, prepare};
+use solidcraft_geom::{Seg2, Vec2, Vec3, brep};
+use solidcraft_kernel::*;
+
+const TOLERANCE: f64 = 1e-3;
+
+/// Points along a drawn piece, a sixteenth of it apart.
+fn points(s: &Seg2) -> Vec<Vec2> {
+    (0..=16).map(|i| s.at(f64::from(i) / 16.0)).collect()
+}
+
+/// The lines of `body`, drawn looking along `d`, with the record they number into.
+fn draw(body: &Body, d: Vec3, up: Vec3) -> (Vec<Line>, brep::Brep) {
+    let record = native::brep_record(body).unwrap();
+    let o = Options::default();
+    let def = prepare(&record, &o).unwrap();
+    let placed = Placement { body: 0, rotation: [Vec3::new(1., 0., 0.), Vec3::new(0., 1., 0.), Vec3::new(0., 0., 1.)], translation: Vec3::ZERO };
+    let view = View::new(d, up).unwrap();
+    let drawing = hidden_lines(std::slice::from_ref(&def), &[placed], &view, &o, &|| false).unwrap();
+    (drawing.lines, record)
+}
+
+/// Whether a line is an edge whose record is a circle (`want`), or another edge.
+fn is_circle(l: &Line, record: &brep::Brep, want: bool) -> bool {
+    matches!(l.source, Source::Edge { edge, .. } if matches!(record.edges.get(edge).map(|e| &e.curve), Some(brep::Curve::Circle { .. })) == want)
+}
+
+/// Total drawn length of `lines` of a visibility.
+fn total<'a>(lines: impl Iterator<Item = &'a Line>, v: Visibility) -> f64 {
+    lines.flat_map(|l| &l.pieces).filter(|p| p.visibility == v).flat_map(|p| &p.geometry).map(Seg2::length).sum()
+}
+
+fn near(got: f64, want: f64, what: &str) {
+    assert!((got - want).abs() <= TOLERANCE, "{what}: {got} against {want}");
+}
+
+/// A cylinder of radius 10 cut by the plane x = 5: rims from -60 to 60 degrees. Seen along z each
+/// rim is one arc of a third of a turn on the side x >= 5, the near one visible, the far one
+/// behind it; a rim on the wrong angles would leave that side.
+#[test]
+fn a_cut_cylinders_rims_are_drawn_on_their_own_arc() {
+    let round = cylinder(Vec3::ZERO, Vec3::new(0., 0., 1.), 10., 20.).unwrap();
+    let keep = box_solid(Vec3::new(5., -20., -5.), Vec3::new(25., 20., 25.)).unwrap();
+    let body = boolean(&round, &keep, BoolOp::Intersect).unwrap().unwrap();
+    let (lines, record) = draw(&body, Vec3::new(0., 0., -1.), Vec3::new(0., 1., 0.));
+    let rim = 10. * 2. * FRAC_PI_3;
+    let rims = || lines.iter().filter(|l| is_circle(l, &record, true));
+    near(total(rims(), Visibility::Visible), rim, "near rim");
+    near(total(rims(), Visibility::Hidden), rim, "far rim");
+    let drawn: Vec<Vec2> = rims().flat_map(|l| &l.pieces).flat_map(|p| &p.geometry).flat_map(points).collect();
+    assert!(drawn.len() > 8, "{lines:?}");
+    // Whichever way the view's x runs, the rims keep to one side of the axis, at the radius.
+    let side = drawn.iter().map(|p| p.x).fold(0., |a: f64, x| if x.abs() > a.abs() { x } else { a }).signum();
+    for p in &drawn {
+        assert!((p.x.hypot(p.y) - 10.).abs() <= TOLERANCE && p.x * side >= 5. - TOLERANCE, "{p:?}");
+    }
+    // The chord closes them: 2 * 10 * sin 60 on each end of the flat.
+    let flats = || lines.iter().filter(|l| is_circle(l, &record, false));
+    let chord = 20. * FRAC_PI_3.sin();
+    near(total(flats(), Visibility::Visible), chord, "near chord");
+    near(total(flats(), Visibility::Hidden), chord, "far chord");
+    // From the flat side, along x, the rims are level lines as wide as the chord.
+    let (lines, record) = draw(&body, Vec3::new(1., 0., 0.), Vec3::new(0., 0., 1.));
+    let wide: Vec<Vec2> = lines.iter().filter(|l| is_circle(l, &record, true)).flat_map(|l| &l.pieces).flat_map(|p| &p.geometry).flat_map(points).collect();
+    let reach = wide.iter().map(|p| p.x.abs()).fold(0., f64::max);
+    near(reach, chord / 2., "rims from the flat side");
+}
+
+/// A hole of radius 5 through a 40 by 40 by 10 plate. Along z its near rim is a whole circle in
+/// sight, the far one behind the plate; from the side, its outline is two hidden lines the
+/// plate's thickness long and 5 either side of the axis, and nothing of it shows.
+#[test]
+fn a_holes_rims_and_outline_are_drawn_where_they_are() {
+    let plate = box_solid(Vec3::new(-20., -20., 0.), Vec3::new(20., 20., 10.)).unwrap();
+    let bore = cylinder(Vec3::new(0., 0., -5.), Vec3::new(0., 0., 1.), 5., 20.).unwrap();
+    let body = boolean(&plate, &bore, BoolOp::Cut).unwrap().unwrap();
+    let (lines, record) = draw(&body, Vec3::new(0., 0., -1.), Vec3::new(0., 1., 0.));
+    let rims = || lines.iter().filter(|l| is_circle(l, &record, true));
+    near(total(rims(), Visibility::Visible), 2. * PI * 5., "near rim");
+    near(total(rims(), Visibility::Hidden), 2. * PI * 5., "far rim");
+    let (lines, _) = draw(&body, Vec3::new(0., 1., 0.), Vec3::new(0., 0., 1.));
+    let outline = || lines.iter().filter(|l| matches!(l.source, Source::Silhouette { .. }));
+    near(total(outline(), Visibility::Hidden), 2. * 10., "the outline's hidden length");
+    near(total(outline(), Visibility::Visible), 0., "the outline in sight");
+    for p in outline().flat_map(|l| &l.pieces).flat_map(|p| &p.geometry).flat_map(points) {
+        assert!((p.x.abs() - 5.).abs() <= TOLERANCE, "{p:?}");
+    }
+}