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::{Point3f, Result, TriangleMesh, 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 {
178            plane: None,
179            front: None,
180            back: None,
181            polygons: Vec::new(),
182        }
183    }
184
185    fn invert(&mut self) {
186        for p in &mut self.polygons {
187            p.flip();
188        }
189        if let Some(ref mut plane) = self.plane {
190            plane.flip();
191        }
192        if let Some(ref mut f) = self.front {
193            f.invert();
194        }
195        if let Some(ref mut b) = self.back {
196            b.invert();
197        }
198        std::mem::swap(&mut self.front, &mut self.back);
199    }
200
201    /// Recursively discard polygons inside this BSP tree, keeping those outside.
202    fn clip_polygons(&self, polygons: Vec<Polygon>) -> Vec<Polygon> {
203        let plane = match &self.plane {
204            Some(p) => p,
205            None => return polygons,
206        };
207
208        let mut f: Vec<Polygon> = Vec::new();
209        let mut b: Vec<Polygon> = Vec::new();
210
211        for poly in polygons {
212            let mut cf = Vec::new();
213            let mut cb = Vec::new();
214            let mut pf = Vec::new();
215            let mut pb = Vec::new();
216            plane.split_polygon(&poly, &mut cf, &mut cb, &mut pf, &mut pb);
217            f.extend(cf);
218            f.extend(pf);
219            b.extend(cb);
220            b.extend(pb);
221        }
222
223        let mut result = match &self.front {
224            Some(node) => node.clip_polygons(f),
225            None => f,
226        };
227        let back_result = match &self.back {
228            Some(node) => node.clip_polygons(b),
229            None => Vec::new(), // discard: these are inside the solid
230        };
231        result.extend(back_result);
232        result
233    }
234
235    /// Remove all polygons in this tree that are inside `other`.
236    fn clip_to(&mut self, other: &BspNode) {
237        self.polygons = other.clip_polygons(std::mem::take(&mut self.polygons));
238        if let Some(ref mut f) = self.front {
239            f.clip_to(other);
240        }
241        if let Some(ref mut b) = self.back {
242            b.clip_to(other);
243        }
244    }
245
246    /// Collect all polygons stored in this tree.
247    fn all_polygons(&self) -> Vec<Polygon> {
248        let mut out = self.polygons.clone();
249        if let Some(ref f) = self.front {
250            out.extend(f.all_polygons());
251        }
252        if let Some(ref b) = self.back {
253            out.extend(b.all_polygons());
254        }
255        out
256    }
257
258    /// Insert `polygons` into the BSP tree (building the tree as needed).
259    fn build(&mut self, polygons: Vec<Polygon>) {
260        if polygons.is_empty() {
261            return;
262        }
263
264        if self.plane.is_none() {
265            self.plane = Some(polygons[0].plane.clone());
266        }
267        let plane = self.plane.clone().unwrap();
268
269        let mut front = Vec::new();
270        let mut back = Vec::new();
271
272        for poly in polygons {
273            let mut cf = Vec::new();
274            let mut cb = Vec::new();
275            let mut pf = Vec::new();
276            let mut pb = Vec::new();
277            plane.split_polygon(&poly, &mut cf, &mut cb, &mut pf, &mut pb);
278            self.polygons.extend(cf);
279            self.polygons.extend(cb);
280            front.extend(pf);
281            back.extend(pb);
282        }
283
284        if !front.is_empty() {
285            if self.front.is_none() {
286                self.front = Some(Box::new(BspNode::new()));
287            }
288            self.front.as_mut().unwrap().build(front);
289        }
290        if !back.is_empty() {
291            if self.back.is_none() {
292                self.back = Some(Box::new(BspNode::new()));
293            }
294            self.back.as_mut().unwrap().build(back);
295        }
296    }
297}
298
299// ---------------------------------------------------------------------------
300// Mesh ↔ polygon conversion helpers
301// ---------------------------------------------------------------------------
302
303fn mesh_to_polygons(mesh: &TriangleMesh) -> Vec<Polygon> {
304    mesh.faces
305        .iter()
306        .filter_map(|face| {
307            Polygon::new(vec![
308                mesh.vertices[face[0]],
309                mesh.vertices[face[1]],
310                mesh.vertices[face[2]],
311            ])
312        })
313        .collect()
314}
315
316fn polygons_to_mesh(polygons: Vec<Polygon>) -> TriangleMesh {
317    // Each Polygon is already a triangle (fan-triangulation keeps them as triangles),
318    // so we can flatten directly.
319    let capacity = polygons
320        .iter()
321        .map(|p| p.vertices.len().saturating_sub(2))
322        .sum::<usize>();
323    let mut vertices = Vec::with_capacity(capacity * 3);
324    let mut faces = Vec::with_capacity(capacity);
325
326    for poly in &polygons {
327        let n = poly.vertices.len();
328        if n < 3 {
329            continue;
330        }
331        // Fan triangulation
332        for i in 1..n - 1 {
333            let base = vertices.len();
334            vertices.push(poly.vertices[0]);
335            vertices.push(poly.vertices[i]);
336            vertices.push(poly.vertices[i + 1]);
337            faces.push([base, base + 1, base + 2]);
338        }
339    }
340
341    TriangleMesh::from_vertices_and_faces(vertices, faces)
342}
343
344fn build_bsp(mesh: &TriangleMesh) -> BspNode {
345    let mut node = BspNode::new();
346    node.build(mesh_to_polygons(mesh));
347    node
348}
349
350// ---------------------------------------------------------------------------
351// Public API
352// ---------------------------------------------------------------------------
353
354/// Which boolean operation to perform.
355#[derive(Debug, Clone, Copy, PartialEq, Eq)]
356pub enum BooleanOp {
357    /// All regions inside A **or** B.
358    Union,
359    /// All regions inside both A **and** B.
360    Intersection,
361    /// All regions inside A but **not** B.
362    Difference,
363}
364
365/// Perform a mesh boolean operation on two closed triangle meshes.
366///
367/// # Arguments
368/// * `a` - First operand mesh (must be closed / watertight)
369/// * `b` - Second operand mesh (must be closed / watertight)
370/// * `op` - The boolean operation to apply
371///
372/// # Returns
373/// A new `TriangleMesh` representing the result of the operation.
374pub fn mesh_boolean(a: &TriangleMesh, b: &TriangleMesh, op: BooleanOp) -> Result<TriangleMesh> {
375    match op {
376        BooleanOp::Union => mesh_union(a, b),
377        BooleanOp::Intersection => mesh_intersection(a, b),
378        BooleanOp::Difference => mesh_difference(a, b),
379    }
380}
381
382/// Compute the boolean **union** of two closed triangle meshes.
383///
384/// Returns a mesh containing all regions inside either `a` or `b`.
385///
386/// ```text
387///     A.union(B)
388///
389///     +-------+            +-------+
390///     |       |            |       |
391///     |   A   |            |       |
392///     |    +--+----+   =   |       +----+
393///     +----+--+    |       +----+       |
394///          |   B   |            |       |
395///          |       |            |       |
396///          +-------+            +-------+
397/// ```
398pub fn mesh_union(a: &TriangleMesh, b: &TriangleMesh) -> Result<TriangleMesh> {
399    if a.is_empty() {
400        return Ok(b.clone());
401    }
402    if b.is_empty() {
403        return Ok(a.clone());
404    }
405
406    let mut a_bsp = build_bsp(a);
407    let mut b_bsp = build_bsp(b);
408
409    a_bsp.clip_to(&b_bsp);
410    b_bsp.clip_to(&a_bsp);
411    b_bsp.invert();
412    b_bsp.clip_to(&a_bsp);
413    b_bsp.invert();
414    a_bsp.build(b_bsp.all_polygons());
415
416    Ok(polygons_to_mesh(a_bsp.all_polygons()))
417}
418
419/// Compute the boolean **intersection** of two closed triangle meshes.
420///
421/// Returns a mesh containing only the regions inside both `a` and `b`.
422///
423/// ```text
424///     A.intersection(B)
425///
426///     +-------+
427///     |       |
428///     |   A   |
429///     |    +--+----+   =   +--+
430///     +----+--+    |       +--+
431///          |   B   |
432///          |       |
433///          +-------+
434/// ```
435pub fn mesh_intersection(a: &TriangleMesh, b: &TriangleMesh) -> Result<TriangleMesh> {
436    if a.is_empty() || b.is_empty() {
437        return Ok(TriangleMesh::new());
438    }
439
440    let mut a_bsp = build_bsp(a);
441    let mut b_bsp = build_bsp(b);
442
443    a_bsp.invert();
444    b_bsp.clip_to(&a_bsp);
445    b_bsp.invert();
446    a_bsp.clip_to(&b_bsp);
447    b_bsp.clip_to(&a_bsp);
448    a_bsp.build(b_bsp.all_polygons());
449    a_bsp.invert();
450
451    Ok(polygons_to_mesh(a_bsp.all_polygons()))
452}
453
454/// Compute the boolean **difference** of two closed triangle meshes (A minus B).
455///
456/// Returns a mesh containing only the regions inside `a` but not `b`.
457///
458/// ```text
459///     A.difference(B)
460///
461///     +-------+            +-------+
462///     |       |            |       |
463///     |   A   |            |   A   |
464///     |    +--+----+   =   |    +--+
465///     +----+--+    |       +----+
466///          |   B   |
467///          |       |
468///          +-------+
469/// ```
470pub fn mesh_difference(a: &TriangleMesh, b: &TriangleMesh) -> Result<TriangleMesh> {
471    if a.is_empty() {
472        return Ok(TriangleMesh::new());
473    }
474    if b.is_empty() {
475        return Ok(a.clone());
476    }
477
478    let mut a_bsp = build_bsp(a);
479    let mut b_bsp = build_bsp(b);
480
481    a_bsp.invert();
482    a_bsp.clip_to(&b_bsp);
483    b_bsp.clip_to(&a_bsp);
484    b_bsp.invert();
485    b_bsp.clip_to(&a_bsp);
486    b_bsp.invert();
487    a_bsp.build(b_bsp.all_polygons());
488    a_bsp.invert();
489
490    Ok(polygons_to_mesh(a_bsp.all_polygons()))
491}
492
493// ---------------------------------------------------------------------------
494// Tests
495// ---------------------------------------------------------------------------
496
497#[cfg(test)]
498mod tests {
499    use super::*;
500
501    /// Build a closed box mesh with outward-facing normals.
502    /// Vertices are the 8 corners of [min, max]^3.
503    fn make_box(min: Point3f, max: Point3f) -> TriangleMesh {
504        let (x0, y0, z0) = (min.x, min.y, min.z);
505        let (x1, y1, z1) = (max.x, max.y, max.z);
506        let verts = vec![
507            Point3f::new(x0, y0, z0), // 0
508            Point3f::new(x1, y0, z0), // 1
509            Point3f::new(x1, y1, z0), // 2
510            Point3f::new(x0, y1, z0), // 3
511            Point3f::new(x0, y0, z1), // 4
512            Point3f::new(x1, y0, z1), // 5
513            Point3f::new(x1, y1, z1), // 6
514            Point3f::new(x0, y1, z1), // 7
515        ];
516        // Each face is CCW when viewed from outside → outward normal
517        let faces = vec![
518            // Bottom (z-): normal (0,0,-1)
519            [0, 2, 1],
520            [0, 3, 2],
521            // Top (z+): normal (0,0,+1)
522            [4, 5, 6],
523            [4, 6, 7],
524            // Front (y-): normal (0,-1,0)
525            [0, 1, 5],
526            [0, 5, 4],
527            // Back (y+): normal (0,+1,0)
528            [3, 7, 6],
529            [3, 6, 2],
530            // Left (x-): normal (-1,0,0)
531            [0, 4, 7],
532            [0, 7, 3],
533            // Right (x+): normal (+1,0,0)
534            [1, 2, 6],
535            [1, 6, 5],
536        ];
537        TriangleMesh::from_vertices_and_faces(verts, faces)
538    }
539
540    // ---- empty-mesh edge cases ----
541
542    #[test]
543    fn test_union_empty_a() {
544        let empty = TriangleMesh::new();
545        let b = make_box(Point3f::new(0.0, 0.0, 0.0), Point3f::new(1.0, 1.0, 1.0));
546        let result = mesh_union(&empty, &b).unwrap();
547        assert!(!result.is_empty());
548    }
549
550    #[test]
551    fn test_union_empty_b() {
552        let a = make_box(Point3f::new(0.0, 0.0, 0.0), Point3f::new(1.0, 1.0, 1.0));
553        let empty = TriangleMesh::new();
554        let result = mesh_union(&a, &empty).unwrap();
555        assert!(!result.is_empty());
556    }
557
558    #[test]
559    fn test_intersection_empty() {
560        let empty = TriangleMesh::new();
561        let b = make_box(Point3f::new(0.0, 0.0, 0.0), Point3f::new(1.0, 1.0, 1.0));
562        let result = mesh_intersection(&empty, &b).unwrap();
563        assert!(result.is_empty());
564    }
565
566    #[test]
567    fn test_difference_empty_a() {
568        let empty = TriangleMesh::new();
569        let b = make_box(Point3f::new(0.0, 0.0, 0.0), Point3f::new(1.0, 1.0, 1.0));
570        let result = mesh_difference(&empty, &b).unwrap();
571        assert!(result.is_empty());
572    }
573
574    #[test]
575    fn test_difference_empty_b() {
576        let a = make_box(Point3f::new(0.0, 0.0, 0.0), Point3f::new(1.0, 1.0, 1.0));
577        let empty = TriangleMesh::new();
578        let result = mesh_difference(&a, &empty).unwrap();
579        assert!(!result.is_empty());
580    }
581
582    // ---- non-overlapping meshes ----
583
584    #[test]
585    fn test_union_non_overlapping() {
586        // Two unit cubes separated along X
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_union(&a, &b).unwrap();
590        // Both boxes should be fully preserved
591        assert!(!result.is_empty());
592        assert!(
593            result.face_count() >= 24,
594            "expected ≥24 faces, got {}",
595            result.face_count()
596        );
597    }
598
599    #[test]
600    fn test_intersection_non_overlapping() {
601        let a = make_box(Point3f::new(0.0, 0.0, 0.0), Point3f::new(1.0, 1.0, 1.0));
602        let b = make_box(Point3f::new(5.0, 0.0, 0.0), Point3f::new(6.0, 1.0, 1.0));
603        let result = mesh_intersection(&a, &b).unwrap();
604        // No overlap → empty result
605        assert!(
606            result.is_empty(),
607            "non-overlapping intersection should be empty, got {} faces",
608            result.face_count()
609        );
610    }
611
612    #[test]
613    fn test_difference_non_overlapping() {
614        let a = make_box(Point3f::new(0.0, 0.0, 0.0), Point3f::new(1.0, 1.0, 1.0));
615        let b = make_box(Point3f::new(5.0, 0.0, 0.0), Point3f::new(6.0, 1.0, 1.0));
616        let result = mesh_difference(&a, &b).unwrap();
617        // B doesn't touch A → result should be A unchanged
618        assert!(!result.is_empty());
619        assert!(
620            result.face_count() >= 12,
621            "expected ≥12 faces, got {}",
622            result.face_count()
623        );
624    }
625
626    // ---- overlapping meshes ----
627
628    #[test]
629    fn test_union_overlapping() {
630        // Two unit cubes that partially overlap
631        let a = make_box(Point3f::new(0.0, 0.0, 0.0), Point3f::new(2.0, 2.0, 2.0));
632        let b = make_box(Point3f::new(1.0, 1.0, 1.0), Point3f::new(3.0, 3.0, 3.0));
633        let result = mesh_union(&a, &b).unwrap();
634        assert!(
635            !result.is_empty(),
636            "union of overlapping boxes must not be empty"
637        );
638    }
639
640    #[test]
641    fn test_intersection_overlapping() {
642        let a = make_box(Point3f::new(0.0, 0.0, 0.0), Point3f::new(2.0, 2.0, 2.0));
643        let b = make_box(Point3f::new(1.0, 1.0, 1.0), Point3f::new(3.0, 3.0, 3.0));
644        let result = mesh_intersection(&a, &b).unwrap();
645        assert!(
646            !result.is_empty(),
647            "intersection of overlapping boxes must not be empty"
648        );
649        // The intersection is the unit cube [1,2]^3 → should have at most as many faces as either input
650        assert!(result.face_count() <= a.face_count() + b.face_count() + 50);
651    }
652
653    #[test]
654    fn test_difference_overlapping() {
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        let result = mesh_difference(&a, &b).unwrap();
658        assert!(
659            !result.is_empty(),
660            "A minus partially-overlapping B must not be empty"
661        );
662    }
663
664    #[test]
665    fn test_difference_a_fully_contains_b() {
666        // B is entirely inside A → difference cuts a hole
667        let a = make_box(Point3f::new(0.0, 0.0, 0.0), Point3f::new(4.0, 4.0, 4.0));
668        let b = make_box(Point3f::new(1.0, 1.0, 1.0), Point3f::new(2.0, 2.0, 2.0));
669        let result = mesh_difference(&a, &b).unwrap();
670        assert!(
671            !result.is_empty(),
672            "difference with interior hole must not be empty"
673        );
674    }
675
676    #[test]
677    fn test_intersection_a_fully_contains_b() {
678        // B fully inside A → intersection == B
679        let a = make_box(Point3f::new(0.0, 0.0, 0.0), Point3f::new(4.0, 4.0, 4.0));
680        let b = make_box(Point3f::new(1.0, 1.0, 1.0), Point3f::new(2.0, 2.0, 2.0));
681        let result = mesh_intersection(&a, &b).unwrap();
682        assert!(
683            !result.is_empty(),
684            "intersection when B inside A must not be empty"
685        );
686    }
687
688    // ---- BooleanOp dispatch ----
689
690    #[test]
691    fn test_mesh_boolean_dispatch() {
692        let a = make_box(Point3f::new(0.0, 0.0, 0.0), Point3f::new(2.0, 2.0, 2.0));
693        let b = make_box(Point3f::new(1.0, 1.0, 1.0), Point3f::new(3.0, 3.0, 3.0));
694
695        let u = mesh_boolean(&a, &b, BooleanOp::Union).unwrap();
696        let i = mesh_boolean(&a, &b, BooleanOp::Intersection).unwrap();
697        let d = mesh_boolean(&a, &b, BooleanOp::Difference).unwrap();
698
699        assert!(!u.is_empty());
700        assert!(!i.is_empty());
701        assert!(!d.is_empty());
702    }
703}
threecrate_algorithms/mesh_boolean.rs

threecrate_algorithms/
mesh_boolean.rs
