Skip to content

Commit 543c4c8

Browse files
asinghvi17asinghvi17
authored
2D Voronoi tessellation from DelaunayTriangulation.jl (#355)
2 parents 7650849 + 2ede894 commit 543c4c8

File tree

4 files changed

+344
-0
lines changed

4 files changed

+344
-0
lines changed

‎src/GeometryOps.jl‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -33,6 +33,7 @@ import SortTileRecursiveTree
3333
import SortTileRecursiveTree: STRtree
3434

3535
const GI = GeoInterface
36+
const DelTri = DelaunayTriangulation
3637

3738
const TuplePoint{T} = Tuple{T, T} where T <: AbstractFloat
3839
const Edge{T} = Tuple{TuplePoint{T},TuplePoint{T}} where T
@@ -85,6 +86,7 @@ include("methods/geom_relations/within.jl")
8586
include("methods/geom_relations/common.jl")
8687
include("methods/orientation.jl")
8788
include("methods/polygonize.jl")
89+
include("methods/voronoi.jl")
8890

8991
include("transformations/extent.jl")
9092
include("transformations/flip.jl")

‎src/methods/voronoi.jl‎

Lines changed: 210 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,210 @@
1+
#=
2+
# Voronoi Tessellation
3+
4+
The [_Voronoi tessellation_](https://en.wikipedia.org/wiki/Voronoi_diagram) of a set of points is a partitioning of the plane into regions based on distance to points.
5+
Each region contains all points closer to one generator point than to any other.
6+
7+
GeometryOps.jl provides a method for computing the Voronoi tessellation of a set of points,
8+
using the [DelaunayTriangulation.jl](https://github.com/JuliaGeometry/DelaunayTriangulation.jl) package.
9+
10+
Right now, the GeometryOps.jl method can only provide clipped voronoi tesselations, as the function returns a list of GeoInterface polygons.
11+
If you need an unbounded tessellation, open an issue and we can discuss the best way to represent unbounded polygons within GeometryOps.
12+
13+
## Example
14+
15+
### Simple tessellation
16+
```@example simple
17+
import GeometryOps as GO, GeoInterface as GI
18+
using CairoMakie # to plot
19+
20+
points = tuple.(randn(20), randn(20))
21+
polygons = GO.voronoi(points)
22+
f, a, p = plot(polygons[1]; label = "Voronoi cell 1")
23+
for (i, poly) in enumerate(polygons[2:end])
24+
plot!(a, poly; label = "Voronoi cell $(i+1)")
25+
end
26+
scatter!(a, points; color = :black, markersize = 10, label = "Generators")
27+
axislegend(a)
28+
f
29+
```
30+
31+
## Implementation
32+
33+
This implementation mainly just preforms some assertion checks before passing the Arguments
34+
to the DelaunayTriangulation package. We always set the argument `clip` to the DelaunayTriangulation
35+
`voronoi` call to `True` such that we can return a list of valid polygons. The default clipping polygon
36+
is the convex hull of the tessleation, but the user can pass in a bounding polygon with the `clip_polygon`
37+
argument. After the call to `voronoi`, the call then unpacks the voronoi output into GeoInterface
38+
polygons, whose point match the float type input by the user.
39+
=#
40+
41+
struct __NoCRSProvided end
42+
43+
"""
44+
voronoi(geometries, [T = Float64]; clip_polygon = nothing, kwargs...)
45+
46+
Compute the Voronoi tessellation of the points in `geometries`.
47+
Returns a vector of `GI.Polygon` objects representing the Voronoi cells,
48+
in the same order as the input points.
49+
50+
## Arguments
51+
- `geometries`: Any GeoInterface-compatible geometry or collection of geometries
52+
that can be decomposed into points
53+
- `T`: Float-type for returned polygons points (default: Float64)
54+
55+
## Keyword Arguments
56+
- `clip_polygon`: what bounding shape should the Voronoi cells be clipped to? (default: nothing -> clipped to the convex hull)
57+
clip_polygon can of several types: (1) a GeoInterface polygon, (2) a two-element tuple where the first element is a list of tuple points
58+
and the second element is a list of integer indices to indicate the order of the provided points, or (3) a a two-element tuple where the
59+
first element is a tuple of tuple points and the second element is a tuple of integer indices to indicate the order of the provided points
60+
- $CRS_KEYWORD
61+
- `rng`: random number generator to generating the voronoi tesselation
62+
63+
!!! warning
64+
This interface only computes the 2-dimensional Voronoi tessellation!
65+
Only clipped voronoi tesselations can be created!
66+
Only `T = Float64` or `Float32` are guaranteed good results by the underlying package DelaunayTriangulation.
67+
68+
!!! note
69+
The polygons are returned in the same order as the input points after flattening.
70+
Each polygon corresponds to the Voronoi cell of the point at the same index.
71+
72+
## Examples
73+
An example with default clipping to the convex hull.
74+
75+
```jldoctest voronoi
76+
import GeometryOps as GO
77+
import GeoInterface as GI
78+
using Random
79+
80+
rng = Xoshiro(0)
81+
points = [(rand(rng), rand(rng)) .* 5 for i in range(1, 3)]
82+
GO.voronoi(points; rng = rng)
83+
# output
84+
3-element Vector{GeoInterface.Wrappers.Polygon{false, false, Vector{GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Nothing, Nothing}}, Nothing, Nothing}}:
85+
GeoInterface.Wrappers.Polygon{false, false}([GeoInterface.Wrappers.LinearRing([(4.310704285977424, 0.42985432929210976), … (2) … , (4.310704285977424, 0.42985432929210976)])])
86+
GeoInterface.Wrappers.Polygon{false, false}([GeoInterface.Wrappers.LinearRing([(3.7949144210695653, 0.4101636087384888), … (4) … , (3.7949144210695653, 0.4101636087384888)])])
87+
GeoInterface.Wrappers.Polygon{false, false}([GeoInterface.Wrappers.LinearRing([(2.685897788908803, 0.3678259474564151), … (2) … , (2.685897788908803, 0.3678259474564151)])])
88+
```
89+
90+
An example with clipping to a GeoInterface polygon.
91+
```jldoctest voronoi
92+
clip_points = ((0.0,0.0), (5.0,0.0), (5.0,5.0), (0.0,5.0), (0.0,0.0))
93+
clip_order = (1, 2, 3, 4, 1)
94+
clip_poly1 = GI.Polygon([collect(clip_points)])
95+
GO.voronoi(points; clip_polygon = clip_poly1, rng = rng)
96+
# output
97+
3-element Vector{GeoInterface.Wrappers.Polygon{false, false, Vector{GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Nothing, Nothing}}, Nothing, Nothing}}:
98+
GeoInterface.Wrappers.Polygon{false, false}([GeoInterface.Wrappers.LinearRing([(5.0, 0.0), … (3) … , (5.0, 0.0)])])
99+
GeoInterface.Wrappers.Polygon{false, false}([GeoInterface.Wrappers.LinearRing([(3.7328227614527916, 0.0), … (3) … , (3.7328227614527916, 0.0)])])
100+
GeoInterface.Wrappers.Polygon{false, false}([GeoInterface.Wrappers.LinearRing([(0.0, 5.0), … (3) … , (0.0, 5.0)])])
101+
```
102+
103+
An example with clipping to a tuple of tuples.
104+
```jldoctest voronoi
105+
clip_poly2 = (clip_points, clip_order) # tuples
106+
GO.voronoi(points; clip_polygon = clip_poly2, rng = rng)
107+
# output
108+
3-element Vector{GeoInterface.Wrappers.Polygon{false, false, Vector{GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Nothing, Nothing}}, Nothing, Nothing}}:
109+
GeoInterface.Wrappers.Polygon{false, false}([GeoInterface.Wrappers.LinearRing([(5.0, 0.0), … (3) … , (5.0, 0.0)])])
110+
GeoInterface.Wrappers.Polygon{false, false}([GeoInterface.Wrappers.LinearRing([(3.7328227614527916, 0.0), … (3) … , (3.7328227614527916, 0.0)])])
111+
GeoInterface.Wrappers.Polygon{false, false}([GeoInterface.Wrappers.LinearRing([(0.0, 5.0), … (3) … , (0.0, 5.0)])])
112+
```
113+
114+
An example with clipping to a tuple of vectors.
115+
```jldoctest voronoi
116+
clip_poly3 = (collect(clip_points), collect(clip_order)) # vectors
117+
GO.voronoi(points; clip_polygon = clip_poly3, rng = rng)
118+
# output
119+
3-element Vector{GeoInterface.Wrappers.Polygon{false, false, Vector{GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Nothing, Nothing}}, Nothing, Nothing}}:
120+
GeoInterface.Wrappers.Polygon{false, false}([GeoInterface.Wrappers.LinearRing([(5.0, 0.0), … (3) … , (5.0, 0.0)])])
121+
GeoInterface.Wrappers.Polygon{false, false}([GeoInterface.Wrappers.LinearRing([(3.7328227614527916, 0.0), … (3) … , (3.7328227614527916, 0.0)])])
122+
GeoInterface.Wrappers.Polygon{false, false}([GeoInterface.Wrappers.LinearRing([(0.0, 5.0), … (3) … , (0.0, 5.0)])])
123+
```
124+
125+
"""
126+
function voronoi(geometries, ::Type{T} = Float64; kwargs...) where T
127+
return voronoi(Planar(), geometries, T; kwargs...)
128+
end
129+
130+
function voronoi(::Planar, geometries, ::Type{T} = Float64; clip_polygon = nothing, crs = __NoCRSProvided(), kwargs...) where T
131+
# Extract all points as tuples using GO.flatten
132+
# This handles any GeoInterface-compatible input
133+
points_iter = collect(flatten(tuples, GI.PointTrait, geometries))
134+
if crs isa __NoCRSProvided
135+
crs = GI.crs(geometries)
136+
end
137+
# if we have not figured it out yet, we can't do anything
138+
if crs isa __NoCRSProvided
139+
error("This code should be unreachable; please file an issue at https://github.com/JuliaGeometry/GeometryOps.jl/issues with the stacktrace and a reproducible example.")
140+
end
141+
142+
# Handle edge case of too few points
143+
if length(points_iter) < 3
144+
throw(ArgumentError("Voronoi tessellation requires at least 3 points, got $(length(points_iter))"))
145+
end
146+
147+
# Compute Delaunay triangulation
148+
tri = DelTri.triangulate(points_iter; kwargs...)
149+
150+
# Compute Voronoi tessellation from the triangulation
151+
_clip_polygon = if isnothing(clip_polygon)
152+
nothing
153+
elseif GI.geomtrait(clip_polygon) isa GI.PolygonTrait
154+
_clean_voronoi_clip_polygon_inputs(clip_polygon)
155+
else
156+
_clean_voronoi_clip_point_inputs(clip_polygon)
157+
end
158+
# if isclockwise(clip_polygon)
159+
vorn = DelTri.voronoi(tri; clip = true, clip_polygon = _clip_polygon)
160+
161+
polygons = GeoInterface.Wrappers.Polygon{false, false, Vector{GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{T, T}}, Nothing, Nothing}}, Nothing, typeof(crs)}[]
162+
sizehint!(polygons, DelTri.num_polygons(vorn))
163+
# Implementation below copied from Makie.jl
164+
# see https://github.com/MakieOrg/Makie.jl/blob/687c4466ce00154714297e36a7f610443c6ad5be/Makie/src/basic_recipes/voronoiplot.jl#L101-L110
165+
for i in DelTri.each_generator(vorn)
166+
!DelTri.has_polygon(vorn, i) && continue
167+
polygon_coords = DelTri.getxy.(DelTri.get_polygon_coordinates(vorn, i))
168+
push!(polygons, GI.Polygon([GI.LinearRing(polygon_coords)], crs = crs))
169+
# The code below gets the generator point, but we don't need it
170+
# gp = DelTri.getxy(DelTri.get_generator(vorn, i))
171+
# !isempty(polygon_coords) && push!(generators, gp)
172+
end
173+
174+
return polygons
175+
end
176+
177+
function _clean_voronoi_clip_polygon_inputs(clip_polygon)
178+
@assert GI.nhole(clip_polygon) == 0
179+
points = collect(flatten(tuples, GI.PointTrait, clip_polygon))
180+
npoints = GI.npoint(clip_polygon)
181+
if points[1] == points[end]
182+
npoints -= 1
183+
points = points[1:npoints]
184+
end
185+
point_order = collect(1:npoints)
186+
return _clean_voronoi_clip_point_inputs((points, point_order))
187+
end
188+
189+
function _clean_voronoi_clip_point_inputs((points, point_order)::Tuple{Vector{<:Tuple{<:Any, <:Any}}, Vector{<:Integer}})
190+
combined_data = collect(zip(points, point_order))
191+
sort!(combined_data, by = last)
192+
unique!(combined_data)
193+
194+
points, point_order = first.(combined_data), last.(combined_data)
195+
push!(points, points[1])
196+
push!(point_order, 1)
197+
198+
if isclockwise(GI.LineString(points))
199+
reverse!(points)
200+
end
201+
return points, point_order
202+
end
203+
204+
_clean_voronoi_clip_point_inputs((points, point_order)::Tuple{NTuple{<:Any, <:Tuple{<:Any, <:Any}}, NTuple{<:Any, <:Integer}}) =
205+
_clean_voronoi_clip_point_inputs((collect(points), collect(point_order)))
206+
207+
function _clean_voronoi_clip_point_inputs(clip_polygon)
208+
error("Clip polygon must be a polygon or other recognizable form, see the docstring for `DelaunayTriangulation.voronoi` for the recognizable form. Was neither, got $(typeof(clip_polygon))")
209+
return
210+
end

‎test/methods/voronoi.jl‎

Lines changed: 131 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,131 @@
1+
using Test
2+
using DelaunayTriangulation
3+
4+
import GeometryOps as GO
5+
import GeoInterface as GI
6+
using Random
7+
8+
@testset "Voronoi" begin
9+
@testset "Basic voronoi tessellation" begin
10+
# Test with simple points
11+
points = [(0.0, 0.0), (1.0, 0.0), (0.5, 1.0)]
12+
polygons = GO.voronoi(points)
13+
14+
# Should return 3 polygons, one for each point
15+
@test length(polygons) == 3
16+
17+
# Each should be a valid polygon
18+
for poly in polygons
19+
@test GI.isgeometry(poly)
20+
@test GI.geomtrait(poly) isa GI.PolygonTrait
21+
end
22+
end
23+
24+
@testset "Voronoi with various input types" begin
25+
# Test with GeoInterface points
26+
points = [GI.Point(0.0, 0.0), GI.Point(1.0, 0.0), GI.Point(0.5, 1.0)]
27+
polygons = GO.voronoi(points)
28+
@test length(polygons) == 3
29+
30+
# Test with mixed geometry collection
31+
geoms = [
32+
GI.Point(0.0, 0.0),
33+
GI.LineString([(1.0, 0.0), (1.5, 0.5)]), # Will extract endpoints
34+
GI.Point(0.5, 1.0)
35+
]
36+
polygons = GO.voronoi(geoms)
37+
@test length(polygons) == 4 # 1 + 2 + 1 points
38+
end
39+
40+
@testset "Voronoi with custom boundary" begin
41+
rng = Xoshiro(0)
42+
43+
points = [(0.25, 0.25), (0.75, 0.25), (0.75, 0.75), (0.25, 0.75), (0.5, 0.5)]
44+
boundary1 = (((0.0, 0.0), (0.0, 1.0), (1.0, 1.0), (1.0, 0.0), (0.0, 0.0)), (1, 2, 3, 4, 1))
45+
boundary2 = ([(0.0, 0.0), (1.0, 0.0), (1.0, 1.0), (0.0, 1.0), (0.0, 0.0)], [1, 2, 3, 4, 1])
46+
boundary3 = GI.Polygon([[(0.0, 0.0), (1.0, 0.0), (1.0, 1.0), (0.0, 1.0), (0.0, 0.0)]])
47+
48+
vorn0 = GO.voronoi(points) # clipped to convex hull
49+
vorn1 = GO.voronoi(points, clip_polygon = boundary1, rng = Xoshiro(0))
50+
vorn2 = GO.voronoi(points, clip_polygon = boundary2, rng = Xoshiro(0))
51+
vorn3 = GO.voronoi(points, clip_polygon = boundary3, rng = Xoshiro(0))
52+
53+
for vorn in (vorn0, vorn1, vorn2, vorn3)
54+
@test length(vorn) == 5
55+
# All polygons should be valid
56+
for poly in vorn
57+
@test GI.isgeometry(poly)
58+
@test GI.geomtrait(poly) isa GI.PolygonTrait
59+
end
60+
end
61+
end
62+
63+
@testset "Grid of points" begin
64+
# Create a regular grid
65+
xs = 0.0:0.5:2.0
66+
ys = 0.0:0.5:2.0
67+
points = [(x, y) for x in xs for y in ys]
68+
69+
polygons = GO.voronoi(points)
70+
@test length(polygons) == length(points)
71+
72+
# Each polygon should be valid
73+
for poly in polygons
74+
@test GI.isgeometry(poly)
75+
76+
# Get the exterior ring to check it's closed
77+
ring = GI.getexterior(poly)
78+
coords = GI.coordinates(ring)
79+
@test first(coords) ≈ last(coords) # Ring should be closed
80+
end
81+
end
82+
83+
@testset "Error handling" begin
84+
# Too few points
85+
@test_throws ArgumentError GO.voronoi([(0.0, 0.0)])
86+
@test_throws ArgumentError GO.voronoi([(0.0, 0.0), (1.0, 0.0)])
87+
88+
# Empty input
89+
@test_throws ArgumentError GO.voronoi([])
90+
end
91+
92+
@testset "Random points stress test" begin
93+
# Test with more points
94+
n = 50
95+
points = [(rand(), rand()) for _ in 1:n]
96+
97+
polygons = GO.voronoi(points)
98+
@test length(polygons) == n
99+
100+
# All should be valid polygons
101+
for poly in polygons
102+
@test GI.isgeometry(poly)
103+
@test GI.geomtrait(poly) isa GI.PolygonTrait
104+
end
105+
end
106+
107+
@testset "Clean clipping input" begin
108+
points = ((0.0, 0.0), (1.0, 0.0), (1.0, 1.0), (0.0, 1.0), (0.0, 0.0))
109+
order = (1, 2, 3, 4, 1)
110+
new_points, new_order = GO._clean_voronoi_clip_point_inputs((points, order))
111+
@test all(points .== new_points)
112+
@test all(order .== new_order)
113+
114+
reverse_points = reverse(points)
115+
new_points, new_order = GO._clean_voronoi_clip_point_inputs((reverse_points, order))
116+
@test all(points .== new_points)
117+
@test all(order .== new_order)
118+
119+
short_points = points[1:end-1]
120+
short_order = order[1:end-1]
121+
new_points, new_order = GO._clean_voronoi_clip_point_inputs((short_points, short_order))
122+
@test all(points .== new_points)
123+
@test all(order .== new_order)
124+
125+
shuffled_combos = shuffle(Xoshiro(0), collect(zip(points, order)))
126+
shuffled_points, shuffled_order = first.(shuffled_combos), last.(shuffled_combos)
127+
new_points, new_order = GO._clean_voronoi_clip_point_inputs((shuffled_points, shuffled_order))
128+
@test all(points .== new_points)
129+
@test all(order .== new_order)
130+
end
131+
end

‎test/runtests.jl‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -27,6 +27,7 @@ end
2727
@safetestset "Orientation" begin include("methods/orientation.jl") end
2828
@safetestset "Centroid" begin include("methods/centroid.jl") end
2929
@safetestset "Convex Hull" begin include("methods/convex_hull.jl") end
30+
@safetestset "Voronoi" begin include("methods/voronoi.jl") end
3031
@safetestset "DE-9IM Geom Relations" begin include("methods/geom_relations.jl") end
3132
@safetestset "Distance" begin include("methods/distance.jl") end
3233
@safetestset "Equals" begin include("methods/equals.jl") end

0 commit comments

Comments
 (0)