Implement weighted Delaunay triangulations and power diagrams (#180)

NEWS.md

@@ -1,5 +1,18 @@
 # Changelog
 
+## 1.3.0
+
+This release finally introduces weighted triangulations and power diagrams, and also allows for users to provide a generic convex polygon to for clipping a Voronoi tessellation instead of only the convex hull.
+
+- Weighted triangulations have now been implemented, as have power diagrams. The weights are also no longer restricted to `Float64` type. See [#180](https://github.com/JuliaGeometry/DelaunayTriangulation.jl/pull/180).
+- `intersection_of_edge_and_bisector_ray` now accepts a `project` keyword argument. See [#180](https://github.com/JuliaGeometry/DelaunayTriangulation.jl/pull/180).
+- `get_weight(w, i)` now returns, when `i` is not an integer, either `i[3]` if it represents a point in space or `0`. See [#180](https://github.com/JuliaGeometry/DelaunayTriangulation.jl/pull/180).
+- Define `project_onto_line(p, q, r)` for projecting a point `r` onto the line defined by `p` and `q`. See [#180](https://github.com/JuliaGeometry/DelaunayTriangulation.jl/pull/180).
+- Fixed a bug with clipping Voronoi tessellations in cases where there are no intersections of any Voronoi polygon with the convex hull. See [#180](https://github.com/JuliaGeometry/DelaunayTriangulation.jl/pull/180).
+- `voronoi` now accepts an optional `clip_polygon` keyword argument, defaulting to `nothing` (corresponding to the convex hull), allowing for a convex clip polygon to be used instead of the convex hull. The `clip_polygon` should be a `Tuple` of the form `(points, boundary_nodes)` where the `boundary_nodes` give vertices of `points` adhering to the usual convention. Note that this could be used as an alternative to looping over `get_polygon_coordinates` for clipping to a rectangle. See [#180](https://github.com/JuliaGeometry/DelaunayTriangulation.jl/pull/180).
+- `centroidal_smooth` now accepts `clip_points` and `clip_vertices` as keyword arguments, defaulting to `nothing` (corresponding to the convex hull), to accommodate the new `clip_polygon` keyword argument in `voronoi`. See [#180](https://github.com/JuliaGeometry/DelaunayTriangulation.jl/pull/180).
+- `has_multiple_curves`, `has_multiple_sections`, and `num_boundary_edges` now have methods for `Tuple`s of integers. A bug was also fixed with `number_type` of a `Tuple` of `Tuple`s of coordinates returning the `Tuple` type instead of the coordinate type. See [#180](https://github.com/JuliaGeometry/DelaunayTriangulation.jl/pull/180).
+
 ## v1.2.0
 
 - Warnings are now thrown when you try and triangulate point sets not in the plane. The `is_planar` function has been introduced for this. See [#178](https://github.com/JuliaGeometry/DelaunayTriangulation.jl/pull/178).

Project.toml

@@ -1,7 +1,7 @@
 name = "DelaunayTriangulation"
 uuid = "927a84f5-c5f4-47a5-9785-b46e178433df"
 authors = ["Daniel VandenHeuvel <[email protected]>"]
-version = "1.2.0"
+version = "1.3.0"
 
 [deps]
 AdaptivePredicates = "35492f91-a3bd-45ad-95db-fcad7dcfedb7"

README.md

@@ -6,7 +6,7 @@
 [![DOI](https://zenodo.org/badge/540660309.svg)](https://zenodo.org/badge/latestdoi/540660309)
 [![Aqua QA](https://raw.githubusercontent.com/JuliaTesting/Aqua.jl/master/badge.svg)](https://github.com/JuliaTesting/Aqua.jl)
 
-This is a package for constructing Delaunay triangulations and Voronoi tessellations of planar point sets. Supports unconstrained and constrained triangulations, mesh refinement, triangulation of curve bounded domains, Voronoi tessellations, and clipped and centroidal Voronoi tessellations. To install the package, do
+This is a package for constructing Delaunay triangulations and Voronoi tessellations of planar point sets. Supports unconstrained and constrained triangulations, weighted triangulations, mesh refinement, triangulation of curve bounded domains, Voronoi tessellations, power diagrams, and clipped and centroidal Voronoi tessellations. To install the package, do
 
 ```julia
 julia>] add DelaunayTriangulation
@@ -15,7 +15,9 @@ julia>] add DelaunayTriangulation
 Many features are available, some of these being:
 
 - Unconstrained and constrained Delaunay triangulations, supporting many types of domains.
-- Computation of Voronoi tessellations, clipped Voronoi tessellations where the Voronoi tiles get clipped to the convex hull, and centroidal Voronoi tessellations where each Voronoi tile's generator is the tile's centroid.
+- Weighted Delaunay triangulations.
+- Computation of Voronoi tessellations, clipped Voronoi tessellations where the Voronoi tiles get clipped to a convex polygon, and centroidal Voronoi tessellations where each Voronoi tile's generator is the tile's centroid.
+- Power diagrams.
 - Mesh refinement, with support for custom angle and area constraints, as well as refinement of curve-bounded domains.
 - Dynamic point insertion, point deletion, and segment insertion, amongst many other operations.
 - Computation of convex hulls.
@@ -49,7 +51,7 @@ vorn2 = voronoi(tri1)
 vorn3 = voronoi(tri1, clip=true, predicates = ExactKernel())
 
 # Smoothed Voronoi 
-vorn4 = centroidal_smooth(vorn3, predicates = FastKernel())
+vorn4 = centroidal_smooth(vorn3; predicates = FastKernel()) # or do voronoi(tri1, clip = true, smooth = true)
 
 # Constrained example with refinement 
 boundary_points = [(0.0, 0.0), (1.0, 0.0), (1.0, 0.3), (0.5, 0.3),
@@ -88,16 +90,36 @@ boundary_nodes = [[square], [[ellipse]], [[catmull_spl]], [[circle]], [[circle2]
 tri8 = triangulate(points; boundary_nodes)
 refine!(tri8; max_area=1e-2get_area(tri8)) # could also use find_polygon to help define a custom refinement function for each shape
 
+# Weighted triangulation example
+points = tuple.(rand(20), rand(20))
+weights = 3randn(20)
+tri9 = triangulate(points; weights)
+
+# Power diagram example 
+vorn10 = voronoi(tri9) # can also use clip/smooth here 
+
+# Clipped Voronoi example with a generic convex polygon
+points = 10randn(2, 100)
+weights = rand(100)
+circ = CircularArc((0.0, 2.0), (0.0, 2.0), (0.0, 0.0)) # clip to a circle 
+clip_points = [circ(t) for t in LinRange(0, 1, 100)]
+clip_vertices = [1:(length(clip_points)-1); 1]
+tri11 = triangulate(points; weights)
+vorn11 = voronoi(tri11, clip=true, clip_polygon=(clip_points, clip_vertices))
+
 # Plotting 
-fig = Figure(fontsize = 42, size = (2800, 1480))
-ax = Axis(fig[1, 1], title="Unconstrained", width=600,height=600);            triplot!(ax, tri1)
-ax = Axis(fig[1, 2], title="Voronoi", width=600,height=600);                  voronoiplot!(ax, vorn2)
-ax = Axis(fig[1, 3], title="Clipped Voronoi", width=600,height=600);          voronoiplot!(ax, vorn3)
-ax = Axis(fig[1, 4], title="Centroidal Voronoi", width=600,height=600);       voronoiplot!(ax, vorn4)
-ax = Axis(fig[2, 1], title="Constrained", width=600,height=600);              triplot!(ax, tri5)
-ax = Axis(fig[2, 2], title="Disjoint Constrained", width=600,height=600);     triplot!(ax, tri6)
-ax = Axis(fig[2, 3], title="Curve-Bounded", width=600,height=600);            triplot!(ax, tri7)
-ax = Axis(fig[2, 4], title="Disjoint Curve-Bounded", width=600,height=600);   triplot!(ax, tri8)
+fig = Figure(fontsize = 42, size = (2800, 2200))
+ax = Axis(fig[1, 1], title = "Unconstrained", width = 600, height = 600);            triplot!(ax, tri1)
+ax = Axis(fig[1, 2], title = "Voronoi", width = 600, height = 600);                  voronoiplot!(ax, vorn2)
+ax = Axis(fig[1, 3], title = "Clipped Voronoi", width = 600, height = 600);          voronoiplot!(ax, vorn3)
+ax = Axis(fig[1, 4], title = "Centroidal Voronoi", width = 600, height = 600);       voronoiplot!(ax, vorn4)
+ax = Axis(fig[2, 1], title = "Constrained", width = 600, height = 600);              triplot!(ax, tri5)
+ax = Axis(fig[2, 2], title = "Disjoint Constrained", width = 600, height = 600);     triplot!(ax, tri6)
+ax = Axis(fig[2, 3], title = "Curve-Bounded", width = 600, height = 600);            triplot!(ax, tri7)
+ax = Axis(fig[2, 4], title = "Disjoint Curve-Bounded", width = 600, height = 600);   triplot!(ax, tri8)
+ax = Axis(fig[3, 1], title = "Weighted", width = 600, height = 600);                 triplot!(ax, tri9)
+ax = Axis(fig[3, 2], title = "Power Diagram", width = 600, height = 600);            voronoiplot!(ax, vorn10)
+ax = Axis(fig[3, 3], title = "Clipped Voronoi", width = 600, height = 600);          voronoiplot!(ax, vorn11)
 ```
 
 ![](readme.png)

docs/Project.toml

@@ -6,6 +6,7 @@ DelaunayTriangulation = "927a84f5-c5f4-47a5-9785-b46e178433df"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 DocumenterTools = "35a29f4d-8980-5a13-9543-d66fff28ecb8"
 ElasticArrays = "fdbdab4c-e67f-52f5-8c3f-e7b388dad3d4"
+GeometryBasics = "5c1252a2-5f33-56bf-86c9-59e7332b4326"
 InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"

docs/make.jl

@@ -126,12 +126,15 @@ const _PAGES = [
         "Triangulating Rectangular Regions" => "tutorials/lattice.md",
         "Triangulating Convex Polygons" => "tutorials/convex.md",
         "Triangulating Curve-Bounded Domains" => "tutorials/curve_bounded.md",
+        "Weighted Triangulations" => "tutorials/weighted.md",
         "Voronoi Tessellations" => "tutorials/voronoi.md",
         "Clipped Voronoi Tessellations" => [
             "Clipping to the Convex Hull" => "tutorials/clipped.md",
             "Clipping to a Rectangle" => "tutorials/clipped_rectangle.md",
+            "Clipping to a Generic Convex Polygon" => "tutorials/clipped_polygon.md",
         ],
         "Centroidal Voronoi Tessellations" => "tutorials/centroidal.md",
+        "Power Diagrams" => "tutorials/power.md",
         "Point Location" => "tutorials/point_location.md",
         "Nearest Neighbour Queries" => "tutorials/nearest.md",
         "Convex Hulls" => "tutorials/convex_hull.md",
@@ -179,9 +182,11 @@ const _PAGES = [
         "Mesh Refinement" => "math/refinement.md",
         "Curves" => "math/curves.md",
         "Triangulating Curve-Bounded Domains" => "math/curve_bounded.md",
+        "Weighted Delaunay Triangulations" => "math/weighted.md",
         "Voronoi Tessellations" => "math/voronoi.md",
         "Clipped Voronoi Tessellations" => "math/clipped.md",
         "Centroidal Voronoi Tessellations" => "math/centroidal.md",
+        "Power Diagrams" => "math/power.md",
     ],
     "Example Applications" => [
         "Overview" => "applications/overview.md",

docs/src/api/statistics.md

@@ -19,6 +19,8 @@ triangle_circumcenter
 triangle_offcenter
 triangle_edge_midpoints
 triangle_sink
+triangle_orthocenter
+triangle_orthoradius_squared
 TriangulationStatistics
 statistics
 num_boundary_segments

docs/src/extended/algorithms.md

@@ -89,8 +89,6 @@ Pages = ["src/data_structures/triangulation/methods/segments.jl"]
 
 ## Weighted Triangulations
 
-While weighted triangulations are not yet completely implemented in the package, there are some functions related to them in preparation.
-
 ```@autodocs 
 Modules = [DelaunayTriangulation]
 Pages = ["src/data_structures/triangulation/methods/weights.jl"]

docs/src/index.md

@@ -9,7 +9,9 @@ This is the documentation for DelaunayTriangulation.jl. [Click here to go back t
 This is a package for computing Delaunay triangulations and Voronoi tessellations of points in two dimensions, amongst many other features:
 
 - Unconstrained and constrained Delaunay triangulations, supporting many types of domains.
-- Computation of Voronoi tessellations, clipped Voronoi tessellations where the Voronoi tiles get clipped to the convex hull, and centroidal Voronoi tessellations where each Voronoi tile's generator is the tile's centroid.
+- Weighted Delaunay triangulations.
+- Computation of Voronoi tessellations, clipped Voronoi tessellations where the Voronoi tiles get clipped to a convex polygon, and centroidal Voronoi tessellations where each Voronoi tile's generator is the tile's centroid.
+- Power diagrams.
 - Mesh refinement, with support for custom angle and area constraints, as well as refinement of curve-bounded domains.
 - Dynamic point insertion, point deletion, and segment insertion, amongst many other operations.
 - Computation of convex hulls.

docs/src/literate_tutorials/clipped.jl

@@ -4,8 +4,8 @@
 
 # One issue that may arise when dealing with Voronoi tessellations is the 
 # presence of unbounded polygons occurring on the boundary. One way to deal with this 
-# is to clip polygons to the convex hull of the tessellation. (Arbitrary clipping boundaries 
-# are on the to-do list, but they are not yet implemented.)
+# is to clip polygons to the convex hull of the tessellation. We describe how to also clip
+# the tessellation to a generic convex polygon, instead of just the convex hull, in [this tutorial](clipped_polygon.md).
 
 # In the example below, we clip the tessellation to the convex hull of the point set by using `clip=true`
 # in the keyword arguments.
@@ -36,4 +36,4 @@ fig
 @test_reference joinpath(fig_path, "voronoi_ex_2.png") fig #src
 
 # As you can see, the unbounded polygons, and any polygons that included points 
-# outside of the convex hull, have now been clipped to the convex hull.
+# outside of the convex hull, have now been clipped to the convex hull.

docs/src/literate_tutorials/clipped_polygon.jl

@@ -0,0 +1,73 @@
+# # Clipped Voronoi Tessellations 
+# ## Clipping to a Generic Convex Polygon
+#
+
+# In this tutorial we show how to clip a Voronoi tessellation to 
+# more generic convex polygons (non-convex polygons are not currently supported) than just a convex hull or a rectangle. This is
+# done by using the `clip_polygon` keyword argument in `voronoi`.
+
+# We start by clipping the tessellation to a rectangle, showing an alternative
+# to the [previous tutorial](clipped_rectangle.md). To start, we load in the packages 
+# we need and generate some data.
+using DelaunayTriangulation
+using DelaunayTriangulation: EllipticalArc
+using CairoMakie
+using StableRNGs
+using ReferenceTests #src
+using Test #src
+fig_path = joinpath(@__DIR__, "../figures") #src
+
+rng = StableRNG(123)
+points = randn(rng, 2, 50)
+tri = triangulate(points; rng)
+vorn = voronoi(tri)
+
+# To define the polygon, we define the points and vertices just as we would, for example,
+# the boundary of a triangulation.
+xmin, xmax, ymin, ymax = -1/2, 1/2, -1.0, 1.0
+clip_points = ((xmin, ymin), (xmax, ymin), (xmax, ymax), (xmin, ymax))
+clip_vertices = (1, 2, 3, 4, 1)
+clip_polygon = (clip_points, clip_vertices)
+
+# Now we simply pass the polygon into `voronoi`.
+clipped_vorn = voronoi(tri, clip = true, clip_polygon = clip_polygon)
+
+# Now let's look at the results.
+fig = Figure()
+ax1 = Axis(fig[1, 1], title = "Unclipped", width = 600, height = 400)
+ax2 = Axis(fig[1, 2], title = "Clipped", width = 600, height = 400)
+voronoiplot!(ax1, vorn, show_generators = false, colormap = :matter, strokewidth = 4)
+xlims!(ax1, -2, 2)
+ylims!(ax1, -2, 2)
+lines!(ax1, [clip_points..., clip_points[begin]], color = :black, linewidth = 4, linestyle = :dash)
+voronoiplot!(ax2, clipped_vorn, show_generators = false, colormap = :matter, strokewidth = 4)
+xlims!(ax2, -2, 2)
+ylims!(ax2, -2, 2)
+resize_to_layout!(fig)
+fig
+@test_reference joinpath(fig_path, "generic_clipped_voronoi_ex_1.png") fig #src
+
+# We can clip to any convex polygon that we want to. For example, below we clip to an elliptical boundary.
+rng = StableRNG(123333)
+points = randn(rng, 2, 50)
+tri = triangulate(points; rng)
+vorn = voronoi(tri)
+ellip = EllipticalArc((1/2, 0.0), (1/2, 0.0), (0.0, 0.0), 1/2, 1.0, 0.0)
+t = LinRange(0, 1, 50)
+clip_points = ellip.(t)
+clip_vertices = [1:(length(clip_points)-1); 1]
+clip_polygon = (clip_points, clip_vertices)
+clipped_vorn = voronoi(tri, clip = true, clip_polygon = clip_polygon)
+fig = Figure()
+ax1 = Axis(fig[1, 1], title = "Unclipped", width = 600, height = 400)
+ax2 = Axis(fig[1, 2], title = "Clipped", width = 600, height = 400)
+voronoiplot!(ax1, vorn, show_generators = false, colormap = :matter, strokewidth = 4)
+xlims!(ax1, -2, 2)
+ylims!(ax1, -2, 2)
+lines!(ax1, [clip_points..., clip_points[begin]], color = :black, linewidth = 4, linestyle = :dash)
+voronoiplot!(ax2, clipped_vorn, show_generators = false, colormap = :matter, strokewidth = 4)
+xlims!(ax2, -2, 2)
+ylims!(ax2, -2, 2)
+resize_to_layout!(fig)
+fig
+@test_reference joinpath(fig_path, "generic_clipped_voronoi_ex_2.png") fig #src

docs/src/literate_tutorials/clipped_rectangle.jl

@@ -5,9 +5,9 @@
 # In the previous tutorial, we demonstrated how we can clip to
 # the convex hull of the point set. However, it is often useful to clip
 # to a rectangle, for example if you want to clip to a region of interest
-# in a simulation. We do not yet support this within [`voronoi`](@ref) itself, 
-# but we provide the function [`get_polygon_coordinates`](@ref) for this (this is what 
-# `voronoiplot` uses to plot inside a bounding box).
+# in a simulation. Here we obtain the coordinates just by looping over all 
+# the polygons, but you can also use this by providing a `clip_polygon` into 
+# `voronoi`, as described in [this tutorial](clipped_polygon.md).
 
 # Let us now demonstrate. First, we construct a tessellation of 
 # some example point set.

docs/src/literate_tutorials/curve_bounded.jl

@@ -188,7 +188,7 @@ tri = triangulate(copy(points); boundary_nodes = curve, rng) # copying so that w
 refine!(tri; max_area = 1.0e-2)
 fig, ax, sc = triplot(tri)
 fig
-@test_reference joinpath(fig_path, "triangulate_curve_bounded_ex_4.png") fig by = psnr_equality(15) #src
+@test_reference joinpath(fig_path, "triangulate_curve_bounded_ex_4.png") fig by = psnr_equality(7) #src
 
 # ### Using Custom Constraints to Control Refinement 
 # Let's give another example of using custom constraints to better control the refinement within different domains. Referencing the 
@@ -221,7 +221,7 @@ tri = triangulate(points; boundary_nodes = curve, rng)
 refine!(tri; custom_constraint = poly_constraint, rng)
 fig, ax, sc = triplot(tri)
 fig
-@test_reference joinpath(fig_path, "triangulate_curve_bounded_ex_5.png") fig #src
+@test_reference joinpath(fig_path, "triangulate_curve_bounded_ex_5.png") fig by=psnr_equality(7) #src
 
 # ## Defining a New Parametric Curve 
 # Let us now give an example where we define a domain by a parametric curve that is not provided natively by this package. For this 

docs/src/literate_tutorials/nearest.jl

@@ -7,7 +7,11 @@
 # finding a nearest neighbour is the same as a point location problem, meaning 
 # given a point `p` find the Voronoi tile `P` containing it. Here we give an example 
 # of how we can use triangulations or tessellations to find the nearest neighbour 
-# in the point set to a given point. First, we load in the packages we need.
+# in the point set to a given point. We note that these same ideas 
+# could be applied for power diagrams, except that the metric used for defining distances
+# is based on the power distance instead of the Euclidean distance (see the [power diagram tutorial](power.md)
+# for more details); this is not demonstrated
+# in this tutorial. First, we load in the packages we need.
 using DelaunayTriangulation
 using CairoMakie
 using ReferenceTests #src