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 if !self.delaunay.can_add_constraint(handle1, handle2) {
53 let _ = self.delaunay.add_constraint(handle1, handle2);
54 }
55 }
56}
57
58impl TriMesh {
59 pub fn canonical_split(&self, axis: usize, bias: Real, epsilon: Real) -> SplitResult<Self> {
70 self.local_split(&Vector::ith_axis(axis), bias, epsilon)
72 }
73
74 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 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 let vertices = self.vertices();
104 let indices = self.indices();
105 let mut colors = vec![0u8; self.vertices().len()];
106
107 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 if !found_negative {
125 return SplitResult::Positive;
126 }
127
128 if !found_positive {
129 return SplitResult::Negative;
130 }
131
132 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 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 (0, _) => FeatureId::Vertex(ia),
150 _ => continue,
151 };
152
153 if intersection_features.0 == FeatureId::Unknown {
154 intersection_features.0 = fid;
155 } else {
156 intersection_features.1 = fid;
159 }
160 }
161
162 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 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 match intersection_features {
188 (_, FeatureId::Unknown) => {
189 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 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 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 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 if e2 != (e1 + 1) % 3 {
234 core::mem::swap(&mut e1, &mut e2);
235 }
236
237 let ia = e2; 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 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]; 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 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 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); indices_rhs.push(idx2);
330 }
331 }
332 }
333
334 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 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 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 pub fn intersection_with_local_plane(
384 &self,
385 local_axis: &UnitVector<Real>,
386 bias: Real,
387 epsilon: Real,
388 ) -> IntersectResult<Polyline> {
389 let vertices = self.vertices();
391 let indices = self.indices();
392 let mut colors = vec![0u8; self.vertices().len()];
393
394 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 if !found_negative {
412 return IntersectResult::Positive;
413 }
414
415 if !found_positive {
416 return IntersectResult::Negative;
417 }
418
419 let mut index_adjacencies: Vec<Vec<usize>> = Vec::new(); 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 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 (0, _) => FeatureId::Vertex(ia),
459 _ => continue,
460 };
461
462 if intersection_features.0 == FeatureId::Unknown {
463 intersection_features.0 = fid;
464 } else {
465 intersection_features.1 = fid;
468 }
469 }
470
471 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 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 match intersection_features {
495 (_, FeatureId::Unknown) => {
496 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 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 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 if e2 != (e1 + 1) % 3 {
552 core::mem::swap(&mut e1, &mut e2);
553 }
554
555 let ia = e2; 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 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 let first = idx;
581 let mut prev = first;
582 let mut next = neighbors.first(); '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 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 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 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 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 .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 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 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 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}