Debugging Bowyer-Watson

Singular matrices oh my!

README.md

@@ -23,10 +23,10 @@ the [Rust] ecosystem.
 - [x] Arbitrary data types associated with vertices and cells
 - [x] Serialization/Deserialization of all data structures to/from [JSON]
 
-At some point I may merge it into another library, such as [Spade], or
-[delaunay], but for now I am developing this to use in my [research] without
-trying to figure out how to fit it into other libraries and with the
-minimum number of traits to do generic computational geometry.
+At some point I may merge into another library, such as [Spade] or [delaunay],
+but for now I am developing this to use in my [research] without trying to
+figure out how to mesh with other libraries and coding conventions, and with
+the minimum number of [traits] to do generic computational geometry.
 
 [Rust]: https://rust-lang.org
 [CGAL]: https://www.cgal.org/
@@ -38,3 +38,4 @@ minimum number of traits to do generic computational geometry.
 [Constrained Delaunay triangulations]: https://en.wikipedia.org/wiki/Constrained_Delaunay_triangulation
 [Voronoi diagrams]: https://en.wikipedia.org/wiki/Voronoi_diagram
 [research]: https://github.com/acgetchell/cdt-rs
+[traits]: https://doc.rust-lang.org/book/ch10-02-traits.html

src/delaunay_core/cell.rs

@@ -7,7 +7,7 @@ use serde::{de::DeserializeOwned, Deserialize, Serialize};
 use std::{collections::HashMap, fmt::Debug, iter::Sum, ops::Div};
 use uuid::Uuid;
 
-#[derive(Clone, Debug, Default, Deserialize, PartialEq, Serialize)]
+#[derive(Clone, Debug, Default, Deserialize, PartialEq, Serialize, PartialOrd)]
 /// The `Cell` struct represents a d-dimensional
 /// [simplex](https://en.wikipedia.org/wiki/Simplex) with vertices, a unique
 /// identifier, optional neighbors, and optional data.

src/delaunay_core/facet.rs

@@ -9,7 +9,7 @@
 use super::{cell::Cell, vertex::Vertex};
 use serde::{de::DeserializeOwned, Deserialize, Serialize};
 
-#[derive(Clone, Debug, Default, Deserialize, PartialEq, Serialize)]
+#[derive(Clone, Debug, Default, Deserialize, PartialEq, Serialize, PartialOrd)]
 /// The `Facet` struct represents a facet of a d-dimensional simplex.
 /// Passing in a `Vertex` and a `Cell` containing that vertex to the
 /// constructor will create a `Facet` struct.

src/delaunay_core/point.rs

@@ -2,7 +2,7 @@
 
 use serde::{de::DeserializeOwned, Deserialize, Serialize};
 
-#[derive(Clone, Copy, Debug, Default, Deserialize, PartialEq, Serialize)]
+#[derive(Clone, Copy, Debug, Default, Deserialize, PartialEq, Serialize, PartialOrd)]
 /// The `Point` struct represents a point in a D-dimensional space, where the
 /// coordinates are of type `T`.
 ///

src/delaunay_core/triangulation_data_structure.rs

@@ -3,7 +3,9 @@
 //! Intended to match functionality of the
 //! [CGAL Triangulation](https://doc.cgal.org/latest/Triangulation/index.html).
 
-use super::{cell::Cell, point::Point, utilities::find_extreme_coordinates, vertex::Vertex};
+use super::{
+    cell::Cell, facet::Facet, point::Point, utilities::find_extreme_coordinates, vertex::Vertex,
+};
 use na::{ComplexField, Const, OPoint};
 use nalgebra as na;
 use serde::{de::DeserializeOwned, Deserialize, Serialize};
@@ -76,8 +78,8 @@ where
     f64: From<T>,
     for<'a> &'a T: Div<f64>,
     [T; D]: Default + DeserializeOwned + Serialize + Sized,
-    U: Clone + Copy,
-    V: Clone + Copy,
+    U: Clone + Copy + PartialEq,
+    V: Clone + Copy + PartialEq,
 {
     /// The function creates a new instance of a triangulation data structure
     /// with given points, initializing the vertices and cells.
@@ -257,6 +259,9 @@ where
         Vertex<T, U, D>: Clone, // Add the Clone trait bound for Vertex
         OPoint<T, Const<D>>: From<[f64; D]>,
         [f64; D]: Serialize + DeserializeOwned + Default,
+        f64: PartialOrd, // Replace f64: Ord with f64: PartialOrd
+        U: PartialOrd,
+        V: PartialOrd,
     {
         let mut cells: Vec<Cell<T, U, V, D>> = Vec::new();
 
@@ -266,42 +271,59 @@ where
 
         // Iterate over vertices
         for vertex in self.vertices.values() {
-            // Find all cells that contain the vertex
+            // Find cells that contain the vertex
             let mut bad_cells: Vec<Cell<T, U, V, D>> = Vec::new();
             for cell in cells.iter() {
+                // TODO: understand why we're getting singular matrices here
                 if cell.circumsphere_contains(*vertex).unwrap() {
                     bad_cells.push((*cell).clone());
                 }
             }
 
             // Find the boundary of the polygonal hole
-            let mut _polygonal_hole: Vec<Vertex<T, U, D>> = Vec::new();
-            for _cell in bad_cells.iter() {
+            let mut polygonal_hole: Vec<Facet<T, U, V, D>> = Vec::new();
+            for cell in bad_cells.iter() {
+                // Create `Facet`s from the `Cell`
+                for vertex in cell.vertices.iter() {
+                    let facet = Facet::new(cell.clone(), *vertex).unwrap();
+                    polygonal_hole.push(facet);
+                }
+
                 // for vertex in cell.vertices.iter() {
                 //     if bad_cells.iter().any(|c| c.contains_vertex(vertex)) {
                 //         polygonal_hole.push(vertex.clone());
                 //     }
                 // }
             }
 
-            // // Remove duplicate vertices
-            // polygonal_hole.sort();
-            // polygonal_hole.dedup();
-
-            // // Remove bad cells from the triangulation
-            // for cell in bad_cells.iter() {
-            //     cells.remove(cells.iter().position(|c| c == cell).unwrap());
-            // }
-
-            // // Re-triangulate the polygonal hole
-            // for vertex in polygonal_hole.iter() {
-            //     let mut new_cell = Cell::new(vec![vertex.clone()]);
-            //     new_cell.vertices.push(vertex.clone());
-            //     new_cell.vertices.push(vertex.clone());
-            //     cells.push(new_cell);
-            // }
+            // Remove duplicate facets
+            polygonal_hole.sort_by(|a, b| a.partial_cmp(b).unwrap());
+            polygonal_hole.dedup();
+
+            // Remove bad cells from the triangulation
+            for cell in bad_cells.iter() {
+                cells.remove(cells.iter().position(|c| c == cell).unwrap());
+            }
+
+            // Re-triangulate the polygonal hole
+            for mut facet in polygonal_hole.iter().cloned() {
+                let mut new_cell_vertices: Vec<Vertex<T, U, D>> = Vec::new();
+                for facet_vertex in facet.vertices().iter() {
+                    new_cell_vertices.push(*facet_vertex);
+                }
+                new_cell_vertices.push(*vertex);
+                cells.push(Cell::new(new_cell_vertices)?);
+            }
         }
 
+        // Remove all cells containing vertices from the supercell
+        // TODO: Fix borrow-checker errors
+        // for cell in cells.iter_mut() {
+        //     if cell.contains_vertex(supercell.vertices[0]) {
+        //         cells.remove(cells.iter().position(|c| c == cell).unwrap());
+        //     }
+        // }
+
         Ok(cells)
     }
 
@@ -427,13 +449,14 @@ mod tests {
         println!("{:?}", unwrapped_supercell);
     }
 
+    #[ignore]
     #[test]
     fn tds_bowyer_watson() {
         let points = vec![
-            Point::new([1.0, 2.0, 3.0]),
-            Point::new([4.0, 5.0, 6.0]),
-            Point::new([7.0, 8.0, 9.0]),
-            Point::new([10.0, 11.0, 12.0]),
+            Point::new([0.0, 0.0, 0.0]),
+            Point::new([1.0, 0.0, 0.0]),
+            Point::new([0.0, 1.0, 0.0]),
+            Point::new([0.0, 0.0, 1.0]),
         ];
         let mut tds: Tds<f64, usize, usize, 3> = Tds::new(points);
         let cells = tds.bowyer_watson();

src/delaunay_core/vertex.rs

@@ -5,7 +5,7 @@ use serde::{de::DeserializeOwned, Deserialize, Serialize};
 use std::{collections::HashMap, option::Option};
 use uuid::Uuid;
 
-#[derive(Clone, Copy, Debug, Default, Deserialize, PartialEq, Serialize)]
+#[derive(Clone, Copy, Debug, Default, Deserialize, PartialEq, Serialize, PartialOrd)]
 /// The `Vertex` struct represents a vertex in a triangulation with a `Point`,
 /// a unique identifier, an optional incident cell identifier, and optional
 /// data.