# Voronoi Tessellations

We provide support for constructing Voronoi tessellations, the dual graph to a Delaunay triangulation (strictly true only for unconstrained triangulations). The main function for this is `voronoi`

:

`DelaunayTriangulation.voronoi`

— Function`voronoi(tri::Triangulation, clip=false) -> VoronoiTessellation`

Construct the Voronoi tessellation of the points in `tri`

. If `clip`

is `true`

then the Voronoi tessellation will be clipped to the convex hull of the points in `tri`

.

If you are interested instead in clipping the tessellation to a rectangular bounding box, see `get_polygon_coordinates`

which allows for a bounding box to be applied to an unclipped tessellation, returning a vector of the coordinates of polygons, clipping to the bounding box where applicable. `polygon_bounds`

can be used to obtain good default bounding boxes. Note that if you are not worried about this clipping, and you know that your polygon is not unbounded (which would mean it is not in the `unbounded_polygons`

field of the `VoronoiTessellation`

output), then you should instead use `get_polygon(vorn, i)`

to get the indices of the points in `vorn`

that define the polygon, and then use `get_polygon_point`

to get the coordinates.

Exact predicates are used only for classifying intersections, but no special methods are used for computing the intersections themselves.

Clipping is not yet guaranteed to work for constrained triangulations.

The algorithm is really quite simple:

- Construct the Delaunay triangulation.
- For each point in the triangulation, connect the circumcenters of all triangles with that point as a vertex.

(Omitting, of course, some details like duplicate circumcenters from cocircular susets – as with many things in computational geometry, anything is "simple" until you think about reality.) Note that extra care is taken for the unbounded edges, which we simply represent using negative indices. The negative indices we use have no real meaning, but they are unique so that we can uniquely identify for any unbounded ray the associated ghost triangle (shown by example later).

## Example

Let us give an example. The first step in any tessellation is the construction of the triangulation.

```
using DelaunayTriangulation, CairoMakie
pts = rand(2, 25)
tri = triangulate(pts)
```

To now construct the tessellation, use `voronoi`

.

`vorn = voronoi(tri)`

We can visualise the diagram using `voronoiplot`

. This function uses `get_polygon_coordinates`

to get a finite representation of the unbounded edges, clipping the polygons to some bounding box determined by `DelaunayTriangulation.polygon_bounds`

.

```
fig, ax, sc = voronoiplot(vorn, strokecolor=:red, markersize=9)
triplot!(ax, tri, strokewidth=1, strokecolor=(:black, 0.4))
```

You can set the bounding box using the `bounding_box`

keyword if you want it determined manually rather than programatically. The function that handles chopping to a box is `get_polygon_coordinates`

, which primarily relies on `intersection_of_ray_with_bounding_box`

:

`DelaunayTriangulation.intersection_of_ray_with_bounding_box`

— Function`intersection_of_ray_with_bounding_box(p, q, a, b, c, d)`

Given a ray starting at `p`

and in the direction of `q`

, finds the intersection of the ray with the bounding box `[a, b] × [c, d]`

. It is assumed that `p`

is inside the bounding box, but `q`

can be inside or outside.

We will see later how to clip these polygons to the convex hull show in the red dashed line.

Let us now actually show the data structure. There are several fields:

```
julia> vorn.
adjacent circumcenter_to_triangle generators polygons triangulation
boundary_polygons cocircular_circumcenters polygon_points triangle_to_circumcenter unbounded_polygon
```

The `adjacent`

field is similar to the one we deal with for triangulations, except edges now map to their associated polygon. For example, we have

```
julia> get_polygon(vorn, 1)
7-element Vector{Int}:
23
17
-8
-6
32
34
23
julia> get_adjacent(vorn, 23, 17) == get_adjacent(vorn, 17, -8) == 1
true
```

The `circumcenter_to_triangle`

field allows us to go from the index of a circumcenter to its associated triangle. For instance, if we look at the vertices of the first polygon above, we may want to know where this vertex `34`

came from:

```
julia> DelaunayTriangulation.get_circumcenter_to_triangle(vorn, 34)
(10, 4, 1)
```

Thus, it came from the triangle `(10, 4, 1)`

. Note that the triangles used in this structure are always sorted so that the minimum index is last (while preserving the counter-clockwise orientation).

The `generators`

field is a repackaged version of `tri.points`

:

```
julia> DelaunayTriangulation.get_generators(vorn)
Dict{Int, Tuple{Float64, Float64}} with 25 entries:
5 => (0.684537, 0.618674)
16 => (0.585988, 0.554157)
20 => (0.43467, 0.161223)
12 => (0.00191717, 0.70299)
24 => (0.562014, 0.246805)
8 => (0.446736, 0.717801)
17 => (0.00141105, 0.639474)
23 => (0.953932, 0.352213)
1 => (0.234991, 0.0601079)
22 => (0.559084, 0.538079)
19 => (0.650726, 0.439868)
6 => (0.966998, 0.627038)
11 => (0.464044, 0.594492)
9 => (0.512324, 0.0257408)
14 => (0.559852, 0.55035)
3 => (0.139261, 0.472539)
7 => (0.560707, 0.767171)
25 => (0.0731719, 0.791413)
4 => (0.389615, 0.253273)
15 => (0.597794, 0.761001)
⋮ =>
```

Here, we need to use a `Dict`

to represent the points since the original triangulation may not include all points. With this, we always know that `get_generator(vorn, i) == get_point(tri, i)`

. If we want the coordinates of the `i`

th generator, we could use

```
julia> get_generator(vorn, 16)
(0.5859881889826812, 0.5541565595408635)
julia> get_generator(vorn, 7, 22, 1)
((0.560707474714347, 0.7671706685373769), (0.5590841028546631, 0.5380792727814472), (0.2349910552466512, 0.060107935323076456))
```

The `polygons`

field is a `Dict`

that maps the indices of the generators to their polygon:

```
julia> DelaunayTriangulation.get_polygons(vorn)
Dict{Int, Vector{Int}} with 25 entries:
5 => [24, 29, 15, 8, 24]
16 => [24, 37, 26, 20, 16, 29, 24]
20 => [30, 38, 14, 6, 30]
12 => [3, -3, -9, 39, 3]
24 => [33, 6, 14, 9, 18, 2, 33]
8 => [35, 12, 5, 21, 31, 10, 35]
17 => [39, -9, -8, 17, 22, 39]
23 => [11, 18, 9, -2, -5, 11]
1 => [23, 17, -8, -6, 32, 34, 23]
22 => [33, 2, 16, 20, 1, 4, 33]
19 => [29, 16, 2, 18, 11, 15, 29]
6 => [15, 11, -5, -7, 8, 15]
11 => [1, 7, 10, 31, 36, 19, 4, 1]
9 => [38, 32, -6, -2, 9, 14, 38]
14 => [20, 26, 7, 1, 20]
3 => [19, 36, 22, 17, 23, 19]
7 => [13, 25, 12, 35, 13]
25 => [27, 28, -1, -3, 3, 27]
4 => [4, 19, 23, 34, 30, 6, 33, 4]
15 => [8, -7, -4, 25, 13, 37, 24, 8]
⋮ => ⋮
```

So, if we wanted the polygon for the fifth generator, we would use:

```
julia> get_polygon(vorn, 5)
5-element Vector{Int}:
24
29
15
8
24
```

These polygons are always sorted to be counter-clockwise, with the first element equalling the last also.

The `triangulation`

field is simply `tri`

. The `boundary_polygons`

field is essentially the same as `unbounded_polygons`

in most cases, and is only relevant later on when we consider clipped tessellations (for unbounded tessellations, this field is always empty). The `cocircular_circumcenters`

field is just for keeping track of cocircular circumcenters to avoid adding duplicate points into the tessellation; it is an implementation detail, you shouldn't need to use it.

The `polygon_points`

field stores the coordinates of the points used for the polygons:

```
julia> DelaunayTriangulation.get_polygon_points(vorn)
39-element Vector{Tuple{Float64, Float64}}:
(0.5006477732564246, 0.5478946133538716)
(0.5008095733914651, 0.3918411792218543)
(0.10565601252653685, 0.692314380288574)
(0.4270143529620006, 0.4238420900557963)
(0.33079352748966667, 0.7624750560470457)
(0.4753703522300822, 0.23819561748527615)
(0.5226957825549574, 0.5957485305429776)
(0.8205472522054272, 0.799163351762504)
(0.8189082811407147, 0.07294487509120647)
(0.4026965184932325, 0.6487501595269026)
⋮
(0.330822898041236, 0.16744471621902937)
(0.3146257041841115, 0.6363882426877013)
(0.34474111887983666, -0.1904242188538559)
(0.4811272594567467, 0.391643188260307)
(0.2982325567607496, 0.16795358951510841)
(0.5199378327863162, 0.7050503627303082)
(0.2840418508918076, 0.5804155201166311)
(0.5862985419779859, 0.657897964228626)
(0.4650876739497739, 0.08866174194893103)
(0.10972835538971261, 0.6703707990569139)
```

If we wanted the coordinates for the sixth vertex, we would do:

```
julia> get_polygon_point(vorn, 6)
(0.4753703522300822, 0.23819561748527615)
```

The `triangle_to_circumcenter`

field is similar to `circumcenter_to_triangle`

, except this takes triangles to the index of their associated circumcenter:

```
julia> DelaunayTriangulation.get_triangle_to_circumcenter(vorn)
Dict{Tuple{Int, Int, Int}, Int} with 48 entries:
(22, 14, 11) => 1
(8, 13, 7) => 35
(22, 19, 16) => 16
(25, 18, -1) => -1
(3, 17, 1) => 17
(17, 21, 12) => 39
(14, 13, 11) => 7
(23, 9, -1) => -2
(6, 15, 5) => 8
(23, 24, 9) => 9
(19, 23, 6) => 11
(12, 25, -1) => -3
(18, 15, -1) => -4
(16, 19, 5) => 29
(13, 15, 7) => 13
(4, 3, 1) => 23
(19, 6, 5) => 15
(21, 11, 8) => 31
(6, 23, -1) => -5
(18, 25, 2) => 28
⋮ => ⋮
```

Lastly, the `unbounded_polygons`

field stores the indices of all generators whose associated polygon is unbounded:

```
julia> DelaunayTriangulation.get_unbounded_polygons(vorn)
Set{Int} with 9 elements:
6
15
25
18
9
12
17
23
1
```

## Relevant Docstrings

Here are some relevant docstrings for the construction of the tessellation.

`DelaunayTriangulation.initialise_voronoi_tessellation`

— Function`initialise_voronoi_tessellation(tri::Triangulation)`

Initialise a `VoronoiTessellation`

from the triangulation `tri`

.

`DelaunayTriangulation.prepare_add_voronoi_polygon`

— Function`prepare_add_voronoi_polygon(vorn::VoronoiTessellation, i)`

Prepare the arrays for adding the Voronoi polygon for the point `i`

to the `VoronoiTessellation`

`vorn`

, returning:

`S`

: The vertices surrounding the generator`i`

.`B`

: An empty array for storing vertices.

`DelaunayTriangulation.get_next_triangle_for_voronoi_polygon`

— Function`get_next_triangle_for_voronoi_polygon(vorn::VoronoiTessellation, i, k, S, m)`

Get the next triangle for the Voronoi polygon for the point `i`

in the `VoronoiTessellation`

`vorn`

, returning:

`ci`

: The index for the circumcenter of the next triangle.`k`

: The index of the next vertex in`S`

.

`DelaunayTriangulation.connect_circumcenters!`

— Function`connect_circumcenters!(B, ci)`

Add the circumcenter `ci`

to the array `B`

.

`DelaunayTriangulation.add_edge_to_voronoi_polygon!`

— Function`add_edge_to_voronoi_polygon!(B, vorn::VoronoiTessellation, i, k, S, m, encountered_duplicate_circumcenter)`

Add the next edge to the Voronoi polygon for the point `i`

in the `VoronoiTessellation`

`vorn`

, returning:

`ci`

: The index for the circumcenter of the triangle considered.`encountered_duplicate_circumcenter`

: Whether or not a duplicate circumcenter has been encountered.`k`

: The index of the latest vertex in`S`

.

`DelaunayTriangulation.close_voronoi_polygon!`

— Function`close_voronoi_polygon!(vorn::VoronoiTessellation, B, i, encountered_duplicate_circumcenter, prev_ci)`

Close the Voronoi polygon for the point `i`

in the `VoronoiTessellation`

`vorn`

.

`DelaunayTriangulation.add_voronoi_polygon!`

— Function`add_voronoi_polygon!(vorn::VoronoiTessellation, i)`

Add the Voronoi polygon for the point `i`

to the `VoronoiTessellation`

`vorn`

.