fix(brep): standard loop winding, enforced by validate()

main/opengeometry/src/brep/builder/brep_builder.rs

@@ -193,6 +193,7 @@ impl Builder {
             outer_shell: 0,
             cavity_shells: Vec::new(),
         });
+        self.brep.orient_loops()?;
         self.brep.validate()?;
         Ok(self.brep)
     }

main/opengeometry/src/brep/topology/envelope.rs

@@ -61,8 +61,10 @@ impl BrepEnvelope {
                 "BRep JSON exceeds 64 MiB".into(),
             ));
         }
-        let brep: Self = serde_json::from_str(json)
+        let mut brep: Self = serde_json::from_str(json)
             .map_err(|e| GeometryError::InvalidGeometry(format!("invalid BRep JSON: {e}")))?;
+        brep.validate_structure()?;
+        brep.orient_loops()?;
         brep.validate()?;
         Ok(brep)
     }

main/opengeometry/src/brep/validate/mod.rs

@@ -1,5 +1,8 @@
 mod edges;
 mod face_use;
+#[cfg(test)]
+mod tests;
+mod winding;
 
 use super::error::{invalid, missing, positive, GeometryError};
 use super::geometry::{CurveGeometry, Surface};
@@ -10,6 +13,11 @@ use std::collections::HashSet;
 
 impl BrepEnvelope {
     pub fn validate(&self) -> Result<(), GeometryError> {
+        self.validate_structure()?;
+        self.check_loop_winding()
+    }
+
+    pub(crate) fn validate_structure(&self) -> Result<(), GeometryError> {
         self.check_header()?;
         self.check_geometry()?;
         self.check_pcurves()?;

main/opengeometry/src/brep/validate/tests.rs

@@ -0,0 +1,96 @@
+use crate::brep::{
+    coarse_accuracy, BrepEnvelope, CurveGeometry, EdgeGeometry, Frame3, GeometryError, Orientation,
+};
+use crate::primitives;
+
+fn tube() -> BrepEnvelope {
+    primitives::annular_cylinder(
+        "tube".into(),
+        Frame3::IDENTITY,
+        0.4,
+        1.0,
+        2.0,
+        coarse_accuracy(1e-6),
+    )
+    .unwrap()
+}
+
+fn tube_with_old_winding() -> BrepEnvelope {
+    let mut brep = tube();
+    let mut loops = Vec::new();
+    for face in &brep.topology.faces {
+        loops.extend(&face.trim.holes);
+        if face.key == "inner" {
+            loops.push(face.trim.outer);
+        }
+    }
+    brep.reverse_loops(&loops).unwrap();
+    brep
+}
+
+fn circular_turn(brep: &BrepEnvelope, face: u32, loop_id: u32) -> f64 {
+    let use_ =
+        &brep.topology.halfedges[brep.topology.loops[loop_id as usize].start_halfedge as usize];
+    assert_eq!(use_.next, Some(use_.id), "one edge per cap loop");
+    let edge = &brep.topology.edges[use_.edge as usize];
+    let EdgeGeometry::Curve { curve, .. } = edge.geometry else {
+        panic!("cap loop edge is a curve")
+    };
+    let CurveGeometry::Circle { frame, .. } = brep.geometry.curves[curve as usize] else {
+        panic!("cap loop edge is a circle")
+    };
+    let point = brep.topology.vertices[use_.from as usize].position;
+    let normal = brep.face_normal_at(face, point).unwrap();
+    let along: f64 = (0..3).map(|i| frame.z[i] * normal[i]).sum();
+    along.signum() * use_.geometry_use.sense.multiplier()
+}
+
+#[test]
+fn cap_holes_run_clockwise_and_outer_loops_counter_clockwise() {
+    let brep = tube();
+    brep.validate().unwrap();
+    for face in brep
+        .topology
+        .faces
+        .iter()
+        .filter(|f| f.key.ends_with("cap"))
+    {
+        assert_eq!(circular_turn(&brep, face.id, face.trim.outer), 1.0);
+        for &hole in &face.trim.holes {
+            assert_eq!(circular_turn(&brep, face.id, hole), -1.0, "{}", face.key);
+        }
+    }
+}
+
+#[test]
+fn validate_rejects_loops_against_the_material_normal() {
+    let err = tube_with_old_winding().validate().unwrap_err();
+    assert!(
+        matches!(&err, GeometryError::InvalidTopology(m) if m.contains("material normal")),
+        "{err:?}"
+    );
+}
+
+#[test]
+fn a_flipped_face_sense_with_its_loop_unturned_is_rejected() {
+    let mut brep = primitives::cuboid(
+        "box".into(),
+        Frame3::IDENTITY,
+        [1.0, 1.0, 1.0],
+        coarse_accuracy(1e-6),
+    )
+    .unwrap();
+    for face in &mut brep.topology.faces {
+        face.sense = Orientation::Reverse;
+    }
+    assert!(brep.validate().is_err());
+    brep.orient_loops().unwrap();
+    brep.validate().unwrap();
+}
+
+#[test]
+fn bodies_with_the_old_winding_load_with_their_loops_turned() {
+    let old = serde_json::to_string(&tube_with_old_winding()).unwrap();
+    let loaded = BrepEnvelope::from_json(&old).unwrap();
+    assert_eq!(loaded.topology.halfedges, tube().topology.halfedges);
+}

main/opengeometry/src/brep/validate/winding.rs

@@ -0,0 +1,215 @@
+use crate::brep::error::{invalid, GeometryError};
+use crate::brep::frame::UV;
+use crate::brep::geometry::Surface;
+use crate::brep::pcurve::{period_lifted_uv, PcurveGeometry};
+use crate::brep::topology::{BrepEnvelope, Face, Orientation};
+use crate::math::{cross, dot};
+
+const SAMPLES_PER_USE: usize = 16;
+
+struct Winding {
+    id: u32,
+    hole: bool,
+    area: f64,
+    delta: UV,
+    mean: UV,
+}
+
+impl BrepEnvelope {
+    pub(crate) fn orient_loops(&mut self) -> Result<(), GeometryError> {
+        let mut flips = Vec::new();
+        for face in &self.topology.faces {
+            flips.extend(self.misoriented_loops(face)?);
+        }
+        if flips.is_empty() {
+            return Ok(());
+        }
+        self.reverse_loops(&flips)
+    }
+
+    pub(super) fn reverse_loops(&mut self, loops: &[u32]) -> Result<(), GeometryError> {
+        for &loop_id in loops {
+            let start = self.topology.loops[loop_id as usize].start_halfedge;
+            let mut h = start;
+            for _ in 0..=self.topology.halfedges.len() {
+                let use_ = &mut self.topology.halfedges[h as usize];
+                use_.reverse();
+                h = use_.prev.ok_or_else(|| invalid("nonterminating loop"))?;
+                if h == start {
+                    break;
+                }
+            }
+        }
+        let mut outgoing = vec![None; self.topology.vertices.len()];
+        for use_ in &self.topology.halfedges {
+            outgoing[use_.from as usize].get_or_insert(use_.id);
+        }
+        let halfedges = &self.topology.halfedges;
+        for vertex in &mut self.topology.vertices {
+            if let Some(h) = vertex.outgoing_halfedge {
+                if halfedges[h as usize].from != vertex.id {
+                    vertex.outgoing_halfedge = outgoing[vertex.id as usize];
+                }
+            }
+        }
+        Ok(())
+    }
+
+    pub(super) fn check_loop_winding(&self) -> Result<(), GeometryError> {
+        for face in &self.topology.faces {
+            if !self.misoriented_loops(face)?.is_empty() {
+                return Err(invalid(format!(
+                    "face {} ('{}') has loops running against its material normal",
+                    face.id, face.key
+                )));
+            }
+        }
+        Ok(())
+    }
+
+    fn misoriented_loops(&self, face: &Face) -> Result<Vec<u32>, GeometryError> {
+        let Some(material) = self.material_sign(face)? else {
+            return Ok(Vec::new());
+        };
+        let periods = self.face_periods(face)?;
+        let mut windings = Vec::new();
+        for &loop_id in std::iter::once(&face.trim.outer).chain(&face.trim.holes) {
+            windings.push(self.winding(loop_id, periods)?);
+        }
+        let wraps = |w: &Winding, axis: usize| {
+            periods[axis].is_some_and(|period| w.delta[axis].abs() > 0.5 * period)
+        };
+        let mut flips = Vec::new();
+        let mut wrapping = Vec::new();
+        for w in &windings {
+            if wraps(w, 0) || wraps(w, 1) {
+                wrapping.push(w);
+            } else if (w.area * material > 0.0) == w.hole {
+                flips.push(w.id);
+            }
+        }
+        if let [a, b] = wrapping[..] {
+            let axis = match (wraps(a, 0) && wraps(b, 0), wraps(a, 1) && wraps(b, 1)) {
+                (true, false) => Some(0),
+                (false, true) => Some(1),
+                _ => None,
+            };
+            if let Some(axis) = axis {
+                flips.extend(band_flips(a, b, axis, material));
+            }
+        }
+        Ok(flips)
+    }
+
+    fn material_sign(&self, face: &Face) -> Result<Option<f64>, GeometryError> {
+        let surface = self.geometry.surface(face.surface)?;
+        let [u, v] = face.trim.uv_bounds;
+        let centre = [(u.lo + u.hi) * 0.5, (v.lo + v.hi) * 0.5];
+        let Ok(normal) = surface.normal_at(centre) else {
+            return Ok(None);
+        };
+        let jet = surface.derivatives(centre)?;
+        let chart = dot(cross(jet.du, jet.dv), normal);
+        if chart == 0.0 {
+            return Ok(None);
+        }
+        Ok(Some(chart.signum() * face.sense.multiplier()))
+    }
+
+    fn face_periods(&self, face: &Face) -> Result<[Option<f64>; 2], GeometryError> {
+        let surface = self.geometry.surface(face.surface)?;
+        Ok(surface.charts()[face.trim.chart as usize].periods)
+    }
+
+    fn winding(&self, loop_id: u32, periods: [Option<f64>; 2]) -> Result<Winding, GeometryError> {
+        let points = self.loop_samples(loop_id, periods)?;
+        let mut area = 0.0;
+        let mut mean = [0.0; 2];
+        let count = (points.len() - 1) as f64;
+        for pair in points.windows(2) {
+            area -= 0.5 * (pair[0][1] + pair[1][1]) * (pair[1][0] - pair[0][0]);
+            for axis in 0..2 {
+                mean[axis] += pair[0][axis] / count;
+            }
+        }
+        let (first, last) = (points[0], points[points.len() - 1]);
+        Ok(Winding {
+            id: loop_id,
+            hole: self.topology.loops[loop_id as usize].is_hole,
+            area,
+            delta: [last[0] - first[0], last[1] - first[1]],
+            mean,
+        })
+    }
+
+    fn loop_samples(
+        &self,
+        loop_id: u32,
+        periods: [Option<f64>; 2],
+    ) -> Result<Vec<UV>, GeometryError> {
+        let start = self.topology.loops[loop_id as usize].start_halfedge;
+        let mut points: Vec<UV> = Vec::new();
+        let mut h = start;
+        for _ in 0..=self.topology.halfedges.len() {
+            let use_ = &self.topology.halfedges[h as usize];
+            let range = self.topology.edges[use_.edge as usize].geometry.range();
+            let pcurve = use_
+                .geometry_use
+                .pcurve
+                .ok_or_else(|| invalid("missing face pcurve"))?;
+            let straight = matches!(
+                self.geometry.pcurves[pcurve as usize],
+                PcurveGeometry::Line2 { .. }
+            );
+            let samples = if straight && periods == [None, None] {
+                1
+            } else {
+                SAMPLES_PER_USE
+            };
+            for i in 0..=samples {
+                let s = i as f64 / samples as f64;
+                let s = match use_.geometry_use.sense {
+                    Orientation::Forward => s,
+                    Orientation::Reverse => 1.0 - s,
+                };
+                let uv = self
+                    .geometry
+                    .pcurve_at(pcurve, range.lo + range.width() * s)?;
+                let mut uv = period_lifted_uv(uv, periods, use_.geometry_use.periodic_lift);
+                if let Some(previous) = points.last() {
+                    for axis in 0..2 {
+                        if let Some(period) = periods[axis] {
+                            uv[axis] += ((previous[axis] - uv[axis]) / period).round() * period;
+                        }
+                    }
+                }
+                if i > 0 || points.is_empty() {
+                    points.push(uv);
+                }
+            }
+            h = use_.next.ok_or_else(|| invalid("nonterminating loop"))?;
+            if h == start {
+                return Ok(points);
+            }
+        }
+        Err(invalid("nonterminating loop"))
+    }
+}
+
+fn band_flips(a: &Winding, b: &Winding, axis: usize, material: f64) -> Vec<u32> {
+    let other = 1 - axis;
+    let (low, high) = if a.mean[other] <= b.mean[other] {
+        (a, b)
+    } else {
+        (b, a)
+    };
+    let forward = if axis == 0 { 1.0 } else { -1.0 };
+    let mut flips = Vec::new();
+    if low.delta[axis] * material * forward < 0.0 {
+        flips.push(low.id);
+    }
+    if high.delta[axis] * material * forward > 0.0 {
+        flips.push(high.id);
+    }
+    flips
+}