threecrate_algorithms/
mesh_boolean.rs

1//! Mesh boolean operations: union, intersection, and difference
2//!
3//! Implements CSG (Constructive Solid Geometry) using a BSP-tree approach.
4//! Each closed `TriangleMesh` is converted into a set of `Polygon`s, which are
5//! inserted into a BSP tree.  The classic three-step clip/invert/build sequences
6//! then produce the desired solid.
7//!
8//! **Limitations**
9//! - Meshes must be closed (watertight) and have outward-facing normals.
10//! - Very thin or near-degenerate triangles may be dropped (area < `EPSILON²`).
11//! - Shared planar faces are handled via coplanar classification but exact
12//!   floating-point coincidence can produce small artefacts at seam edges.
13
14use threecrate_core::{TriangleMesh, Result, Point3f, Vector3f};
15
16/// Numeric tolerance used for plane-side classification and degenerate checks.
17const EPSILON: f32 = 1e-5;
18
19// ---------------------------------------------------------------------------
20// Internal geometry types
21// ---------------------------------------------------------------------------
22
23/// A plane in the form `normal · x = w`.
24#[derive(Debug, Clone)]
25struct Plane {
26    normal: Vector3f,
27    w: f32,
28}
29
30impl Plane {
31    fn from_points(a: &Point3f, b: &Point3f, c: &Point3f) -> Option<Self> {
32        let n = (b - a).cross(&(c - a));
33        let len = n.magnitude();
34        if len < EPSILON * EPSILON {
35            return None; // degenerate triangle
36        }
37        let normal = n / len;
38        let w = normal.dot(&a.coords);
39        Some(Plane { normal, w })
40    }
41
42    fn flip(&mut self) {
43        self.normal = -self.normal;
44        self.w = -self.w;
45    }
46
47    #[inline]
48    fn signed_distance(&self, p: &Point3f) -> f32 {
49        self.normal.dot(&p.coords) - self.w
50    }
51
52    /// Split `poly` into up to four output lists: coplanar-front, coplanar-back,
53    /// front, and back.  Spanning polygons are clipped and triangulated.
54    fn split_polygon(
55        &self,
56        poly: &Polygon,
57        coplanar_front: &mut Vec<Polygon>,
58        coplanar_back: &mut Vec<Polygon>,
59        front: &mut Vec<Polygon>,
60        back: &mut Vec<Polygon>,
61    ) {
62        let verts = &poly.vertices;
63        let n = verts.len();
64
65        // Classify every vertex
66        let dists: Vec<f32> = verts.iter().map(|v| self.signed_distance(v)).collect();
67        let mut poly_mask: u8 = 0;
68        for &d in &dists {
69            if d > EPSILON {
70                poly_mask |= 1; // front
71            } else if d < -EPSILON {
72                poly_mask |= 2; // back
73            }
74            // else on-plane → 0
75        }
76
77        match poly_mask {
78            0 => {
79                // Coplanar – classify by normal alignment
80                if self.normal.dot(&poly.plane.normal) > 0.0 {
81                    coplanar_front.push(poly.clone());
82                } else {
83                    coplanar_back.push(poly.clone());
84                }
85            }
86            1 => front.push(poly.clone()),
87            2 => back.push(poly.clone()),
88            _ => {
89                // Spanning – clip and re-triangulate
90                let mut f_verts: Vec<Point3f> = Vec::new();
91                let mut b_verts: Vec<Point3f> = Vec::new();
92
93                for i in 0..n {
94                    let j = (i + 1) % n;
95                    let di = dists[i];
96                    let dj = dists[j];
97                    let vi = verts[i];
98                    let vj = verts[j];
99
100                    let i_front = di > EPSILON;
101                    let i_back = di < -EPSILON;
102                    let j_front = dj > EPSILON;
103                    let j_back = dj < -EPSILON;
104
105                    if !i_back {
106                        f_verts.push(vi);
107                    }
108                    if !i_front {
109                        b_verts.push(vi);
110                    }
111
112                    // Edge crosses the plane?
113                    if (i_front && j_back) || (i_back && j_front) {
114                        let t = di / (di - dj);
115                        let mid = vi + (vj - vi) * t;
116                        f_verts.push(mid);
117                        b_verts.push(mid);
118                    }
119                }
120
121                front.extend(Polygon::fan_triangulate(f_verts));
122                back.extend(Polygon::fan_triangulate(b_verts));
123            }
124        }
125    }
126}
127
128/// A convex polygon stored as an ordered vertex list plus its precomputed plane.
129#[derive(Debug, Clone)]
130struct Polygon {
131    vertices: Vec<Point3f>,
132    plane: Plane,
133}
134
135impl Polygon {
136    fn new(vertices: Vec<Point3f>) -> Option<Self> {
137        if vertices.len() < 3 {
138            return None;
139        }
140        let plane = Plane::from_points(&vertices[0], &vertices[1], &vertices[2])?;
141        Some(Polygon { vertices, plane })
142    }
143
144    fn flip(&mut self) {
145        self.vertices.reverse();
146        self.plane.flip();
147    }
148
149    /// Fan-triangulate a convex vertex list into `Polygon` triangles.
150    fn fan_triangulate(verts: Vec<Point3f>) -> Vec<Polygon> {
151        if verts.len() < 3 {
152            return Vec::new();
153        }
154        let mut out = Vec::with_capacity(verts.len() - 2);
155        for i in 1..verts.len() - 1 {
156            if let Some(p) = Polygon::new(vec![verts[0], verts[i], verts[i + 1]]) {
157                out.push(p);
158            }
159        }
160        out
161    }
162}
163
164// ---------------------------------------------------------------------------
165// BSP tree
166// ---------------------------------------------------------------------------
167
168struct BspNode {
169    plane: Option<Plane>,
170    front: Option<Box<BspNode>>,
171    back: Option<Box<BspNode>>,
172    polygons: Vec<Polygon>,
173}
174
175impl BspNode {
176    fn new() -> Self {
177        BspNode { plane: None, front: None, back: None, polygons: Vec::new() }
178    }
179
180    fn invert(&mut self) {
181        for p in &mut self.polygons {
182            p.flip();
183        }
184        if let Some(ref mut plane) = self.plane {
185            plane.flip();
186        }
187        if let Some(ref mut f) = self.front {
188            f.invert();
189        }
190        if let Some(ref mut b) = self.back {
191            b.invert();
192        }
193        std::mem::swap(&mut self.front, &mut self.back);
194    }
195
196    /// Recursively discard polygons inside this BSP tree, keeping those outside.
197    fn clip_polygons(&self, polygons: Vec<Polygon>) -> Vec<Polygon> {
198        let plane = match &self.plane {
199            Some(p) => p,
200            None => return polygons,
201        };
202
203        let mut f: Vec<Polygon> = Vec::new();
204        let mut b: Vec<Polygon> = Vec::new();
205
206        for poly in polygons {
207            let mut cf = Vec::new();
208            let mut cb = Vec::new();
209            let mut pf = Vec::new();
210            let mut pb = Vec::new();
211            plane.split_polygon(&poly, &mut cf, &mut cb, &mut pf, &mut pb);
212            f.extend(cf);
213            f.extend(pf);
214            b.extend(cb);
215            b.extend(pb);
216        }
217
218        let mut result = match &self.front {
219            Some(node) => node.clip_polygons(f),
220            None => f,
221        };
222        let back_result = match &self.back {
223            Some(node) => node.clip_polygons(b),
224            None => Vec::new(), // discard: these are inside the solid
225        };
226        result.extend(back_result);
227        result
228    }
229
230    /// Remove all polygons in this tree that are inside `other`.
231    fn clip_to(&mut self, other: &BspNode) {
232        self.polygons = other.clip_polygons(std::mem::take(&mut self.polygons));
233        if let Some(ref mut f) = self.front {
234            f.clip_to(other);
235        }
236        if let Some(ref mut b) = self.back {
237            b.clip_to(other);
238        }
239    }
240
241    /// Collect all polygons stored in this tree.
242    fn all_polygons(&self) -> Vec<Polygon> {
243        let mut out = self.polygons.clone();
244        if let Some(ref f) = self.front {
245            out.extend(f.all_polygons());
246        }
247        if let Some(ref b) = self.back {
248            out.extend(b.all_polygons());
249        }
250        out
251    }
252
253    /// Insert `polygons` into the BSP tree (building the tree as needed).
254    fn build(&mut self, polygons: Vec<Polygon>) {
255        if polygons.is_empty() {
256            return;
257        }
258
259        if self.plane.is_none() {
260            self.plane = Some(polygons[0].plane.clone());
261        }
262        let plane = self.plane.clone().unwrap();
263
264        let mut front = Vec::new();
265        let mut back = Vec::new();
266
267        for poly in polygons {
268            let mut cf = Vec::new();
269            let mut cb = Vec::new();
270            let mut pf = Vec::new();
271            let mut pb = Vec::new();
272            plane.split_polygon(&poly, &mut cf, &mut cb, &mut pf, &mut pb);
273            self.polygons.extend(cf);
274            self.polygons.extend(cb);
275            front.extend(pf);
276            back.extend(pb);
277        }
278
279        if !front.is_empty() {
280            if self.front.is_none() {
281                self.front = Some(Box::new(BspNode::new()));
282            }
283            self.front.as_mut().unwrap().build(front);
284        }
285        if !back.is_empty() {
286            if self.back.is_none() {
287                self.back = Some(Box::new(BspNode::new()));
288            }
289            self.back.as_mut().unwrap().build(back);
290        }
291    }
292}
293
294// ---------------------------------------------------------------------------
295// Mesh ↔ polygon conversion helpers
296// ---------------------------------------------------------------------------
297
298fn mesh_to_polygons(mesh: &TriangleMesh) -> Vec<Polygon> {
299    mesh.faces
300        .iter()
301        .filter_map(|face| {
302            Polygon::new(vec![
303                mesh.vertices[face[0]],
304                mesh.vertices[face[1]],
305                mesh.vertices[face[2]],
306            ])
307        })
308        .collect()
309}
310
311fn polygons_to_mesh(polygons: Vec<Polygon>) -> TriangleMesh {
312    // Each Polygon is already a triangle (fan-triangulation keeps them as triangles),
313    // so we can flatten directly.
314    let capacity = polygons.iter().map(|p| p.vertices.len().saturating_sub(2)).sum::<usize>();
315    let mut vertices = Vec::with_capacity(capacity * 3);
316    let mut faces = Vec::with_capacity(capacity);
317
318    for poly in &polygons {
319        let n = poly.vertices.len();
320        if n < 3 {
321            continue;
322        }
323        // Fan triangulation
324        for i in 1..n - 1 {
325            let base = vertices.len();
326            vertices.push(poly.vertices[0]);
327            vertices.push(poly.vertices[i]);
328            vertices.push(poly.vertices[i + 1]);
329            faces.push([base, base + 1, base + 2]);
330        }
331    }
332
333    TriangleMesh::from_vertices_and_faces(vertices, faces)
334}
335
336fn build_bsp(mesh: &TriangleMesh) -> BspNode {
337    let mut node = BspNode::new();
338    node.build(mesh_to_polygons(mesh));
339    node
340}
341
342// ---------------------------------------------------------------------------
343// Public API
344// ---------------------------------------------------------------------------
345
346/// Which boolean operation to perform.
347#[derive(Debug, Clone, Copy, PartialEq, Eq)]
348pub enum BooleanOp {
349    /// All regions inside A **or** B.
350    Union,
351    /// All regions inside both A **and** B.
352    Intersection,
353    /// All regions inside A but **not** B.
354    Difference,
355}
356
357/// Perform a mesh boolean operation on two closed triangle meshes.
358///
359/// # Arguments
360/// * `a` - First operand mesh (must be closed / watertight)
361/// * `b` - Second operand mesh (must be closed / watertight)
362/// * `op` - The boolean operation to apply
363///
364/// # Returns
365/// A new `TriangleMesh` representing the result of the operation.
366pub fn mesh_boolean(
367    a: &TriangleMesh,
368    b: &TriangleMesh,
369    op: BooleanOp,
370) -> Result<TriangleMesh> {
371    match op {
372        BooleanOp::Union => mesh_union(a, b),
373        BooleanOp::Intersection => mesh_intersection(a, b),
374        BooleanOp::Difference => mesh_difference(a, b),
375    }
376}
377
378/// Compute the boolean **union** of two closed triangle meshes.
379///
380/// Returns a mesh containing all regions inside either `a` or `b`.
381///
382/// ```text
383///     A.union(B)
384///
385///     +-------+            +-------+
386///     |       |            |       |
387///     |   A   |            |       |
388///     |    +--+----+   =   |       +----+
389///     +----+--+    |       +----+       |
390///          |   B   |            |       |
391///          |       |            |       |
392///          +-------+            +-------+
393/// ```
394pub fn mesh_union(a: &TriangleMesh, b: &TriangleMesh) -> Result<TriangleMesh> {
395    if a.is_empty() {
396        return Ok(b.clone());
397    }
398    if b.is_empty() {
399        return Ok(a.clone());
400    }
401
402    let mut a_bsp = build_bsp(a);
403    let mut b_bsp = build_bsp(b);
404
405    a_bsp.clip_to(&b_bsp);
406    b_bsp.clip_to(&a_bsp);
407    b_bsp.invert();
408    b_bsp.clip_to(&a_bsp);
409    b_bsp.invert();
410    a_bsp.build(b_bsp.all_polygons());
411
412    Ok(polygons_to_mesh(a_bsp.all_polygons()))
413}
414
415/// Compute the boolean **intersection** of two closed triangle meshes.
416///
417/// Returns a mesh containing only the regions inside both `a` and `b`.
418///
419/// ```text
420///     A.intersection(B)
421///
422///     +-------+
423///     |       |
424///     |   A   |
425///     |    +--+----+   =   +--+
426///     +----+--+    |       +--+
427///          |   B   |
428///          |       |
429///          +-------+
430/// ```
431pub fn mesh_intersection(a: &TriangleMesh, b: &TriangleMesh) -> Result<TriangleMesh> {
432    if a.is_empty() || b.is_empty() {
433        return Ok(TriangleMesh::new());
434    }
435
436    let mut a_bsp = build_bsp(a);
437    let mut b_bsp = build_bsp(b);
438
439    a_bsp.invert();
440    b_bsp.clip_to(&a_bsp);
441    b_bsp.invert();
442    a_bsp.clip_to(&b_bsp);
443    b_bsp.clip_to(&a_bsp);
444    a_bsp.build(b_bsp.all_polygons());
445    a_bsp.invert();
446
447    Ok(polygons_to_mesh(a_bsp.all_polygons()))
448}
449
450/// Compute the boolean **difference** of two closed triangle meshes (A minus B).
451///
452/// Returns a mesh containing only the regions inside `a` but not `b`.
453///
454/// ```text
455///     A.difference(B)
456///
457///     +-------+            +-------+
458///     |       |            |       |
459///     |   A   |            |   A   |
460///     |    +--+----+   =   |    +--+
461///     +----+--+    |       +----+
462///          |   B   |
463///          |       |
464///          +-------+
465/// ```
466pub fn mesh_difference(a: &TriangleMesh, b: &TriangleMesh) -> Result<TriangleMesh> {
467    if a.is_empty() {
468        return Ok(TriangleMesh::new());
469    }
470    if b.is_empty() {
471        return Ok(a.clone());
472    }
473
474    let mut a_bsp = build_bsp(a);
475    let mut b_bsp = build_bsp(b);
476
477    a_bsp.invert();
478    a_bsp.clip_to(&b_bsp);
479    b_bsp.clip_to(&a_bsp);
480    b_bsp.invert();
481    b_bsp.clip_to(&a_bsp);
482    b_bsp.invert();
483    a_bsp.build(b_bsp.all_polygons());
484    a_bsp.invert();
485
486    Ok(polygons_to_mesh(a_bsp.all_polygons()))
487}
488
489// ---------------------------------------------------------------------------
490// Tests
491// ---------------------------------------------------------------------------
492
493#[cfg(test)]
494mod tests {
495    use super::*;
496
497    /// Build a closed box mesh with outward-facing normals.
498    /// Vertices are the 8 corners of [min, max]^3.
499    fn make_box(min: Point3f, max: Point3f) -> TriangleMesh {
500        let (x0, y0, z0) = (min.x, min.y, min.z);
501        let (x1, y1, z1) = (max.x, max.y, max.z);
502        let verts = vec![
503            Point3f::new(x0, y0, z0), // 0
504            Point3f::new(x1, y0, z0), // 1
505            Point3f::new(x1, y1, z0), // 2
506            Point3f::new(x0, y1, z0), // 3
507            Point3f::new(x0, y0, z1), // 4
508            Point3f::new(x1, y0, z1), // 5
509            Point3f::new(x1, y1, z1), // 6
510            Point3f::new(x0, y1, z1), // 7
511        ];
512        // Each face is CCW when viewed from outside → outward normal
513        let faces = vec![
514            // Bottom (z-): normal (0,0,-1)
515            [0, 2, 1], [0, 3, 2],
516            // Top (z+): normal (0,0,+1)
517            [4, 5, 6], [4, 6, 7],
518            // Front (y-): normal (0,-1,0)
519            [0, 1, 5], [0, 5, 4],
520            // Back (y+): normal (0,+1,0)
521            [3, 7, 6], [3, 6, 2],
522            // Left (x-): normal (-1,0,0)
523            [0, 4, 7], [0, 7, 3],
524            // Right (x+): normal (+1,0,0)
525            [1, 2, 6], [1, 6, 5],
526        ];
527        TriangleMesh::from_vertices_and_faces(verts, faces)
528    }
529
530    // ---- empty-mesh edge cases ----
531
532    #[test]
533    fn test_union_empty_a() {
534        let empty = TriangleMesh::new();
535        let b = make_box(Point3f::new(0.0, 0.0, 0.0), Point3f::new(1.0, 1.0, 1.0));
536        let result = mesh_union(&empty, &b).unwrap();
537        assert!(!result.is_empty());
538    }
539
540    #[test]
541    fn test_union_empty_b() {
542        let a = make_box(Point3f::new(0.0, 0.0, 0.0), Point3f::new(1.0, 1.0, 1.0));
543        let empty = TriangleMesh::new();
544        let result = mesh_union(&a, &empty).unwrap();
545        assert!(!result.is_empty());
546    }
547
548    #[test]
549    fn test_intersection_empty() {
550        let empty = TriangleMesh::new();
551        let b = make_box(Point3f::new(0.0, 0.0, 0.0), Point3f::new(1.0, 1.0, 1.0));
552        let result = mesh_intersection(&empty, &b).unwrap();
553        assert!(result.is_empty());
554    }
555
556    #[test]
557    fn test_difference_empty_a() {
558        let empty = TriangleMesh::new();
559        let b = make_box(Point3f::new(0.0, 0.0, 0.0), Point3f::new(1.0, 1.0, 1.0));
560        let result = mesh_difference(&empty, &b).unwrap();
561        assert!(result.is_empty());
562    }
563
564    #[test]
565    fn test_difference_empty_b() {
566        let a = make_box(Point3f::new(0.0, 0.0, 0.0), Point3f::new(1.0, 1.0, 1.0));
567        let empty = TriangleMesh::new();
568        let result = mesh_difference(&a, &empty).unwrap();
569        assert!(!result.is_empty());
570    }
571
572    // ---- non-overlapping meshes ----
573
574    #[test]
575    fn test_union_non_overlapping() {
576        // Two unit cubes separated along X
577        let a = make_box(Point3f::new(0.0, 0.0, 0.0), Point3f::new(1.0, 1.0, 1.0));
578        let b = make_box(Point3f::new(5.0, 0.0, 0.0), Point3f::new(6.0, 1.0, 1.0));
579        let result = mesh_union(&a, &b).unwrap();
580        // Both boxes should be fully preserved
581        assert!(!result.is_empty());
582        assert!(result.face_count() >= 24, "expected ≥24 faces, got {}", result.face_count());
583    }
584
585    #[test]
586    fn test_intersection_non_overlapping() {
587        let a = make_box(Point3f::new(0.0, 0.0, 0.0), Point3f::new(1.0, 1.0, 1.0));
588        let b = make_box(Point3f::new(5.0, 0.0, 0.0), Point3f::new(6.0, 1.0, 1.0));
589        let result = mesh_intersection(&a, &b).unwrap();
590        // No overlap → empty result
591        assert!(result.is_empty(), "non-overlapping intersection should be empty, got {} faces", result.face_count());
592    }
593
594    #[test]
595    fn test_difference_non_overlapping() {
596        let a = make_box(Point3f::new(0.0, 0.0, 0.0), Point3f::new(1.0, 1.0, 1.0));
597        let b = make_box(Point3f::new(5.0, 0.0, 0.0), Point3f::new(6.0, 1.0, 1.0));
598        let result = mesh_difference(&a, &b).unwrap();
599        // B doesn't touch A → result should be A unchanged
600        assert!(!result.is_empty());
601        assert!(result.face_count() >= 12, "expected ≥12 faces, got {}", result.face_count());
602    }
603
604    // ---- overlapping meshes ----
605
606    #[test]
607    fn test_union_overlapping() {
608        // Two unit cubes that partially overlap
609        let a = make_box(Point3f::new(0.0, 0.0, 0.0), Point3f::new(2.0, 2.0, 2.0));
610        let b = make_box(Point3f::new(1.0, 1.0, 1.0), Point3f::new(3.0, 3.0, 3.0));
611        let result = mesh_union(&a, &b).unwrap();
612        assert!(!result.is_empty(), "union of overlapping boxes must not be empty");
613    }
614
615    #[test]
616    fn test_intersection_overlapping() {
617        let a = make_box(Point3f::new(0.0, 0.0, 0.0), Point3f::new(2.0, 2.0, 2.0));
618        let b = make_box(Point3f::new(1.0, 1.0, 1.0), Point3f::new(3.0, 3.0, 3.0));
619        let result = mesh_intersection(&a, &b).unwrap();
620        assert!(!result.is_empty(), "intersection of overlapping boxes must not be empty");
621        // The intersection is the unit cube [1,2]^3 → should have at most as many faces as either input
622        assert!(result.face_count() <= a.face_count() + b.face_count() + 50);
623    }
624
625    #[test]
626    fn test_difference_overlapping() {
627        let a = make_box(Point3f::new(0.0, 0.0, 0.0), Point3f::new(2.0, 2.0, 2.0));
628        let b = make_box(Point3f::new(1.0, 1.0, 1.0), Point3f::new(3.0, 3.0, 3.0));
629        let result = mesh_difference(&a, &b).unwrap();
630        assert!(!result.is_empty(), "A minus partially-overlapping B must not be empty");
631    }
632
633    #[test]
634    fn test_difference_a_fully_contains_b() {
635        // B is entirely inside A → difference cuts a hole
636        let a = make_box(Point3f::new(0.0, 0.0, 0.0), Point3f::new(4.0, 4.0, 4.0));
637        let b = make_box(Point3f::new(1.0, 1.0, 1.0), Point3f::new(2.0, 2.0, 2.0));
638        let result = mesh_difference(&a, &b).unwrap();
639        assert!(!result.is_empty(), "difference with interior hole must not be empty");
640    }
641
642    #[test]
643    fn test_intersection_a_fully_contains_b() {
644        // B fully inside A → intersection == B
645        let a = make_box(Point3f::new(0.0, 0.0, 0.0), Point3f::new(4.0, 4.0, 4.0));
646        let b = make_box(Point3f::new(1.0, 1.0, 1.0), Point3f::new(2.0, 2.0, 2.0));
647        let result = mesh_intersection(&a, &b).unwrap();
648        assert!(!result.is_empty(), "intersection when B inside A must not be empty");
649    }
650
651    // ---- BooleanOp dispatch ----
652
653    #[test]
654    fn test_mesh_boolean_dispatch() {
655        let a = make_box(Point3f::new(0.0, 0.0, 0.0), Point3f::new(2.0, 2.0, 2.0));
656        let b = make_box(Point3f::new(1.0, 1.0, 1.0), Point3f::new(3.0, 3.0, 3.0));
657
658        let u = mesh_boolean(&a, &b, BooleanOp::Union).unwrap();
659        let i = mesh_boolean(&a, &b, BooleanOp::Intersection).unwrap();
660        let d = mesh_boolean(&a, &b, BooleanOp::Difference).unwrap();
661
662        assert!(!u.is_empty());
663        assert!(!i.is_empty());
664        assert!(!d.is_empty());
665    }
666}
threecrate_algorithms/mesh_boolean.rs

threecrate_algorithms/
mesh_boolean.rs
