Drawings hidden lines: read develop's B-rep records

Edges come from the record's edge list, cut to their interval (C3::of_edge), and
faces from its face list; a cylinder, cone, sphere or torus chart takes its bounds
from where its NURBS form lies on its label (S3::of_face). A rigidly placed
revolution has main's revolution chart; other placements stay refused. Moving a body
to its centre goes through Brep::moved. Loops are read against the natural normal of
the face's parameters where a chart wraps.

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

crates/geom/src/hlr/curve.rs

@@ -1,9 +1,12 @@
-//! Evaluation of portable B-rep curves: lines, circle and ellipse arcs, polynomial and rational
-//! splines (fitted intersection records are ordinary splines with a stated fitting tolerance).
+//! Evaluation of portable B-rep edges: lines, circle and ellipse arcs, polynomial and rational
+//! splines (a traced intersection whose leader is a spline is that spline, within its tolerance).
+//! Records' curves are unbounded: every curve here is cut to its edge's interval, so `domain` and
+//! `breaks` are the edge's.
 
 use super::HlrError;
 use crate::Vec3;
-use crate::brep::{Curve, Frame};
+use crate::brep::eval::CurveEval;
+use crate::brep::{Curve, Edge, Frame};
 use crate::section::CurveAccuracy;
 
 #[derive(Clone, Debug)]
@@ -42,29 +45,82 @@ pub(crate) fn frame_axes(f: &Frame) -> (Vec3, Vec3, Vec3) {
 }
 
 impl C3 {
-    pub fn new(c: &Curve) -> Result<(C3, CurveAccuracy), HlrError> {
+    /// 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
+    /// from [`Edge::curve_interval`].
+    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) => {
+                // 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.
+                let fitted = match (form, leader.as_ref()) {
+                    (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 }))
+            }
+            (Curve::Line { .. }, Some(form)) => {
+                let eval = CurveEval::new(form).map_err(|_| HlrError::Invalid("edge form"))?;
+                let end = |t: f64| eval.point(t).filter(|p| p.is_finite()).ok_or(HlrError::Invalid("edge form"));
+                Ok(exact(C3::Line { a: end(e.interval[0])?, b: end(e.interval[1])? }))
+            }
+            (Curve::Circle { .. } | Curve::Ellipse { .. }, _) => {
+                let angles = e.curve_interval().ok_or(HlrError::Invalid("edge form"))?;
+                Ok(exact(C3::bounded(&e.curve, angles)?))
+            }
+            (curve, _) => Ok(exact(C3::bounded(curve, e.interval)?)),
+        }
+    }
+
+    /// `c` between the parameters `t` on its own parameter (a line then runs over [0, 1], a
+    /// polyline over [0, its segments]).
+    pub fn bounded(c: &Curve, t: [f64; 2]) -> Result<C3, HlrError> {
+        if !(t[0].is_finite() && t[1].is_finite() && t[0] < t[1]) {
+            return Err(HlrError::Invalid("curve interval"));
+        }
         Ok(match c {
-            Curve::Line { start, end } => (C3::Line { a: v3(*start), b: v3(*end) }, CurveAccuracy::Analytic),
-            Curve::CircleArc { frame, radius, angle } => {
-                let (o, x, y) = frame_axes(frame);
-                (C3::Arc { o, x, y, r: [*radius, *radius], range: *angle }, CurveAccuracy::Analytic)
+            Curve::Line { start, end } => {
+                let (a, d) = (v3(*start), v3(*end) - v3(*start));
+                let (p, q) = (a + d * t[0], a + d * t[1]);
+                if !(p.is_finite() && q.is_finite()) {
+                    return Err(HlrError::Invalid("curve interval"));
+                }
+                C3::Line { a: p, b: q }
             }
-            Curve::EllipseArc { frame, radii, angle } => {
+            Curve::Circle { frame, radius } => {
                 let (o, x, y) = frame_axes(frame);
-                (C3::Arc { o, x, y, r: *radii, range: *angle }, CurveAccuracy::Analytic)
+                C3::Arc { o, x, y, r: [*radius, *radius], range: t }
             }
-            Curve::Polyline { points } => (C3::Polyline(points.iter().copied().map(v3).collect()), CurveAccuracy::Analytic),
-            Curve::Spline { knots, points } => (spline(knots, points.iter().map(|p| [p[0], p[1], p[2], 1.]).collect())?, CurveAccuracy::Analytic),
-            Curve::RationalSpline { knots, homogeneous } => (spline(knots, homogeneous.clone())?, CurveAccuracy::Analytic),
-            Curve::Intersection { knots, points, fitting_tolerance, .. } => (
-                spline(knots, points.iter().map(|p| [p[0], p[1], p[2], 1.]).collect())?,
-                CurveAccuracy::Bounded { max_deviation_mm: *fitting_tolerance },
-            ),
-            // Procedural: the curve is traced on demand from its leader, which carries no bound.
+            Curve::Ellipse { frame, radii } => {
+                let (o, x, y) = frame_axes(frame);
+                C3::Arc { o, x, y, r: *radii, range: t }
+            }
+            Curve::Polyline { points } => C3::Polyline(polyline_part(points, t)?),
+            Curve::Spline { knots, points } => part(knots, points.iter().map(|p| [p[0], p[1], p[2], 1.]).collect(), t)?,
+            Curve::RationalSpline { knots, homogeneous } => part(knots, homogeneous.clone(), t)?,
             Curve::SurfaceIntersection { .. } => return Err(HlrError::Unsupported("surface intersection curve")),
         })
     }
 
+    /// The curve under a placement: `point` maps positions, `direction` directions (a rigid motion,
+    /// which keeps the arcs' parameter and a spline's weights).
+    pub fn placed(&self, point: &dyn Fn(Vec3) -> Vec3, direction: &dyn Fn(Vec3) -> Vec3) -> Result<C3, HlrError> {
+        Ok(match self {
+            C3::Line { a, b } => C3::Line { a: point(*a), b: point(*b) },
+            C3::Arc { o, x, y, r, range } => C3::Arc { o: point(*o), x: direction(*x), y: direction(*y), r: *r, range: *range },
+            C3::Polyline(p) => C3::Polyline(p.iter().map(|q| point(*q)).collect()),
+            C3::Spline { knots, points, .. } => {
+                let moved = points.iter().map(|h| {
+                    let q = point(Vec3::new(h[0] / h[3], h[1] / h[3], h[2] / h[3]));
+                    [q.x * h[3], q.y * h[3], q.z * h[3], h[3]]
+                });
+                spline(knots, moved.collect())?
+            }
+        })
+    }
+
     pub fn domain(&self) -> [f64; 2] {
         match self {
             C3::Line { .. } => [0., 1.],
@@ -139,6 +195,57 @@ fn spline(knots: &[f64], points: Vec<[f64; 4]>) -> Result<C3, HlrError> {
     Ok(C3::Spline { degree, knots: knots.to_vec(), points, pieces })
 }
 
+/// The polyline's points between the parameters `t` (segment i spans [i, i + 1]).
+fn polyline_part(points: &[[f64; 3]], t: [f64; 2]) -> Result<Vec<Vec3>, HlrError> {
+    let last = points.len().saturating_sub(1) as f64;
+    let t = [t[0].clamp(0., last), t[1].clamp(0., last)];
+    let at = |x: f64| {
+        let i = (x.floor().max(0.) as usize).min(points.len().saturating_sub(2));
+        match (points.get(i), points.get(i + 1)) {
+            (Some(a), Some(b)) => Ok(v3(*a).lerp(v3(*b), x - i as f64)),
+            _ => Err(HlrError::Invalid("polyline")),
+        }
+    };
+    let inner = (t[0].floor() as i64 + 1..t[1].ceil() as i64).filter_map(|i| usize::try_from(i).ok()).filter_map(|i| points.get(i).copied()).map(v3);
+    Ok(std::iter::once(at(t[0])?).chain(inner).chain(std::iter::once(at(t[1])?)).collect())
+}
+
+/// The spline's part between the parameters `t`, itself a clamped spline on them: the whole
+/// spline as it stands when `t` is its whole active interval.
+fn part(knots: &[f64], points: Vec<[f64; 4]>, t: [f64; 2]) -> Result<C3, HlrError> {
+    use crate::bspline::{HPoint, clamp_row, insert_knot};
+    let refused = HlrError::Unsupported("spline edge interval");
+    let p = knots.len().checked_sub(points.len() + 1).ok_or(HlrError::Unsupported("spline knot count"))?;
+    let (Some(&lo), Some(&hi)) = (knots.get(p), knots.get(points.len())) else { return Err(HlrError::Unsupported("spline knot count")) };
+    if t[0] <= lo && hi <= t[1] {
+        return spline(knots, points);
+    }
+    let row: Vec<HPoint> = points.iter().map(|h| HPoint { x: h[0], y: h[1], z: h[2], w: h[3] }).collect();
+    let Ok((mut k, mut c)) = clamp_row(p, knots, &row) else { return Err(refused) };
+    let t = [t[0].max(lo), t[1].min(hi)];
+    if !(t[0] < t[1]) {
+        return Err(refused);
+    }
+    for u in t {
+        if u > lo && u < hi {
+            // Each end a knot of full multiplicity: the control net splits there.
+            for _ in k.iter().filter(|x| **x == u).count()..p {
+                let Ok((nk, nc)) = insert_knot(p, &k, &c, u) else { return Err(refused) };
+                (k, c) = (nk, nc);
+            }
+        }
+    }
+    let last_at_start = k.partition_point(|x| *x <= t[0]).checked_sub(1).ok_or(HlrError::Unsupported("spline edge interval"))?;
+    let first_at_end = k.partition_point(|x| *x < t[1]);
+    let (Some(inner), Some(net)) = (k.get(last_at_start + 1..first_at_end), last_at_start.checked_sub(p).and_then(|s| c.get(s..first_at_end))) else {
+        return Err(refused);
+    };
+    let mut cut = vec![t[0]; p + 1];
+    cut.extend_from_slice(inner);
+    cut.extend(std::iter::repeat_n(t[1], p + 1));
+    spline(&cut, net.iter().map(|h| [h.x, h.y, h.z, h.w]).collect())
+}
+
 /// A spline split into Bézier pieces: clamped, then every inner knot raised to the degree
 /// (NURBS Book §5.2); empty when that fails.
 fn bezier_pieces(p: usize, knots: &[f64], points: &[[f64; 4]]) -> Vec<([f64; 2], Vec<[f64; 4]>)> {

crates/geom/src/hlr/mod.rs

@@ -4,8 +4,9 @@
 //! trimmed faces. Every piece is visible, hidden or explicitly uncertain; nothing is guessed.
 //!
 //! Supported: planes, and cylinders, cones, spheres and tori whether stored as analytic records,
-//! revolutions of lines and circles or spline patches (recognised to 1e-9 relative and then
-//! trimmed and intersected in their own chart). Other spline and revolution faces are freeform:
+//! rigidly placed revolutions of lines and circles or spline patches (recognised to 1e-9 relative
+//! and then trimmed and intersected in their own chart). Other spline and rigidly placed revolution
+//! faces are freeform:
 //! unless [`Options::freeform`] is set they are listed in [`Drawing::unsupported`], their
 //! silhouettes are not drawn and anything they could hide is reported uncertain.
 //!
@@ -200,149 +201,147 @@ pub fn prepare_cancellable(b: &Brep, opts: &Options, cancel: &dyn Fn() -> bool)
     let mut edges = Vec::new();
     let mut faces = Vec::new();
     let mut report = PrepareReport::default();
-    for shell in &b.shells {
-        let base_face = faces.len();
-        let mut curves = Vec::with_capacity(shell.edges.len());
-        for (ei, e) in shell.edges.iter().enumerate() {
-            if ei % 256 == 0 && cancel() {
-                return Err(HlrError::Cancelled);
-            }
-            let (c, acc) = C3::new(&e.curve)?;
-            let record = match &e.curve {
-                crate::brep::Curve::Line { .. } => "line",
-                crate::brep::Curve::CircleArc { .. } => "circle_arc",
-                crate::brep::Curve::EllipseArc { .. } => "ellipse_arc",
-                crate::brep::Curve::Polyline { .. } => "polyline",
-                crate::brep::Curve::Spline { .. } => "spline",
-                crate::brep::Curve::RationalSpline { .. } => "rational_spline",
-                crate::brep::Curve::Intersection { .. } => "intersection",
-                crate::brep::Curve::SurfaceIntersection { .. } => "surface_intersection",
-            };
-            let shape = if matches!(c, C3::Line { .. }) {
-                "line"
-            } else if as_circle(&c, 1e-9).is_some() {
-                "circle"
-            } else {
-                "other"
-            };
-            *report.curves.entry(format!("{record}/{shape}")).or_insert(0) += 1;
-            curves.push((c, acc));
+    let mut curves = Vec::with_capacity(b.edges.len());
+    for (ei, e) in b.edges.iter().enumerate() {
+        if ei % 256 == 0 && cancel() {
+            return Err(HlrError::Cancelled);
         }
-        let local_curves: Vec<C3> = curves.iter().map(|(c, _)| c.clone()).collect();
-        let mut users: Vec<Vec<usize>> = vec![Vec::new(); shell.edges.len()];
-        for (fi, f) in shell.faces.iter().enumerate() {
-            if cancel() {
-                return Err(HlrError::Cancelled);
-            }
-            let mut surface = S3::new(&f.surface)?;
-            // A flat spline patch is a plane: trimmed and intersected as one.
-            if let Some(plane) = surface.flat(1e-9) {
-                surface = plane;
-            }
-            let analytic = surface.analytic(1e-9);
-            // Trim and intersect recognised spline surfaces in their own analytic chart.
-            if matches!(surface, S3::Nurbs(_)) {
-                let axis = match analytic {
-                    Analytic::Cylinder { p, a, .. } => Some((p, a)),
-                    Analytic::Cone { apex, a, .. } => Some((apex, a)),
-                    _ => None,
-                };
-                let mut h = [0., 0.];
-                if let Some((p, a)) = axis {
-                    let along = f.loops.iter().flatten().filter_map(|u| local_curves.get(u.edge as usize)).flat_map(|c| {
-                        let [t0, t1] = c.domain();
-                        (0..=64).map(move |i| c.at(t0 + (t1 - t0) * i as f64 / 64.))
-                    });
-                    let (lo, hi) = along.fold((f64::MAX, f64::MIN), |(lo, hi), q| ((q - p).dot(a).min(lo), (q - p).dot(a).max(hi)));
-                    // Generous: the chart only has to contain the face; the trim bounds it exactly.
-                    let pad = 0.01 * (hi - lo).abs() + 1e-6;
-                    h = [lo - pad, hi + pad];
-                }
-                let toward = match analytic {
-                    Analytic::Sphere { c, .. } => {
-                        let pts: Vec<Vec3> = f
-                            .loops
-                            .iter()
-                            .flatten()
-                            .filter_map(|u| local_curves.get(u.edge as usize))
-                            .map(|c| c.at(0.5 * (c.domain()[0] + c.domain()[1])))
-                            .collect();
-                        let n = pts.len().max(1) as f64;
-                        pts.iter().fold(Vec3::ZERO, |acc, p| acc + *p) / n - c
-                    }
-                    _ => Vec3::ZERO,
-                };
-                if let Some(chart) = S3::from_analytic(&analytic, h, toward) {
-                    surface = chart;
-                }
-            }
-            *report.surfaces.entry(analytic_name(&analytic).to_string()).or_insert(0) += 1;
-            let trim = trim::build(f, &local_curves, &surface, opts.trim_samples);
-            if let Some(why) = trim.unsupported {
-                report.unsupported_trims.push((base_face + fi, why.to_string()));
-            }
-            report.max_trim_chord = report.max_trim_chord.max(trim.chord_error);
-            let (mut lo, mut hi) = (trim.lo, trim.hi);
-            // Curved faces can bulge past their boundary: bound them by the whole surface's
-            // guaranteed box (spline hull, turned profile), clipped to nothing tighter.
-            if !matches!(surface, S3::Plane { .. })
-                && let Some((blo, bhi)) = surface.bound()
-            {
-                lo = lo.min(blo);
-                hi = hi.max(bhi);
-            }
-            // Every face is padded, planes too: a plane record and its edges may disagree by
-            // rounding, and a box flat along the normal would cull the plane's own hit. Rounding
-            // grows with where the body was modelled (its offset), not with the face's size.
-            let pad = (hi - lo).len() * 1e-6 + 1e-12 * offset.x.abs().max(offset.y.abs()).max(offset.z.abs());
-            lo -= Vec3::new(pad, pad, pad);
-            hi += Vec3::new(pad, pad, pad);
-            for u in f.loops.iter().flatten() {
-                if let Some(list) = users.get_mut(u.edge as usize) {
-                    list.push(base_face + fi);
-                }
-            }
-            let freeform = matches!(analytic, Analytic::Other);
-            if freeform {
-                report.freeform_faces.push(base_face + fi);
-            }
-            faces.push(FaceData { surface, analytic, group: Vec::new(), trim, lo, hi, freeform });
+        let (c, acc) = C3::of_edge(e)?;
+        let record = match &e.curve {
+            crate::brep::Curve::Line { .. } => "line",
+            crate::brep::Curve::Circle { .. } => "circle_arc",
+            crate::brep::Curve::Ellipse { .. } => "ellipse_arc",
+            crate::brep::Curve::Polyline { .. } => "polyline",
+            crate::brep::Curve::Spline { .. } => "spline",
+            crate::brep::Curve::RationalSpline { .. } => "rational_spline",
+            // Only one with a spline to follow gets here: what was an `intersection` record.
+            crate::brep::Curve::SurfaceIntersection { .. } => "intersection",
+        };
+        let shape = if matches!(c, C3::Line { .. }) {
+            "line"
+        } else if as_circle(&c, 1e-9).is_some() {
+            "circle"
+        } else {
+            "other"
+        };
+        *report.curves.entry(format!("{record}/{shape}")).or_insert(0) += 1;
+        curves.push((c, acc));
+    }
+    let local_curves: Vec<C3> = curves.iter().map(|(c, _)| c.clone()).collect();
+    let mut users: Vec<Vec<usize>> = vec![Vec::new(); b.edges.len()];
+    // Faces and edges by their numbers in the record. Validation puts every face in a shell; its
+    // loops are stored the same whichever side a shell uses, and hidden lines look at both, so the
+    // shells' orientation does not come into it (a face two shells share is one face). An edge no
+    // face uses (a wire) is drawn like any other.
+    for (fi, f) in b.faces.iter().enumerate() {
+        if cancel() {
+            return Err(HlrError::Cancelled);
         }
-        for (ei, (c, acc)) in curves.into_iter().enumerate() {
-            let list = users.get(ei).cloned().unwrap_or_default();
-            let mut distinct = list.clone();
-            distinct.sort_unstable();
-            distinct.dedup();
-            let kind = if list.len() >= 2 && distinct.len() == 1 {
-                EdgeKind::Seam
-            } else if let [f1, f2] = distinct.as_slice() {
-                let [t0, t1] = c.domain();
-                let normals_parallel = [0.25, 0.5, 0.75].iter().all(|k| {
-                    let m = c.at(t0 + (t1 - t0) * k);
-                    let n = |f: usize| faces.get(f).and_then(|fd: &FaceData| fd.surface.inverse(m, None).and_then(|uv| fd.surface.normal(uv)));
-                    n(*f1).zip(n(*f2)).is_some_and(|(a, b)| a.cross(b).len() < 1e-6)
+        let mut surface = S3::of_face(f)?;
+        // A flat spline patch is a plane: trimmed and intersected as one.
+        if let Some(plane) = surface.flat(1e-9) {
+            surface = plane;
+        }
+        let analytic = surface.analytic(1e-9);
+        // Trim and intersect recognised spline surfaces in their own analytic chart.
+        if matches!(surface, S3::Nurbs(_)) {
+            let axis = match analytic {
+                Analytic::Cylinder { p, a, .. } => Some((p, a)),
+                Analytic::Cone { apex, a, .. } => Some((apex, a)),
+                _ => None,
+            };
+            let mut h = [0., 0.];
+            if let Some((p, a)) = axis {
+                let along = trim::uses(f).filter_map(|u| local_curves.get(u.edge as usize)).flat_map(|c| {
+                    let [t0, t1] = c.domain();
+                    (0..=64).map(move |i| c.at(t0 + (t1 - t0) * i as f64 / 64.))
                 });
-                let same_surface = match (faces.get(*f1).map(|f: &FaceData| f.analytic), faces.get(*f2).map(|f: &FaceData| f.analytic)) {
-                    (Some(a), Some(b)) => same_analytic(&a, &b, 1e-9),
-                    _ => false,
-                };
-                if same_surface {
-                    EdgeKind::Seam
-                } else if normals_parallel {
-                    EdgeKind::Tangent
-                } else {
-                    EdgeKind::Sharp
+                let (lo, hi) = along.fold((f64::MAX, f64::MIN), |(lo, hi), q| ((q - p).dot(a).min(lo), (q - p).dot(a).max(hi)));
+                // Generous: the chart only has to contain the face; the trim bounds it exactly.
+                let pad = 0.01 * (hi - lo).abs() + 1e-6;
+                h = [lo - pad, hi + pad];
+            }
+            let toward = match analytic {
+                Analytic::Sphere { c, .. } => {
+                    let pts: Vec<Vec3> = trim::uses(f)
+                        .filter_map(|u| local_curves.get(u.edge as usize))
+                        .map(|c| c.at(0.5 * (c.domain()[0] + c.domain()[1])))
+                        .collect();
+                    let n = pts.len().max(1) as f64;
+                    pts.iter().fold(Vec3::ZERO, |acc, p| acc + *p) / n - c
                 }
+                _ => Vec3::ZERO,
+            };
+            if let Some(chart) = S3::from_analytic(&analytic, h, toward) {
+                surface = chart;
+            }
+        }
+        *report.surfaces.entry(analytic_name(&analytic).to_string()).or_insert(0) += 1;
+        let trim = trim::build(f, &local_curves, &surface, opts.trim_samples);
+        if let Some(why) = trim.unsupported {
+            report.unsupported_trims.push((fi, why.to_string()));
+        }
+        report.max_trim_chord = report.max_trim_chord.max(trim.chord_error);
+        let (mut lo, mut hi) = (trim.lo, trim.hi);
+        // Curved faces can bulge past their boundary: bound them by the whole surface's
+        // guaranteed box (spline hull, turned profile), clipped to nothing tighter.
+        if !matches!(surface, S3::Plane { .. })
+            && let Some((blo, bhi)) = surface.bound()
+        {
+            lo = lo.min(blo);
+            hi = hi.max(bhi);
+        }
+        // Every face is padded, planes too: a plane record and its edges may disagree by
+        // rounding, and a box flat along the normal would cull the plane's own hit. Rounding
+        // grows with where the body was modelled (its offset), not with the face's size.
+        let pad = (hi - lo).len() * 1e-6 + 1e-12 * offset.x.abs().max(offset.y.abs()).max(offset.z.abs());
+        lo -= Vec3::new(pad, pad, pad);
+        hi += Vec3::new(pad, pad, pad);
+        for u in trim::uses(f) {
+            if let Some(list) = users.get_mut(u.edge as usize) {
+                list.push(fi);
+            }
+        }
+        let freeform = matches!(analytic, Analytic::Other);
+        if freeform {
+            report.freeform_faces.push(fi);
+        }
+        faces.push(FaceData { surface, analytic, group: Vec::new(), trim, lo, hi, freeform });
+    }
+    for (ei, (c, acc)) in curves.into_iter().enumerate() {
+        let list = users.get(ei).cloned().unwrap_or_default();
+        let mut distinct = list.clone();
+        distinct.sort_unstable();
+        distinct.dedup();
+        let kind = if list.len() >= 2 && distinct.len() == 1 {
+            EdgeKind::Seam
+        } else if let [f1, f2] = distinct.as_slice() {
+            let [t0, t1] = c.domain();
+            let normals_parallel = [0.25, 0.5, 0.75].iter().all(|k| {
+                let m = c.at(t0 + (t1 - t0) * k);
+                let n = |f: usize| faces.get(f).and_then(|fd: &FaceData| fd.surface.inverse(m, None).and_then(|uv| fd.surface.normal(uv)));
+                n(*f1).zip(n(*f2)).is_some_and(|(a, b)| a.cross(b).len() < 1e-6)
+            });
+            let same_surface = match (faces.get(*f1).map(|f: &FaceData| f.analytic), faces.get(*f2).map(|f: &FaceData| f.analytic)) {
+                (Some(a), Some(b)) => same_analytic(&a, &b, 1e-9),
+                _ => false,
+            };
+            if same_surface {
+                EdgeKind::Seam
+            } else if normals_parallel {
+                EdgeKind::Tangent
             } else {
                 EdgeKind::Sharp
-            };
-            match kind {
-                EdgeKind::Tangent => report.tangent_edges += 1,
-                EdgeKind::Seam => report.seam_edges += 1,
-                EdgeKind::Sharp => {}
             }
-            edges.push(EdgeData { curve: c, acc, faces: list, kind });
+        } else {
+            EdgeKind::Sharp
+        };
+        match kind {
+            EdgeKind::Tangent => report.tangent_edges += 1,
+            EdgeKind::Seam => report.seam_edges += 1,
+            EdgeKind::Sharp => {}
         }
+        edges.push(EdgeData { curve: c, acc, faces: list, kind });
     }
     // Faces on one quadric or torus. Only surfaces of one kind whose size (radius, half angle,
     // major radius) agrees to the tolerance can match, so each face is compared with the faces
@@ -394,93 +393,20 @@ pub fn prepare_cancellable(b: &Brep, opts: &Options, cancel: &dyn Fn() -> bool)
     Ok(Prepared { edges, faces, offset, lo, hi, report })
 }
 
-/// The centre of the body's vertices' box (of its faces' reference points when it has no
-/// vertices, a sphere without edges say).
+/// The centre of the body's vertices' box (a validated body has vertices: every face has an edge
+/// loop, every other shell an edge or a vertex).
 fn anchor(b: &Brep) -> Vec3 {
-    use crate::brep::Surface;
-    let v3 = |p: &[f64; 3]| Vec3::new(p[0], p[1], p[2]);
-    let mut pts: Vec<Vec3> = b.shells.iter().flat_map(|s| s.vertices.iter().map(v3)).collect();
-    if pts.is_empty() {
-        pts = b
-            .shells
-            .iter()
-            .flat_map(|s| s.faces.iter())
-            .filter_map(|f| match &f.surface {
-                Surface::Plane { origin, .. } | Surface::PlaneThrough { origin, .. } | Surface::Revolution { origin, .. } => Some(v3(origin)),
-                Surface::PlacedRevolution { origin, matrix: [x, y, z, w], .. } => {
-                    let column = |c: &[f64; 4]| Vec3::new(c[0], c[1], c[2]);
-                    Some(column(x) * origin[0] + column(y) * origin[1] + column(z) * origin[2] + column(w))
-                }
-                Surface::Cylinder { frame, .. } | Surface::Cone { frame, .. } | Surface::Sphere { frame, .. } | Surface::Torus { frame, .. } => {
-                    Some(v3(&frame.origin))
-                }
-                Surface::Spline { points, .. } => points.first().and_then(|r| r.first()).map(v3),
-                Surface::RationalSpline { homogeneous, .. } => {
-                    homogeneous.first().and_then(|r| r.first()).filter(|h| h[3] != 0.).map(|h| Vec3::new(h[0] / h[3], h[1] / h[3], h[2] / h[3]))
-                }
-            })
-            .collect();
-    }
-    let Some(first) = pts.first().copied() else { return Vec3::ZERO };
-    let (lo, hi) = pts.iter().fold((first, first), |(lo, hi), p| (lo.min(*p), hi.max(*p)));
+    let mut pts = b.vertices.iter().map(|v| Vec3::new(v.point[0], v.point[1], v.point[2]));
+    let Some(first) = pts.next() else { return Vec3::ZERO };
+    let (lo, hi) = pts.fold((first, first), |(lo, hi), p| (lo.min(p), hi.max(p)));
     let c = (lo + hi) * 0.5;
     if c.is_finite() { c } else { Vec3::ZERO }
 }
 
-/// The body's records with every position moved by `-d`.
+/// The body's records with every position moved by `-d`: points, control points and NURBS forms,
+/// a placed revolution through its matrix; regions, parameters and pcurves as they are.
 fn shifted(b: &Brep, d: Vec3) -> Brep {
-    use crate::brep::{Curve, Surface};
-    fn point(p: &mut [f64; 3], d: Vec3) {
-        p[0] -= d.x;
-        p[1] -= d.y;
-        p[2] -= d.z;
-    }
-    fn weighted(h: &mut [f64; 4], d: Vec3) {
-        h[0] -= d.x * h[3];
-        h[1] -= d.y * h[3];
-        h[2] -= d.z * h[3];
-    }
-    fn curve(c: &mut Curve, d: Vec3) {
-        match c {
-            Curve::Line { start, end } => {
-                point(start, d);
-                point(end, d);
-            }
-            Curve::CircleArc { frame, .. } | Curve::EllipseArc { frame, .. } => point(&mut frame.origin, d),
-            Curve::Polyline { points } | Curve::Spline { points, .. } | Curve::Intersection { points, .. } => {
-                points.iter_mut().for_each(|p| point(p, d))
-            }
-            Curve::RationalSpline { homogeneous, .. } => homogeneous.iter_mut().for_each(|h| weighted(h, d)),
-            Curve::SurfaceIntersection { surfaces, leader, .. } => {
-                surfaces.iter_mut().for_each(|s| surface(s, d));
-                curve(leader, d);
-            }
-        }
-    }
-    fn surface(s: &mut Surface, d: Vec3) {
-        match s {
-            Surface::Plane { origin, .. } => point(origin, d),
-            Surface::PlaneThrough { origin, u_point, v_point } => [origin, u_point, v_point].into_iter().for_each(|p| point(p, d)),
-            Surface::Cylinder { frame, .. } | Surface::Cone { frame, .. } | Surface::Sphere { frame, .. } | Surface::Torus { frame, .. } => {
-                point(&mut frame.origin, d)
-            }
-            Surface::Spline { points, .. } => points.iter_mut().flatten().for_each(|p| point(p, d)),
-            Surface::RationalSpline { homogeneous, .. } => homogeneous.iter_mut().flatten().for_each(|h| weighted(h, d)),
-            Surface::Revolution { profile, origin, .. } => {
-                curve(profile, d);
-                point(origin, d);
-            }
-            // The placement's last column carries its translation (and the homogeneous weight).
-            Surface::PlacedRevolution { matrix: [.., w], .. } => weighted(w, d),
-        }
-    }
-    let mut out = b.clone();
-    for s in &mut out.shells {
-        s.vertices.iter_mut().for_each(|p| point(p, d));
-        s.edges.iter_mut().for_each(|e| curve(&mut e.curve, d));
-        s.faces.iter_mut().for_each(|f| surface(&mut f.surface, d));
-    }
-    out
+    b.moved(&crate::brep::placement::Rigid { rotation: [[1., 0., 0.], [0., 1., 0.], [0., 0., 1.]], translation: [-d.x, -d.y, -d.z] })
 }
 
 #[derive(Clone, Copy, Debug)]

crates/geom/src/hlr/surface.rs

@@ -1,11 +1,13 @@
 //! Surfaces of portable B-rep faces with a parameter chart, inversion, ray intersection and
 //! analytic recognition. Planes are closed form; revolutions are reduced to their 2D profile;
-//! spline surfaces use Newton iterations seeded from a sample grid.
+//! spline surfaces use Newton iterations seeded from a sample grid. Records' surfaces are unbounded:
+//! a cylinder, cone, sphere, torus or revolution gets its chart's bounds from its face.
 
 use super::HlrError;
 use super::curve::{C3, as_circle, frame_axes, v3};
 use crate::Vec3;
-use crate::brep::{Frame, Surface};
+use crate::brep::eval::SurfaceEval;
+use crate::brep::{Curve, Face, Frame, Surface, form, placement::Rigid};
 use crate::bspline::{HPoint, NurbsSurface};
 
 const TAU: f64 = std::f64::consts::TAU;
@@ -116,8 +118,11 @@ fn perp(a: Vec3) -> Vec3 {
 }
 
 impl S3 {
-    pub fn new(s: &Surface) -> Result<S3, HlrError> {
-        Ok(match s {
+    /// The chart of a face: its surface, in its own parameters where it is a plane or a spline
+    /// and in angle and height (angle and angle) where it is a cylinder, cone, sphere or torus,
+    /// over the part its NURBS form spans ([`Face::domain`]).
+    pub fn of_face(face: &Face) -> Result<S3, HlrError> {
+        Ok(match &face.surface {
             Surface::Plane { origin, u_axis, v_axis } => {
                 let (u, v) = (v3(*u_axis), v3(*v_axis));
                 let n = u.cross(v).normalized().ok_or(HlrError::Unsupported("degenerate plane"))?;
@@ -130,39 +135,45 @@ impl S3 {
                 let n = u.cross(v).normalized().ok_or(HlrError::Unsupported("degenerate plane"))?;
                 S3::Plane { o, u, v, n }
             }
-            // A revolution mapped by an affine placement, which may scale it unevenly.
-            Surface::PlacedRevolution { .. } => return Err(HlrError::Unsupported("placed revolution surface")),
             Surface::Spline { u_knots, v_knots, points } => {
                 let h = points.iter().map(|r| r.iter().map(|p| [p[0], p[1], p[2], 1.]).collect()).collect::<Vec<Vec<_>>>();
                 S3::Nurbs(Box::new(Nurbs::new(u_knots, v_knots, &h)?))
             }
             Surface::RationalSpline { u_knots, v_knots, homogeneous } => S3::Nurbs(Box::new(Nurbs::new(u_knots, v_knots, homogeneous)?)),
-            Surface::Revolution { profile, origin, axis, angle } => {
-                let (c, _) = C3::new(profile)?;
-                S3::Rev(Box::new(Rev::new(v3(*origin), v3(*axis), None, c, *angle)?))
-            }
-            Surface::Cylinder { frame, radius, angle, height } => {
+            Surface::PlacedRevolution { .. } => placed_revolution(face)?,
+            label => S3::bounded(label, form_box(face, label)?)?,
+        })
+    }
+
+    /// A cylinder, cone, sphere or torus over `[angle, along]`: `along` the height on the frame's
+    /// axis of a cylinder or cone, the latitude of a sphere, a torus's tube angle from the frame
+    /// (its tube start counted in).
+    pub fn bounded(s: &Surface, [angle, along]: [[f64; 2]; 2]) -> Result<S3, HlrError> {
+        let rev = match s {
+            Surface::Cylinder { frame, radius } => {
                 let (o, x, _, z) = axes(frame);
-                let line = C3::Line { a: o + x * *radius + z * height[0], b: o + x * *radius + z * height[1] };
-                S3::Rev(Box::new(Rev::new(o, z, Some(x), line, *angle)?))
+                let line = C3::Line { a: o + x * *radius + z * along[0], b: o + x * *radius + z * along[1] };
+                Rev::new(o, z, Some(x), line, angle)?
             }
-            Surface::Cone { frame, radius_at_origin, half_angle, angle, height } => {
+            Surface::Cone { frame, radius_at_origin, half_angle } => {
                 let (o, x, _, z) = axes(frame);
                 let r = |h: f64| radius_at_origin + h * half_angle.tan();
-                let line = C3::Line { a: o + x * r(height[0]) + z * height[0], b: o + x * r(height[1]) + z * height[1] };
-                S3::Rev(Box::new(Rev::new(o, z, Some(x), line, *angle)?))
+                let line = C3::Line { a: o + x * r(along[0]) + z * along[0], b: o + x * r(along[1]) + z * along[1] };
+                Rev::new(o, z, Some(x), line, angle)?
             }
-            Surface::Sphere { frame, radius, longitude, latitude } => {
+            Surface::Sphere { frame, radius } => {
                 let (o, x, _, z) = axes(frame);
-                let arc = C3::Arc { o, x, y: z, r: [*radius, *radius], range: *latitude };
-                S3::Rev(Box::new(Rev::new(o, z, Some(x), arc, *longitude)?))
+                let arc = C3::Arc { o, x, y: z, r: [*radius, *radius], range: along };
+                Rev::new(o, z, Some(x), arc, angle)?
             }
-            Surface::Torus { frame, major_radius, minor_radius, longitude, meridian } => {
+            Surface::Torus { frame, major_radius, minor_radius, .. } => {
                 let (o, x, _, z) = axes(frame);
-                let arc = C3::Arc { o: o + x * *major_radius, x, y: z, r: [*minor_radius, *minor_radius], range: *meridian };
-                S3::Rev(Box::new(Rev::new(o, z, Some(x), arc, *longitude)?))
+                let arc = C3::Arc { o: o + x * *major_radius, x, y: z, r: [*minor_radius, *minor_radius], range: along };
+                Rev::new(o, z, Some(x), arc, angle)?
             }
-        })
+            _ => return Err(HlrError::Unsupported("surface chart")),
+        };
+        Ok(S3::Rev(Box::new(rev)))
     }
 
     /// The recognised quadric or torus as a revolution chart: closed-form inversion and exact
@@ -525,6 +536,123 @@ pub(crate) fn torus_roots(w: Vec3, d: Vec3, a: Vec3, big: f64, small: f64, range
     out
 }
 
+/// Steps per direction when a face's NURBS form is read back onto its label: along a direction
+/// that turns (at most a turn, in rational quarter arcs) each step stays well within half a turn.
+const FORM_STEPS: usize = 16;
+
+/// The angle `x` on the turn nearest `near`.
+fn unwrapped(x: f64, near: f64) -> f64 {
+    x + TAU * ((near - x) / TAU).round()
+}
+
+/// A cylinder, cone, sphere or torus face's box on its label ([`S3::bounded`]): where its form's
+/// points over its domain lie on the label's frame, whichever way round the form runs.
+fn form_box(face: &Face, label: &Surface) -> Result<[[f64; 2]; 2], HlrError> {
+    let frame = match label {
+        Surface::Cylinder { frame, .. } | Surface::Cone { frame, .. } | Surface::Sphere { frame, .. } | Surface::Torus { frame, .. } => frame,
+        _ => return Err(HlrError::Unsupported("surface chart")),
+    };
+    let (o, x, y, z) = axes(frame);
+    let along = |w: Vec3, rho: f64| match label {
+        Surface::Sphere { .. } => w.dot(z).atan2(rho),
+        Surface::Torus { major_radius, .. } => w.dot(z).atan2(rho - major_radius),
+        _ => w.dot(z),
+    };
+    let turns = matches!(label, Surface::Torus { .. });
+    let (angle, mut second) = sampled_box(face, turns, &|p| {
+        let w = p - o;
+        let (a, b) = (w.dot(x), w.dot(y));
+        let rho = a.hypot(b);
+        [rho, b.atan2(a), along(w, rho)]
+    })?;
+    if !(second[0] < second[1]) {
+        return Err(HlrError::Invalid("surface form domain"));
+    }
+    if turns && second[1] - second[0] >= TAU - 1e-9 {
+        second[1] = second[0] + TAU;
+    }
+    Ok([angle, second])
+}
+
+/// The range of the angle about the axis and of a second coordinate over a form's points on a grid
+/// of the face's domain; `at` gives a point's distance from the axis, angle and second coordinate,
+/// the second an angle too when `turns`. Angles are unwrapped from point to point, row after row, so
+/// a face that goes round comes to one turn. On the axis (a pole, an apex) there is no angle.
+fn sampled_box(face: &Face, turns: bool, at: &dyn Fn(Vec3) -> [f64; 3]) -> Result<([f64; 2], [f64; 2]), HlrError> {
+    let form = face.nurbs.as_ref().ok_or(HlrError::Invalid("analytic surface without its form"))?;
+    let eval = SurfaceEval::new(form).map_err(|_| HlrError::Invalid("surface form"))?;
+    let [[u0, u1], [v0, v1]] = face.domain;
+    let step = |a: f64, b: f64, k: usize| a + (b - a) * k as f64 / FORM_STEPS as f64;
+    let mut samples = Vec::with_capacity((FORM_STEPS + 1) * (FORM_STEPS + 1));
+    for j in 0..=FORM_STEPS {
+        for i in 0..=FORM_STEPS {
+            let p = eval.point([step(u0, u1, i), step(v0, v1, j)]).filter(|p| p.is_finite()).ok_or(HlrError::Invalid("surface form"))?;
+            samples.push(at(p));
+        }
+    }
+    let widest = samples.iter().fold(0f64, |m, s| m.max(s[0]));
+    let (mut angle, mut second) = ([f64::INFINITY, f64::NEG_INFINITY], [f64::INFINITY, f64::NEG_INFINITY]);
+    // Each row starts next to the start of the one before.
+    let (mut row_a, mut row_s): (Option<f64>, Option<f64>) = (None, None);
+    for row in samples.chunks(FORM_STEPS + 1) {
+        let (mut last_a, mut last_s) = (row_a, row_s);
+        for (i, &[rho, a, s]) in row.iter().enumerate() {
+            if rho > 1e-6 * widest {
+                let a = last_a.map_or(a, |l| unwrapped(a, l));
+                angle = [angle[0].min(a), angle[1].max(a)];
+                last_a = Some(a);
+                if i == 0 || row_a.is_none() {
+                    row_a = Some(a);
+                }
+            }
+            let s = if turns { last_s.map_or(s, |l| unwrapped(s, l)) } else { s };
+            second = [second[0].min(s), second[1].max(s)];
+            last_s = Some(s);
+            if i == 0 {
+                row_s = Some(s);
+            }
+        }
+    }
+    if !(angle[0] < angle[1]) {
+        return Err(HlrError::Invalid("surface form domain"));
+    }
+    Ok((angle, second))
+}
+
+/// A placed revolution's chart. Its placement must be rigid: then it is a revolution of the placed
+/// profile about the placed axis. Any other (which may scale it unevenly) is refused. The profile's
+/// range is `domain` along it: its own parameter, a conic's on its form's quarter grid.
+fn placed_revolution(face: &Face) -> Result<S3, HlrError> {
+    let Surface::PlacedRevolution { profile, origin, axis, matrix, forward } = &face.surface else {
+        return Err(HlrError::Unsupported("surface chart"));
+    };
+    let Some(rigid) = Rigid::of_matrix(matrix) else { return Err(HlrError::Unsupported("placed revolution surface")) };
+    let along = if *forward { face.domain[0] } else { face.domain[1] };
+    let range = match profile.as_ref() {
+        Curve::Circle { .. } | Curve::Ellipse { .. } => [form::angle(along[0]), form::angle(along[1])],
+        _ => along,
+    };
+    if matches!(profile.as_ref(), Curve::Circle { .. } | Curve::Ellipse { .. }) && !(range[1] - range[0] <= TAU * (1. + 1e-9)) {
+        return Err(HlrError::Invalid("revolution profile range"));
+    }
+    let rows = rigid.rotation.map(v3);
+    let direction = |d: Vec3| Vec3::new(rows[0].dot(d), rows[1].dot(d), rows[2].dot(d));
+    let point = |p: Vec3| direction(p) + v3(rigid.translation);
+    let curve = C3::bounded(profile, range)?.placed(&point, &direction)?;
+    let (o, a) = (point(v3(*origin)), direction(v3(*axis)).normalized().ok_or(HlrError::Unsupported("zero revolution axis"))?);
+    let [t0, _] = curve.domain();
+    let start = curve.at(t0) - o;
+    let e1 = (start - a * start.dot(a)).normalized().unwrap_or_else(|| perp(a));
+    let e2 = a.cross(e1);
+    let (angle, _) = sampled_box(face, false, &|p| {
+        let w = p - o;
+        let h = w.dot(a);
+        let radial = w - a * h;
+        [radial.len(), radial.dot(e2).atan2(radial.dot(e1)), h]
+    })?;
+    Ok(S3::Rev(Box::new(Rev::new(o, a, Some(e1), curve, angle)?)))
+}
+
 fn axes(f: &Frame) -> (Vec3, Vec3, Vec3, Vec3) {
     let (o, x, y) = frame_axes(f);
     (o, x, y, x.cross(y))
@@ -828,12 +956,9 @@ impl Rev {
         let r = radial.len();
         let theta = if r > 1e-12 * (1. + self.rmax) { radial.dot(self.e2).atan2(radial.dot(self.e1)) } else { hint.map_or(0., |u| u[0]) };
         let (s, _) = self.closest([r, h]);
-        let mut t = theta;
-        if self.theta[1] - self.theta[0] >= TAU - 1e-9 {
-            t = t.rem_euclid(TAU);
-        } else if t < self.theta[0] {
-            t += TAU;
-        }
+        // Of a part of a turn, the copy nearest its middle: an end rounds either way.
+        let t =
+            if self.theta[1] - self.theta[0] >= TAU - 1e-9 { theta.rem_euclid(TAU) } else { unwrapped(theta, 0.5 * (self.theta[0] + self.theta[1])) };
         Some([t, s])
     }
 

crates/geom/src/hlr/trim.rs

@@ -6,7 +6,8 @@
 use super::curve::C3;
 use super::surface::S3;
 use crate::Vec3;
-use crate::brep::Face;
+use crate::brep::eval::SurfaceEval;
+use crate::brep::{Face, Loop, OrientedEdge};
 
 #[derive(Clone, Debug)]
 pub(crate) struct Trim {
@@ -159,6 +160,27 @@ fn index(segs: &[([f64; 2], [f64; 2])], period: [Option<f64>; 2], r: [f64; 2]) -
     Index { copies, vcross, buckets, v0, dv }
 }
 
+/// Every use of an edge in a face's loops.
+pub(crate) fn uses(face: &Face) -> impl Iterator<Item = &OrientedEdge> {
+    face.loops.iter().flat_map(|l| match l {
+        Loop::Edges { edges } => edges.as_slice(),
+        Loop::Vertex { .. } => &[],
+    })
+}
+
+/// +1 where the chart's normal is the natural normal of the face's parameters, which its loops
+/// are stored against, and -1 where it is the other side (a form turned against its label, a
+/// recognised patch's own chart): the side the face is on, walking its loops in the chart.
+fn loops_left(face: &Face, surface: &S3) -> f64 {
+    let [[u0, u1], [v0, v1]] = face.domain;
+    let natural = SurfaceEval::new(face.parameters()).ok().and_then(|s| s.derivs(0.5 * (u0 + u1), 0.5 * (v0 + v1)));
+    let along = natural.and_then(|(p, su, sv)| {
+        let chart = surface.inverse(p, None).and_then(|uv| surface.normal(uv))?;
+        Some(su.cross(sv).dot(chart))
+    });
+    if along.is_some_and(|a| a < 0.) { -1. } else { 1. }
+}
+
 fn lift(x: f64, near: f64, p: Option<f64>) -> f64 {
     match p {
         Some(p) if p > 0. => x + p * ((near - x) / p).round(),
@@ -173,15 +195,17 @@ pub(crate) fn build(face: &Face, edges: &[C3], surface: &S3, per_break: usize) -
     let mut chords = Vec::new();
     let mut chord_error = 0f64;
     let (mut lo, mut hi) = (Vec3::new(f64::MAX, f64::MAX, f64::MAX), Vec3::new(f64::MIN, f64::MIN, f64::MIN));
-    let mut uses = std::collections::BTreeMap::new();
-    for u in face.loops.iter().flatten() {
-        *uses.entry(u.edge).or_insert(0usize) += 1;
+    let mut counts = std::collections::BTreeMap::new();
+    for u in self::uses(face) {
+        *counts.entry(u.edge).or_insert(0usize) += 1;
     }
     let mut area_u = 0.; // -loop integral of v du
     let mut area_v = 0.; // loop integral of u dv
     let mut wraps = [false; 2];
     let mut unsupported = None;
     for lp in &face.loops {
+        // A vertex loop (a pole) bounds nothing in the chart.
+        let Loop::Edges { edges: lp } = lp else { continue };
         let mut pts: Vec<([f64; 2], bool)> = Vec::new();
         let mut hint: Option<[f64; 2]> = None;
         for u in lp {
@@ -189,7 +213,7 @@ pub(crate) fn build(face: &Face, edges: &[C3], surface: &S3, per_break: usize) -
                 unsupported = Some("edge reference");
                 continue;
             };
-            let is_seam = uses.get(&u.edge).copied().unwrap_or(0) > 1;
+            let is_seam = counts.get(&u.edge).copied().unwrap_or(0) > 1;
             let mut ts = Vec::new();
             let b = c.breaks();
             for w in b.windows(2) {
@@ -253,11 +277,11 @@ pub(crate) fn build(face: &Face, edges: &[C3], surface: &S3, per_break: usize) -
     };
     // Loops that wind once around a periodic direction without coming back (a hemisphere's
     // equator) close against the chart's own edge, so their area depends on where the chart
-    // starts: there the face is taken on the loops' left (the B-rep convention), and below the
-    // side within the surface's domain wins either way.
+    // starts: there the face is taken on the loops' left seen from the natural normal (the B-rep
+    // convention), and below the side within the surface's domain wins either way.
     let net = segs.iter().fold([0f64; 2], |n, (a, b)| [n[0] + b[0] - a[0], n[1] + b[1] - a[1]]);
     let open = (0..2).any(|k| period[k].is_some_and(|p| net[k].abs() > 0.5 * p));
-    let sigma = if open { 1. } else { area.signum() };
+    let sigma = if open { loops_left(face, surface) } else { area.signum() };
     let dom = surface.domain();
     // Inside point: beside the longest non-seam segment, on the face's side.
     let mut inside = None;