parry3d/query/split/
split_trimesh.rs

1use crate::bounding_volume::Aabb;
2use crate::math::{Isometry, Point, Real, UnitVector, Vector};
3use crate::query::{IntersectResult, PointQuery, SplitResult};
4use crate::shape::{Cuboid, FeatureId, Polyline, Segment, Shape, TriMesh, TriMeshFlags, Triangle};
5use crate::transformation::{intersect_meshes, MeshIntersectionError};
6use crate::utils::{hashmap::HashMap, SortedPair, WBasis};
7use alloc::{vec, vec::Vec};
8use core::cmp::Ordering;
9use spade::{handles::FixedVertexHandle, ConstrainedDelaunayTriangulation, Triangulation as _};
10
11struct Triangulation {
12    delaunay: ConstrainedDelaunayTriangulation<spade::Point2<Real>>,
13    basis: [Vector<Real>; 2],
14    basis_origin: Point<Real>,
15    spade2index: HashMap<FixedVertexHandle, u32>,
16    index2spade: HashMap<u32, FixedVertexHandle>,
17}
18
19impl Triangulation {
20    fn new(axis: UnitVector<Real>, basis_origin: Point<Real>) -> Self {
21        Triangulation {
22            delaunay: ConstrainedDelaunayTriangulation::new(),
23            basis: axis.orthonormal_basis(),
24            basis_origin,
25            spade2index: HashMap::default(),
26            index2spade: HashMap::default(),
27        }
28    }
29
30    fn project(&self, pt: Point<Real>) -> spade::Point2<Real> {
31        let dpt = pt - self.basis_origin;
32        spade::Point2::new(dpt.dot(&self.basis[0]), dpt.dot(&self.basis[1]))
33    }
34
35    fn add_edge(&mut self, id1: u32, id2: u32, points: &[Point<Real>]) {
36        let proj1 = self.project(points[id1 as usize]);
37        let proj2 = self.project(points[id2 as usize]);
38
39        let handle1 = *self.index2spade.entry(id1).or_insert_with(|| {
40            let h = self.delaunay.insert(proj1).unwrap();
41            let _ = self.spade2index.insert(h, id1);
42            h
43        });
44
45        let handle2 = *self.index2spade.entry(id2).or_insert_with(|| {
46            let h = self.delaunay.insert(proj2).unwrap();
47            let _ = self.spade2index.insert(h, id2);
48            h
49        });
50
51        // NOTE: the naming of the `ConstrainedDelaunayTriangulation::can_add_constraint` method is misleading.
52        if !self.delaunay.can_add_constraint(handle1, handle2) {
53            let _ = self.delaunay.add_constraint(handle1, handle2);
54        }
55    }
56}
57
58impl TriMesh {
59    /// Splits this `TriMesh` along the given canonical axis.
60    ///
61    /// This will split the Aabb by a plane with a normal with it’s `axis`-th component set to 1.
62    /// The splitting plane is shifted wrt. the origin by the `bias` (i.e. it passes through the point
63    /// equal to `normal * bias`).
64    ///
65    /// # Result
66    /// Returns the result of the split. The first mesh returned is the piece lying on the negative
67    /// half-space delimited by the splitting plane. The second mesh returned is the piece lying on the
68    /// positive half-space delimited by the splitting plane.
69    pub fn canonical_split(&self, axis: usize, bias: Real, epsilon: Real) -> SplitResult<Self> {
70        // TODO: optimize this.
71        self.local_split(&Vector::ith_axis(axis), bias, epsilon)
72    }
73
74    /// Splits this mesh, transformed by `position` by a plane identified by its normal `local_axis`
75    /// and the `bias` (i.e. the plane passes through the point equal to `normal * bias`).
76    pub fn split(
77        &self,
78        position: &Isometry<Real>,
79        axis: &UnitVector<Real>,
80        bias: Real,
81        epsilon: Real,
82    ) -> SplitResult<Self> {
83        let local_axis = position.inverse_transform_unit_vector(axis);
84        let added_bias = -position.translation.vector.dot(axis);
85        self.local_split(&local_axis, bias + added_bias, epsilon)
86    }
87
88    /// Splits this mesh by a plane identified by its normal `local_axis`
89    /// and the `bias` (i.e. the plane passes through the point equal to `normal * bias`).
90    pub fn local_split(
91        &self,
92        local_axis: &UnitVector<Real>,
93        bias: Real,
94        epsilon: Real,
95    ) -> SplitResult<Self> {
96        let mut triangulation = if self.pseudo_normals_if_oriented().is_some() {
97            Some(Triangulation::new(*local_axis, self.vertices()[0]))
98        } else {
99            None
100        };
101
102        // 1. Partition the vertices.
103        let vertices = self.vertices();
104        let indices = self.indices();
105        let mut colors = vec![0u8; self.vertices().len()];
106
107        // Color 0 = on plane.
108        //       1 = on negative half-space.
109        //       2 = on positive half-space.
110        let mut found_negative = false;
111        let mut found_positive = false;
112        for (i, pt) in vertices.iter().enumerate() {
113            let dist_to_plane = pt.coords.dot(local_axis) - bias;
114            if dist_to_plane < -epsilon {
115                found_negative = true;
116                colors[i] = 1;
117            } else if dist_to_plane > epsilon {
118                found_positive = true;
119                colors[i] = 2;
120            }
121        }
122
123        // Exit early if `self` isn’t crossed by the plane.
124        if !found_negative {
125            return SplitResult::Positive;
126        }
127
128        if !found_positive {
129            return SplitResult::Negative;
130        }
131
132        // 2. Split the triangles.
133        let mut intersections_found = HashMap::default();
134        let mut new_indices = indices.to_vec();
135        let mut new_vertices = vertices.to_vec();
136
137        for (tri_id, idx) in indices.iter().enumerate() {
138            let mut intersection_features = (FeatureId::Unknown, FeatureId::Unknown);
139
140            // First, find where the plane intersects the triangle.
141            for ia in 0..3 {
142                let ib = (ia + 1) % 3;
143                let idx_a = idx[ia as usize];
144                let idx_b = idx[ib as usize];
145
146                let fid = match (colors[idx_a as usize], colors[idx_b as usize]) {
147                    (1, 2) | (2, 1) => FeatureId::Edge(ia),
148                    // NOTE: the case (_, 0) will be dealt with in the next loop iteration.
149                    (0, _) => FeatureId::Vertex(ia),
150                    _ => continue,
151                };
152
153                if intersection_features.0 == FeatureId::Unknown {
154                    intersection_features.0 = fid;
155                } else {
156                    // TODO: this assertion may fire if the triangle is coplanar with the edge?
157                    // assert_eq!(intersection_features.1, FeatureId::Unknown);
158                    intersection_features.1 = fid;
159                }
160            }
161
162            // Helper that intersects an edge with the plane.
163            let mut intersect_edge = |idx_a, idx_b| {
164                *intersections_found
165                    .entry(SortedPair::new(idx_a, idx_b))
166                    .or_insert_with(|| {
167                        let segment = Segment::new(
168                            new_vertices[idx_a as usize],
169                            new_vertices[idx_b as usize],
170                        );
171                        // Intersect the segment with the plane.
172                        if let Some((intersection, _)) = segment
173                            .local_split_and_get_intersection(local_axis, bias, epsilon)
174                            .1
175                        {
176                            new_vertices.push(intersection);
177                            colors.push(0);
178                            (new_vertices.len() - 1) as u32
179                        } else {
180                            unreachable!()
181                        }
182                    })
183            };
184
185            // Perform the intersection, push new triangles, and update
186            // triangulation constraints if needed.
187            match intersection_features {
188                (_, FeatureId::Unknown) => {
189                    // The plane doesn’t intersect the triangle, or intersects it at
190                    // a single vertex, so we don’t have anything to do.
191                    assert!(
192                        matches!(intersection_features.0, FeatureId::Unknown)
193                            || matches!(intersection_features.0, FeatureId::Vertex(_))
194                    );
195                }
196                (FeatureId::Vertex(v1), FeatureId::Vertex(v2)) => {
197                    // The plane intersects the triangle along one of its edge.
198                    // We don’t have to split the triangle, but we need to add
199                    // a constraint to the triangulation.
200                    if let Some(triangulation) = &mut triangulation {
201                        let id1 = idx[v1 as usize];
202                        let id2 = idx[v2 as usize];
203                        triangulation.add_edge(id1, id2, &new_vertices);
204                    }
205                }
206                (FeatureId::Vertex(iv), FeatureId::Edge(ie))
207                | (FeatureId::Edge(ie), FeatureId::Vertex(iv)) => {
208                    // The plane splits the triangle into exactly two triangles.
209                    let ia = ie;
210                    let ib = (ie + 1) % 3;
211                    let ic = (ie + 2) % 3;
212                    let idx_a = idx[ia as usize];
213                    let idx_b = idx[ib as usize];
214                    let idx_c = idx[ic as usize];
215                    assert_eq!(iv, ic);
216
217                    let intersection_idx = intersect_edge(idx_a, idx_b);
218
219                    // Compute the indices of the two triangles.
220                    let new_tri_a = [idx_c, idx_a, intersection_idx];
221                    let new_tri_b = [idx_b, idx_c, intersection_idx];
222
223                    new_indices[tri_id] = new_tri_a;
224                    new_indices.push(new_tri_b);
225
226                    if let Some(triangulation) = &mut triangulation {
227                        triangulation.add_edge(intersection_idx, idx_c, &new_vertices);
228                    }
229                }
230                (FeatureId::Edge(mut e1), FeatureId::Edge(mut e2)) => {
231                    // The plane splits the triangle into 1 + 2 triangles.
232                    // First, make sure the edge indices are consecutive.
233                    if e2 != (e1 + 1) % 3 {
234                        core::mem::swap(&mut e1, &mut e2);
235                    }
236
237                    let ia = e2; // The first point of the second edge is the vertex shared by both edges.
238                    let ib = (e2 + 1) % 3;
239                    let ic = (e2 + 2) % 3;
240                    let idx_a = idx[ia as usize];
241                    let idx_b = idx[ib as usize];
242                    let idx_c = idx[ic as usize];
243
244                    let intersection1 = intersect_edge(idx_c, idx_a);
245                    let intersection2 = intersect_edge(idx_a, idx_b);
246
247                    let new_tri1 = [idx_a, intersection2, intersection1];
248                    let new_tri2 = [intersection2, idx_b, idx_c];
249                    let new_tri3 = [intersection2, idx_c, intersection1];
250                    new_indices[tri_id] = new_tri1;
251                    new_indices.push(new_tri2);
252                    new_indices.push(new_tri3);
253
254                    if let Some(triangulation) = &mut triangulation {
255                        triangulation.add_edge(intersection1, intersection2, &new_vertices);
256                    }
257                }
258                _ => unreachable!(),
259            }
260        }
261
262        // 3. Partition the new triangles into two trimeshes.
263        let mut vertices_lhs = vec![];
264        let mut vertices_rhs = vec![];
265        let mut indices_lhs = vec![];
266        let mut indices_rhs = vec![];
267        let mut remap = vec![];
268
269        for i in 0..new_vertices.len() {
270            match colors[i] {
271                0 => {
272                    remap.push((vertices_lhs.len() as u32, vertices_rhs.len() as u32));
273                    vertices_lhs.push(new_vertices[i]);
274                    vertices_rhs.push(new_vertices[i]);
275                }
276                1 => {
277                    remap.push((vertices_lhs.len() as u32, u32::MAX));
278                    vertices_lhs.push(new_vertices[i]);
279                }
280                2 => {
281                    remap.push((u32::MAX, vertices_rhs.len() as u32));
282                    vertices_rhs.push(new_vertices[i]);
283                }
284                _ => unreachable!(),
285            }
286        }
287
288        for idx in new_indices {
289            let idx = [idx[0] as usize, idx[1] as usize, idx[2] as usize]; // Convert to usize.
290            let colors = [colors[idx[0]], colors[idx[1]], colors[idx[2]]];
291            let remap = [remap[idx[0]], remap[idx[1]], remap[idx[2]]];
292
293            if colors[0] == 1 || colors[1] == 1 || colors[2] == 1 {
294                assert!(colors[0] != 2 && colors[1] != 2 && colors[2] != 2);
295                indices_lhs.push([remap[0].0, remap[1].0, remap[2].0]);
296            } else if colors[0] == 2 || colors[1] == 2 || colors[2] == 2 {
297                assert!(colors[0] != 1 && colors[1] != 1 && colors[2] != 1);
298                indices_rhs.push([remap[0].1, remap[1].1, remap[2].1]);
299            } else {
300                // The colors are all 0, so push into both trimeshes.
301                indices_lhs.push([remap[0].0, remap[1].0, remap[2].0]);
302                indices_rhs.push([remap[0].1, remap[1].1, remap[2].1]);
303            }
304        }
305
306        // Push the triangulation if there is one.
307        if let Some(triangulation) = triangulation {
308            for face in triangulation.delaunay.inner_faces() {
309                let vtx = face.vertices();
310                let mut idx1 = [0; 3];
311                let mut idx2 = [0; 3];
312                for k in 0..3 {
313                    let vid = triangulation.spade2index[&vtx[k].fix()];
314                    assert_eq!(colors[vid as usize], 0);
315                    idx1[k] = remap[vid as usize].0;
316                    idx2[k] = remap[vid as usize].1;
317                }
318
319                let tri = Triangle::new(
320                    vertices_lhs[idx1[0] as usize],
321                    vertices_lhs[idx1[1] as usize],
322                    vertices_lhs[idx1[2] as usize],
323                );
324
325                if self.contains_local_point(&tri.center()) {
326                    indices_lhs.push(idx1);
327
328                    idx2.swap(1, 2); // Flip orientation for the second half of the split.
329                    indices_rhs.push(idx2);
330                }
331            }
332        }
333
334        // TODO: none of the index buffers should be empty at this point unless perhaps
335        //       because of some rounding errors?
336        //       Should we just panic if they are empty?
337        if indices_rhs.is_empty() {
338            SplitResult::Negative
339        } else if indices_lhs.is_empty() {
340            SplitResult::Positive
341        } else {
342            let mesh_lhs = TriMesh::new(vertices_lhs, indices_lhs).unwrap();
343            let mesh_rhs = TriMesh::new(vertices_rhs, indices_rhs).unwrap();
344            SplitResult::Pair(mesh_lhs, mesh_rhs)
345        }
346    }
347
348    /// Computes the intersection [`Polyline`]s between this mesh and the plane identified by
349    /// the given canonical axis.
350    ///
351    /// This will intersect the mesh by a plane with a normal with it’s `axis`-th component set to 1.
352    /// The splitting plane is shifted wrt. the origin by the `bias` (i.e. it passes through the point
353    /// equal to `normal * bias`).
354    ///
355    /// Note that the resultant polyline may have multiple connected components
356    pub fn canonical_intersection_with_plane(
357        &self,
358        axis: usize,
359        bias: Real,
360        epsilon: Real,
361    ) -> IntersectResult<Polyline> {
362        self.intersection_with_local_plane(&Vector::ith_axis(axis), bias, epsilon)
363    }
364
365    /// Computes the intersection [`Polyline`]s between this mesh, transformed by `position`,
366    /// and a plane identified by its normal `axis` and the `bias`
367    /// (i.e. the plane passes through the point equal to `normal * bias`).
368    pub fn intersection_with_plane(
369        &self,
370        position: &Isometry<Real>,
371        axis: &UnitVector<Real>,
372        bias: Real,
373        epsilon: Real,
374    ) -> IntersectResult<Polyline> {
375        let local_axis = position.inverse_transform_unit_vector(axis);
376        let added_bias = -position.translation.vector.dot(axis);
377        self.intersection_with_local_plane(&local_axis, bias + added_bias, epsilon)
378    }
379
380    /// Computes the intersection [`Polyline`]s between this mesh
381    /// and a plane identified by its normal `local_axis`
382    /// and the `bias` (i.e. the plane passes through the point equal to `normal * bias`).
383    pub fn intersection_with_local_plane(
384        &self,
385        local_axis: &UnitVector<Real>,
386        bias: Real,
387        epsilon: Real,
388    ) -> IntersectResult<Polyline> {
389        // 1. Partition the vertices.
390        let vertices = self.vertices();
391        let indices = self.indices();
392        let mut colors = vec![0u8; self.vertices().len()];
393
394        // Color 0 = on plane.
395        //       1 = on negative half-space.
396        //       2 = on positive half-space.
397        let mut found_negative = false;
398        let mut found_positive = false;
399        for (i, pt) in vertices.iter().enumerate() {
400            let dist_to_plane = pt.coords.dot(local_axis) - bias;
401            if dist_to_plane < -epsilon {
402                found_negative = true;
403                colors[i] = 1;
404            } else if dist_to_plane > epsilon {
405                found_positive = true;
406                colors[i] = 2;
407            }
408        }
409
410        // Exit early if `self` isn’t crossed by the plane.
411        if !found_negative {
412            return IntersectResult::Positive;
413        }
414
415        if !found_positive {
416            return IntersectResult::Negative;
417        }
418
419        // 2. Split the triangles.
420        let mut index_adjacencies: Vec<Vec<usize>> = Vec::new(); // Adjacency list of indices
421
422        // Helper functions for adding polyline segments to the adjacency list
423        let mut add_segment_adjacencies = |idx_a: usize, idx_b| {
424            assert!(idx_a <= index_adjacencies.len());
425
426            match idx_a.cmp(&index_adjacencies.len()) {
427                Ordering::Less => index_adjacencies[idx_a].push(idx_b),
428                Ordering::Equal => index_adjacencies.push(vec![idx_b]),
429                Ordering::Greater => {}
430            }
431        };
432        let mut add_segment_adjacencies_symmetric = |idx_a: usize, idx_b| {
433            if idx_a < idx_b {
434                add_segment_adjacencies(idx_a, idx_b);
435                add_segment_adjacencies(idx_b, idx_a);
436            } else {
437                add_segment_adjacencies(idx_b, idx_a);
438                add_segment_adjacencies(idx_a, idx_b);
439            }
440        };
441
442        let mut intersections_found = HashMap::default();
443        let mut existing_vertices_found = HashMap::default();
444        let mut new_vertices = Vec::new();
445
446        for idx in indices.iter() {
447            let mut intersection_features = (FeatureId::Unknown, FeatureId::Unknown);
448
449            // First, find where the plane intersects the triangle.
450            for ia in 0..3 {
451                let ib = (ia + 1) % 3;
452                let idx_a = idx[ia as usize];
453                let idx_b = idx[ib as usize];
454
455                let fid = match (colors[idx_a as usize], colors[idx_b as usize]) {
456                    (1, 2) | (2, 1) => FeatureId::Edge(ia),
457                    // NOTE: the case (_, 0) will be dealt with in the next loop iteration.
458                    (0, _) => FeatureId::Vertex(ia),
459                    _ => continue,
460                };
461
462                if intersection_features.0 == FeatureId::Unknown {
463                    intersection_features.0 = fid;
464                } else {
465                    // TODO: this assertion may fire if the triangle is coplanar with the edge?
466                    // assert_eq!(intersection_features.1, FeatureId::Unknown);
467                    intersection_features.1 = fid;
468                }
469            }
470
471            // Helper that intersects an edge with the plane.
472            let mut intersect_edge = |idx_a, idx_b| {
473                *intersections_found
474                    .entry(SortedPair::new(idx_a, idx_b))
475                    .or_insert_with(|| {
476                        let segment =
477                            Segment::new(vertices[idx_a as usize], vertices[idx_b as usize]);
478                        // Intersect the segment with the plane.
479                        if let Some((intersection, _)) = segment
480                            .local_split_and_get_intersection(local_axis, bias, epsilon)
481                            .1
482                        {
483                            new_vertices.push(intersection);
484                            colors.push(0);
485                            new_vertices.len() - 1
486                        } else {
487                            unreachable!()
488                        }
489                    })
490            };
491
492            // Perform the intersection, push new triangles, and update
493            // triangulation constraints if needed.
494            match intersection_features {
495                (_, FeatureId::Unknown) => {
496                    // The plane doesn’t intersect the triangle, or intersects it at
497                    // a single vertex, so we don’t have anything to do.
498                    assert!(
499                        matches!(intersection_features.0, FeatureId::Unknown)
500                            || matches!(intersection_features.0, FeatureId::Vertex(_))
501                    );
502                }
503                (FeatureId::Vertex(iv1), FeatureId::Vertex(iv2)) => {
504                    // The plane intersects the triangle along one of its edge.
505                    // We don’t have to split the triangle, but we need to add
506                    // the edge to the polyline indices
507
508                    let id1 = idx[iv1 as usize];
509                    let id2 = idx[iv2 as usize];
510
511                    let out_id1 = *existing_vertices_found.entry(id1).or_insert_with(|| {
512                        let v1 = vertices[id1 as usize];
513
514                        new_vertices.push(v1);
515                        new_vertices.len() - 1
516                    });
517                    let out_id2 = *existing_vertices_found.entry(id2).or_insert_with(|| {
518                        let v2 = vertices[id2 as usize];
519
520                        new_vertices.push(v2);
521                        new_vertices.len() - 1
522                    });
523
524                    add_segment_adjacencies_symmetric(out_id1, out_id2);
525                }
526                (FeatureId::Vertex(iv), FeatureId::Edge(ie))
527                | (FeatureId::Edge(ie), FeatureId::Vertex(iv)) => {
528                    // The plane splits the triangle into exactly two triangles.
529                    let ia = ie;
530                    let ib = (ie + 1) % 3;
531                    let ic = (ie + 2) % 3;
532                    let idx_a = idx[ia as usize];
533                    let idx_b = idx[ib as usize];
534                    let idx_c = idx[ic as usize];
535                    assert_eq!(iv, ic);
536
537                    let intersection_idx = intersect_edge(idx_a, idx_b);
538
539                    let out_idx_c = *existing_vertices_found.entry(idx_c).or_insert_with(|| {
540                        let v2 = vertices[idx_c as usize];
541
542                        new_vertices.push(v2);
543                        new_vertices.len() - 1
544                    });
545
546                    add_segment_adjacencies_symmetric(out_idx_c, intersection_idx);
547                }
548                (FeatureId::Edge(mut e1), FeatureId::Edge(mut e2)) => {
549                    // The plane splits the triangle into 1 + 2 triangles.
550                    // First, make sure the edge indices are consecutive.
551                    if e2 != (e1 + 1) % 3 {
552                        core::mem::swap(&mut e1, &mut e2);
553                    }
554
555                    let ia = e2; // The first point of the second edge is the vertex shared by both edges.
556                    let ib = (e2 + 1) % 3;
557                    let ic = (e2 + 2) % 3;
558                    let idx_a = idx[ia as usize];
559                    let idx_b = idx[ib as usize];
560                    let idx_c = idx[ic as usize];
561
562                    let intersection1 = intersect_edge(idx_c, idx_a);
563                    let intersection2 = intersect_edge(idx_a, idx_b);
564
565                    add_segment_adjacencies_symmetric(intersection1, intersection2);
566                }
567                _ => unreachable!(),
568            }
569        }
570
571        // 3. Ensure consistent edge orientation by traversing the adjacency list
572        let mut polyline_indices: Vec<[u32; 2]> = Vec::with_capacity(index_adjacencies.len() + 1);
573
574        let mut seen = vec![false; index_adjacencies.len()];
575        for (idx, neighbors) in index_adjacencies.iter().enumerate() {
576            if !seen[idx] {
577                // Start a new component
578                // Traverse the adjencies until the loop closes
579
580                let first = idx;
581                let mut prev = first;
582                let mut next = neighbors.first(); // Arbitrary neighbor
583
584                'traversal: while let Some(current) = next {
585                    seen[*current] = true;
586                    polyline_indices.push([prev as u32, *current as u32]);
587
588                    for neighbor in index_adjacencies[*current].iter() {
589                        if *neighbor != prev && *neighbor != first {
590                            prev = *current;
591                            next = Some(neighbor);
592                            continue 'traversal;
593                        } else if *neighbor != prev && *neighbor == first {
594                            // If the next index is same as the first, close the polyline and exit
595                            polyline_indices.push([*current as u32, first as u32]);
596                            next = None;
597                            continue 'traversal;
598                        }
599                    }
600                }
601            }
602        }
603
604        IntersectResult::Intersect(Polyline::new(new_vertices, Some(polyline_indices)))
605    }
606
607    /// Computes the intersection mesh between an Aabb and this mesh.
608    pub fn intersection_with_aabb(
609        &self,
610        position: &Isometry<Real>,
611        flip_mesh: bool,
612        aabb: &Aabb,
613        flip_cuboid: bool,
614        epsilon: Real,
615    ) -> Result<Option<Self>, MeshIntersectionError> {
616        let cuboid = Cuboid::new(aabb.half_extents());
617        let cuboid_pos = Isometry::from(aabb.center());
618        self.intersection_with_cuboid(
619            position,
620            flip_mesh,
621            &cuboid,
622            &cuboid_pos,
623            flip_cuboid,
624            epsilon,
625        )
626    }
627
628    /// Computes the intersection mesh between a cuboid and this mesh transformed by `position`.
629    pub fn intersection_with_cuboid(
630        &self,
631        position: &Isometry<Real>,
632        flip_mesh: bool,
633        cuboid: &Cuboid,
634        cuboid_position: &Isometry<Real>,
635        flip_cuboid: bool,
636        epsilon: Real,
637    ) -> Result<Option<Self>, MeshIntersectionError> {
638        self.intersection_with_local_cuboid(
639            flip_mesh,
640            cuboid,
641            &position.inv_mul(cuboid_position),
642            flip_cuboid,
643            epsilon,
644        )
645    }
646
647    /// Computes the intersection mesh between a cuboid and this mesh.
648    pub fn intersection_with_local_cuboid(
649        &self,
650        flip_mesh: bool,
651        cuboid: &Cuboid,
652        cuboid_position: &Isometry<Real>,
653        flip_cuboid: bool,
654        _epsilon: Real,
655    ) -> Result<Option<Self>, MeshIntersectionError> {
656        if self.topology().is_some() && self.pseudo_normals_if_oriented().is_some() {
657            let (cuboid_vtx, cuboid_idx) = cuboid.to_trimesh();
658            let cuboid_trimesh = TriMesh::with_flags(
659                cuboid_vtx,
660                cuboid_idx,
661                TriMeshFlags::HALF_EDGE_TOPOLOGY | TriMeshFlags::ORIENTED,
662            )
663            // `TriMesh::with_flags` can't fail for `cuboid.to_trimesh()`.
664            .unwrap();
665
666            let intersect_meshes = intersect_meshes(
667                &Isometry::identity(),
668                self,
669                flip_mesh,
670                cuboid_position,
671                &cuboid_trimesh,
672                flip_cuboid,
673            );
674            return intersect_meshes;
675        }
676
677        let cuboid_aabb = cuboid.compute_aabb(cuboid_position);
678        let intersecting_tris: Vec<_> = self.bvh().intersect_aabb(&cuboid_aabb).collect();
679
680        if intersecting_tris.is_empty() {
681            return Ok(None);
682        }
683
684        // First, very naive version that outputs a triangle soup without
685        // index buffer (shared vertices are duplicated).
686        let vertices = self.vertices();
687        let indices = self.indices();
688
689        let mut clip_workspace = vec![];
690        let mut new_vertices = vec![];
691        let mut new_indices = vec![];
692        let aabb = cuboid.local_aabb();
693        let inv_pos = cuboid_position.inverse();
694        let mut to_clip = vec![];
695
696        for tri in intersecting_tris {
697            let idx = indices[tri as usize];
698            to_clip.extend_from_slice(&[
699                inv_pos * vertices[idx[0] as usize],
700                inv_pos * vertices[idx[1] as usize],
701                inv_pos * vertices[idx[2] as usize],
702            ]);
703
704            // There is no need to clip if the triangle is fully inside of the Aabb.
705            // Note that we can’t take a shortcut for the case where all the vertices are
706            // outside of the Aabb, because the Aabb can still instersect the edges or face.
707            if !(aabb.contains_local_point(&to_clip[0])
708                && aabb.contains_local_point(&to_clip[1])
709                && aabb.contains_local_point(&to_clip[2]))
710            {
711                aabb.clip_polygon_with_workspace(&mut to_clip, &mut clip_workspace);
712            }
713
714            if to_clip.len() >= 3 {
715                let base_i = new_vertices.len();
716                for i in 1..to_clip.len() - 1 {
717                    new_indices.push([base_i as u32, (base_i + i) as u32, (base_i + i + 1) as u32]);
718                }
719                new_vertices.append(&mut to_clip);
720            }
721        }
722
723        // The clipping outputs points in the local-space of the cuboid.
724        // So we need to transform it back.
725        for pt in &mut new_vertices {
726            *pt = cuboid_position * *pt;
727        }
728
729        Ok(if new_vertices.len() >= 3 {
730            Some(TriMesh::new(new_vertices, new_indices).unwrap())
731        } else {
732            None
733        })
734    }
735}