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

`DelaunayTriangulation.brute_force_search`

— Function`brute_force_search(tri::Triangulation, q; itr = each_triangle(tri)`

Returns the triangle in `tri`

containing `q`

using brute force search.

## The Main Method

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

`DelaunayTriangulation.jump_and_march`

— Function```
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.

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`

.

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

— Function```
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)`

.

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

— Function`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_point`

— Function```
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_node`

— Function`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 `k`

th 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`k`

th 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_edges`

— Function`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_edges`

— Function```
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`

.

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

`DelaunayTriangulation.check_for_intersections_with_interior_edges_adjacent_to_boundary_node`

— Function```
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 `k`

th point of `pts`

to `q`

intersects any of the edges neighbouring the boundary node `k`

.

**Arguments**

`tri`

: The`Triangulation`

.`k`

: The boundary node.`q`

: The point we are searching for.`right_cert`

: A certificate giving the position of`q`

to the right of the`k`

th point. This comes from`check_for_intersections_with_adjacent_boundary_edges`

.`left_cert`

: A certificate giving the position of`q`

to the left of the`k`

th point. This comes from`check_for_intersections_with_adjacent_boundary_edges`

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

.`bnd_idx=I(BoundaryIndex)`

: The boundary index for the boundary that`k`

is on.

**Outputs**

There are several possible forms for the returned values. These are listed below, letting `p`

be the `k`

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

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.PointLocationHistory`

— Type`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.