Point Location

The most expensive step for building a Delaunay triangulation is the point location step, wherein we need to find a triangle that contains a given point. The code for this turns out to be very complicated so that we can correctly handle points outside of the domain, inside interior holes, collinear with other points, on the corner, etc. The main function that handles all of this is jump_and_march, derived from the jump-and-march algorithm of Mücke, Saias, and Zhu (1999). We also provide a method that simply searches over all triangles, given below.

The Main Method

Below we list docstrings for the main jump and march algorithm.

DelaunayTriangulation.jump_and_marchFunction
jump_and_march(tri::Triangulation, q;
    point_indices=each_point_index(tri),
    m=default_num_samples(length(point_indices)),
    try_points=(),
    k=select_initial_point(get_points(tri), q; m, point_indices, try_points),
    check_existence::C=Val(has_multiple_segments(tri)),
    store_history::F=Val(false),
    history=nothing,
    rng::AbstractRNG=Random.default_rng(),
    exterior_curve_index=1,
    maxiters=2 + length(exterior_curve_index) - num_solid_vertices(tri) + num_solid_edges(tri)) where {C,F}

Returns the triangle containing q using the jump-and-march algorithm.

Arguments

  • tri::Triangulation: The triangulation.
  • q: The query point.

Keyword Arguments

  • point_indices=each_point_index(tri): The indices of the points in the triangulation.
  • m=default_num_samples(length(point_indices)): The number of samples to use when sampling the point to start the algorithm at.
  • try_points=(): Additional points to try when determining which point to start at.
  • k=select_initial_point(get_points(tri), q; m, point_indices, try_points): The index of the point to start the algorithm at.
  • store_history::F=Val(false): Whether to record the history of the algorithm. See also PointLocationHistory.
  • history=nothing: The object to store the history into, if is_true(store_history).
  • rng::AbstractRNG=Random.default_rng(): The random number generator to use.
  • exterior_curve_index=1: The curve (or curves) corresponding to the outermost boundary.
  • maxiters = num_triangles(tri): Maximum number of iterations to perform before restarting the algorithm at a new initial point.
  • concavity_protection=false: When your triangulation has concave boundaries, it is possible that a ghost triangle is incorrectly classified as containing the point. By setting this to true, this will be protected against.
Note

You shouldn't ever need maxiters if your triangulation is convex everywhere, as Delaunay triangulations have no problems with jump-and-march, as the sequence of triangles visited is acyclic (H. Edelsbrunner, An acyclicity theorem for cell complexes in d dimensions, Combinatorica 10 (1990) 251– 260.) However, if the triangulation is not convex, e.g. if you have a constrained triangulation with boundaries and excavations, then an infinite loop can be found where we just keep walking in circles. In this case, you can use the maxiters keyword argument to specify the maximum number of iterations to perform before reinitialising the algorithm at a random vertex. When reinitialising, the value of m is doubled each time.

Outputs

Returns V, the triangle in tri containing q.

Warning

If your triangulation does not have ghost triangles, and the point q is outside of the triangulation, this function may fail to terminate. You may like to add ghost triangles in this case (using add_ghost_triangles!), noting that there is no actual triangle that q is inside of when it is outside of the triangulation unless ghost triangles are present.

The variant of the jump and march algorithm for points outside of the triangle is also accessed via the jump_and_march, calling into exterior_jump_and_march:

DelaunayTriangulation.exterior_jump_and_marchFunction
exterior_jump_and_march(
    pts, 
    adj, 
    boundary_index_ranges, 
    representative_point_list,
    boundary_map, 
    k, 
    q, 
    check_existence=Val(has_multiple_segments(boundary_map)),
    bnd_idx=I(BoundaryIndex))

Given a point q outside of the triangulation, finds the ghost triangle containing it.

Arguments

  • pts: The collection of points.
  • adj: The Adjacent map.
  • boundary_index_ranges: A Dict mapping boundary indices to ranges from construct_boundary_index_ranges.
  • representative_point_list: A Dict mapping curve indices to representative points.
  • boundary_map: A boundary map from construct_boundary_map.
  • k: The vertex to start from.
  • q: The point outside of the triangulation.
  • check_existence=Val(has_multiple_segments(boundary_map))): Used to check that the edge exists when using get_adjacent, in case there are multiple segments.
  • bnd_idx=I(BoundaryIndex): A boundary index corresponding to the boundary that k is on.

Ouptut

  • The edge (i, j) such that the ghost triangle (i, j, g) contains q, and g = get_adjacent(adj, i, j).
Warning

The result is meaningless if q is inside of the triangulation.

You should not need to call into this method directly.

The core complexity of the algorithm comes from having to find the direction of the point from an initial search point. The docstrings for some of these initialisers are given below.

DelaunayTriangulation.default_num_samplesFunction
default_num_samples(num_points)

Returns ceil(cbrt(num_points)). This is the default number of samples to use when sampling in the jump-and-march algorithm.

DelaunayTriangulation.select_initial_pointFunction
select_initial_point(tri::Triangulation{P,Ts,I}, q;
    point_indices=each_solid_vertex(tri),
    m=default_num_samples(num_vertices(point_indices)),
    try_points=(),
    rng::AbstractRNG=Random.default_rng()) where {P,Ts,I}

Given a triangulation tri and a point q, select m random points from tri and return the one that is closet to q.

Arguments

  • tri: The Triangulation.
  • q: The coordinates of the query point, or its index in pts.

Keyword Arguments

  • point_indices=each_point_index(pts): The indices of the points to consider.
  • m=default_num_samples(length(point_indices)): The number of points to sample.
  • try_points=(): The indices of points to try in addition to the sampled points.
  • rng::AbstractRNG=Random.default_rng(): The random number generator to use.

Outputs

  • k: The index of the point in pts that is closest to q out of the points sampled.
DelaunayTriangulation.select_initial_triangle_interior_nodeFunction
select_initial_triangle_interior_node(tri, k, q, store_history::F=Val(false), history=nothing, rng::AbstractRNG=Random.default_rng()) where {F}

Selects an initial triangle for the jump-and-march algorithm, starting from a point with index k and finding the triangle such that the line from the kth point to q intersects it. It is assumed that k is a point that is not on the boundary.

Arguments

  • tri: The Triangulation.
  • k: The index of the point in pts that we are starting at.
  • q: The point being searched for.
  • boundary_index_ranges: The Dict handling the mapping of boundary indices to the range of boundary indices belonging to the same curve. See construct_boundary_index_ranges.
  • store_history=Val(false): Whether to store visited triangles. Exterior ghost triangles will not be stored.
  • history=nothing: The history. This should be a PointLocationHistory type if store_history is true.
  • rng::AbstractRNG = Random.default_rng(): The random number generator.

Outputs

  • p: The kth point in tri.
  • i, j: These are indices defining the edge of a triangle including the point p, such that i is to the left of the line pq and j is to the right of pq.
  • pᵢ, pⱼ: The points in tri corresponding to the indices in i and j, respectively.
DelaunayTriangulation.check_for_intersections_with_adjacent_boundary_edgesFunction
check_for_intersections_with_adjacent_boundary_edges(tri, k, q, bnd_idx=I(BoundaryIndex))

Assuming that k is on the outer boundary, this function searches down the boundary edges adjacent to k to try and locate a triangle or edge containing q.

See also search_down_adjacent_boundary_edges, which uses this function to determine an initial direction to search along a straight boundary in case q is collinear with it.

Arguments

  • tri: The Triangulation.
  • k: The outer boundary point k to start searching from.
  • q: The query point q.
  • bnd_idx=I(BoundaryIndex): The boundary index for the boundary that k is on.

Outputs

The output is a 5-tuple, with the first three elements having several possible forms:

  • (Certificate.Outside, Certificate.Outside, k): The point q is not collinear with either of the adjacent boundary edges.
  • (Certificate.Right, C, r), where C is either Certificate.On or Certificate.Right and r is the vertex right of k: The point q is collinear with the edge to the right of k. If C is Certificate.On, then q is on the edge, whereas C being Certificate.Right means it is right of the edge.
  • (Certificate.Left, C, ℓ), where C is either Certificate.On or Certificate.Left and is the vertex left of k: The point q is collinear with the edge to the left of k. If C is Certificate.On, then q is on the edge, whereas C being Certificate.Left means it is left of the edge.

In the latter two outputs above, C could also be Certificate.Degenerate, which means that q is get_point(pts, r) or get_point(pts, ℓ), respectively.

The latter two elements of the tuple are:

  • right_cert: The position of q relative to the edge (p, p_right), where p_right is the point on the boundary to the right of p.
  • left_cert: The position of q relative to the edge (p, p_left), where p_left is the point on the boundary to the left of p.

These returned values are useful in case we need to go to check_for_intersections_with_interior_edges_adjacent_to_boundary_node, since we can reuse these certificates.

DelaunayTriangulation.search_down_adjacent_boundary_edgesFunction
search_down_adjacent_boundary_edges(
    tri,
    k, 
    q,
    direction, 
    q_pos, 
    next_vertex,
    store_history=Val(false),
    history=nothing,
    bnd_idx=I(BoundaryIndex))

Starting at the outer boundary node k, walks down the boundary in the direction of q until finding q or finding that it is outside of the triangulation.

Arguments

  • tri: The Triangulation.
  • k: The outer boundary index.
  • q: The point being searched for.
  • direction: The direction of q from get_point(pts, k).
  • q_pos: The certificate for the position of q from this point k.
  • next_vertex: The next vertex in the direction of q (this argument comes from check_for_intersections_with_adjacent_boundary_edges).
  • store_history=Val(false): Whether to store the history of the search.
  • history=nothing: The history of the search.
  • bnd_idx=I(BoundaryIndex): The boundary index for the boundary that k is on.

Outputs

The output takes the form (cert, u, v, w), where:

  • cert: This is Certificate.On if q is on the edge (u, v), and Certificate.Outside if q is outside of the triangulation.
  • (u, v, w): If is_on(cert), then this is a positively oriented triangle with q on the edge (u, v). Otherwise, (u, v, w) is a ghost triangle close to q.
Warning

This function relies on the assumption that the geometry is convex.

DelaunayTriangulation.check_for_intersections_with_interior_edges_adjacent_to_boundary_nodeFunction
check_for_intersections_with_interior_edges_adjacent_to_boundary_node(
    tri,
    k, 
    q, 
    right_cert, 
    left_cert, 
    store_history=Val(false), 
    history=nothing,
    bnd_idx=I(BoundaryIndex))

Checks if the line connecting the kth point of pts to q intersects any of the edges neighbouring the boundary node k.

Arguments

Outputs

There are several possible forms for the returned values. These are listed below, letting p be the kth point, pᵢ the point corresponding to the index i, and pⱼ the point corresponding to the index j:

  • (i, j, Certificate.Single, Certificate.Outside)

The line pq intersects the edge pᵢpⱼ, and (j, i, k) is a positively oriented triangle. In particular, pᵢ is left of pq and pⱼ is right of pq.

  • (i, j, Certificate.None, Certificate.Inside)

The point q is inside the positively oriented triangle (i, j, k).

  • (zero(I), zero(I), Cert.None, Cert.Outside)

The point q is outside of the triangulation. Note that I is the integer type.

  • (i, j, Cert.On, Cert.Inside)

The point q is on the edge pᵢpⱼ, and so is inside the positively oriented triangle (i, j, k).

  • (i, j, Cert.Right, Cert.Outside)

The point q is collinear with the edge pᵢpⱼ, but is off of it and further into the triangulation.

Warning

This function should only be used after check_for_intersections_with_adjacent_boundary_edges, and currently is only guaranteed to work on convex geometries.

History

If you need to, you can also store the history of the algorithm. This is primarily only used for (and only tested for) constrained triangulations, as it helps us locate which triangles to delete. The struct we use for storing history is given below.

DelaunayTriangulation.PointLocationHistoryType
PointLocationHistory{T,E,I}

History from using jump_and_march.

Fields

  • triangles::Vector{T}: The visited triangles.
  • collinear_segments::Vector{E}: Segments collinear with the original line pq using to jump.
  • collinear_point_indices::Vector{I}: This field contains indices to segments in collinear_segments that refer to points that were on the original segment, but there is no valid segment for them. We use manually fix this after the fact. For example, we could add an edge (1, 14), when really we mean something like (7, 14) which isn't a valid edge.
  • left_vertices::Vector{I}: Vertices from the visited triangles to the left of pq.
  • right_verices::Vector{I}: Vertices from the visited triangles to the right of pq.

Basic Description of the Algorithm

Let us give a basic description of what jump_and_march does.

First, using select_initial_point, an initial point to start the algorithm is selected. This function samples some number $m$ of points, and then selects the point that is closest to the query point $q$.

Let the initially selected point be $p_k$. We break the discussion into two cases, where $p_k$ is an interior point and $p_k$ is a point on the boundary.

If $p_k$ is not a point on the boundary, then it is possible to completely rotate around the point searching for a triangle such that line $\overrightarrow{p_kq}$ intersects an edge of the triangle. This will give us an edge $e_{ij}$ that $\overrightarrow{p_kq}$ intersects, and we will put $p_i$ to the left of $\overrightarrow{p_kq}$ and $p_j$ to the right. The function that handles this selection is select_initial_triangle_interior_node, which starts by handling the case of collinear edges and then rotates around.

Now suppose that $p_k$ is a point on the outer boundary. There are several possibilities in this case. First, $\overrightarrow{p_kq}$ might intersect a neighbouring boundary edge; secondly, $\overrightarrow{p_kq}$ could intersect a neighbouring interior edge; thirdly, $\overrightarrow{p_kq}$ could point away from the boundary, meaning $q$ is outside of the boundary. The function starts by checking the neighbouring boundary edges, done via check_for_intersections_with_adjacent_boundary_edges, making use of get_right_boundary_node and get_left_boundary_node to obtain the neighbouring boundary nodes. If we find that $\overrightarrow{p_kq}$ does intersect a neighbouring boundary edge, then we can search down adjacent boundary edges via search_down_adjacent_boundary_edges until we either find an edge that $q$ is on, or until we identify that $q$ is outside of the triangulation – this function assumes the domain is convex, as triangulations being built are. If $q$ is outside of the triangulation, then we can use the exterior variant of the jump and march algorithm, exterior_jump_and_march, to find the ghost triangle containing $q$. This functon simply rotates around the boundary until we find two ghost edges enclosing the point $q$. Now let us assume we did not find a neighbouring boundary edge that intersects $\overrightarrow{p_kq}$. When this happens, we need to check the neighbouring interior edges, done via check_for_intersections_with_interior_edges_adjacent_to_boundary_node, which just rotates around the edges until we find an intersection. If we have still not found any intersection then, again assuming convexity, $q$ must be outside of the triangulation and so we use exterior_jump_and_march.

Now, if the algorithm is still going, then we need to start marching along the triangulation from $p_k$ towards $q$. The idea is to keep marching, keeping $p_i$ and $p_j$ to the left and right of $\overrightarrow{p_kq}$, respectively, until we find a case where $p_ip_jq$ is no longer a positively oriented triangle. When this happens, it must mean that we have passed $q$, and so we have found the triangle.

Let us describe how this marching is done in more detail. First, we need to be careful of boundary indices, first checking for an outer boundary index. If we have found an outer boundary index, then we have marched into the boundary, and so $q$ will be outside of the domain, meaning we go into exterior_jump_and_march. If this is not the case, then we march using get_adjacent to step onto the next triangle from a given edge $e_{ij}$. Then, checking the positions of the new points relative to $\overrightarrow{p_kq}$ and rearranging accordingly, we can step forward. We keep doing this until we get a negatively oriented triangle $p_ip_jq$. Unfortunately, this loop could terminate even if $q$ is not in the found triangle, which can occasionally happen if $p_ip_jq$ is a degenerate triangle. In this case, we just restart the jump and march algorithm. This latter worry is a very rare concern and does not alter the runtime in any significant manner.