Voronoi tessellations (#47)

.gitignore

@@ -16,3 +16,4 @@ docs/src/triangulations/waste.jl
 docs/animation.jl
 test/animation.jl
 test/tassy.txt
+lloyd_animation.mp4

Project.toml

@@ -1,7 +1,7 @@
 name = "DelaunayTriangulation"
 uuid = "927a84f5-c5f4-47a5-9785-b46e178433df"
 authors = ["Daniel VandenHeuvel <[email protected]>"]
-version = "0.5.2"
+version = "0.6.0"
 
 [deps]
 DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"

README.md

@@ -18,6 +18,8 @@ This is a package for constructing Delaunay triangulations of planar point sets.
 - [Fully customisable interface](https://danielvandh.github.io/DelaunayTriangulation.jl/dev/interface/interface/) for defining geometric primitives.
 - [Simple iteration over mesh elements, including points, edges, or triangles](https://danielvandh.github.io/DelaunayTriangulation.jl/dev/data_structures/triangulation/).
 - Computation of [statistics](https://danielvandh.github.io/DelaunayTriangulation.jl/dev/data_structures/statistics/) over individual triangular elements and over a complete triangulation.
+- [Computation of Voronoi tessellations](https://danielvandh.github.io/DelaunayTriangulation.jl/dev/tessellations/voronoi.md), including [clipping of polygons to the convex hull]((https://danielvandh.github.io/DelaunayTriangulation.jl/dev/tessellations/clipped.md)). I hope to get this working for constrained triangulations, but it's difficult.
+- Computation of [centroidal Voronoi tessellations](https://danielvandh.github.io/DelaunayTriangulation.jl/dev/tessellations/lloyd.md) using Lloyd's algorithm.
 
 Much of the work in this package is derived from the book *Delaunay Mesh Generation* by Cheng, Dey, and Shewchuk (2013). 
 
@@ -27,5 +29,4 @@ Just as a nice demonstration of the incremental behaviour of the algorithms in t
 
 https://user-images.githubusercontent.com/95613936/232210266-615e08bd-d4f2-4a40-849e-e832778a4b71.mp4
 
-
-
+Here is also a nice animation showing the computation of a centroidal Voronoi tessellation.

docs/make.jl

@@ -24,14 +24,21 @@ makedocs(;
             "Plotting" => "triangulations/plotting.md",
             "Convex Polygons" => "triangulations/convex.md"
         ],
+        "Voronoi Tessellations" => [
+            "Voronoi Tessellations" => "tessellations/voronoi.md",
+            "Clipped Voronoi Tessellations" => "tessellations/clipped.md",
+            "Centroidal Voronoi Tessellation" => "tessellations/lloyd.md",
+            "Plotting" => "tessellations/plotting.md"
+        ],
         "Boundary Handling" => "boundary_handling.md",
         "Data Structures" => [
             "Adjacent" => "data_structures/adjacent.md",
             "Adjacent2Vertex" => "data_structures/adjacent2vertex.md",
             "Graph" => "data_structures/graph.md",
             "Convex Hull" => "data_structures/convex_hull.md",
             "Triangulation" => "data_structures/triangulation.md",
-            "Statistics" => "data_structures/statistics.md"
+            "Statistics" => "data_structures/statistics.md",
+            "Voronoi Tessellation" => "data_structures/voronoi.md"
         ],
         "Operations" => "operations.md",
         "Other Features" => [

docs/src/data_structures/statistics.md

@@ -77,6 +77,7 @@ Delaunay Triangulation Statistics.
 All the relevant docstrings for working with these structs are below.
 
 ```@docs 
+get_all_stat
 num_vertices(::TriangulationStatistics)
 num_solid_vertices(::TriangulationStatistics)
 num_ghost_vertices(::TriangulationStatistics)

docs/src/data_structures/voronoi.md

@@ -0,0 +1,104 @@
+```@meta
+CurrentModule = DelaunayTriangulation
+```
+
+# Voronoi Tessellation 
+
+The data structure for a Voronoi tessellation is reasonably simple. The data structure here is probably not as exhaustive as it could be, but it is sufficient.
+
+```@docs 
+VoronoiTessellation 
+```
+
+Each field has its own accessor:
+
+```@docs 
+get_triangulation(::VoronoiTessellation) 
+get_generators(::VoronoiTessellation) 
+get_polygon_points(::VoronoiTessellation)
+get_polygons(::VoronoiTessellation) 
+get_circumcenter_to_triangle(::VoronoiTessellation)
+get_triangle_to_circumcenter(::VoronoiTessellation) 
+get_unbounded_polygons(::VoronoiTessellation) 
+get_cocircular_circumcenters(::VoronoiTessellation) 
+get_adjacent(::VoronoiTessellation) 
+get_boundary_polygons(::VoronoiTessellation) 
+```
+
+There are several useful methods available for working with this data structure. We list some of these below; for functions that actually construct the tessellation, see the dedicated tessellation section in the sidebar.
+
+## Type queries 
+
+```@docs 
+edge_type(::VoronoiTessellation{Tr,P,I,T,S,E}) where {Tr,P,I,T,S,E}
+number_type(::VoronoiTessellation{Tr,P}) where {Tr,P}
+integer_type(::VoronoiTessellation{Tr,P,I}) where {Tr,P,I}
+triangle_type(::VoronoiTessellation{Tr,P,I,T}) where {Tr,P,I,T}
+```
+
+## Getters 
+
+```@docs 
+get_generator(::VoronoiTessellation, ::Any)
+get_polygon_point(::VoronoiTessellation, ::Any)
+get_polygon(::VoronoiTessellation, ::Any)
+get_circumcenter_to_triangle(::VoronoiTessellation, ::Any)
+get_triangle_to_circumcenter(::VoronoiTessellation, ::Any)
+get_polygon_coordinates 
+get_neighbouring_boundary_edges 
+```
+
+## Nums 
+
+```@docs 
+num_polygons 
+num_polygon_vertices 
+num_generators 
+``` 
+
+## Adders 
+
+```@docs
+add_polygon! 
+push_polygon_point! 
+add_unbounded_polygon! 
+delete_unbounded_polygon! 
+add_boundary_polygon! 
+```
+
+## Iterators 
+
+```@docs 
+each_generator 
+each_polygon_vertex 
+each_unbounded_polygon 
+each_polygon 
+each_polygon_index 
+```
+
+## Adjacent 
+
+```@docs 
+get_adjacent(::VoronoiTessellation, ::Any)
+add_adjacent!(::VoronoiTessellation, ::Any, ::Any)
+delete_adjacent!(::VoronoiTessellation, ::Any, ::Any)
+delete_polygon_adjacent!
+add_polygon_adjacent!
+```
+
+## Features 
+
+```@docs 
+polygon_features(::VoronoiTessellation, ::Any)
+get_area
+get_centroid 
+polygon_bounds(::VoronoiTessellation, ::Any)
+jump_and_march(::VoronoiTessellation, ::Any)
+```
+
+## Utilities 
+
+```@docs 
+get_surrounding_polygon(::VoronoiTessellation, ::Any)
+convert_to_boundary_edge(::VoronoiTessellation, :Any)
+```

docs/src/index.md

@@ -6,7 +6,9 @@ CurrentModule = DelaunayTriangulation
 
 Documentation for [DelaunayTriangulation](https://github.com/DanielVandH/DelaunayTriangulation.jl).
 
-This is a package for computing Delaunay triangulations of planar point sets. We support both unconstrained and constrained Delaunay triangulations, as well as mesh refinement with Rupper's algorithm. An interface for computing  constrained Delaunay triangulations with Gmsh is also available if needed; see the Gmsh discussion in the sidebar. Unconstrained Delaunay triangulations are computed with the Bowyer-Watson algorithm, and constrained Delaunay triangulations are computed with the incremental algorithm given by [https://doi.org/10.1016/j.comgeo.2015.04.006](https://doi.org/10.1016/j.comgeo.2015.04.006).
+This is a package for computing Delaunay triangulations of planar point sets. We support both unconstrained and constrained Delaunay triangulations, as well as mesh refinement with Rupper's algorithm. We also support Voronoi tessellations, clipped Voronoi tessellations, and central Voronoi tessellations; for these latter two cases, the triangulation is treated as constrained with the convex hull edges, but we do not support general boundary constraints for tessellations - if you do know about this, get in touch!
+
+An interface for computing  constrained Delaunay triangulations with Gmsh is also available if needed; see the Gmsh discussion in the sidebar. Unconstrained Delaunay triangulations are computed with the Bowyer-Watson algorithm, and constrained Delaunay triangulations are computed with the incremental algorithm given by [https://doi.org/10.1016/j.comgeo.2015.04.006](https://doi.org/10.1016/j.comgeo.2015.04.006).
 
 To ensure that the triangulations are robust to degeneracies, we use ExactPredicates.jl for all geometrical predicates. The results from these predicates are handled through a Certificate module, as outlined in the predicates section in the sidebar.
 

docs/src/interface/edges.md

@@ -19,6 +19,7 @@ terminal
 ```@docs 
 edge_indices 
 reverse_edge
+compare_unoriented_edge
 ```
 
 ## Collection of Edges 

docs/src/interface/points.md

@@ -33,12 +33,14 @@ num_points
 number_type 
 push_point!
 pop_point!
+set_point!
 ```
 
 ### Generic Methods 
 
 ```@docs 
 get_point 
 points_are_unique 
-lexicographic_order 
+lexicographic_order
+mean_points 
 ```

docs/src/other_features/pole_of_inaccessibility.md

@@ -17,10 +17,11 @@ Below we also list some other relevant docstrings.
 ```@docs 
 Cell 
 CellQueue 
-polygon_features 
+polygon_features(::Any, ::Any)
 squared_distance_to_segment 
 distance_to_polygon 
-polygon_bounds 
+polygon_bounds(::Any, ::Any)
+sort_convex_polygon!
 ```
 
 `distance_to_polygon` is also useful for point location, as shown in the examples below.

docs/src/tessellations/clipped.md

@@ -0,0 +1,174 @@
+```@meta
+CurrentModule = DelaunayTriangulation
+```
+
+# Clipped Voronoi Tessellations
+
+Often, it is useful to chop off the unbounded polygons in a tessellation, truncating them to some boundary. Usually this is just a box, but that's not so useful in general. We currently provide support for chopping to a convex hull, but the aim is to eventually support chopping to any boundary (the same algorithm we use should apply here, but there are some special cases that are really quite difficult). 
+
+The same function `voronoi` is used for this, just setting the second argument to be `true`. The algorithm is described at the end of this section.
+
+## Example
+
+Let us give an example.
+
+```julia
+using DelaunayTriangulation, CairoMakie
+pts = randn(2, 250)
+tri = triangulate(pts)
+vorn = voronoi(tri)
+vorn_clip = voronoi(tri, true)
+cmap = Makie.cgrad(:jet)
+colors = get_polygon_colors(vorn, cmap)
+fig = Figure()
+ax = Axis(fig[1, 1], aspect=1)
+voronoiplot!(ax, vorn, strokecolor=:red, strokewidth=0.2, polygon_color=colors, show_generators=false)
+triplot!(ax, tri, strokewidth=0.0, strokecolor=(:black, 0.4), show_all_points=false)
+xlims!(ax, -5, 5)
+ylims!(ax, -5, 5)
+ax = Axis(fig[1, 2], aspect=1)
+voronoiplot!(ax, vorn_clip, strokecolor=:red, strokewidth=0.2, polygon_color=colors, show_generators=false)
+triplot!(ax, tri, strokewidth=0.0, strokecolor=(:black, 0.4), show_all_points=false)
+xlims!(ax, -5, 5)
+ylims!(ax, -5, 5)
+```
+
+```@raw html
+<figure>
+    <img src='../figs/bounded.png', alt='Clipped Voronoi Tessellation'><br>
+</figure>
+```
+
+As we can see, all the polygons have now been chopped so that the entire tessellation fits into the original boundary of the dual triangulation. Also, the `unbounded_polygons` and `boundary_polygons` fields have been updated:
+
+```julia-repl
+julia> DelaunayTriangulation.get_unbounded_polygons(vorn_clip)
+Set{Int64}()
+
+julia> DelaunayTriangulation.get_boundary_polygons(vorn_clip)
+Set{Int64} with 31 elements:
+  169
+  56
+  200
+  195
+  72
+  180
+  221
+  8
+  37
+  187
+  32
+  6
+  171
+  190
+  69
+  219
+  73
+  ⋮
+```
+
+## Algorithm
+
+Since the code for clipping is quite involved, most likely more involved than it should be, it is probably worthwhile describing the actual implementation here. I want to one day refine this approach, and allow it to work on any domain. The approach is mostly based on the ideas in the paper "Efficient Computation of Clipped Voronoi Diagram for Mesh Generation"
+by Yan, Wang, Lévy, and Liu. I imagine if I had access to their code, everything would be a lot nicer (that seems to be true in many computational geometry papers).
+
+The function that performs the clipping is `clip_voronoi_tessellation!`:
+
+```@docs 
+clip_voronoi_tessellation!
+```
+
+This function starts off by finding all the intersections and then adding them into the point set, then using these new points to clip the polygons. We describe these steps below.
+
+### Finding all intersections 
+
+The function `find_all_intersections` finds the intersections of the polygons with the boundary.
+
+```@docs 
+find_all_intersections
+initialise_clipping_arrays
+```
+
+The idea is to use a first-in first-out (FIFO) queue to prioritise the edges and polygons to consider. In particular, we pick some random edge and then use its midpoint to find what polygon that midpoint is in (via `jump_and_march`). This edge and the polygon are then added into the queue to be processed. Starting with this edge-polygon pair, we keep looking for intersections until we have processed all boundary edges and there are no more pairs in the queue. The function that does this processing of queued pairs is `dequeue_and_process!`:
+
+```@docs 
+dequeue_and_process!
+enqueue_new_edge!
+```
+
+Essentially, this function first checks if the pair has already been processed and, if so, returns. Otherwise, the first step is to go into `process_polygon!` to find the intersections of the polygons with the nearby boundary:
+
+```@docs 
+process_polygon!
+```
+
+The goal of `process_polygon!` is to iterate over the edges of the current polygon to look for intersections, searching (1) the current edge, (2) the edge to the left of the current edge, and (3) the edge to the right of the current edge. Suppose we have an edge $e_{uv}$ of the polygon being considered. There are four possibilities:
+
+1. First, both $u$ and $v$ could correspond to ghost vertices, meaning no intersection with the boundary is possible.
+2. Alternatively, $e_{uv}$ could be a ray going inwards, meaning $u$ is a ghost vertex while $v$ corresponds to an actual circumcenter. In this case, `process_ray_intersection!` is used to look for an intersection of $e_{uv}$ with the edge of the ghost triangle corresponding to $u$. If an intersection is found, we do need to be careful of rays that intersect multiple boundary edges (as is common when we have very few triangles in the underlying triangulation). This check is done with `process_ray_intersection_with_other_edges!`, which simply looks at the left and right edges along with the current edge. 
+3. Same as the above, except $e_{uv}$ could be a ray going inwards, meaning $u$ now corresponds to an actual cicrcumcenter while $v$ is a ghost vertex.
+4. Lastly, $e_{uv}$ could be a fniite segment, in which case we can use simple intersection formula to test for intersections with the left and right edges along with the current edge. This check is done via `process_segment_intersection!`.
+
+In these steps, when there are no intersections we also check if it was because a segment lies entirely outside of the domain. We keep track of these points as these will need to be deleted later.
+
+Some of the relevant docstrings for working with `process_polygon!` are below.
+
+```@docs 
+is_segment_between_two_ghosts
+is_ray_going_in
+is_ray_going_out
+is_finite_segment
+process_segment_intersection!
+process_ray_intersection!
+add_to_intersected_edge_cache!
+add_segment_intersection!
+segment_intersection_coordinates
+intersection_of_edge_and_bisector_ray
+classify_and_compute_segment_intersection
+process_ray_intersection_with_other_edges!
+```
+
+Once we have gone through `process_polygon!`, we need to assign the identified intersections to the current edge, the left edge, or to the right edge. This is done with `classify_intersections!`.
+
+```@docs 
+classify_intersections!
+```
+
+Once this is done, we need to actually consider the intersection points. The two goals here are: (1) Check if a corner point inside the polygon has to be added, and (2) look for other incident polygons to process. This processing is done via `process_intersection_points!`:
+
+```@docs 
+process_intersection_points!
+```
+
+Basically, this function first looks at the left and right intersections. If we have intersections on the left and on the current edge, then we will also have the vertex shared by the two edges included as a point we need to chop to, unless the current polygon being considered corresponds to a generator not on a boundary. We check this for each intersection and update the intersections with `add_segment_intersection!` accordingly. In such a case, we add the current edge, together with the indices of the left edge (corresponding to polygons) to the queue, ready for the next iteration. The same is done for the right edge. After handling this case, we check all the intersections without worry for corner points. For each intersection, we take the polygon on the other side of the intersecting segment and add it to the queue together with the boundary ede that was intersected. 
+
+We do all the above steps until we run out of edges to process and the queue is empty, adding all the intersection points into the point list using `add_intersection_points!`:
+
+```@docs 
+add_intersection_points!
+```
+
+### Clipping the polygons
+
+The next step is to clip the polygons to the boundary using the computed intersections. This is handled via `clip_all_polygons!`:
+
+```@docs 
+clip_all_polygons!
+```
+
+This function iterates over all the identified boundary polygons and calls `clip_polygon!` on each:
+
+```@docs 
+clip_polygon!
+```
+
+Consider a single polygon. The `clip_polygon!` function starts by deleting the adjacencies for this polygon, ready for updating later. We then remove all ghost vertices from the vertex list, and all circumcenters identified from the previous step that lie outside of the domain. We then add all the vertices. With this processing, the vertex list is no longer sorted, so `sort_convex_polygon!` is called to get a counter-clockwise representation. Adding in the adjacencies via `add_polygon_adjacent!`, this completes the clipping of this polygon. 
+
+With all polygons processed, we add in all that we processed into the `boundary_polygons` field via `add_all_boundary_polygons!`:
+
+```@docs 
+add_all_boundary_polygons!
+```
+
+With this last step, the clipping is complete.
+