1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
use super::delaunay_core::Dcel;
use crate::{
    delaunay_core::bulk_load, handles::VertexHandle, HasPosition, HintGenerator, InsertionError,
    LastUsedVertexHintGenerator, NaturalNeighbor, Point2, Triangulation, TriangulationExt,
};

use alloc::vec::Vec;
use num_traits::Float;

#[cfg(feature = "serde")]
use serde::{Deserialize, Serialize};

/// A two dimensional [Delaunay triangulation](https://en.wikipedia.org/wiki/Delaunay_triangulation).
///
/// A Delaunay triangulation  a triangulation that fulfills the *Delaunay Property*: No
/// vertex of the triangulation is contained in the
/// [circumcircle](https://en.wikipedia.org/wiki/Circumscribed_circle) of any triangle.
/// As a consequence, Delaunay triangulations are well suited to support interpolation
/// algorithms and nearest neighbor searches. It is often constructed in order to retrieve its dual
/// graph, the [Voronoi diagram](#voronoi-diagram).
///
#[doc = include_str!("../images/circumcircle.svg")]
///
/// *An example triangulation with a few circumcircles drawn. No circumcircle contains more than three vertices.*
///
/// Most methods on this type require the [Triangulation] trait. Refer to its documentation
/// for more details on how to use `DelaunayTriangulation`.
///
/// # Basic Usage
/// Vertices need to implement the [HasPosition] trait. Spade bundles
/// the [Point2](crate::Point2) struct for basic use cases.
///
/// ## Basic example
///  ```
/// use spade::{DelaunayTriangulation, Triangulation, Point2, InsertionError};
///
/// fn main() -> Result<(), InsertionError> {
///
///     let mut triangulation: DelaunayTriangulation<_> = DelaunayTriangulation::new();
///
///     // Insert three vertices that span one triangle (face)
///     triangulation.insert(Point2::new(0.0, 1.0))?;
///     triangulation.insert(Point2::new(1.0, 1.0))?;
///     triangulation.insert(Point2::new(0.5, -1.0))?;
///
///     assert_eq!(triangulation.num_vertices(), 3);
///     assert_eq!(triangulation.num_inner_faces(), 1);
///     assert_eq!(triangulation.num_undirected_edges(), 3);
///     Ok(())
/// }
/// ```
/// ## Right handed and left handed coordinate systems
/// For simplicity, all method names and their documentation assume that the underlying coordinate system
/// is right handed (e.g. x axis points to the right, y axis points upwards). If a left handed system
/// (lhs) is used, any term related to orientation needs to be reversed:
///  - "left" becomes "right" (example: the face of a directed edge is on the right side for a lhs
///  - "counter clock wise" becomes "clockwise" (example: the vertices of a face are returned in clock wise order for a lhs)
///  
/// <table>
/// <tr><th>left handed system</th><th>right handed system</th></tr>
/// <tr><td>
#[doc = concat!(include_str!("../images/lhs.svg"), "</td><td>",include_str!("../images/rhs.svg"), " </td></tr></table>")]
/// # Extracting geometry information
///
/// Spade uses [handles](crate::handles) to extract the triangulation's geometry.
/// Handles are usually retrieved by inserting a vertex or by iterating.
///  
/// ## Example
/// ```
///  fn main() -> Result<(), spade::InsertionError> {
/// use crate::spade::{DelaunayTriangulation, Triangulation, Point2};
///
/// let mut triangulation: DelaunayTriangulation<Point2<f64>> = DelaunayTriangulation::new();
///
/// triangulation.insert(Point2::new(0.0, 1.0))?;
/// triangulation.insert(Point2::new(1.0, 1.0))?;
/// triangulation.insert(Point2::new(0.5, -1.0))?;
///
/// for face in triangulation.inner_faces() {
///   // face is a FaceHandle
///   // edges is an array containing 3 directed edge handles
///   let edges = face.adjacent_edges();
///   for edge in &edges {
///     let from = edge.from();
/// let to = edge.to();
///     // from and to are vertex handles
///     println!("found an edge: {:?} -> {:?}", from, to);
///   }
///
///   // vertices is an array containing 3 vertex handles
///   let vertices = face.vertices();
///   for vertex in &vertices {
///     println!("Found vertex with position {:?}", vertex.position());
///   }
/// }
/// # Ok(()) }
/// ```
///
/// # Type parameters
/// The triangulation's vertices, edges and faces can contain custom data.
/// By default, the edge and face types are set to `()`. The vertex type must
/// be specified.
///
///  * `V: HasPosition` The vertex type
///  * `DE: Default` The directed edge type.
///  * `UE: Default` The undirected edge type.
///  * `F: Default` The face type.
///
/// Only vertices can be inserted directly. Faces and edges are create via `Default::default()`.
/// Usually, edge and face data will need to be modified in a separate pass.
///
/// Setting any custom data works by calling [vertex_data_mut](Triangulation::vertex_data_mut),
/// [directed_edge_data_mut](Triangulation::directed_edge_data_mut),
/// [undirected_edge_data_mut](Triangulation::undirected_edge_data_mut) and
/// [face_data_mut](Triangulation::face_data_mut).
///  
/// ## Example
/// ```
/// fn main() -> Result<(), spade::InsertionError> {
/// use crate::spade::{DelaunayTriangulation, Triangulation, Point2};
///
/// // A custom undirected edge type used to cache the length of an edge
/// #[derive(Default)]
/// struct EdgeWithLength { length: f64 }
///
/// // Creates a new triangulation with a custom undirected edge type
/// let mut triangulation: DelaunayTriangulation<Point2<f64>, (), EdgeWithLength>
///                          = DelaunayTriangulation::new();
///
/// triangulation.insert(Point2::new(0.0, 1.0))?;
/// triangulation.insert(Point2::new(1.0, 1.0))?;
/// triangulation.insert(Point2::new(0.5, -1.0))?;
///
/// for edge in triangulation.fixed_undirected_edges() {
///   let positions = triangulation.undirected_edge(edge).positions();
///   let length = positions[0].distance_2(positions[1]).sqrt();
///   // Write length into the edge data
///   triangulation.undirected_edge_data_mut(edge).length = length;
/// }
///
/// for edge in triangulation.undirected_edges() {
///    let length = edge.data().length;
///    assert!(length > 0.0);
/// }
/// # Ok(()) }
/// ```
///
/// # Outer face
/// Every triangulation contains an *outer face* that is adjacent to all edges of the
/// triangulation's convex hull. The outer face is even present for a triangulation without
/// vertices. It is either referenced by calling [Triangulation::outer_face()] or
/// [handles::OUTER_FACE](crate::handles::OUTER_FACE)
///
#[doc = include_str!("../images/outer_faces.svg")]
///
/// # Voronoi Diagram
///
/// the dual graph of the Delaunay triangulation is the *Voronoi diagram*. The Voronoi diagram
/// subdivides the plane into several *Voronoi cells* (called `VoronoiFace` by Spade for consistency)
/// which contain all points in the plane that share a common nearest neighbor.
///
/// Keep in mind that "faces" and "vertices" are swapped - an (inner) Voronoi *vertex*
/// corresponds to a single Delaunay *face*.
/// The position of an inner voronoi vertex is the *circumcenter* of its dual Delaunay
/// face.
///
#[doc = include_str!("../images/basic_voronoi.svg")]
///
/// *A Delaunay triangulation (black lines) and its dual graph, the Voronoi diagram (blue lines)*
///
/// ## Extracting the Voronoi Diagram
/// Spade defines various functions to extract information about the Voronoi diagram:
///
/// **Types**
///  * [DirectedVoronoiEdge](crate::handles::DirectedVoronoiEdge)
///  * [UndirectedVoronoiEdge](crate::handles::UndirectedVoronoiEdge)
///  * [VoronoiVertex](crate::handles::VoronoiVertex)
///  * [VoronoiFace](crate::handles::VoronoiFace)
///
/// **Iterators**
///  * [Triangulation::directed_voronoi_edges()]
///  * [Triangulation::undirected_voronoi_edges()]
///
/// **Conversion**
///  * [DirectedVoronoiEdge::as_undirected()](crate::handles::DirectedVoronoiEdge::as_undirected())
///  * [UndirectedVoronoiEdge::as_directed()](crate::handles::UndirectedVoronoiEdge::as_directed())
///  * [DirectedEdgeHandle::as_voronoi_edge()](crate::handles::DirectedEdgeHandle::as_voronoi_edge())
///  * [DirectedVoronoiEdge::as_delaunay_edge()](crate::handles::DirectedVoronoiEdge::as_delaunay_edge())
///  * [UndirectedEdgeHandle::as_voronoi_edge()](crate::handles::UndirectedEdgeHandle::as_voronoi_edge())
///  * [UndirectedVoronoiEdge::as_delaunay_edge()](crate::handles::UndirectedVoronoiEdge::as_delaunay_edge())
///
/// ## Extracting the Voronoi Diagram (Example)
/// Extracting the geometry of the voronoi diagram can be slightly tricky as some of the voronoi
/// extend into infinity (see the dashed lines in the example above).
///
/// ```
/// use spade::handles::{VoronoiVertex::Inner, VoronoiVertex::Outer};
/// use spade::{DelaunayTriangulation, Point2, Triangulation};
///
/// // Prints out the location of all voronoi edges in a triangulation
/// fn log_voronoi_diagram(triangulation: &DelaunayTriangulation<Point2<f64>>) {
///     for edge in triangulation.undirected_voronoi_edges() {
///         match edge.vertices() {
///             [Inner(from), Inner(to)] => {
///                 // "from" and "to" are inner faces of the Delaunay triangulation
///                 println!(
///                     "Found voronoi edge between {:?} and {:?}",
///                     from.circumcenter(),
///                     to.circumcenter()
///                 );
///             }
///             [Inner(from), Outer(edge)] | [Outer(edge), Inner(from)] => {
///                 // Some lines don't have a finite end and extend into infinity.
///                 println!(
///                     "Found infinite voronoi edge going out of {:?} into the direction {:?}",
///                     from.circumcenter(),
///                     edge.direction_vector()
///                 );
///             }
///             [Outer(_), Outer(_)] => {
///                 // This case only happens if all vertices of the triangulation lie on the
///                 // same line and can probably be ignored.
///             }
///         }
///     }
/// }
/// ```
///
/// # Performance tuning
///
/// Fine-tuning a Delaunay triangulation can be more tricky from time to time. However, some will *nearly always* be
/// the right thing to do:
///
/// - Measure, don't guess. Execution times are hard to predict.
/// - If you plan to perform several random access queries (e.g. looking up the point at an arbitrary position):
///   Consider using `[HierarchyHintGenerator](crate::HierarchyHintGenerator)
/// - For data sets with uniformly distributed vertices: Use [HierarchyHintGenerator](crate::HierarchyHintGenerator) if
///   bulk loading is not applicable.
/// - For data sets where vertices are inserted in close local proximity (each vertex is not too far away from the
///   previously inserted vertex): Consider using [LastUsedVertexHintGenerator](crate::LastUsedVertexHintGenerator).
/// - Try to avoid large custom data types for edges, vertices and faces.
/// - Using `f64` and `f32` as scalar type will usually end up roughly having the same run time performance.
/// - Prefer using [bulk_load](Triangulation::bulk_load) over [insert](Triangulation::insert).
/// - The run time of all vertex operations (insertion, removal and lookup) is roughly the same for larger triangulations.
///   
/// ## Complexity classes
///
/// This table display the average and amortized cost for inserting a vertex into a triangulation with `n` vertices.
///
/// |                             | Uniformly distributed vertices | Insertion of vertices with local proximity |
/// |-----------------------------|--------------------------------|--------------------------------------------|
/// | LastUsedVertexHintGenerator |        O(sqrt(n)) (worst case) |                  O(1) (best case), fastest |
/// | HierarchyHintGenerator      |       O(log(n)) (Average case) |                           O(1) (best case) |
///
/// # See also
/// Delaunay triangulations are closely related to [constrained Delaunay triangulations](crate::ConstrainedDelaunayTriangulation)
#[doc(alias = "Voronoi")]
#[doc(alias = "Voronoi diagram")]
#[doc(alias = "Delaunay")]
#[derive(Debug, Clone)]
#[cfg_attr(
    feature = "serde",
    derive(Serialize, Deserialize),
    serde(crate = "serde")
)]
pub struct DelaunayTriangulation<V, DE = (), UE = (), F = (), L = LastUsedVertexHintGenerator>
where
    V: HasPosition,
    DE: Default,
    UE: Default,
    F: Default,
    L: HintGenerator<<V as HasPosition>::Scalar>,
{
    pub(crate) dcel: Dcel<V, DE, UE, F>,
    pub(crate) hint_generator: L,
}

impl<V, DE, UE, F, L> DelaunayTriangulation<V, DE, UE, F, L>
where
    V: HasPosition,
    DE: Default,
    UE: Default,
    F: Default,
    L: HintGenerator<<V as HasPosition>::Scalar>,
{
    /// Returns the nearest neighbor of a given input vertex.
    ///
    /// Returns `None` if the triangulation is empty.
    ///
    /// # Runtime
    /// This method takes `O(sqrt(n))` on average where n is the number of vertices.
    pub fn nearest_neighbor(
        &self,
        position: Point2<<V as HasPosition>::Scalar>,
    ) -> Option<VertexHandle<V, DE, UE, F>> {
        if self.num_vertices() == 0 {
            return None;
        }

        let hint = self.hint_generator().get_hint(position);
        let hint = self.validate_vertex_handle(hint);

        let vertex = self.walk_to_nearest_neighbor(hint, position);
        self.hint_generator().notify_vertex_lookup(vertex.fix());
        Some(vertex)
    }

    /// Creates a new delaunay triangulation with an efficient bulk loading strategy.
    ///
    /// In contrast to [Triangulation::bulk_load], this method will create a triangulation with
    /// vertices returned *in the same order* as the input vertices.
    ///
    /// # Duplicate handling
    ///
    /// If two vertices have the same position, only one of them will be included in the final
    /// triangulation. It is undefined which of them is discarded.
    ///
    /// For example, if the input vertices are [v0, v1, v2, v1] (where v1 is duplicated), the
    /// resulting triangulation will be either
    /// [v0, v1, v2] or [v0, v2, v1]
    ///
    /// Consider checking the triangulation's size after calling this method to ensure that no
    /// duplicates were present.
    ///
    /// # Example
    /// ```
    /// # use spade::InsertionError;
    /// use spade::{DelaunayTriangulation, Point2, Triangulation};
    ///
    /// # fn main() -> Result<(), InsertionError> {
    /// let vertices = vec![
    ///      Point2::new(0.5, 1.0),
    ///      Point2::new(-0.5, 2.0),
    ///      Point2::new(0.2, 3.0),
    ///      Point2::new(0.0, 4.0),
    ///      Point2::new(-0.2, 5.0)
    /// ];
    /// let triangulation = DelaunayTriangulation::<Point2<f64>>::bulk_load_stable(vertices.clone())?;
    /// // This assert will not hold for regular bulk loading!
    /// assert_eq!(triangulation.vertices().map(|v| *v.data()).collect::<Vec<_>>(), vertices);
    ///
    /// // This is how you would check for duplicates:
    /// let duplicates_found = triangulation.num_vertices() < vertices.len();
    /// assert!(!duplicates_found);
    /// # Ok(()) }
    /// ```
    pub fn bulk_load_stable(elements: Vec<V>) -> Result<Self, InsertionError> {
        let mut result: Self = crate::delaunay_core::bulk_load_stable::<
            _,
            _,
            DelaunayTriangulation<_, _, _, _, _>,
            _,
        >(bulk_load, elements)?;
        *result.hint_generator_mut() = L::initialize_from_triangulation(&result);
        Ok(result)
    }
}

impl<V, DE, UE, F, L> Default for DelaunayTriangulation<V, DE, UE, F, L>
where
    V: HasPosition,
    DE: Default,
    UE: Default,
    F: Default,
    L: HintGenerator<<V as HasPosition>::Scalar>,
{
    fn default() -> Self {
        Self {
            dcel: Default::default(),
            hint_generator: Default::default(),
        }
    }
}

impl<V, DE, UE, F, L> DelaunayTriangulation<V, DE, UE, F, L>
where
    V: HasPosition,
    DE: Default,
    UE: Default,
    F: Default,
    V::Scalar: Float,
    L: HintGenerator<<V as HasPosition>::Scalar>,
{
    /// Allows using natural neighbor interpolation on this triangulation. Refer to the documentation
    /// of [NaturalNeighbor] for more information.
    pub fn natural_neighbor(&self) -> NaturalNeighbor<Self> {
        NaturalNeighbor::new(self)
    }
}

impl<V, DE, UE, F, L> Triangulation for DelaunayTriangulation<V, DE, UE, F, L>
where
    V: HasPosition,
    DE: Default,
    UE: Default,
    F: Default,
    L: HintGenerator<<V as HasPosition>::Scalar>,
{
    type Vertex = V;
    type DirectedEdge = DE;
    type UndirectedEdge = UE;
    type Face = F;
    type HintGenerator = L;

    fn s(&self) -> &Dcel<V, DE, UE, F> {
        &self.dcel
    }

    fn s_mut(&mut self) -> &mut Dcel<V, DE, UE, F> {
        &mut self.dcel
    }

    fn hint_generator(&self) -> &Self::HintGenerator {
        &self.hint_generator
    }

    fn hint_generator_mut(&mut self) -> &mut Self::HintGenerator {
        &mut self.hint_generator
    }

    fn from_parts(
        dcel: Dcel<Self::Vertex, Self::DirectedEdge, Self::UndirectedEdge, Self::Face>,
        hint_generator: Self::HintGenerator,
        num_constraints: usize,
    ) -> Self {
        assert_eq!(num_constraints, 0);
        Self {
            dcel,
            hint_generator,
        }
    }

    fn into_parts(
        self,
    ) -> (
        Dcel<Self::Vertex, Self::DirectedEdge, Self::UndirectedEdge, Self::Face>,
        Self::HintGenerator,
        usize,
    ) {
        (self.dcel, self.hint_generator, 0)
    }
}

#[cfg(test)]
mod test {
    use crate::test_utilities::{random_points_with_seed, SEED};

    use crate::{DelaunayTriangulation, InsertionError, Point2, Triangulation};

    #[allow(unused)]
    #[cfg(feature = "serde")]
    // Just needs to compile
    fn check_serde() {
        use serde::{Deserialize, Serialize};

        use crate::{HierarchyHintGenerator, LastUsedVertexHintGenerator, Point2};

        fn requires_serde<'de, T: Serialize + Deserialize<'de>>() {}

        type DT<L> = super::DelaunayTriangulation<Point2<f64>, (), (), (), L>;

        requires_serde::<DT<LastUsedVertexHintGenerator>>();
        requires_serde::<DT<HierarchyHintGenerator<f64>>>();
    }

    #[test]
    fn test_nearest_neighbor() -> Result<(), InsertionError> {
        const SIZE: usize = 54;
        let points = random_points_with_seed(SIZE, SEED);

        let d = DelaunayTriangulation::<_>::bulk_load(points.clone())?;

        let sample_points = random_points_with_seed(SIZE * 3, SEED);
        for p in sample_points {
            let nn_delaunay = d.nearest_neighbor(p);
            let nn_linear_search = points.iter().min_by(|l, r| {
                let d1 = l.distance_2(p);
                let d2 = r.distance_2(p);
                d1.partial_cmp(&d2).unwrap()
            });
            assert_eq!(nn_delaunay.map(|p| p.position()), nn_linear_search.cloned());
        }
        Ok(())
    }

    #[test]
    fn test_nearest_neighbor_small() -> Result<(), InsertionError> {
        let mut d = DelaunayTriangulation::<_>::new();
        assert_eq!(None, d.nearest_neighbor(Point2::new(0.0, 0.0)));

        d.insert(Point2::new(0.0, 0.0))?;
        assert!(d.nearest_neighbor(Point2::new(0.0, 1.0)).is_some());
        Ok(())
    }

    #[test]
    #[allow(clippy::redundant_clone)]
    #[allow(unused_must_use)]
    fn test_clone_is_implemented() {
        // Just needs to compile
        DelaunayTriangulation::<Point2<f64>>::new().clone();
    }
}