geop_core_topology/
mesh.rs

1//! A body of triangles ([`BodySpec::of_triangles`]): what a mesh read from
2//! a file (STL) is as a B-rep.
3//!
4//! The topology comes from the mesh's own indexing, not from stitching
5//! geometry: two triangles share an edge where they share its two corners.
6//! Each triangle is a face on its plane. The plane is a patch over a domain
7//! the triangle does not decide — the box of every corner of the mesh seen
8//! along the axis the plane faces most — so triangles in one plane lie on
9//! the very same patch, and [`BodySpec::merged`] makes them one face at
10//! once.
11
12use geop_core_geometry::nurb_surface::NurbSurface3D;
13use geop_core_math::{
14    geop_error::{GeopError, GeopResult},
15    scalars::Scalar,
16    vector::{Vector2, Vector3},
17};
18
19use crate::{
20    Curve2, Curve3, Sense,
21    build::{BodySpec, CoedgeOn, CoedgeSpec, EdgeSpec, FaceSpec},
22};
23
24impl BodySpec {
25    /// The body the triangles `triangles` make, each three indices into
26    /// `points`, counter-clockwise seen from the side its face points to:
27    /// a solid where every edge is shared by two triangles, a sheet where
28    /// some are by one. A triangle whose corners repeat an index has no
29    /// area and no edge of its own, and is left out. Refused, naming the
30    /// triangle: corners in one line, and an edge shared by more than two
31    /// triangles or by two running it the same way.
32    pub fn of_triangles(points: &[Vector3], triangles: &[[usize; 3]]) -> GeopResult<BodySpec> {
33        let triangles: Vec<(usize, [usize; 3])> = (triangles.iter().copied().enumerate())
34            .filter(|(_, [a, b, c])| a != b && b != c && c != a)
35            .collect();
36        if triangles.is_empty() {
37            return Err(GeopError::new("the mesh has no triangles"));
38        }
39        if let Some((t, _)) =
40            (triangles.iter()).find(|(_, corners)| corners.iter().any(|&p| p >= points.len()))
41        {
42            return Err(GeopError::new(format!(
43                "triangle {t} names a corner the mesh does not have"
44            )));
45        }
46        // The box of every corner, which each plane's patch spans.
47        let lo = Vector3::from_array(std::array::from_fn(|k| {
48            points
49                .iter()
50                .map(|p| p[k].lower())
51                .fold(points[0][k].lower(), |a, b| {
52                    if b.definitely_less(a) { b } else { a }
53                })
54        }));
55        let hi = Vector3::from_array(std::array::from_fn(|k| {
56            points
57                .iter()
58                .map(|p| p[k].upper())
59                .fold(points[0][k].upper(), |a, b| {
60                    if b.definitely_greater(a) { b } else { a }
61                })
62        }));
63        let mut spec = BodySpec::new(false);
64        spec.vertices = points.to_vec();
65        let mut edge_of: std::collections::HashMap<(usize, usize), usize> = Default::default();
66        let mut senses: Vec<Vec<Sense>> = Vec::new();
67        for &(t, [a, b, c]) in &triangles {
68            let corners = [points[a], points[b], points[c]];
69            let normal = corners[1]
70                .sub(&corners[0])
71                .prod_cross(&corners[2].sub(&corners[0]));
72            if normal.norm_sq().could_be_equal(Scalar::ZERO) {
73                return Err(GeopError::new(format!(
74                    "triangle {t} has its corners in one line: it has no plane"
75                )));
76            }
77            // Seen along the axis the plane faces most — a free choice of
78            // parametrization; any it faces at all would do.
79            let k = (0..3)
80                .max_by(|&i, &j| {
81                    normal[i]
82                        .abs()
83                        .to_f64()
84                        .total_cmp(&normal[j].abs().to_f64())
85                })
86                .expect("three axes");
87            let (mut i, mut j) = ((k + 1) % 3, (k + 2) % 3);
88            if normal[k].definitely_less(Scalar::ZERO) {
89                std::mem::swap(&mut i, &mut j);
90            }
91            // The plane: `x_k = (n . p0 - n_i u - n_j v) / n_k`.
92            let d = normal.prod_dot(&corners[0]);
93            let point = |u: Scalar, v: Scalar| -> GeopResult<Vector3> {
94                let mut p = Vector3::zero();
95                p[i] = u;
96                p[j] = v;
97                p[k] = d
98                    .sub(normal[i].mul(u))
99                    .sub(normal[j].mul(v))
100                    .div(normal[k])?;
101                Ok(p)
102            };
103            let surface = NurbSurface3D::bilinear(
104                point(lo[i], lo[j])?,
105                point(hi[i], lo[j])?,
106                point(hi[i], hi[j])?,
107                point(lo[i], hi[j])?,
108            )?;
109            // Its `(u, v)` run over `[0, 1]`, the box's.
110            let uv = |p: &Vector3| -> GeopResult<Vector2> {
111                Ok(Vector2::from_array([
112                    p[i].sub(lo[i]).div(hi[i].sub(lo[i]))?,
113                    p[j].sub(lo[j]).div(hi[j].sub(lo[j]))?,
114                ]))
115            };
116            let mut outer = Vec::with_capacity(3);
117            for (from, to) in [(a, b), (b, c), (c, a)] {
118                let key = (from.min(to), from.max(to));
119                let e = *edge_of.entry(key).or_insert_with(|| {
120                    spec.edges.push(EdgeSpec {
121                        curve: Curve3::line(points[key.0], points[key.1])
122                            .expect("two corners apart"),
123                        start: key.0,
124                        end: key.1,
125                    });
126                    senses.push(Vec::new());
127                    spec.edges.len() - 1
128                });
129                let sense = if from == key.0 {
130                    Sense::Forward
131                } else {
132                    Sense::Reversed
133                };
134                if senses[e].len() >= 2 || senses[e].contains(&sense) {
135                    return Err(GeopError::new(format!(
136                        "triangle {t} shares its edge from corner {from} to {to} with {} other triangle(s), {}: the mesh is no surface there",
137                        senses[e].len(),
138                        if senses[e].contains(&sense) {
139                            "running it the same way"
140                        } else {
141                            "more than one"
142                        }
143                    )));
144                }
145                senses[e].push(sense);
146                outer.push(CoedgeSpec {
147                    on: CoedgeOn::Edge(e, sense),
148                    pcurve: Curve2::line(uv(&points[from])?, uv(&points[to])?)?,
149                });
150            }
151            spec.faces.push(FaceSpec {
152                surface,
153                outer,
154                holes: Vec::new(),
155            });
156        }
157        spec.solid = senses.iter().all(|s| s.len() == 2);
158        spec.shells = vec![(0..spec.faces.len()).collect()];
159        // Corners no triangle uses are no vertices.
160        let mut used = vec![false; points.len()];
161        for edge in &spec.edges {
162            used[edge.start] = true;
163            used[edge.end] = true;
164        }
165        let mut index = vec![usize::MAX; points.len()];
166        let mut vertices = Vec::new();
167        for (p, &u) in used.iter().enumerate() {
168            if u {
169                index[p] = vertices.len();
170                vertices.push(points[p]);
171            }
172        }
173        for edge in &mut spec.edges {
174            edge.start = index[edge.start];
175            edge.end = index[edge.end];
176        }
177        spec.vertices = vertices;
178        spec.check()?;
179        Ok(spec)
180    }
181}
182
183#[cfg(test)]
184mod tests {
185    use geop_core_math::{scalars::Scalar, vector::Vector3};
186
187    use crate::{
188        Model,
189        build::BodySpec,
190        validation::{ValidationParameters, validate},
191    };
192
193    /// The unit cube as twelve triangles, two per side, outward.
194    fn cube() -> (Vec<Vector3>, Vec<[usize; 3]>) {
195        let points = (0..8)
196            .map(|k| {
197                Vector3::from_array(
198                    [k & 1, (k >> 1) & 1, (k >> 2) & 1].map(|c| Scalar::from_f64(c as f64)),
199                )
200            })
201            .collect();
202        let quads = [
203            [0, 2, 3, 1],
204            [4, 5, 7, 6],
205            [0, 1, 5, 4],
206            [2, 6, 7, 3],
207            [0, 4, 6, 2],
208            [1, 3, 7, 5],
209        ];
210        let triangles = quads
211            .iter()
212            .flat_map(|[a, b, c, d]| [[*a, *b, *c], [*a, *c, *d]])
213            .collect();
214        (points, triangles)
215    }
216
217    /// Twelve triangles of a cube: a valid solid, and merged, six faces.
218    #[test]
219    fn a_cube_of_triangles_is_a_box() {
220        let (points, triangles) = cube();
221        let spec = BodySpec::of_triangles(&points, &triangles).unwrap();
222        assert!(spec.solid);
223        assert_eq!(spec.faces.len(), 12);
224        let mut model = Model::new();
225        model.build_body(spec.clone()).unwrap();
226        validate(&ValidationParameters::default(), &model).unwrap();
227        let (merged, _) = spec.merged(|_| true, |_| true).unwrap();
228        assert_eq!(merged.faces.len(), 6);
229        assert_eq!(merged.edges.len(), 12);
230        let mut model = Model::new();
231        model.build_body(merged).unwrap();
232        validate(&ValidationParameters::default(), &model).unwrap();
233    }
234
235    /// A triangle left out: an open box, a sheet.
236    #[test]
237    fn a_cube_missing_a_triangle_is_a_sheet() {
238        let (points, mut triangles) = cube();
239        triangles.pop();
240        let spec = BodySpec::of_triangles(&points, &triangles).unwrap();
241        assert!(!spec.solid);
242    }
243
244    /// A triangle turned the other way shares edges running them the same
245    /// way as its neighbours: refused, naming it.
246    #[test]
247    fn a_triangle_turned_round_is_refused() {
248        let (points, mut triangles) = cube();
249        triangles[11].swap(1, 2);
250        let error = BodySpec::of_triangles(&points, &triangles).unwrap_err();
251        assert!(format!("{error}").contains("same way"), "{error}");
252    }
253}
geop_core_topology/mesh.rs

geop_core_topology/
mesh.rs
