Implements CDT bulk loading.

Fixes some edge cases for regular DT bulk loading.

CHANGELOG.md

@@ -4,6 +4,12 @@ All notable changes to this project will be documented in this file.
 The format is based on [Keep a Changelog](http://keepachangelog.com/) 
 and this project adheres to [Semantic Versioning](http://semver.org/).
 
+## [2.7.0] - 2024-MM-DD
+
+### Added
+ - Implements `ConstrainedDelaunayTriangulation::bulk_load_cdt`
+ - Implements `ConstrainedDelaunayTriangulation::bulk_load_cdt_stable`
+
 ## [2.6.0] - 2024-01-14
 
 ### Added

fuzz/Cargo.toml

@@ -31,3 +31,8 @@ path = "fuzz_targets/bulk_load_int_fuzz.rs"
 test = false
 doc = false
 
+[[bin]]
+name = "bulk_load_cdt_fuzz"
+path = "fuzz_targets/bulk_load_cdt_fuzz.rs"
+test = false
+doc = false

fuzz/fuzz_targets/bulk_load_cdt_fuzz.rs

@@ -0,0 +1,53 @@
+#![no_main]
+mod fuzz_shared;
+use fuzz_shared::FuzzPoint;
+use libfuzzer_sys::fuzz_target;
+
+use spade::{ConstrainedDelaunayTriangulation, HasPosition, Triangulation};
+
+fuzz_target!(|input: (Vec<FuzzPoint>, Vec<[usize; 2]>)| {
+    let (data, mut edges) = input;
+    for p in &data {
+        if spade::validate_coordinate(p.x).is_err() || spade::validate_coordinate(p.y).is_err() {
+            return;
+        }
+        if p.x.abs() > 20.0 || p.y.abs() > 20.0 {
+            return;
+        }
+    }
+
+    for &[from, to] in &edges {
+        if from >= data.len() || to >= data.len() || from == to {
+            return;
+        }
+    }
+
+    let mut reference_cdt =
+        ConstrainedDelaunayTriangulation::<FuzzPoint>::bulk_load(data.clone()).unwrap();
+
+    let mut last_index = 0;
+    for (index, [from, to]) in edges.iter().copied().enumerate() {
+        let from = reference_cdt
+            .locate_vertex(data[from].position())
+            .unwrap()
+            .fix();
+        let to = reference_cdt
+            .locate_vertex(data[to].position())
+            .unwrap()
+            .fix();
+
+        if reference_cdt.can_add_constraint(from, to) {
+            reference_cdt.add_constraint(from, to);
+        } else {
+            last_index = index;
+            break;
+        }
+    }
+
+    edges.truncate(last_index);
+
+    let bulk_loaded =
+        ConstrainedDelaunayTriangulation::<FuzzPoint>::bulk_load_cdt(data, edges).unwrap();
+
+    bulk_loaded.cdt_sanity_check();
+});

fuzz/fuzz_targets/bulk_load_fuzz.rs

@@ -1,7 +1,8 @@
 #![no_main]
+mod fuzz_shared;
+use fuzz_shared::FuzzPoint;
 use libfuzzer_sys::fuzz_target;
-
-use spade::{DelaunayTriangulation, Point2, Triangulation, TriangulationExt};
+use spade::{DelaunayTriangulation, Triangulation, TriangulationExt};
 
 fuzz_target!(|data: Vec<FuzzPoint>| {
     for p in &data {
@@ -12,23 +13,6 @@ fuzz_target!(|data: Vec<FuzzPoint>| {
             return;
         }
     }
-    let triangulation = DelaunayTriangulation::<_>::bulk_load(
-        data.iter()
-            .map(|point| Point2::new(point.x as f64, point.y as f64))
-            .collect(),
-    )
-    .unwrap();
+    let triangulation = DelaunayTriangulation::<_>::bulk_load(data.clone()).unwrap();
     triangulation.sanity_check();
 });
-
-#[derive(Clone, Copy, arbitrary::Arbitrary)]
-pub struct FuzzPoint {
-    x: f64,
-    y: f64,
-}
-
-impl core::fmt::Debug for FuzzPoint {
-    fn fmt(&self, f: &mut core::fmt::Formatter<'_>) -> core::fmt::Result {
-        f.write_fmt(format_args!("Point2::new({:?}, {:?})", self.x, self.y))
-    }
-}

fuzz/fuzz_targets/fuzz_shared.rs

@@ -0,0 +1,20 @@
+use spade::{HasPosition, Point2};
+
+#[derive(Clone, Copy, arbitrary::Arbitrary)]
+pub struct FuzzPoint {
+    pub x: f64,
+    pub y: f64,
+}
+
+impl HasPosition for FuzzPoint {
+    type Scalar = f64;
+    fn position(&self) -> Point2<f64> {
+        Point2::new(self.x, self.y)
+    }
+}
+
+impl core::fmt::Debug for FuzzPoint {
+    fn fmt(&self, f: &mut core::fmt::Formatter<'_>) -> core::fmt::Result {
+        f.write_fmt(format_args!("Point2::new({:?}, {:?})", self.x, self.y))
+    }
+}

src/cdt.rs

@@ -1,9 +1,11 @@
+use crate::delaunay_core::{bulk_load_cdt, bulk_load_stable};
 use crate::{delaunay_core::Dcel, intersection_iterator::LineIntersectionIterator};
 use crate::{handles::*, intersection_iterator::Intersection};
 use crate::{
     DelaunayTriangulation, HasPosition, HintGenerator, InsertionError, LastUsedVertexHintGenerator,
     Point2, Triangulation, TriangulationExt,
 };
+
 #[cfg(feature = "serde")]
 use serde::{Deserialize, Serialize};
 
@@ -235,8 +237,9 @@ where
     ) -> (
         Dcel<Self::Vertex, Self::DirectedEdge, Self::UndirectedEdge, Self::Face>,
         Self::HintGenerator,
+        usize,
     ) {
-        todo!()
+        (self.dcel, self.hint_generator, self.num_constraints)
     }
 }
 
@@ -270,6 +273,114 @@ where
     F: Default,
     L: HintGenerator<<V as HasPosition>::Scalar>,
 {
+    /// Efficient bulk loading of a constraint delaunay triangulation, including both vertices and constraint edges.
+    ///
+    /// The edges are given as pairs of vertex indices.
+    ///
+    /// Note that the vertex order is not preserved by this function - iterating through all vertices will not result in
+    /// the same sequence as the input vertices. Use [ConstrainedDelaunayTriangulation::bulk_load_cdt_stable] for a
+    /// slower but order preserving variant.
+    ///
+    /// Input vertices may have the same position. However, only one vertex for each position will be kept. Edges
+    /// that go to a discarded vertex are rerouted and still inserted.
+    /// It is arbitrary which duplicated vertex remains.
+    ///
+    /// # Example
+    /// ```
+    /// # fn main() -> Result<(), spade::InsertionError> {
+    /// use spade::{ConstrainedDelaunayTriangulation, Point2, Triangulation};
+    /// let mut vertices = vec![
+    ///     Point2::new(0.0, 1.0),
+    ///     Point2::new(1.0, 2.0),
+    ///     Point2::new(3.0, -3.0),
+    ///     Point2::new(-1.0, -2.0),
+    ///     Point2::new(-4.0, -5.0),
+    /// ];
+    /// let mut edges = vec![[0, 1], [1, 2], [2, 3], [3, 4]];
+    /// let cdt = ConstrainedDelaunayTriangulation::<_>::bulk_load_cdt(vertices.clone(), edges)?;
+    ///
+    /// assert_eq!(cdt.num_vertices(), 5);
+    /// assert_eq!(cdt.num_constraints(), 4);
+    /// // The order will usually change
+    /// assert_ne!(cdt.vertices().map(|v| v.position()).collect::<Vec<_>>(), vertices);
+    /// # Ok(())
+    /// # }
+    /// ```
+    ///
+    /// # Panics
+    /// Panics if any constraint edges overlap. Panics if the edges contain an invalid index (out of range).
+    pub fn bulk_load_cdt(vertices: Vec<V>, edges: Vec<[usize; 2]>) -> Result<Self, InsertionError> {
+        let mut result = bulk_load_cdt(vertices, edges)?;
+        *result.hint_generator_mut() = L::initialize_from_triangulation(&result);
+        Ok(result)
+    }
+
+    /// Stable bulk load variant that preserves the input vertex order
+    ///
+    /// The resulting vertex set will be equal to the input vertex set if their positions are all distinct.
+    /// See [ConstrainedDelaunayTriangulation::bulk_load_cdt] for additional details like panic behavior and duplicate
+    /// handling.
+    ///
+    /// # Example
+    /// ```
+    /// # fn main() -> Result<(), spade::InsertionError> {
+    /// use spade::{ConstrainedDelaunayTriangulation, Point2, Triangulation};
+    /// let mut vertices = vec![
+    ///     Point2::new(0.0, 1.0),
+    ///     Point2::new(1.0, 2.0),
+    ///     Point2::new(3.0, -3.0),
+    ///     Point2::new(-1.0, -2.0),
+    ///     Point2::new(-4.0, -5.0),
+    /// ];
+    /// let mut edges = vec![[0, 1], [1, 2], [2, 3], [3, 4]];
+    /// let cdt = ConstrainedDelaunayTriangulation::<_>::bulk_load_cdt_stable(vertices.clone(), edges)?;
+    ///
+    /// // The ordered will be preserved:
+    /// assert_eq!(cdt.vertices().map(|v| v.position()).collect::<Vec<_>>(), vertices);
+    /// # Ok(())
+    /// # }
+    /// ```
+    ///
+    /// It is fine to include vertex positions multiple times. The resulting order will be the same as if
+    /// the duplicates were removed prior to insertion. However, it is unclear *which* duplicates are
+    /// removed - e.g. do not assume that always the first duplicated vertex remains.
+    ///
+    /// ```
+    /// # fn main() -> Result<(), spade::InsertionError> {
+    /// use spade::{ConstrainedDelaunayTriangulation, Point2, Triangulation};
+    /// let mut vertices = vec![
+    ///     Point2::new(0.0, 1.0),
+    ///     Point2::new(1.0, 2.0), // Duplicate
+    ///     Point2::new(1.0, 2.0),
+    ///     Point2::new(3.0, -3.0),
+    ///     Point2::new(3.0, -3.0), // Duplicate
+    ///     Point2::new(-4.0, -5.0),
+    /// ];
+    /// let mut edges = vec![[0, 1], [2, 3], [4, 5]];
+    /// let cdt = ConstrainedDelaunayTriangulation::<_>::bulk_load_cdt_stable(vertices.clone(), edges)?;
+    ///
+    /// // The choice of deduplicated vertices is arbitrary. In this example, dedup[1] and dedup[2] could
+    /// // have been swapped
+    /// let dedup = [
+    ///     Point2::new(0.0, 1.0),
+    ///     Point2::new(1.0, 2.0),
+    ///     Point2::new(3.0, -3.0),
+    ///     Point2::new(-4.0, -5.0),
+    /// ];
+    /// assert_eq!(cdt.vertices().map(|v| v.position()).collect::<Vec<_>>(), dedup);
+    /// # Ok(())
+    /// # }
+    /// ```
+    pub fn bulk_load_cdt_stable(
+        vertices: Vec<V>,
+        edges: Vec<[usize; 2]>,
+    ) -> Result<Self, InsertionError> {
+        let mut result: Self =
+            bulk_load_stable(move |vertices| bulk_load_cdt(vertices, edges), vertices)?;
+        *result.hint_generator_mut() = L::initialize_from_triangulation(&result);
+        Ok(result)
+    }
+
     /// Removes a vertex from the triangulation.
     ///
     /// This operation runs in O(n²), where n is the degree of the
@@ -622,21 +733,32 @@ where
         }
     }
 
-    #[cfg(test)]
+    #[cfg(any(test, fuzzing))]
     pub fn cdt_sanity_check(&self) {
+        self.cdt_sanity_check_with_params(true);
+    }
+
+    #[cfg(any(test, fuzzing))]
+    pub fn cdt_sanity_check_with_params(&self, check_convexity: bool) {
         let num_undirected_edges = self
             .dcel
             .undirected_edges()
             .filter(|e| e.is_constraint_edge())
             .count();
 
         assert_eq!(num_undirected_edges, self.num_constraints());
-        self.basic_sanity_check();
+
+        if self.num_constraints() == 0 && check_convexity {
+            self.sanity_check();
+        } else {
+            self.basic_sanity_check(check_convexity);
+        }
     }
 }
 
 #[cfg(test)]
 mod test {
+
     use super::ConstrainedDelaunayTriangulation;
     use crate::test_utilities::*;
     use crate::{DelaunayTriangulation, InsertionError, Point2, Triangulation};