1use fenris_traits::Real;
2use nalgebra::allocator::Allocator;
3use nalgebra::{
4 distance_squared, DefaultAllocator, DimName, OPoint, OVector, Point2, Point3, Scalar, Unit, Vector3, U2, U3,
5};
6use numeric_literals::replace_float_literals;
7use serde::{Deserialize, Serialize};
8use std::convert::TryFrom;
9use std::fmt::Debug;
10
11mod polygon;
12mod polytope;
13mod primitives;
14
15use crate::util::index_set_nth_power_iter;
16pub use polygon::*;
17pub use polytope::*;
18pub use primitives::*;
19
20pub mod polymesh;
21pub mod predicates;
22pub mod sdf;
23pub mod util;
24
25#[cfg(feature = "vtkio")]
26pub mod vtkio;
27
28#[cfg(feature = "proptest-support")]
29pub mod proptest;
30
31pub trait BoundedGeometry<T>
32where
33 T: Scalar,
34 DefaultAllocator: Allocator<T, Self::Dimension>,
35{
36 type Dimension: DimName;
37
38 fn bounding_box(&self) -> AxisAlignedBoundingBox<T, Self::Dimension>;
39}
40
41pub trait Distance<T, QueryGeometry>
42where
43 T: Scalar,
44{
45 fn distance_bound(&self, query_geometry: &QueryGeometry) -> [T; 2] {
47 let d = self.distance(query_geometry);
48 [d.clone(), d]
49 }
50
51 fn distance(&self, query_geometry: &QueryGeometry) -> T;
52}
53
54pub trait GeometryCollection<'a> {
55 type Geometry;
56
57 fn num_geometries(&self) -> usize;
58 fn get_geometry(&'a self, index: usize) -> Option<Self::Geometry>;
59}
60
61pub trait DistanceQuery<'a, QueryGeometry>: GeometryCollection<'a> {
62 fn nearest(&'a self, query_geometry: &'a QueryGeometry) -> Option<usize>;
67}
68
69#[derive(Debug, Clone, PartialEq, Eq)]
70pub struct SignedDistanceResult<T, D>
71where
72 T: Scalar,
73 D: DimName,
74 DefaultAllocator: Allocator<T, D>,
75{
76 pub feature_id: usize,
77 pub point: OPoint<T, D>,
78 pub signed_distance: T,
79}
80
81pub trait SignedDistance<T, D>
82where
83 T: Scalar,
84 D: DimName,
85 DefaultAllocator: Allocator<T, D>,
86{
87 fn query_signed_distance(&self, point: &OPoint<T, D>) -> Option<SignedDistanceResult<T, D>>;
88}
89
90#[derive(Debug, Clone, PartialEq, Eq, Serialize, Deserialize)]
91#[serde(bound(
92 serialize = "OPoint<T, D>: Serialize",
93 deserialize = "OPoint<T, D>: Deserialize<'de>"
94))]
95pub struct AxisAlignedBoundingBox<T, D>
96where
97 T: Scalar,
98 D: DimName,
99 DefaultAllocator: Allocator<T, D>,
100{
101 min: OPoint<T, D>,
102 max: OPoint<T, D>,
103}
104
105impl<T, D> Copy for AxisAlignedBoundingBox<T, D>
106where
107 T: Scalar,
108 D: DimName,
109 DefaultAllocator: Allocator<T, D>,
110 OPoint<T, D>: Copy,
111{
112}
113
114pub type AxisAlignedBoundingBox2d<T> = AxisAlignedBoundingBox<T, U2>;
115pub type AxisAlignedBoundingBox3d<T> = AxisAlignedBoundingBox<T, U3>;
116
117impl<T, D> AxisAlignedBoundingBox<T, D>
118where
119 T: Scalar + PartialOrd,
120 D: DimName,
121 DefaultAllocator: Allocator<T, D>,
122{
123 pub fn new(min: OPoint<T, D>, max: OPoint<T, D>) -> Self {
124 for i in 0..D::dim() {
125 assert!(min[i] <= max[i]);
126 }
127 Self { min, max }
128 }
129
130 pub fn min(&self) -> &OPoint<T, D> {
131 &self.min
132 }
133
134 pub fn max(&self) -> &OPoint<T, D> {
135 &self.max
136 }
137}
138
139impl<T, D> From<OPoint<T, D>> for AxisAlignedBoundingBox<T, D>
140where
141 T: Scalar + PartialOrd,
142 D: DimName,
143 DefaultAllocator: Allocator<T, D>,
144{
145 fn from(point: OPoint<T, D>) -> Self {
146 AxisAlignedBoundingBox::new(point.clone(), point)
147 }
148}
149
150impl<T, D> AxisAlignedBoundingBox<T, D>
151where
152 T: Real,
153 D: DimName,
154 DefaultAllocator: Allocator<T, D>,
155{
156 pub fn enclose(&self, other: &AxisAlignedBoundingBox<T, D>) -> Self {
158 let min = self
159 .min
160 .iter()
161 .zip(&other.min.coords)
162 .map(|(a, b)| T::min(*a, *b));
163 let min = OVector::<T, D>::from_iterator(min);
164
165 let max = self
166 .max
167 .iter()
168 .zip(&other.max.coords)
169 .map(|(a, b)| T::max(*a, *b));
170 let max = OVector::<T, D>::from_iterator(max);
171
172 AxisAlignedBoundingBox::new(min.into(), max.into())
173 }
174
175 pub fn from_points<'a>(points: impl IntoIterator<Item = &'a OPoint<T, D>>) -> Option<Self> {
176 let mut points = points.into_iter();
177 points.next().map(|first_point| {
178 points.fold(AxisAlignedBoundingBox::from(first_point.clone()), |aabb, point| {
179 aabb.enclose(&AxisAlignedBoundingBox::from(point.clone()))
180 })
181 })
182 }
183
184 pub fn extents(&self) -> OVector<T, D> {
185 self.max() - self.min()
186 }
187
188 pub fn max_extent(&self) -> T {
189 (self.max() - self.min()).amax()
190 }
191
192 pub fn center(&self) -> OPoint<T, D> {
193 OPoint::from((&self.max().coords + &self.min().coords) / T::from_f64(2.0).unwrap())
194 }
195
196 #[replace_float_literals(T::from_f64(literal).unwrap())]
211 pub fn uniformly_scale(&self, scale: T) -> Self {
212 assert!(scale >= T::zero());
213 let s = scale;
214 let (a, b) = (&self.min, &self.max);
215 let ref c = self.center();
216 Self {
217 min: c + (a - c) * s,
218 max: c + (b - c) * s,
219 }
220 }
221
222 pub fn contains_point(&self, point: &OPoint<T, D>) -> bool {
223 (0..D::dim()).all(|dim| point[dim] >= self.min[dim] && point[dim] <= self.max[dim])
224 }
225
226 pub fn intersects(&self, other: &Self) -> bool {
227 for i in 0..D::dim() {
228 if !intervals_intersect([self.min[i], self.max[i]], [other.min[i], other.max[i]]) {
229 return false;
230 }
231 }
232 true
233 }
234
235 pub fn grow_uniformly(&self, distance: T) -> Self {
249 let min = self.min().map(|b_i| b_i - distance);
250 let max = self.max().map(|b_i| b_i + distance);
251 Self::new(min, max)
252 }
253
254 pub fn corners_iter<'a>(&'a self) -> impl 'a + Iterator<Item = OPoint<T, D>>
256 where
257 DefaultAllocator: Allocator<usize, D>,
258 {
259 index_set_nth_power_iter::<D>(2).map(move |multi_idx| {
262 OVector::<T, D>::from_fn(|idx, _| match multi_idx[idx] {
263 0 => self.min[idx].clone(),
264 1 => self.max[idx].clone(),
265 _ => unreachable!(),
266 })
267 .into()
268 })
269 }
270
271 pub fn closest_point_to(&self, point: &OPoint<T, D>) -> OPoint<T, D> {
273 point
274 .coords
275 .zip_zip_map(&self.min.coords, &self.max.coords, |p_i, a_i, b_i| {
276 if p_i <= a_i {
277 a_i
278 } else if p_i >= b_i {
279 b_i
280 } else {
281 p_i
282 }
283 })
284 .into()
285 }
286
287 pub fn dist_to(&self, point: &OPoint<T, D>) -> T {
289 self.dist2_to(point).sqrt()
290 }
291
292 pub fn dist2_to(&self, point: &OPoint<T, D>) -> T {
294 (self.closest_point_to(point) - point).norm_squared()
295 }
296
297 #[replace_float_literals(T::from_f64(literal).unwrap())]
299 pub fn furthest_point_to(&self, point: &OPoint<T, D>) -> OPoint<T, D>
300 where
301 DefaultAllocator: Allocator<usize, D>,
302 {
303 point
306 .coords
307 .zip_zip_map(&self.min.coords, &self.max.coords, |p_i, a_i, b_i| {
308 let mid = (a_i + b_i) / 2.0;
309 if p_i < mid {
310 b_i
311 } else {
312 a_i
313 }
314 })
315 .into()
316 }
317
318 pub fn max_dist2_to(&self, point: &OPoint<T, D>) -> T
324 where
325 DefaultAllocator: Allocator<usize, D>,
327 {
328 (self.furthest_point_to(point) - point).norm_squared()
329 }
330
331 pub fn max_dist_to(&self, point: &OPoint<T, D>) -> T
337 where
338 DefaultAllocator: Allocator<usize, D>,
340 {
341 self.max_dist2_to(point).sqrt()
342 }
343}
344
345fn intervals_intersect<T: Real>([l1, u1]: [T; 2], [l2, u2]: [T; 2]) -> bool {
346 l2 <= u1 && u2 >= l1
347}
348
349impl<T, D> BoundedGeometry<T> for AxisAlignedBoundingBox<T, D>
350where
351 T: Real,
352 D: DimName,
353 DefaultAllocator: Allocator<T, D>,
354{
355 type Dimension = D;
356
357 fn bounding_box(&self) -> AxisAlignedBoundingBox<T, D> {
358 self.clone()
359 }
360}
361
362#[derive(Copy, Debug, Clone, PartialEq, Eq)]
363pub enum Orientation {
364 Clockwise,
365 Counterclockwise,
366}
367
368#[derive(Debug, Copy, Clone, PartialEq, Eq, Hash)]
369pub enum OrientationTestResult {
370 Positive,
371 Zero,
372 Negative,
373}
374
375#[derive(Debug, Clone, PartialEq)]
376pub struct PolygonClosestPoint<T, D>
377where
378 T: Scalar,
379 D: DimName,
380 DefaultAllocator: Allocator<T, D>,
381{
382 pub closest_point: OPoint<T, D>,
383 pub distance: T,
384}
385
386#[derive(Debug, Copy, Clone, PartialEq)]
387pub struct PolygonPointProjection3d<T: Scalar> {
388 pub projected_point: Point3<T>,
390 pub distance: T,
392 pub signed_plane_distance: T,
396}
397
398pub trait ConvexPolygon3d<'a, T: Scalar>: Debug {
402 fn num_vertices(&self) -> usize;
403 fn get_vertex(&self, index: usize) -> Option<Point3<T>>;
404
405 fn compute_plane(&self) -> Option<Plane<T>>
406 where
407 T: Real,
408 {
409 let normal = self.compute_normal();
410 let point = self.get_vertex(0)?;
411 Some(Plane::from_point_and_normal(point, Unit::new_normalize(normal)))
412 }
413
414 fn compute_half_space(&self) -> Option<HalfSpace<T>>
415 where
416 T: Real,
417 {
418 let normal = self.compute_normal();
419 let point = self.get_vertex(0)?;
420 Some(HalfSpace::from_point_and_normal(point, Unit::new_unchecked(-normal)))
421 }
422
423 #[replace_float_literals(T::from_f64(literal).unwrap())]
429 fn compute_area_vector(&self) -> Vector3<T>
430 where
431 T: Real,
432 {
433 assert!(self.num_vertices() >= 3, "Polygons must have at least 3 vertices.");
434
435 let mut area_vector = Vector3::zeros();
436
437 let a = self.get_vertex(0).unwrap();
438 for i in 1..self.num_vertices() - 1 {
439 let b = self.get_vertex(i).unwrap();
440 let c = self.get_vertex(i + 1).unwrap();
441 let ab = &b.coords - &a.coords;
442 let ac = &c.coords - &a.coords;
443 area_vector += ab.cross(&ac) * 0.5;
444 }
445
446 area_vector
447 }
448
449 fn compute_normal(&self) -> Vector3<T>
451 where
452 T: Real,
453 {
454 self.compute_area_vector().normalize()
465 }
466
467 fn closest_point(&self, point: &Point3<T>) -> PolygonClosestPoint<T, U3>
468 where
469 T: Real,
470 {
471 assert!(self.num_vertices() >= 3, "Polygon must have at least 3 vertices.");
472
473 let n = self.compute_normal();
477
478 let mut inside = true;
479
480 for i in 0..self.num_vertices() {
481 let v1 = self.get_vertex(i).unwrap();
482 let v2 = self.get_vertex((i + 1) % self.num_vertices()).unwrap();
483
484 let e = &v2.coords - &v1.coords;
486
487 let half_space_normal = Unit::new_normalize(e.cross(&n));
492 let half_space = HalfSpace::from_point_and_normal(v1, half_space_normal);
493
494 if !half_space.contains_point(point) {
495 inside = false;
496 break;
497 }
498 }
499
500 if inside {
501 let x0 = self.get_vertex(0).unwrap();
506 let signed_plane_distance = n.dot(&(point - x0));
508 let closest_point = point - &n * signed_plane_distance;
509 PolygonClosestPoint {
510 closest_point,
511 distance: signed_plane_distance.abs(),
512 }
513 } else {
514 let mut closest_dist2 = T::max_value().unwrap();
518 let mut closest_point = Point3::origin();
519
520 for i in 0..self.num_vertices() {
521 let v1 = self.get_vertex(i).unwrap();
522 let v2 = self.get_vertex((i + 1) % self.num_vertices()).unwrap();
523 let segment = LineSegment3d::from_end_points(v1, v2);
524 let projected = segment.closest_point(point);
525 let dist2 = distance_squared(&projected, point);
526
527 if dist2 < closest_dist2 {
528 closest_dist2 = dist2;
529 closest_point = projected;
530 }
531 }
532
533 PolygonClosestPoint {
534 closest_point,
535 distance: closest_dist2.sqrt(),
536 }
537 }
538 }
539}
540
541pub trait ConvexPolyhedron<'a, T: Scalar>: Debug {
542 type Face: ConvexPolygon3d<'a, T>;
543
544 fn num_faces(&self) -> usize;
545 fn get_face(&self, index: usize) -> Option<Self::Face>;
546
547 #[replace_float_literals(T::from_f64(literal).unwrap())]
548 fn compute_signed_distance(&self, point: &Point3<T>) -> SignedDistanceResult<T, U3>
549 where
550 T: Real,
551 {
552 assert!(self.num_faces() >= 4, "Polyhedron must have at least 4 faces.");
553 let mut inside = true;
554 let mut closest_dist = T::max_value().unwrap();
555 let mut closest_point = Point3::origin();
556 let mut closest_face_index = 0;
557
558 for i in 0..self.num_faces() {
559 let face = self.get_face(i).unwrap();
560 let closest_point_result = face.closest_point(point);
561
562 if closest_point_result.distance < closest_dist {
563 closest_dist = closest_point_result.distance;
564 closest_face_index = i;
565 closest_point = closest_point_result.closest_point;
566 }
567
568 let n = face.compute_normal();
569 let x0 = closest_point_result.closest_point;
570 let half_space = HalfSpace::from_point_and_normal(x0, Unit::new_unchecked(-n));
573 if !half_space.contains_point(&point) {
574 inside = false;
575 }
576 }
577
578 let sign = if inside { -1.0 } else { 1.0 };
579
580 debug_assert!(closest_dist >= 0.0);
581 SignedDistanceResult {
582 feature_id: closest_face_index,
583 point: closest_point,
584 signed_distance: sign * closest_dist,
585 }
586 }
587
588 fn compute_volume(&'a self) -> T
589 where
590 T: Real,
591 {
592 let faces = (0..self.num_faces()).map(move |idx| {
593 self.get_face(idx)
594 .expect("Number of faces reported must be correct.")
595 });
596 compute_polyhedron_volume_from_faces(faces)
597 }
598
599 fn contains_point(&'a self, point: &Point3<T>) -> bool
603 where
604 T: Real,
605 {
606 for i in 0..self.num_faces() {
609 let face = self.get_face(i).unwrap();
610 let n = face.compute_normal();
611 let x0 = face
612 .get_vertex(0)
613 .expect("TODO: How to handle empty polygon?");
614 let half_space = HalfSpace::from_point_and_normal(x0, Unit::new_unchecked(-n));
616 if !half_space.contains_point(point) {
617 return false;
618 }
619 }
620
621 true
622 }
623}
624
625#[replace_float_literals(T::from_f64(literal).unwrap())]
626pub fn compute_polyhedron_volume_from_faces<'a, T, F>(boundary_faces: impl 'a + IntoIterator<Item = F>) -> T
627where
628 T: Real,
629 F: ConvexPolygon3d<'a, T>,
630{
631 let mut volume = T::zero();
635 for face in boundary_faces {
636 if face.num_vertices() > 2 {
640 let x0 = face.get_vertex(0).unwrap();
641 let area_vector = face.compute_area_vector();
643 let area = area_vector.magnitude();
644 if area > T::zero() {
645 let normal = area_vector.normalize();
646 volume += (normal.dot(&x0.coords)) * area;
647 }
648 }
649 }
650
651 volume = volume.abs() / 3.0;
654
655 volume
656}
657
658impl<T> From<Triangle2d<T>> for ConvexPolygon<T>
701where
702 T: Scalar,
703{
704 fn from(triangle: Triangle2d<T>) -> Self {
705 ConvexPolygon::from_vertices(triangle.0.iter().map(|v| Point2::from(v.clone())).collect())
707 }
708}
709
710impl<T> TryFrom<Quad2d<T>> for ConvexPolygon<T>
711where
712 T: Real,
713{
714 type Error = ConcavePolygonError;
715
716 fn try_from(value: Quad2d<T>) -> Result<Self, Self::Error> {
717 if value.concave_corner().is_none() {
718 Ok(ConvexPolygon::from_vertices(
720 value.0.iter().map(|v| Point2::from(v.clone())).collect(),
721 ))
722 } else {
723 Err(ConcavePolygonError)
724 }
725 }
726}
727
728#[derive(Debug, Clone, Copy, PartialEq, Eq)]
729#[repr(transparent)]
730pub struct Convex<T>(pub T);
731
732impl<T> Convex<T> {
733 pub fn assume_convex(obj: T) -> Self {
734 Self(obj)
735 }
736}