use spade

Cargo.toml

@@ -47,6 +47,7 @@ bevy-inspector-egui = { version = "0.23.2", optional = true }
 itertools = "0.12.1"
 meshopt = { version = "0.2.0", optional = true }
 rand = "0.8.5"
+spade = "2.6.0"
 
 [dev-dependencies]
 criterion = { version = "0.4", features = ["html_reports"] }

README.md

@@ -66,9 +66,8 @@ For package development, we recommend using the `editor`-subcrate. This example
   - [ ] UV Coordinates
   - [ ] Custom Attributes
 - [ ] Triangulation
-  - [ ] Montone Sweep-Line
-  - [ ] Min Weight
-  - [ ] Constrained Delaunay
+  - [(x)] Montone Sweep-Line
+  - [x] Constrained Delaunay (using Delaunator)
 - [ ] Primitives
   - [x] Polygon, Star
   - [x] Cuboid

editor/src/main.rs

@@ -179,7 +179,7 @@ fn make_2d_shape(_settings: &MeshSettings) -> MeshVec3 {
     p.reverse();
     let mut mesh = MeshVec3::polygon(&p);
     */
-    let mut mesh = MeshVec3::regular_star(2.0, 1.0, 100);
+    let mut mesh = MeshVec3::regular_star(2.0, 1.0, 1000);
     mesh.transform(&Transform::from_translation(Vec3::new(0.0, -0.99, 0.0)));
     mesh
 }

src/math/impls/f32/mod.rs

@@ -21,4 +21,9 @@ impl Scalar for f32 {
     fn acos(self) -> Self {
         f32::acos(self)
     }
+
+    #[inline(always)]
+    fn to_f64(self) -> f64 {
+        self as f64
+    }
 }

src/math/impls/f64/mod.rs

@@ -21,4 +21,9 @@ impl Scalar for f64 {
     fn acos(self) -> Self {
         f64::acos(self)
     }
+
+    #[inline(always)]
+    fn to_f64(self) -> f64 {
+        self
+    }
 }

src/math/scalar.rs

@@ -31,6 +31,9 @@ pub trait Scalar:
     /// Returns whether the scalar is strictly negative.
     fn is_negative(self) -> bool;
 
+    /// Converts the scalar to a 64-bit floating point number.
+    fn to_f64(self) -> f64;
+
     /// Returns the absolute value of the scalar.
     fn abs(self) -> Self {
         if self.is_positive() {

src/math/vector2d.rs

@@ -51,7 +51,7 @@ pub trait Vector2D<S: Scalar>: Vector<S> {
     }
 
     /// Whether the point is inside the circumcircle of the triangle.
-    #[inline(always)]
+    #[deprecated(since="0.1.0", note="use something more robust")]
     fn is_inside_circumcircle(&self, a: Self, b: Self, c: Self, eps: S) -> bool {
         // https://en.wikipedia.org/wiki/Delaunay_triangulation#Algorithms
 

src/representation/face/mod.rs

@@ -74,6 +74,30 @@ impl<E: IndexType, F: IndexType> Face<E, F> {
     pub fn num_triangles<V: IndexType, P: Payload>(&self, mesh: &Mesh<E, V, F, P>) -> usize {
         (self.num_vertices(mesh) - 2) * 3
     }
+
+    /// Whether a triangle is inside or outside of the face
+    pub fn is_inside<V: IndexType, P: Payload>(
+        &self,
+        mesh: &Mesh<E, V, F, P>,
+        v0: V,
+        v1: V,
+        v2: V,
+    ) -> bool {
+        if let Some(e) = mesh.edge_between(v0, v1) {
+            return !e.is_boundary_self();
+        }
+        if let Some(e) = mesh.edge_between(v1, v2) {
+            return !e.is_boundary_self();
+        }
+        if let Some(e) = mesh.edge_between(v2, v0) {
+            return !e.is_boundary_self();
+        }
+
+        //println!("v0: {} {} {}", v0.index(), v1.index(), v2.index());
+
+        return true;
+        //todo!("is_inside: not implemented for faces not touching the boundary");
+    }
 }
 
 impl<E: IndexType, F: IndexType> std::fmt::Display for Face<E, F> {

src/representation/face/tesselate/delaunay/dual.rs

@@ -1,3 +1,6 @@
+//! This whole module is deprecated
+
+
 use crate::{math::IndexType, representation::Face};
 use itertools::Itertools;
 

src/representation/face/tesselate/delaunay/mod.rs

@@ -1,10 +1,18 @@
 use super::{Face, Mesh, Payload};
 use crate::{
-    math::{Scalar, Vector2D, Vector3D},
+    math::{Scalar, Vector, Vector2D, Vector3D},
     representation::IndexType,
 };
-use std::collections::{HashMap, VecDeque};
+use std::collections::VecDeque;
 mod dual;
+use spade::{
+    handles::{
+        FixedDirectedEdgeHandle,
+        VoronoiVertex::{self, Inner, Outer},
+    },
+    AngleLimit, ConstrainedDelaunayTriangulation, FloatTriangulation as _, HasPosition,
+    InsertionError, Point2, RefinementParameters, Triangulation as _,
+};
 
 impl<E, F> Face<E, F>
 where
@@ -14,6 +22,7 @@ where
     /// Flips edges until the delaunay-condition is met.
     /// This is quite slow in the worst case O(n^3) but usually much better than the naive version.
     /// Assumes local indices
+    #[deprecated(since = "0.1.0", note = "please use `delaunator` instead")]
     pub fn delaunayfy<V: IndexType, P: Payload>(
         &self,
         mesh: &Mesh<E, V, F, P>,
@@ -90,78 +99,71 @@ where
         }
     }
 
-    /// Flips edges until the delaunay-condition is met. This is quite slow: O(n^3).
-    pub fn delaunayfy_naive<V: IndexType, P: Payload>(
+    /// Converts the face into a triangle list using the delaunay triangulation.
+    pub fn delaunay_triangulation<V: IndexType, P: Payload>(
         &self,
         mesh: &Mesh<E, V, F, P>,
         indices: &mut Vec<V>,
-        local_indices: bool,
     ) where
         P::Vec: Vector3D<P::S>,
     {
-        let eps = P::S::EPS * P::S::from(10.0); // TODO
-        let vs: Vec<(P::Vec2, V)> = self.vertices_2d::<V, P>(mesh).collect();
-        assert!(vs.len() == self.num_vertices(mesh));
-        let mut vsh: HashMap<V, P::Vec2> = HashMap::new();
-        if local_indices {
-            for (i, (v, p)) in vs.iter().enumerate() {
-                vsh.insert(V::new(i), *v);
-            }
-        } else {
-            for (v, p) in vs {
-                vsh.insert(p, v);
+        debug_assert!(self.may_be_curved() || self.is_planar2(mesh));
+        debug_assert!(!self.has_self_intersections(mesh));
+
+        //let n = indices.len();
+        //self.sweep_line(mesh, indices);
+        //self.delaunayfy(mesh, indices, n);
+
+        // using Delaunator because it is faster than spade for < 100_000 vertices
+        // https://github.com/Stoeoef/spade/tree/master/delaunay_compare
+        /*let vs: Vec<delaunator::Point> = self
+            .vertices_2d::<V, P>(mesh)
+            .map(|(vec2, _)| delaunator::Point {
+                x: vec2.x().to_f64(),
+                y: vec2.y().to_f64(),
+            })
+            .collect();
+        let triangulation = delaunator::triangulate(&vs);
+        indices.extend(triangulation.triangles.iter().map(|i| V::new(*i)));*/
+
+        let mut cdt = ConstrainedDelaunayTriangulation::<Point2<_>>::new();
+        //TODO: faster: ConstrainedDelaunayTriangulation::bulk_load()
+        let mut last = None;
+        let mut first = None;
+        for (i2, (vec2, _)) in self.vertices_2d::<V, P>(mesh).enumerate() {
+            let i = cdt
+                .insert(Point2::new(vec2.x().to_f64(), vec2.y().to_f64()))
+                .unwrap();
+            assert!(i.index() == i2);
+            if let Some(j) = last {
+                assert!(cdt.add_constraint(j, i));
+            } else {
+                first = Some(i);
             }
+            last = Some(i);
         }
-
-        if indices.len() < 3 {
-            return;
-        }
-
-        for _ in 0..indices.len() {
-            let mut changed = false;
-            for i in (0..indices.len()).step_by(3) {
-                for j in ((i + 3)..indices.len()).step_by(3) {
-                    for k in 0..3 {
-                        let a = indices[i + (0 + k) % 3];
-                        let b = indices[i + (1 + k) % 3];
-                        let c = indices[i + (2 + k) % 3];
-                        for l in 0..3 {
-                            let d = indices[j + (0 + l) % 3];
-                            let e = indices[j + (1 + l) % 3];
-                            let f = indices[j + (2 + l) % 3];
-                            if a == e && b == d {
-                                if vsh[&f].is_inside_circumcircle(vsh[&a], vsh[&b], vsh[&c], eps) {
-                                    //if vsh[&a].distance(&vsh[&b]) > vsh[&c].distance(&vsh[&f]) {
-                                    indices[i + (0 + k) % 3] = c;
-                                    indices[i + (1 + k) % 3] = f;
-                                    indices[i + (2 + k) % 3] = b;
-
-                                    indices[j + (0 + l) % 3] = c;
-                                    indices[j + (1 + l) % 3] = a;
-                                    indices[j + (2 + l) % 3] = f;
-
-                                    changed = true;
-
-                                    break;
-                                }
-                            }
-                        }
-                    }
-                }
-            }
-            if !changed {
-                break;
+        assert!(cdt.add_constraint(last.unwrap(), first.unwrap()));
+
+        let i2v = self
+            .vertices_2d::<V, P>(mesh)
+            .map(|(_, i)| i)
+            .collect::<Vec<_>>();
+        let i0 = indices.len();
+        cdt.inner_faces().for_each(|face| {
+            let [p0, p1, p2] = face.vertices();
+            let v0 = V::new(p0.index());
+            let v1 = V::new(p1.index());
+            let v2 = V::new(p2.index());
+
+            //if self.is_inside(mesh, i2v[p0.index()], i2v[p1.index()], i2v[p2.index()]) {
+            if self.is_inside(
+                mesh,
+                V::new(i0 + p0.index()),
+                V::new(i0 + p1.index()),
+                V::new(i0 + p2.index()),
+            ) {
+                indices.extend(&[v0, v1, v2]);
             }
-        }
+        });
     }
-
-    /*/// Converts the face into a triangle list using the delaunay triangulation.
-    pub fn delaunay_triangulation<V: IndexType, P: Payload>(
-        &self,
-        mesh: &Mesh<E, V, F, P>,
-        _indices: &mut Vec<V>,
-    ) {
-        assert!(self.may_be_curved() || self.is_planar2(mesh));
-        // TODO: or at least some other O(n log n) algorithm: https://en.wikipedia.org/wiki/Delaunay_triangulation
-    }*/
 }

src/representation/face/tesselate/mod.rs

@@ -125,13 +125,7 @@ where
             }
             TriangulationAlgorithm::Delaunay => {
                 self.expand_local_indices(mesh, indices, local_indices, |mesh, indices| {
-                    let n = indices.len();
-                    self.sweep_line(mesh, indices);
-                    self.delaunayfy(mesh, indices, n);
-
-                    /*let mut i = indices[n..].to_vec();
-                    self.delaunayfy_naive(mesh, &mut i, true);
-                    indices[n..].copy_from_slice(&i);*/
+                    self.delaunay_triangulation(mesh, indices);
                 });
             }
         }