spade/
delaunay_triangulation.rs

1use super::delaunay_core::Dcel;
2use crate::{
3    handles::VertexHandle, ConstrainedDelaunayTriangulation, HasPosition, HintGenerator,
4    InsertionError, LastUsedVertexHintGenerator, NaturalNeighbor, Point2, Triangulation,
5    TriangulationExt,
6};
7use alloc::vec;
8use alloc::vec::Vec;
9
10use num_traits::Float;
11
12#[cfg(feature = "serde")]
13use serde::{Deserialize, Serialize};
14
15/// A two-dimensional [Delaunay triangulation](https://en.wikipedia.org/wiki/Delaunay_triangulation).
16///
17/// A Delaunay triangulation  a triangulation that fulfills the *Delaunay Property*: No
18/// vertex of the triangulation is contained in the
19/// [circumcircle](https://en.wikipedia.org/wiki/Circumscribed_circle) of any triangle.
20/// As a consequence, Delaunay triangulations are well suited to support interpolation
21/// algorithms and nearest neighbor searches. It is often constructed in order to retrieve its dual
22/// graph, the [Voronoi diagram](#voronoi-diagram).
23///
24#[doc = include_str!("../images/circumcircle.svg")]
25///
26/// *An example triangulation with a few circumcircles drawn. No circumcircle contains more than three vertices.*
27///
28/// Most methods on this type require the [Triangulation] trait. Refer to its documentation
29/// for more details on how to use `DelaunayTriangulation`.
30///
31/// # Basic Usage
32/// Vertices need to implement the [HasPosition] trait. Spade bundles
33/// the [Point2] struct for basic use cases.
34///
35/// ## Basic example
36///  ```
37/// use spade::{DelaunayTriangulation, Triangulation, Point2, InsertionError};
38///
39/// fn main() -> Result<(), InsertionError> {
40///
41///     let mut triangulation: DelaunayTriangulation<_> = DelaunayTriangulation::new();
42///
43///     // Insert three vertices that span one triangle (face)
44///     triangulation.insert(Point2::new(0.0, 1.0))?;
45///     triangulation.insert(Point2::new(1.0, 1.0))?;
46///     triangulation.insert(Point2::new(0.5, -1.0))?;
47///
48///     assert_eq!(triangulation.num_vertices(), 3);
49///     assert_eq!(triangulation.num_inner_faces(), 1);
50///     assert_eq!(triangulation.num_undirected_edges(), 3);
51///     Ok(())
52/// }
53/// ```
54/// ## Right handed and left-handed coordinate systems
55/// For simplicity, all method names and their documentation assume that the underlying coordinate system
56/// is right-handed (e.g. x-axis points to the right, y-axis points upwards). If a left-handed system
57/// (lhs) is used, any term related to orientation needs to be reversed:
58///  - "left" becomes "right" (example: the face of a directed edge is on the right side for a lhs)
59///  - "counter clock wise" becomes "clockwise" (example: the vertices of a face are returned in clock wise order for a lhs)
60///  
61/// <table>
62/// <tr><th>left-handed system</th><th>right-handed system</th></tr>
63/// <tr><td>
64#[doc = concat!(include_str!("../images/lhs.svg"), "</td><td>",include_str!("../images/rhs.svg"), " </td></tr></table>")]
65/// # Extracting geometry information
66///
67/// Spade uses [handles](crate::handles) to extract the triangulation's geometry.
68/// Handles are usually retrieved by inserting a vertex or by iterating.
69///  
70/// ## Example
71/// ```
72///  fn main() -> Result<(), spade::InsertionError> {
73/// use crate::spade::{DelaunayTriangulation, Triangulation, Point2};
74///
75/// let mut triangulation: DelaunayTriangulation<Point2<f64>> = DelaunayTriangulation::new();
76///
77/// triangulation.insert(Point2::new(0.0, 1.0))?;
78/// triangulation.insert(Point2::new(1.0, 1.0))?;
79/// triangulation.insert(Point2::new(0.5, -1.0))?;
80///
81/// for face in triangulation.inner_faces() {
82///   // face is a FaceHandle
83///   // edges is an array containing 3 directed edge handles
84///   let edges = face.adjacent_edges();
85///   for edge in &edges {
86///     let from = edge.from();
87/// let to = edge.to();
88///     // from and to are vertex handles
89///     println!("found an edge: {:?} -> {:?}", from, to);
90///   }
91///
92///   // vertices is an array containing 3 vertex handles
93///   let vertices = face.vertices();
94///   for vertex in &vertices {
95///     println!("Found vertex with position {:?}", vertex.position());
96///   }
97/// }
98/// # Ok(()) }
99/// ```
100///
101/// # Type parameters
102/// The triangulation's vertices, edges and faces can contain custom data.
103/// By default, the edge and face types are set to `()`. The vertex type must
104/// be specified.
105///
106///  * `V: HasPosition` The vertex type
107///  * `DE: Default` The directed edge type.
108///  * `UE: Default` The undirected edge type.
109///  * `F: Default` The face type.
110///
111/// Only vertices can be inserted directly. Faces and edges are create via `Default::default()`.
112/// Usually, edge and face data will need to be modified in a separate pass.
113///
114/// Setting any custom data works by calling [vertex_data_mut](Triangulation::vertex_data_mut),
115/// [directed_edge_data_mut](Triangulation::directed_edge_data_mut),
116/// [undirected_edge_data_mut](Triangulation::undirected_edge_data_mut) and
117/// [face_data_mut](Triangulation::face_data_mut).
118///  
119/// ## Example
120/// ```
121/// fn main() -> Result<(), spade::InsertionError> {
122/// use crate::spade::{DelaunayTriangulation, Triangulation, Point2};
123///
124/// // A custom undirected edge type used to cache the length of an edge
125/// #[derive(Default)]
126/// struct EdgeWithLength { length: f64 }
127///
128/// // Creates a new triangulation with a custom undirected edge type
129/// let mut triangulation: DelaunayTriangulation<Point2<f64>, (), EdgeWithLength>
130///                          = DelaunayTriangulation::new();
131///
132/// triangulation.insert(Point2::new(0.0, 1.0))?;
133/// triangulation.insert(Point2::new(1.0, 1.0))?;
134/// triangulation.insert(Point2::new(0.5, -1.0))?;
135///
136/// for edge in triangulation.fixed_undirected_edges() {
137///   let positions = triangulation.undirected_edge(edge).positions();
138///   let length = positions[0].distance_2(positions[1]).sqrt();
139///   // Write length into the edge data
140///   triangulation.undirected_edge_data_mut(edge).length = length;
141/// }
142///
143/// for edge in triangulation.undirected_edges() {
144///    let length = edge.data().length;
145///    assert!(length > 0.0);
146/// }
147/// # Ok(()) }
148/// ```
149///
150/// # Outer face
151/// Every triangulation contains an *outer face* that is adjacent to all edges of the
152/// triangulation's convex hull. The outer face is even present for a triangulation without
153/// vertices. It is either referenced by calling [Triangulation::outer_face()] or
154/// [handles::OUTER_FACE](crate::handles::OUTER_FACE)
155///
156#[doc = include_str!("../images/outer_faces.svg")]
157///
158/// # Voronoi Diagram
159///
160/// the dual graph of the Delaunay triangulation is the *Voronoi diagram*. The Voronoi diagram
161/// subdivides the plane into several *Voronoi cells* (called `VoronoiFace` by Spade for consistency)
162/// which contain all points in the plane that share a common nearest neighbor.
163///
164/// Keep in mind that "faces" and "vertices" are swapped - an (inner) Voronoi *vertex*
165/// corresponds to a single Delaunay *face*.
166/// The position of an inner voronoi vertex is the *circumcenter* of its dual Delaunay
167/// face.
168///
169#[doc = include_str!("../images/basic_voronoi.svg")]
170///
171/// *A Delaunay triangulation (black lines) and its dual graph, the Voronoi diagram (blue lines)*
172///
173/// ## Extracting the Voronoi Diagram
174/// Spade defines various functions to extract information about the Voronoi diagram:
175///
176/// **Types**
177///  * [DirectedVoronoiEdge](crate::handles::DirectedVoronoiEdge)
178///  * [UndirectedVoronoiEdge](crate::handles::UndirectedVoronoiEdge)
179///  * [VoronoiVertex](crate::handles::VoronoiVertex)
180///  * [VoronoiFace](crate::handles::VoronoiFace)
181///
182/// **Iterators**
183///  * [Triangulation::directed_voronoi_edges()]
184///  * [Triangulation::undirected_voronoi_edges()]
185///
186/// **Conversion**
187///  * [DirectedVoronoiEdge::as_undirected()](crate::handles::DirectedVoronoiEdge::as_undirected())
188///  * [UndirectedVoronoiEdge::as_directed()](crate::handles::UndirectedVoronoiEdge::as_directed())
189///  * [DirectedEdgeHandle::as_voronoi_edge()](crate::handles::DirectedEdgeHandle::as_voronoi_edge())
190///  * [DirectedVoronoiEdge::as_delaunay_edge()](crate::handles::DirectedVoronoiEdge::as_delaunay_edge())
191///  * [UndirectedEdgeHandle::as_voronoi_edge()](crate::handles::UndirectedEdgeHandle::as_voronoi_edge())
192///  * [UndirectedVoronoiEdge::as_delaunay_edge()](crate::handles::UndirectedVoronoiEdge::as_delaunay_edge())
193///
194/// ## Extracting the Voronoi Diagram (Example)
195/// Extracting the geometry of the voronoi diagram can be slightly tricky as some of the voronoi
196/// extend into infinity (see the dashed lines in the example above).
197///
198/// ```
199/// use spade::handles::{VoronoiVertex::Inner, VoronoiVertex::Outer};
200/// use spade::{DelaunayTriangulation, Point2, Triangulation};
201///
202/// // Prints out the location of all voronoi edges in a triangulation
203/// fn log_voronoi_diagram(triangulation: &DelaunayTriangulation<Point2<f64>>) {
204///     for edge in triangulation.undirected_voronoi_edges() {
205///         match edge.vertices() {
206///             [Inner(from), Inner(to)] => {
207///                 // "from" and "to" are inner faces of the Delaunay triangulation
208///                 println!(
209///                     "Found voronoi edge between {:?} and {:?}",
210///                     from.circumcenter(),
211///                     to.circumcenter()
212///                 );
213///             }
214///             [Inner(from), Outer(edge)] | [Outer(edge), Inner(from)] => {
215///                 // Some lines don't have a finite end and extend into infinity.
216///                 println!(
217///                     "Found infinite voronoi edge going out of {:?} into the direction {:?}",
218///                     from.circumcenter(),
219///                     edge.direction_vector()
220///                 );
221///             }
222///             [Outer(_), Outer(_)] => {
223///                 // This case only happens if all vertices of the triangulation lie on the
224///                 // same line and can probably be ignored.
225///             }
226///         }
227///     }
228/// }
229/// ```
230///
231/// # Performance tuning
232///
233/// Fine-tuning a Delaunay triangulation can be more tricky from time to time. However, some will *nearly always* be
234/// the right thing to do:
235///
236/// - Measure, don't guess. Execution times are hard to predict.
237/// - If you plan to perform several random access queries (e.g. looking up the point at an arbitrary position):
238///   Consider using `[HierarchyHintGenerator](crate::HierarchyHintGenerator)
239/// - For data sets with uniformly distributed vertices: Use [HierarchyHintGenerator](crate::HierarchyHintGenerator) if
240///   bulk loading is not applicable.
241/// - For data sets where vertices are inserted in close local proximity (each vertex is not too far away from the
242///   previously inserted vertex): Consider using [LastUsedVertexHintGenerator](LastUsedVertexHintGenerator).
243/// - Try to avoid large custom data types for edges, vertices and faces.
244/// - Using `f64` and `f32` as scalar type will usually end up roughly having the same run time performance.
245/// - Prefer using [bulk_load](Triangulation::bulk_load) over [insert](Triangulation::insert).
246/// - The run time of all vertex operations (insertion, removal and lookup) is roughly the same for larger triangulations.
247///   
248/// ## Complexity classes
249///
250/// This table display the average and amortized cost for inserting a vertex into a triangulation with `n` vertices.
251///
252/// |                             | Uniformly distributed vertices | Insertion of vertices with local proximity |
253/// |-----------------------------|--------------------------------|--------------------------------------------|
254/// | LastUsedVertexHintGenerator |        O(sqrt(n)) (worst case) |                  O(1) (best case), fastest |
255/// | HierarchyHintGenerator      |       O(log(n)) (Average case) |                           O(1) (best case) |
256///
257/// # See also
258/// Delaunay triangulations are closely related to [constrained Delaunay triangulations](crate::ConstrainedDelaunayTriangulation)
259#[doc(alias = "Voronoi")]
260#[doc(alias = "Voronoi diagram")]
261#[doc(alias = "Delaunay")]
262#[derive(Debug, Clone)]
263#[cfg_attr(
264    feature = "serde",
265    derive(Serialize, Deserialize),
266    serde(crate = "serde")
267)]
268pub struct DelaunayTriangulation<V, DE = (), UE = (), F = (), L = LastUsedVertexHintGenerator>
269where
270    V: HasPosition,
271    DE: Default,
272    UE: Default,
273    F: Default,
274    L: HintGenerator<<V as HasPosition>::Scalar>,
275{
276    pub(crate) dcel: Dcel<V, DE, UE, F>,
277    pub(crate) hint_generator: L,
278}
279
280impl<V, DE, UE, F, L> DelaunayTriangulation<V, DE, UE, F, L>
281where
282    V: HasPosition,
283    DE: Default,
284    UE: Default,
285    F: Default,
286    L: HintGenerator<<V as HasPosition>::Scalar>,
287{
288    /// Returns the nearest neighbor of a given input vertex.
289    ///
290    /// Returns `None` if the triangulation is empty.
291    ///
292    /// # Runtime
293    /// This method takes `O(sqrt(n))` on average where n is the number of vertices.
294    pub fn nearest_neighbor(
295        &self,
296        position: Point2<<V as HasPosition>::Scalar>,
297    ) -> Option<VertexHandle<'_, V, DE, UE, F>> {
298        if self.num_vertices() == 0 {
299            return None;
300        }
301
302        let hint = self.hint_generator().get_hint(position);
303        let hint = self.validate_vertex_handle(hint);
304
305        let vertex = self.walk_to_nearest_neighbor(hint, position);
306        self.hint_generator().notify_vertex_lookup(vertex.fix());
307        Some(vertex)
308    }
309
310    /// Creates a new delaunay triangulation with an efficient bulk loading strategy.
311    ///
312    /// In contrast to [Triangulation::bulk_load], this method will create a triangulation with
313    /// vertices returned *in the same order* as the input vertices.
314    ///
315    /// # Duplicate handling
316    ///
317    /// For every set of vertices with identical positions, only the vertex with the lowest index
318    /// is kept.
319    ///
320    /// For example, if the input vertices are `[v0, v1, v2, v1]` (where v1 is duplicated), the
321    /// resulting triangulation will be `[v0, v1, v2]`.
322    ///
323    /// Consider checking the triangulation's size after calling this method to ensure that no
324    /// duplicates were present. Removing duplicates requires additional work and slightly
325    /// increases the run time.
326    ///
327    /// # Run Time
328    ///
329    /// `O(n*log(n))` for `n` input vertices. Slightly (5% - 10%) slower than [Triangulation::bulk_load].
330    ///
331    /// # Example
332    /// ```
333    /// # use spade::InsertionError;
334    /// use spade::{DelaunayTriangulation, Point2, Triangulation};
335    ///
336    /// # fn main() -> Result<(), InsertionError> {
337    /// let vertices = vec![
338    ///      Point2::new(0.5, 1.0),
339    ///      Point2::new(-0.5, 2.0),
340    ///      Point2::new(0.2, 3.0),
341    ///      Point2::new(0.0, 4.0),
342    ///      Point2::new(-0.2, 5.0)
343    /// ];
344    /// let triangulation = DelaunayTriangulation::<Point2<f64>>::bulk_load_stable(vertices.clone())?;
345    /// // This assert will not hold for regular bulk loading!
346    /// assert_eq!(triangulation.vertices().map(|v| *v.data()).collect::<Vec<_>>(), vertices);
347    ///
348    /// // This is how you would check for duplicates:
349    /// let duplicates_found = triangulation.num_vertices() < vertices.len();
350    /// assert!(!duplicates_found);
351    /// # Ok(()) }
352    /// ```
353    pub fn bulk_load_stable(elements: Vec<V>) -> Result<Self, InsertionError> {
354        let cdt = ConstrainedDelaunayTriangulation::bulk_load_cdt(elements, vec![])?;
355        let result = Self::from_cdt(cdt);
356
357        Ok(result)
358    }
359
360    fn from_cdt(
361        cdt: ConstrainedDelaunayTriangulation<V, DE, UE, F, L>,
362    ) -> DelaunayTriangulation<V, DE, UE, F, L> {
363        let (dcel, hint_generator, num_constraints) = cdt.into_parts();
364        let dcel = dcel.map_undirected_edges(|cdt_edge| cdt_edge.deconstruct().1);
365        assert_eq!(num_constraints, 0);
366        DelaunayTriangulation {
367            dcel,
368            hint_generator,
369        }
370    }
371}
372
373impl<V, DE, UE, F, L> Default for DelaunayTriangulation<V, DE, UE, F, L>
374where
375    V: HasPosition,
376    DE: Default,
377    UE: Default,
378    F: Default,
379    L: HintGenerator<<V as HasPosition>::Scalar>,
380{
381    fn default() -> Self {
382        Self {
383            dcel: Default::default(),
384            hint_generator: Default::default(),
385        }
386    }
387}
388
389impl<V, DE, UE, F, L> DelaunayTriangulation<V, DE, UE, F, L>
390where
391    V: HasPosition,
392    DE: Default,
393    UE: Default,
394    F: Default,
395    V::Scalar: Float,
396    L: HintGenerator<<V as HasPosition>::Scalar>,
397{
398    /// Allows using natural neighbor interpolation on this triangulation. Refer to the documentation
399    /// of [NaturalNeighbor] for more information.
400    pub fn natural_neighbor(&self) -> NaturalNeighbor<'_, Self> {
401        NaturalNeighbor::new(self)
402    }
403}
404
405impl<V, DE, UE, F, L> Triangulation for DelaunayTriangulation<V, DE, UE, F, L>
406where
407    V: HasPosition,
408    DE: Default,
409    UE: Default,
410    F: Default,
411    L: HintGenerator<<V as HasPosition>::Scalar>,
412{
413    type Vertex = V;
414    type DirectedEdge = DE;
415    type UndirectedEdge = UE;
416    type Face = F;
417    type HintGenerator = L;
418
419    fn s(&self) -> &Dcel<V, DE, UE, F> {
420        &self.dcel
421    }
422
423    fn s_mut(&mut self) -> &mut Dcel<V, DE, UE, F> {
424        &mut self.dcel
425    }
426
427    fn hint_generator(&self) -> &Self::HintGenerator {
428        &self.hint_generator
429    }
430
431    fn hint_generator_mut(&mut self) -> &mut Self::HintGenerator {
432        &mut self.hint_generator
433    }
434
435    fn from_parts(
436        dcel: Dcel<Self::Vertex, Self::DirectedEdge, Self::UndirectedEdge, Self::Face>,
437        hint_generator: Self::HintGenerator,
438        num_constraints: usize,
439    ) -> Self {
440        assert_eq!(num_constraints, 0);
441        Self {
442            dcel,
443            hint_generator,
444        }
445    }
446
447    fn into_parts(
448        self,
449    ) -> (
450        Dcel<Self::Vertex, Self::DirectedEdge, Self::UndirectedEdge, Self::Face>,
451        Self::HintGenerator,
452        usize,
453    ) {
454        (self.dcel, self.hint_generator, 0)
455    }
456}
457
458#[cfg(test)]
459mod test {
460    use crate::test_utilities::{random_points_with_seed, SEED};
461
462    use crate::{DelaunayTriangulation, InsertionError, Point2, Triangulation};
463
464    #[allow(unused)]
465    #[cfg(feature = "serde")]
466    // Just needs to compile
467    fn check_serde() {
468        use serde::{Deserialize, Serialize};
469
470        use crate::{HierarchyHintGenerator, LastUsedVertexHintGenerator, Point2};
471
472        fn requires_serde<'de, T: Serialize + Deserialize<'de>>() {}
473
474        type DT<L> = DelaunayTriangulation<Point2<f64>, (), (), (), L>;
475
476        requires_serde::<DT<LastUsedVertexHintGenerator>>();
477        requires_serde::<DT<HierarchyHintGenerator<f64>>>();
478    }
479
480    #[test]
481    fn test_nearest_neighbor() -> Result<(), InsertionError> {
482        const SIZE: usize = 54;
483        let points = random_points_with_seed(SIZE, SEED);
484
485        let d = DelaunayTriangulation::<_>::bulk_load(points.clone())?;
486
487        let sample_points = random_points_with_seed(SIZE * 3, SEED);
488        for p in sample_points {
489            let nn_delaunay = d.nearest_neighbor(p);
490            let nn_linear_search = points.iter().min_by(|l, r| {
491                let d1 = l.distance_2(p);
492                let d2 = r.distance_2(p);
493                d1.partial_cmp(&d2).unwrap()
494            });
495            assert_eq!(nn_delaunay.map(|p| p.position()), nn_linear_search.cloned());
496        }
497        Ok(())
498    }
499
500    #[test]
501    fn test_nearest_neighbor_small() -> Result<(), InsertionError> {
502        let mut d = DelaunayTriangulation::<_>::new();
503        assert_eq!(None, d.nearest_neighbor(Point2::new(0.0, 0.0)));
504
505        d.insert(Point2::new(0.0, 0.0))?;
506        assert!(d.nearest_neighbor(Point2::new(0.0, 1.0)).is_some());
507        Ok(())
508    }
509
510    #[test]
511    #[allow(clippy::redundant_clone)]
512    #[allow(unused_must_use)]
513    fn test_clone_is_implemented() {
514        // Just needs to compile
515        DelaunayTriangulation::<Point2<f64>>::new().clone();
516    }
517}
spade/delaunay_triangulation.rs

spade/
delaunay_triangulation.rs
