scirs2_spatial/mesh/
mod.rs

1//! Triangle mesh operations for spatial computing
2//!
3//! This module provides a comprehensive triangle mesh representation and operations
4//! including simplification, smoothing, normal computation, quality metrics, and I/O.
5//!
6//! # Features
7//!
8//! * **Mesh Representation**: Indexed face set with vertex/face adjacency
9//! * **Simplification**: Edge collapse with quadric error metrics (QEM)
10//! * **Smoothing**: Laplacian and Taubin smoothing
11//! * **Normals**: Vertex and face normal computation
12//! * **Quality Metrics**: Aspect ratio, minimum angle per triangle
13//! * **I/O**: ASCII STL import/export
14
15mod quality;
16mod simplification;
17mod smoothing;
18
19use crate::error::{SpatialError, SpatialResult};
20use std::collections::{HashMap, HashSet};
21use std::fmt;
22
23/// A 3D vertex
24#[derive(Debug, Clone, Copy, PartialEq)]
25pub struct Vertex {
26    pub x: f64,
27    pub y: f64,
28    pub z: f64,
29}
30
31impl Vertex {
32    /// Create a new vertex
33    pub fn new(x: f64, y: f64, z: f64) -> Self {
34        Self { x, y, z }
35    }
36
37    /// Compute the squared distance to another vertex
38    pub fn distance_sq(&self, other: &Vertex) -> f64 {
39        let dx = self.x - other.x;
40        let dy = self.y - other.y;
41        let dz = self.z - other.z;
42        dx * dx + dy * dy + dz * dz
43    }
44
45    /// Compute the distance to another vertex
46    pub fn distance(&self, other: &Vertex) -> f64 {
47        self.distance_sq(other).sqrt()
48    }
49
50    /// Subtract another vertex, returning a vector
51    pub fn sub(&self, other: &Vertex) -> [f64; 3] {
52        [self.x - other.x, self.y - other.y, self.z - other.z]
53    }
54
55    /// Midpoint with another vertex
56    pub fn midpoint(&self, other: &Vertex) -> Vertex {
57        Vertex {
58            x: (self.x + other.x) * 0.5,
59            y: (self.y + other.y) * 0.5,
60            z: (self.z + other.z) * 0.5,
61        }
62    }
63}
64
65impl fmt::Display for Vertex {
66    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
67        write!(f, "({}, {}, {})", self.x, self.y, self.z)
68    }
69}
70
71/// A triangular face referencing 3 vertex indices
72#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash)]
73pub struct Face {
74    pub v0: usize,
75    pub v1: usize,
76    pub v2: usize,
77}
78
79impl Face {
80    /// Create a new face
81    pub fn new(v0: usize, v1: usize, v2: usize) -> Self {
82        Self { v0, v1, v2 }
83    }
84
85    /// Get the indices as an array
86    pub fn indices(&self) -> [usize; 3] {
87        [self.v0, self.v1, self.v2]
88    }
89
90    /// Check if this face contains a vertex index
91    pub fn contains_vertex(&self, v: usize) -> bool {
92        self.v0 == v || self.v1 == v || self.v2 == v
93    }
94
95    /// Replace a vertex index in this face
96    pub fn replace_vertex(&mut self, old: usize, new: usize) {
97        if self.v0 == old {
98            self.v0 = new;
99        }
100        if self.v1 == old {
101            self.v1 = new;
102        }
103        if self.v2 == old {
104            self.v2 = new;
105        }
106    }
107
108    /// Check if this face is degenerate (has duplicate vertices)
109    pub fn is_degenerate(&self) -> bool {
110        self.v0 == self.v1 || self.v1 == self.v2 || self.v0 == self.v2
111    }
112}
113
114/// An edge represented as a pair of vertex indices (ordered so that a < b)
115#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash)]
116pub struct Edge {
117    pub a: usize,
118    pub b: usize,
119}
120
121impl Edge {
122    /// Create a new edge (automatically ordered)
123    pub fn new(v0: usize, v1: usize) -> Self {
124        if v0 <= v1 {
125            Self { a: v0, b: v1 }
126        } else {
127            Self { a: v1, b: v0 }
128        }
129    }
130}
131
132/// A triangle mesh with indexed face set representation
133///
134/// # Examples
135///
136/// ```
137/// use scirs2_spatial::mesh::{TriangleMesh, Vertex, Face};
138///
139/// let vertices = vec![
140///     Vertex::new(0.0, 0.0, 0.0),
141///     Vertex::new(1.0, 0.0, 0.0),
142///     Vertex::new(0.5, 1.0, 0.0),
143///     Vertex::new(0.5, 0.5, 1.0),
144/// ];
145///
146/// let faces = vec![
147///     Face::new(0, 1, 2),
148///     Face::new(0, 1, 3),
149///     Face::new(1, 2, 3),
150///     Face::new(0, 2, 3),
151/// ];
152///
153/// let mesh = TriangleMesh::new(vertices, faces).expect("Operation failed");
154/// assert_eq!(mesh.num_vertices(), 4);
155/// assert_eq!(mesh.num_faces(), 4);
156/// ```
157#[derive(Debug, Clone)]
158pub struct TriangleMesh {
159    /// Vertex positions
160    pub vertices: Vec<Vertex>,
161    /// Triangle faces (indices into vertices)
162    pub faces: Vec<Face>,
163}
164
165impl TriangleMesh {
166    /// Create a new triangle mesh
167    ///
168    /// # Arguments
169    ///
170    /// * `vertices` - List of vertex positions
171    /// * `faces` - List of triangular faces
172    ///
173    /// # Returns
174    ///
175    /// A validated triangle mesh
176    pub fn new(vertices: Vec<Vertex>, faces: Vec<Face>) -> SpatialResult<Self> {
177        // Validate face indices
178        let n = vertices.len();
179        for (i, face) in faces.iter().enumerate() {
180            for &idx in &[face.v0, face.v1, face.v2] {
181                if idx >= n {
182                    return Err(SpatialError::ValueError(format!(
183                        "Face {} references vertex index {} but only {} vertices exist",
184                        i, idx, n
185                    )));
186                }
187            }
188            if face.is_degenerate() {
189                return Err(SpatialError::ValueError(format!(
190                    "Face {} is degenerate (has duplicate vertices)",
191                    i
192                )));
193            }
194        }
195
196        Ok(Self { vertices, faces })
197    }
198
199    /// Create an empty mesh
200    pub fn empty() -> Self {
201        Self {
202            vertices: Vec::new(),
203            faces: Vec::new(),
204        }
205    }
206
207    /// Number of vertices
208    pub fn num_vertices(&self) -> usize {
209        self.vertices.len()
210    }
211
212    /// Number of faces
213    pub fn num_faces(&self) -> usize {
214        self.faces.len()
215    }
216
217    /// Number of edges (computed from face connectivity)
218    pub fn num_edges(&self) -> usize {
219        self.edges().len()
220    }
221
222    /// Get all unique edges in the mesh
223    pub fn edges(&self) -> HashSet<Edge> {
224        let mut edges = HashSet::new();
225        for face in &self.faces {
226            edges.insert(Edge::new(face.v0, face.v1));
227            edges.insert(Edge::new(face.v1, face.v2));
228            edges.insert(Edge::new(face.v0, face.v2));
229        }
230        edges
231    }
232
233    /// Build vertex-to-face adjacency map
234    pub fn vertex_faces(&self) -> HashMap<usize, Vec<usize>> {
235        let mut vf: HashMap<usize, Vec<usize>> = HashMap::new();
236        for (fi, face) in self.faces.iter().enumerate() {
237            for &vi in &[face.v0, face.v1, face.v2] {
238                vf.entry(vi).or_default().push(fi);
239            }
240        }
241        vf
242    }
243
244    /// Build vertex-to-vertex adjacency (neighbors)
245    pub fn vertex_neighbors(&self) -> HashMap<usize, HashSet<usize>> {
246        let mut vn: HashMap<usize, HashSet<usize>> = HashMap::new();
247        for face in &self.faces {
248            let indices = face.indices();
249            for i in 0..3 {
250                let vi = indices[i];
251                for j in 0..3 {
252                    if i != j {
253                        vn.entry(vi).or_default().insert(indices[j]);
254                    }
255                }
256            }
257        }
258        vn
259    }
260
261    /// Compute face normal for a given face index (unnormalized)
262    pub fn face_normal_raw(&self, face_idx: usize) -> SpatialResult<[f64; 3]> {
263        if face_idx >= self.faces.len() {
264            return Err(SpatialError::ValueError(format!(
265                "Face index {} out of range (num faces = {})",
266                face_idx,
267                self.faces.len()
268            )));
269        }
270
271        let face = &self.faces[face_idx];
272        let v0 = &self.vertices[face.v0];
273        let v1 = &self.vertices[face.v1];
274        let v2 = &self.vertices[face.v2];
275
276        let e1 = v1.sub(v0);
277        let e2 = v2.sub(v0);
278
279        // Cross product e1 x e2
280        let nx = e1[1] * e2[2] - e1[2] * e2[1];
281        let ny = e1[2] * e2[0] - e1[0] * e2[2];
282        let nz = e1[0] * e2[1] - e1[1] * e2[0];
283
284        Ok([nx, ny, nz])
285    }
286
287    /// Compute normalized face normal for a given face index
288    pub fn face_normal(&self, face_idx: usize) -> SpatialResult<[f64; 3]> {
289        let n = self.face_normal_raw(face_idx)?;
290        let len = (n[0] * n[0] + n[1] * n[1] + n[2] * n[2]).sqrt();
291        if len < 1e-15 {
292            return Ok([0.0, 0.0, 0.0]);
293        }
294        Ok([n[0] / len, n[1] / len, n[2] / len])
295    }
296
297    /// Compute all face normals (normalized)
298    pub fn compute_face_normals(&self) -> SpatialResult<Vec<[f64; 3]>> {
299        let mut normals = Vec::with_capacity(self.faces.len());
300        for i in 0..self.faces.len() {
301            normals.push(self.face_normal(i)?);
302        }
303        Ok(normals)
304    }
305
306    /// Compute vertex normals by averaging adjacent face normals
307    /// weighted by face area
308    pub fn compute_vertex_normals(&self) -> SpatialResult<Vec<[f64; 3]>> {
309        let n = self.vertices.len();
310        let mut normals = vec![[0.0_f64; 3]; n];
311
312        for (fi, face) in self.faces.iter().enumerate() {
313            let fn_raw = self.face_normal_raw(fi)?;
314            // The magnitude of the raw normal = 2 * face area, so weighting
315            // by face area is automatic with unnormalized normals
316            for &vi in &[face.v0, face.v1, face.v2] {
317                normals[vi][0] += fn_raw[0];
318                normals[vi][1] += fn_raw[1];
319                normals[vi][2] += fn_raw[2];
320            }
321        }
322
323        // Normalize
324        for normal in &mut normals {
325            let len =
326                (normal[0] * normal[0] + normal[1] * normal[1] + normal[2] * normal[2]).sqrt();
327            if len > 1e-15 {
328                normal[0] /= len;
329                normal[1] /= len;
330                normal[2] /= len;
331            }
332        }
333
334        Ok(normals)
335    }
336
337    /// Compute the area of a face
338    pub fn face_area(&self, face_idx: usize) -> SpatialResult<f64> {
339        let n = self.face_normal_raw(face_idx)?;
340        let len = (n[0] * n[0] + n[1] * n[1] + n[2] * n[2]).sqrt();
341        Ok(len * 0.5)
342    }
343
344    /// Compute the total surface area of the mesh
345    pub fn surface_area(&self) -> SpatialResult<f64> {
346        let mut total = 0.0;
347        for i in 0..self.faces.len() {
348            total += self.face_area(i)?;
349        }
350        Ok(total)
351    }
352
353    /// Compute the bounding box of the mesh
354    pub fn bounding_box(&self) -> Option<(Vertex, Vertex)> {
355        if self.vertices.is_empty() {
356            return None;
357        }
358
359        let mut min_v = self.vertices[0];
360        let mut max_v = self.vertices[0];
361
362        for v in &self.vertices {
363            min_v.x = min_v.x.min(v.x);
364            min_v.y = min_v.y.min(v.y);
365            min_v.z = min_v.z.min(v.z);
366            max_v.x = max_v.x.max(v.x);
367            max_v.y = max_v.y.max(v.y);
368            max_v.z = max_v.z.max(v.z);
369        }
370
371        Some((min_v, max_v))
372    }
373
374    /// Export to ASCII STL format string
375    ///
376    /// # Arguments
377    ///
378    /// * `name` - Name for the solid in the STL file
379    pub fn to_stl_ascii(&self, name: &str) -> SpatialResult<String> {
380        let mut out = format!("solid {}\n", name);
381
382        for (fi, face) in self.faces.iter().enumerate() {
383            let n = self.face_normal(fi)?;
384            out.push_str(&format!("  facet normal {} {} {}\n", n[0], n[1], n[2]));
385            out.push_str("    outer loop\n");
386
387            let v0 = &self.vertices[face.v0];
388            let v1 = &self.vertices[face.v1];
389            let v2 = &self.vertices[face.v2];
390
391            out.push_str(&format!("      vertex {} {} {}\n", v0.x, v0.y, v0.z));
392            out.push_str(&format!("      vertex {} {} {}\n", v1.x, v1.y, v1.z));
393            out.push_str(&format!("      vertex {} {} {}\n", v2.x, v2.y, v2.z));
394
395            out.push_str("    endloop\n");
396            out.push_str("  endfacet\n");
397        }
398
399        out.push_str(&format!("endsolid {}\n", name));
400        Ok(out)
401    }
402
403    /// Import from ASCII STL format string
404    ///
405    /// # Arguments
406    ///
407    /// * `stl_data` - The ASCII STL content as a string
408    pub fn from_stl_ascii(stl_data: &str) -> SpatialResult<Self> {
409        let mut vertices: Vec<Vertex> = Vec::new();
410        let mut faces: Vec<Face> = Vec::new();
411        let mut vertex_map: HashMap<[u64; 3], usize> = HashMap::new();
412        let mut current_face_verts: Vec<usize> = Vec::new();
413
414        for line in stl_data.lines() {
415            let trimmed = line.trim();
416
417            if let Some(rest) = trimmed.strip_prefix("vertex") {
418                let parts: Vec<&str> = rest.split_whitespace().collect();
419                if parts.len() < 3 {
420                    return Err(SpatialError::ValueError(
421                        "Invalid vertex line in STL".to_string(),
422                    ));
423                }
424
425                let x: f64 = parts[0].parse().map_err(|e| {
426                    SpatialError::ValueError(format!("Failed to parse vertex x: {}", e))
427                })?;
428                let y: f64 = parts[1].parse().map_err(|e| {
429                    SpatialError::ValueError(format!("Failed to parse vertex y: {}", e))
430                })?;
431                let z: f64 = parts[2].parse().map_err(|e| {
432                    SpatialError::ValueError(format!("Failed to parse vertex z: {}", e))
433                })?;
434
435                // Use bit representation for hashing (exact vertex matching)
436                let key = [x.to_bits(), y.to_bits(), z.to_bits()];
437
438                let idx = if let Some(&existing) = vertex_map.get(&key) {
439                    existing
440                } else {
441                    let idx = vertices.len();
442                    vertices.push(Vertex::new(x, y, z));
443                    vertex_map.insert(key, idx);
444                    idx
445                };
446
447                current_face_verts.push(idx);
448            } else if trimmed == "endloop" {
449                if current_face_verts.len() == 3 {
450                    let v0 = current_face_verts[0];
451                    let v1 = current_face_verts[1];
452                    let v2 = current_face_verts[2];
453
454                    // Skip degenerate faces
455                    if v0 != v1 && v1 != v2 && v0 != v2 {
456                        faces.push(Face::new(v0, v1, v2));
457                    }
458                }
459                current_face_verts.clear();
460            }
461        }
462
463        Ok(Self { vertices, faces })
464    }
465
466    /// Check if the mesh is manifold (every edge has exactly 1 or 2 adjacent faces)
467    pub fn is_manifold(&self) -> bool {
468        let mut edge_count: HashMap<Edge, usize> = HashMap::new();
469        for face in &self.faces {
470            *edge_count.entry(Edge::new(face.v0, face.v1)).or_insert(0) += 1;
471            *edge_count.entry(Edge::new(face.v1, face.v2)).or_insert(0) += 1;
472            *edge_count.entry(Edge::new(face.v0, face.v2)).or_insert(0) += 1;
473        }
474        edge_count.values().all(|&c| c == 1 || c == 2)
475    }
476
477    /// Check if the mesh is closed (every edge has exactly 2 adjacent faces)
478    pub fn is_closed(&self) -> bool {
479        let mut edge_count: HashMap<Edge, usize> = HashMap::new();
480        for face in &self.faces {
481            *edge_count.entry(Edge::new(face.v0, face.v1)).or_insert(0) += 1;
482            *edge_count.entry(Edge::new(face.v1, face.v2)).or_insert(0) += 1;
483            *edge_count.entry(Edge::new(face.v0, face.v2)).or_insert(0) += 1;
484        }
485        edge_count.values().all(|&c| c == 2)
486    }
487
488    /// Get boundary edges (edges with only 1 adjacent face)
489    pub fn boundary_edges(&self) -> Vec<Edge> {
490        let mut edge_count: HashMap<Edge, usize> = HashMap::new();
491        for face in &self.faces {
492            *edge_count.entry(Edge::new(face.v0, face.v1)).or_insert(0) += 1;
493            *edge_count.entry(Edge::new(face.v1, face.v2)).or_insert(0) += 1;
494            *edge_count.entry(Edge::new(face.v0, face.v2)).or_insert(0) += 1;
495        }
496        edge_count
497            .into_iter()
498            .filter(|&(_, c)| c == 1)
499            .map(|(e, _)| e)
500            .collect()
501    }
502
503    /// Compute the Euler characteristic: V - E + F
504    pub fn euler_characteristic(&self) -> i64 {
505        let v = self.num_vertices() as i64;
506        let e = self.num_edges() as i64;
507        let f = self.num_faces() as i64;
508        v - e + f
509    }
510}
511
512// Re-exports from submodules
513pub use quality::{face_aspect_ratio, face_min_angle, mesh_quality_stats, QualityStats};
514pub use simplification::simplify_mesh;
515pub use smoothing::{laplacian_smooth, taubin_smooth};
516
517#[cfg(test)]
518mod tests {
519    use super::*;
520
521    fn tetrahedron() -> TriangleMesh {
522        let vertices = vec![
523            Vertex::new(0.0, 0.0, 0.0),
524            Vertex::new(1.0, 0.0, 0.0),
525            Vertex::new(0.5, 1.0, 0.0),
526            Vertex::new(0.5, 0.5, 1.0),
527        ];
528        let faces = vec![
529            Face::new(0, 1, 2),
530            Face::new(0, 1, 3),
531            Face::new(1, 2, 3),
532            Face::new(0, 2, 3),
533        ];
534        TriangleMesh::new(vertices, faces).expect("tetrahedron should be valid")
535    }
536
537    #[test]
538    fn test_mesh_creation() {
539        let mesh = tetrahedron();
540        assert_eq!(mesh.num_vertices(), 4);
541        assert_eq!(mesh.num_faces(), 4);
542        assert_eq!(mesh.num_edges(), 6);
543    }
544
545    #[test]
546    fn test_mesh_invalid_face() {
547        let vertices = vec![Vertex::new(0.0, 0.0, 0.0), Vertex::new(1.0, 0.0, 0.0)];
548        let faces = vec![Face::new(0, 1, 5)]; // index 5 out of range
549        let result = TriangleMesh::new(vertices, faces);
550        assert!(result.is_err());
551    }
552
553    #[test]
554    fn test_degenerate_face() {
555        let vertices = vec![Vertex::new(0.0, 0.0, 0.0), Vertex::new(1.0, 0.0, 0.0)];
556        let faces = vec![Face::new(0, 0, 1)]; // degenerate
557        let result = TriangleMesh::new(vertices, faces);
558        assert!(result.is_err());
559    }
560
561    #[test]
562    fn test_face_normals() {
563        let mesh = tetrahedron();
564        let normals = mesh.compute_face_normals().expect("compute normals failed");
565        assert_eq!(normals.len(), 4);
566
567        // Face 0 (0,1,2) lies in XY plane, normal should be (0,0,+/-1)
568        let n0 = normals[0];
569        assert!((n0[0].abs()) < 1e-10);
570        assert!((n0[1].abs()) < 1e-10);
571        assert!((n0[2].abs() - 1.0) < 1e-10);
572    }
573
574    #[test]
575    fn test_vertex_normals() {
576        let mesh = tetrahedron();
577        let normals = mesh
578            .compute_vertex_normals()
579            .expect("compute normals failed");
580        assert_eq!(normals.len(), 4);
581
582        // Each vertex normal should be unit length
583        for n in &normals {
584            let len = (n[0] * n[0] + n[1] * n[1] + n[2] * n[2]).sqrt();
585            assert!((len - 1.0).abs() < 1e-10, "Normal length: {}", len);
586        }
587    }
588
589    #[test]
590    fn test_surface_area() {
591        let mesh = tetrahedron();
592        let area = mesh.surface_area().expect("area failed");
593        assert!(area > 0.0);
594    }
595
596    #[test]
597    fn test_bounding_box() {
598        let mesh = tetrahedron();
599        let (min_v, max_v) = mesh.bounding_box().expect("non-empty mesh");
600        assert!((min_v.x - 0.0).abs() < 1e-10);
601        assert!((min_v.y - 0.0).abs() < 1e-10);
602        assert!((min_v.z - 0.0).abs() < 1e-10);
603        assert!((max_v.x - 1.0).abs() < 1e-10);
604        assert!((max_v.y - 1.0).abs() < 1e-10);
605        assert!((max_v.z - 1.0).abs() < 1e-10);
606    }
607
608    #[test]
609    fn test_manifold_closed() {
610        let mesh = tetrahedron();
611        assert!(mesh.is_manifold());
612        assert!(mesh.is_closed());
613        assert!(mesh.boundary_edges().is_empty());
614    }
615
616    #[test]
617    fn test_euler_characteristic() {
618        let mesh = tetrahedron();
619        // V - E + F = 4 - 6 + 4 = 2 (sphere topology)
620        assert_eq!(mesh.euler_characteristic(), 2);
621    }
622
623    #[test]
624    fn test_stl_roundtrip() {
625        let mesh = tetrahedron();
626        let stl_str = mesh.to_stl_ascii("test").expect("STL export failed");
627        assert!(stl_str.starts_with("solid test"));
628        assert!(stl_str.contains("facet normal"));
629        assert!(stl_str.contains("vertex"));
630
631        let mesh2 = TriangleMesh::from_stl_ascii(&stl_str).expect("STL import failed");
632        // Vertex count may differ due to deduplication, but face count should match
633        assert_eq!(mesh2.num_faces(), mesh.num_faces());
634    }
635
636    #[test]
637    fn test_vertex_neighbors() {
638        let mesh = tetrahedron();
639        let neighbors = mesh.vertex_neighbors();
640
641        // In a tetrahedron, every vertex connects to every other vertex
642        for i in 0..4 {
643            let nbs = neighbors.get(&i).expect("vertex should have neighbors");
644            assert_eq!(nbs.len(), 3);
645        }
646    }
647
648    #[test]
649    fn test_open_mesh() {
650        // Single triangle - open mesh
651        let vertices = vec![
652            Vertex::new(0.0, 0.0, 0.0),
653            Vertex::new(1.0, 0.0, 0.0),
654            Vertex::new(0.5, 1.0, 0.0),
655        ];
656        let faces = vec![Face::new(0, 1, 2)];
657        let mesh = TriangleMesh::new(vertices, faces).expect("valid");
658
659        assert!(mesh.is_manifold());
660        assert!(!mesh.is_closed());
661        assert_eq!(mesh.boundary_edges().len(), 3);
662    }
663}
scirs2_spatial/mesh/mod.rs

scirs2_spatial/mesh/
mod.rs
