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)
Clipped Voronoi Tessellation

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

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_intersectionsFunction
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_arraysFunction
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.

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:

  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.

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.intersection_of_edge_and_bisector_rayFunction
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_intersectionFunction
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!:

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

With this last step, the clipping is complete.