delaunayfy

editor/src/main.rs

@@ -54,7 +54,7 @@ impl Default for MeshSettings {
     fn default() -> Self {
         MeshSettings {
             tol: -4.0,
-            n: 7,
+            n: 30,
             r: 1.0,
             d1: Vec3::new(0.4, -1.0, 0.0),
             d2: Vec3::new(-1.0, 0.3, -1.0),
@@ -66,9 +66,9 @@ impl Default for MeshSettings {
 fn make_mesh(settings: &MeshSettings) -> MeshVec3 {
     let mut mesh = MeshVec3::regular_polygon(settings.r, settings.n); //cuboid(1.0, 1.0, 2.0);
     mesh.extrude(mesh.edge_between(1, 0).unwrap().id(), settings.d1, true);
-    let fe = mesh.extrude_face(3, settings.d2, true);
+    /*let fe = mesh.extrude_face(3, settings.d2, true);
     mesh.extrude_face(fe, settings.d3, true);
-    println!("{}", mesh);
+    println!("{}", mesh);*/
     mesh
 }
 

src/representation/deletable.rs

@@ -60,6 +60,11 @@ impl<T: Deletable<I> + Default, I: IndexType> DeletableVector<T, I> {
         self.data.len() - self.deleted.len()
     }
 
+    /// Returns the maximum index of the non-deleted elements.
+    pub fn max_ind(&self) -> usize {
+        self.data.len()
+    }
+
     /// Allocates a new element, moves the given to that index, sets the new id, and returns the index.
     pub fn push(&mut self, mut v: T) -> I {
         assert!(

src/representation/face/tesselate/mod.rs

@@ -1,3 +1,5 @@
+use std::collections::HashMap;
+
 use itertools::Itertools;
 
 use super::{Face, Mesh, Payload, Scalar};
@@ -39,6 +41,7 @@ where
     ) where
         P::Vec: Vector3D<P::S>,
     {
+        // TODO: ear clipping is inefficient
         assert!(self.is_planar(mesh, P::S::EPS * 10.0.into()));
 
         let vs: Vec<(<P::Vec as Vector<P::S>>::Vec2D, V)> =
@@ -98,6 +101,62 @@ where
         }
     }
 
+    /// Flip edges to make the face delaunay.
+    pub fn delaunayfy<V: IndexType, P: Payload>(&self, mesh: &Mesh<E, V, F, P>, indices: &mut Vec<V>)
+    where
+        P::Vec: Vector3D<P::S>,
+    {
+        let vs: Vec<(<P::Vec as Vector<P::S>>::Vec2D, V)> =
+            self.vertices_2d::<V, P>(mesh).collect();
+        let vsl = vs.len();
+        assert!(vs.len() == self.vertices(mesh).count());
+        let mut vsh: HashMap<V, <P::Vec as Vector<P::S>>::Vec2D> = HashMap::new();
+        for (v, p) in vs {
+            vsh.insert(p, v);
+        }
+
+        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]) {
+                                //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;
+            }
+        }
+    }
+
     /*
     /// Converts the face into a triangle list
     pub fn tesselate<V: IndexType, P: Payload>(
@@ -113,7 +172,15 @@ where
     where
         P::Vec: Vector3D<P::S>,
     {
-        self.ear_clipping(mesh, indices);
+        let mut local_indices = Vec::new();
+        self.ear_clipping(mesh, &mut local_indices);
+        self.delaunayfy(mesh, &mut local_indices);
+        indices.extend(local_indices);
+        assert!(indices.len() % 3 == 0, "{:?}", indices.len());
+        assert!(indices.iter().all(|i| i.index() < mesh.max_vertex_index()));
+
+        // Minimize edge length
+        // TODO: https://en.wikipedia.org/wiki/Minimum-weight_triangulation#Variations
     }
 
     /*/// Converts the face into a triangle list using the delaunay triangulation.

src/representation/mesh/mod.rs

@@ -118,6 +118,11 @@ where
     pub fn num_vertices(&self) -> usize {
         self.vertices.len()
     }
+
+    /// Returns the maximum vertex index in the mesh
+    pub fn max_vertex_index(&self) -> usize {
+        self.vertices.max_ind()
+    }
 
     /// Returns the number of edges in the mesh
     pub fn num_edges(&self) -> usize {

src/representation/vertex/payload/vector.rs

@@ -39,7 +39,19 @@ pub trait Scalar:
     /// Returns the arcus cosine of the scalar.
     fn acos(self) -> Self;
 
-
+    fn det3(
+        a: Self,
+        b: Self,
+        c: Self,
+        d: Self,
+        e: Self,
+        f: Self,
+        g: Self,
+        h: Self,
+        i: Self,
+    ) -> Self {
+        a * (e * i - f * h) - b * (d * i - f * g) + c * (d * h - e * g)
+    }
 }
 
 impl Scalar for f32 {
@@ -197,6 +209,29 @@ pub trait Vector2D<ScalarType: Scalar>: Vector<ScalarType> {
         let inside_cw = bs1.is_negative() && bs2.is_negative() && bs3.is_negative();
         inside_ccw || inside_cw
     }
+
+    /// Whether the point is inside the circumcircle of the triangle.
+    fn is_inside_circumcircle(&self, a: Self, b: Self, c: Self) -> bool {
+        // https://en.wikipedia.org/wiki/Delaunay_triangulation#Algorithms
+
+        let adx = a.x() - self.x();
+        let ady = a.y() - self.y();
+        let bdx = b.x() - self.x();
+        let bdy = b.y() - self.y();
+        let cdx = c.x() - self.x();
+        let cdy = c.y() - self.y();
+        ScalarType::det3(
+            adx,
+            ady,
+            adx * adx + ady * ady,
+            bdx,
+            bdy,
+            bdx * bdx + bdy * bdy,
+            cdx,
+            cdy,
+            cdx * cdx + cdy * cdy,
+        ).is_positive()
+    }
 }
 
 /// Trait for coordinates in 3d space.