# 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. 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.

```
using DelaunayTriangulation, CairoMakie
pts = randn(2, 250)
tri = triangulate(pts)
vorn = voronoi(tri)
vorn_clip = voronoi(tri, true)
fig = Figure()
ax = Axis(fig[1, 1], aspect=1)
voronoiplot!(ax, vorn, strokecolor=:red, strokewidth=0.2, show_generators=false)
triplot!(ax, tri, strokewidth=0.0, strokecolor=(:black, 0.4), show_convex_hull=true)
xlims!(ax, -5, 5)
ylims!(ax, -5, 5)
ax = Axis(fig[1, 2], aspect=1)
voronoiplot!(ax, vorn_clip, strokecolor=:red, strokewidth=0.2, show_generators=false)
triplot!(ax, tri, strokewidth=0.0, strokecolor=(:black, 0.4), show_convex_hull=true)
xlims!(ax, -5, 5)
ylims!(ax, -5, 5)
```

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> DelaunayTriangulation.get_unbounded_polygons(vorn_clip)
Set{Int}()
julia> DelaunayTriangulation.get_boundary_polygons(vorn_clip)
Set{Int} 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!`

:

`DelaunayTriangulation.clip_voronoi_tessellation!`

— Function`clip_voronoi_tessellation!(vorn::VoronoiTessellation)`

Clip the Voronoi tessellation `vorn`

to the convex hull of the points in `vorn`

.

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.

`DelaunayTriangulation.find_all_intersections`

— Function`find_all_intersections(vorn::VoronoiTessellation)`

Find all intersections between the edges of the Voronoi tessellation and the boundary of the polygon.

**Outputs**

`boundary_sites`

: A dictionary of boundary sites.`segment_intersections`

: The intersection points.`exterior_circumcenters`

: The circumcenters that are outside of the domain.`equal_circumcenter_mapping`

: A mapping from the indices of the segment intersections that are equal to the circumcenter of a site to the index of the site.

`DelaunayTriangulation.initialise_clipping_arrays`

— Function`initialise_clipping_arrays(vorn::VoronoiTessellation)`

Initialise the arrays used in the clipping algorithm.

**Outputs**

`edges_to_process`

: The set of edges that are to be processed.`polygon_edge_queue`

: The queue of edges that are to be processed.`boundary_sites`

: A mapping from boundary sites to the indices of the segment intersections that are incident to the boundary site.`segment_intersections`

: The list of segment intersections.`processed_pairs`

: The set of pairs of edges and polygons that have been processed.`intersected_edge_cache`

: The list of intersected edges currently being considered.`exterior_circumcenters`

: The list of circumcenters of sites that are outside the boundary.`left_edge_intersectors`

: The set of sites that intersect the edge to the left of an edge currently being considered.`right_edge_intersectors`

: The set of sites that intersect the edge to the right of an edge currently being considered.`current_edge_intersectors`

: The set of sites that intersect the current edge being considered.`equal_circumcenter_mapping`

: A mapping from the indices of the segment intersections that are equal to the circumcenter of a site to the index of the site.

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!`

:

`DelaunayTriangulation.dequeue_and_process!`

— Function`dequeue_and_process!(vorn, polygon_edge_queue, edges_to_process, intersected_edge_cache, left_edge_intersectors, right_edge_intersectors, current_edge_intersectors, processed_pairs, boundary_sites, segment_intersections, exterior_circumcenters, equal_circumcenter_mapping)`

Dequeue an edge from `polygon_edge_queue`

and process it. If `polygon_edge_queue`

is empty, then we process the first edge in `edges_to_process`

.

**Arguments**

`vorn`

: The Voronoi tessellation.`polygon_edge_queue`

: The queue of edges that need to be processed.`edges_to_process`

: The edges that need to be processed.`intersected_edge_cache`

: A cache of intersected edges.`left_edge_intersectors`

: The intersection points of`left_edge`

with other edges.`right_edge_intersectors`

: The intersection points of`right_edge`

with other edges.`current_edge_intersectors`

: The intersection points of`current_edge`

with other edges.`processed_pairs`

: A set of pairs of edges and polygons that have already been processed.`boundary_sites`

: A dictionary of boundary sites.`segment_intersections`

: A dictionary of segment intersections.`exterior_circumcenters`

: A dictionary of exterior circumcenters.`equal_circumcenter_mapping`

: A mapping from the indices of the segment intersections that are equal to the circumcenter of a site to the index of the site.

`DelaunayTriangulation.enqueue_new_edge!`

— Function`enqueue_new_edge!(polygon_edge_queue, vorn::VoronoiTessellation, e)`

Enqueue the edge `e`

of the boundary to be processed.

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:

`DelaunayTriangulation.process_polygon!`

— Function`process_polygon!(vorn::VoronoiTessellation, e, incident_polygon, boundary_sites, segment_intersections, intersected_edge_cache, exterior_circumcenters, equal_circumcenter_mapping)`

Process the polygon `incident_polygon`

for all of its intersections based on the boundary edge `e`

, updating the caches in-place and returning `(left_edge, right_edge, e)`

, where `left_edge`

is the edge to the left of `e`

on the boundary and `right_edge`

is the edge to the right of `e`

on the boundary.

**Arguments**

`vorn::VoronoiTessellation`

: The Voronoi tessellation.`e`

: The edge on the boundary being considered.`incident_polygon`

: The index of the polygon being considered.`boundary_sites`

: The mapping from the indices of the sites on the boundary to the indices of the edges on the boundary that they intersect.`segment_intersections`

: The list of segment intersections.`intersected_edge_cache`

: A cache of the edges that have been intersected by the ray from`u`

to`v`

.`exterior_circumcenters`

: A list of the circumcenters of the sites that are outside the convex hull of the sites on the boundary.`equal_circumcenter_mapping`

: A mapping from the indices of the segment intersections that are equal to the circumcenter of a site to the index of the site.

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:

- First, both $u$ and $v$ could correspond to ghost vertices, meaning no intersection with the boundary is possible.
- 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. - 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.
- 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.

`DelaunayTriangulation.is_segment_between_two_ghosts`

— Function`is_segment_between_two_ghosts(u, v)`

Check if the segment between the sites `u`

and `v`

is between two ghost sites.

`DelaunayTriangulation.is_ray_going_in`

— Function`is_ray_going_in(u, v)`

Check if the ray from the site `u`

to the site `v`

is going in.

`DelaunayTriangulation.is_ray_going_out`

— Function`is_ray_going_out(u, v)`

Check if the ray from the site `u`

to the site `v`

is going out.

`DelaunayTriangulation.is_finite_segment`

— Function`is_finite_segment(u, v)`

Check if the segment between the sites `u`

and `v`

is finite.

`DelaunayTriangulation.process_segment_intersection!`

— Function```
process_segment_intersection!(
vorn::VoronoiTessellation,
u,
v,
e,
incident_polygon,
intersected_edge_cache,
segment_intersections,
boundary_sites,
exterior_circumcenters,
equal_circumcenter_mapping)
```

Process the intersection of the Voronoi polygon's edge `(u, v)`

with the edge `e`

of the boundary, returning the coordinates of the intersection and updating via `add_segment_intersection`

.

**Arguments**

`vorn`

: The Voronoi tessellation.`u`

: The index of the site`u`

.`v`

: The index of the site`v`

.`e`

: The edge`e`

of the boundary.`incident_polygon`

: The index of the Voronoi polygon currently being considered`(u, v)`

.`intersected_edge_cache`

: The list of intersected edges currently being considered.`segment_intersections`

: The list of segment intersections.`boundary_sites`

: A mapping from boundary sites to the indices of the segment intersections that are incident to the boundary site.`exterior_circumcenters`

: The list of circumcenters of sites that are outside the boundary.`equal_circumcenter_mapping`

: A mapping from the indices of the segment intersections that are equal to the circumcenter of a site to the index of the site.

`DelaunayTriangulation.process_ray_intersection!`

— Function```
process_ray_intersection!(
vorn::VoronoiTessellation,
u,
v,
incident_polygon,
intersected_edge_cache,
segment_intersections,
boundary_sites,
exterior_circumcenters,
equal_circumcenter_mapping)
```

Process the intersection of the Voronoi polygon of the site `u`

with the ray emanating from the circumcenter of the site `v`

, returning the coordinates of the intersection and updating via `add_segment_intersection`

.

**Arguments**

`vorn`

: The Voronoi tessellation.`u`

: The index of the site`u`

, given as a boundary index for the associated ghost triangle.`v`

: The index of the site`v`

.`incident_polygon`

: The index of the Voronoi polygon of the site`u`

that is incident to the ray emanating from the circumcenter of the site`v`

.`intersected_edge_cache`

: The list of intersected edges currently being considered.`segment_intersections`

: The list of segment intersections.`boundary_sites`

: A mapping from boundary sites to the indices of the segment intersections that are incident to the boundary site.`exterior_circumcenters`

: The list of circumcenters of sites that are outside the boundary.`equal_circumcenter_mapping`

: A mapping from the indices of the segment intersections that are equal to the circumcenter of a site to the index of the site.

`DelaunayTriangulation.add_to_intersected_edge_cache!`

— Function`add_to_intersected_edge_cache!(intersected_edge_cache, u, v, a, b)`

Add the edge `uv`

to the list of intersected edges, where `uv`

is the edge of the Voronoi polygon intersecting the edge of the boundary `ab`

.

`DelaunayTriangulation.add_segment_intersection!`

— Function`add_segment_intersection!(segment_intersections, boundary_sites, intersection_point, incident_polygon::I) where {I}`

Add the intersection point to the list of segment intersections and return the index of the intersection point in the list. If the intersection point already exists in the list, then the index of the existing point is returned and used instead.

**Arguments**

`segment_intersections`

: The list of segment intersections.`boundary_sites`

: A mapping from boundary sites to the indices of the segment intersections that are incident to the boundary site.`intersection_point`

: The intersection point to add.`incident_polygon`

: The index of the polygon that is incident to the intersection point.

`DelaunayTriangulation.segment_intersection_coordinates`

— Function`segment_intersection_coordinates(a, b, c, d)`

Finds the coordinates of the intersection of the line segment from `a`

to `b`

with the line segment from `c`

to `d`

.

`DelaunayTriangulation.intersection_of_edge_and_bisector_ray`

— Function`intersection_of_edge_and_bisector_ray(a, b, c)`

Given an edge `(a, b)`

and a ray emanating from `c`

perpendicular with the edge and collinear with its midpoint, tests if `c`

intersects the edge, and if so, returns the `(cert, p)`

, where `p`

is the intersection point (which is the midpoint) and `c`

is the position of `c`

relative to `(a, b)`

. If there is no intersection, `p = (NaN, NaN)`

is returned (together with `cert`

). The ray should be directed to the left of the edge.

`DelaunayTriangulation.classify_and_compute_segment_intersection`

— Function`classify_and_compute_segment_intersection(a, b, c, d)`

Given two line segments `(a, b)`

and `(c, d)`

, classifies the intersection of the two segments, returning `(cert, cert_c, cert_d, p)`

, where `p`

is the intersection point or `(NaN, NaN)`

if there is no intersection. The certificate `cert`

determines the intersection type, giving

`Cert.None`

: No intersections.`Cert.Single`

: There is an intersection point.`Cert.Touching`

: There is an intersection point, and one of`c`

and`d`

is the intersection point.`Cert.Multiple`

: The closed segments meet in one or several points.

The certificates `cert_c`

and `cert_d`

similarly return the positions of `c`

and `d`

relative to `(a, b)`

, respectively.

`DelaunayTriangulation.process_ray_intersection_with_other_edges!`

— Function```
process_ray_intersection_with_other_edges!(vorn::VoronoiTessellation, u, v, e, left_edge, right_edge, r, segment_intersections,
boundary_sites, incident_polygon, equal_circumcenter_mapping, intersected_edge_cache)
```

Process the intersection of the ray from the ghost site `u`

to the site `v`

with the edges `e`

, `left_edge`

and `right_edge`

.

**Arguments**

`vorn::VoronoiTessellation`

: The Voronoi tessellation.`u`

: The index of the ghost site.`v`

: The index of the site`u`

is going to.`e`

: The edge on the boundary being considered.`left_edge`

: The edge to the left of`e`

on the boundary.`right_edge`

: The edge to the right of`e`

on the boundary.`r`

: The coordinates of the intersection of the ray from`u`

to`v`

with some edge. If`any(isnan, r)`

, then the ray does not intersect any edge and we skip.`segment_intersections`

: The list of segment intersections.`boundary_sites`

: The mapping from the indices of the sites on the boundary to the indices of the edges on the boundary that they intersect.`incident_polygon`

: The index of the polygon that contains the intersection of the ray from`u`

to`v`

with the boundary.`equal_circumcenter_mapping`

: A mapping from the indices of the segment intersections that are equal to the circumcenter of a site to the index of the site.`intersected_edge_cache`

: A cache of the edges that have been intersected by the ray from`u`

to`v`

.

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!`

.

`DelaunayTriangulation.classify_intersections!`

— Function`classify_intersections!(intersected_edge_cache, left_edge_intersectors, right_edge_intersectors, current_edge_intersectors, left_edge, right_edge, current_edge)`

Classify the intersections in `intersected_edge_cache`

into `left_edge_intersectors`

, `right_edge_intersectors`

, and `current_edge_intersectors`

based on whether they intersect `left_edge`

, `right_edge`

, or `current_edge`

, respectively.

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!`

:

`DelaunayTriangulation.process_intersection_points!`

— Function`process_intersection_points!(polygon_edge_queue, vorn, current_incident_polygon, left_edge_intersectors, right_edge_intersectors, current_edge_intersectors, left_edge, right_edge, current_edge, processed_pairs, segment_intersections, boundary_sites)`

Process the intersection points in `left_edge_intersectors`

, `right_edge_intersectors`

, and `current_edge_intersectors`

and add the new edges to `polygon_edge_queue`

if necessary. Special care is taken to not miss any corner points.

The rules are based on the paper "Efficient Computation of Clipped Voronoi Diagram for Mesh Generation" by Yan, Wang, Levy, and Liu. Namely, an edge that intersects a boundary edge and one to it has its shared vertex added to the queue together with the current polygon (`current_incident_polygon`

) being considered, and any intersections have the adjacent polygon added to the queue together with the intersecting edge. These are not strictly the rules in the paper, but they are the rules that I was able to implement since they do not share their code.

**Arguments**

`polygon_edge_queue`

: The queue of edges that need to be processed.`vorn`

: The Voronoi tessellation.`current_incident_polygon`

: The index of the current polygon being processed.`left_edge_intersectors`

: The intersection points of`left_edge`

with other edges.`right_edge_intersectors`

: The intersection points of`right_edge`

with other edges.`current_edge_intersectors`

: The intersection points of`current_edge`

with other edges.`left_edge`

: The left edge of the current polygon.`right_edge`

: The right edge of the current polygon.`current_edge`

: The current edge of the current polygon.`processed_pairs`

: A set of pairs of edges and polygons that have already been processed.`segment_intersections`

: A dictionary of segment intersections.`boundary_sites`

: A dictionary of boundary sites.

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!`

:

`DelaunayTriangulation.add_intersection_points!`

— Function`add_intersection_points!(vorn::VoronoiTessellation, segment_intersections)`

Add the intersection points to the polygon vertices.

### 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!`

:

`DelaunayTriangulation.clip_all_polygons!`

— Function`clip_all_polygons!(vorn::VoronoiTessellation, n, boundary_sites, exterior_circumcenters, equal_circumcenter_mapping, is_convex)`

Clip all of the polygons in the Voronoi tessellation.

**Arguments**

`vorn`

: The Voronoi tessellation.`n`

: The number of vertices in the tessellation before clipping.`boundary_sites`

: A dictionary of boundary sites.`exterior_circumcenters`

: Any exterior circumcenters to be filtered out.`equal_circumcenter_mapping`

: A mapping from the indices of the segment intersections that are equal to the circumcenter of a site to the index of the site.`is_convex`

: Whether the polygon is convex or not. Not currently used.

This function iterates over all the identified boundary polygons and calls `clip_polygon!`

on each:

`DelaunayTriangulation.clip_polygon!`

— Function`clip_polygon!(vorn::VoronoiTessellation, n, points, polygon, new_verts, exterior_circumcenters, equal_circumcenter_mapping, is_convex)`

Clip the polygon `polygon`

by removing the vertices that are outside of the domain and adding the new vertices `new_verts`

to the polygon.

**Arguments**

`vorn`

: The Voronoi tessellation.`n`

: The number of vertices in the tessellation before clipping.`points`

: The points of the tessellation.`polygon`

: The index of the polygon to be clipped.`new_verts`

: The indices of the new vertices that are added to the polygon.`exterior_circumcenters`

: Any exterior circumcenters to be filtered out.`equal_circumcenter_mapping`

: A mapping from the indices of the segment intersections that are equal to the circumcenter of a site to the index of the site.`is_convex`

: Whether the polygon is convex or not. Not currently used.

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!`

:

`DelaunayTriangulation.add_all_boundary_polygons!`

— Function`add_all_boundary_polygons!(vorn::VoronoiTessellation, boundary_sites)`

Add all of the boundary polygons to the Voronoi tessellation.

With this last step, the clipping is complete.