How to find the Delaunay neighbours of points in two different sets using scipy's Delaunay?

Question

See https://docs.scipy.org/doc/scipy/reference/generated/scipy.spatial.Delaunay.html

Consider two sets of points. For each point in X_, I would like to find the nearest delaunay neighbours in "points". I think a slow way is to form Delaunay triangulations of points plus a points from X_ one at a time and then do the neighbours lookup somehow. Is there a more efficient way of doing this using scipy (or another tool)?

from pylab import *
from scipy.spatial import Delaunay

np.random.seed(1)

points = np.random.randn(10, 2)
X_ = np.random.randn(4, 2)

tri = Delaunay(points, incremental=True)

plt.triplot(points[:,0], points[:,1], tri.simplices, alpha=0.5)

plot(X_[:, 0], X_[:, 1], 'ro', alpha=0.5)

for i, (x, y) in enumerate(points):
    text(x, y, i)
    
for i, (x, y) in enumerate(X_):
    text(x, y, i, color='g')


tri.points, tri.points.shape

Answer 1

Ideally, the simplest solution would probably look something like this:

Generate the Delaunay triangulation of "points"
For each point X in X_
    Add point X to the triangulation
    Get the neighbors of X in the augmented triangulation
    Remove point X from triangulation

Using scipy.spatial.Delaunay, you have the ability to incrementally add points to the triangulation but isn't a mechanism to remove them. So that approach cannot be applied directly.

We can avoid actually adding the new point to the triangulation and instead just determine which points it would be connected to without actually changing the triangulation. Basically, we need to identify the cavity in the Bowyer-Watson algorithm and collect the vertices of the triangles in the cavity. A Delaunay triangle is part of the cavity if the proposed point lies inside its circumcircle. This can be built incrementally / locally starting from the triangle containing the point and then only testing neighboring triangles.

Here the code to perform this algorithm (for a single test point).

import numpy as np
import matplotlib.pyplot as plt
from scipy.spatial import Delaunay

num_tri_points = 50

rng = np.random.default_rng(1)

# points in the Delaunay triangulation
points = rng.random((num_tri_points, 2))

# test point for identifying neighbors in the combined triangulation
X = rng.random((1, 2))*.8 + .1

tri = Delaunay(points)

# determine if point lies in the circumcircle of simplex of Delaunay triangulation tri...
# note: this code is not robust to numerical errors and could be improved with exact predicates
def in_circumcircle(tri, simplex, point):
    coords = tri.points[tri.simplices[simplex]]
    ax = coords[0][0]
    ay = coords[0][1]
    bx = coords[1][0]
    by = coords[1][1]
    cx = coords[2][0]
    cy = coords[2][1]
    d = 2 * (ax * (by - cy) + bx * (cy - ay) + cx * (ay - by))
    ux = ((ax * ax + ay * ay) * (by - cy) + (bx * bx + by * by) * (cy - ay) + (cx * cx + cy * cy) * (ay - by)) / d
    uy = ((ax * ax + ay * ay) * (cx - bx) + (bx * bx + by * by) * (ax - cx) + (cx * cx + cy * cy) * (bx - ax)) / d

    rad_sq = (ax-ux)*(ax-ux)+(ay-uy)*(ay-uy)
    point_dist_sq = (point[0]-ux)*(point[0]-ux)+(point[1]-uy)*(point[1]-uy)
    
    return point_dist_sq < rad_sq

# find the triangle containing point X
starting_tri = tri.find_simplex(X)

# remember triangles that have been tested so we don't repeat
tris_tested = set()

# collect the set of neighboring vertices in the potential combined triangulation
neighbor_vertices = set()

# queue triangles for performing the incircle check
tris_to_process = [ starting_tri[0] ]

while len(tris_to_process):        
    # get the next triangle
    tri_to_process = tris_to_process.pop()

    # remember that we have checked this triangle
    tris_tested.add(tri_to_process)

    # is the proposed point inside the circumcircle of the triangle
    if in_circumcircle(tri, tri_to_process, X[0,:]):
        # if so, add the vertices of this triangle as neighbors in the combined triangulation
        neighbor_vertices.update(tri.simplices[tri_to_process].flatten())
        # queue the neighboring triangles for processing
        for nbr in tri.neighbors[tri_to_process]:
            if nbr not in tris_tested:
                tris_to_process.append(nbr)

# plot the results    
plt.triplot(points[:,0], points[:,1], tri.simplices, alpha=0.7)
plt.plot(X[0, 0], X[0, 1], 'ro', alpha=0.5)
plt.plot(tri.points[list(neighbor_vertices)][:,0], tri.points[list(neighbor_vertices)][:,1], 'bo', alpha=0.7)
plt.show()

Running the code gives the following plot containing the original triangulation, the test point (red) and the identified neighboring vertices (blue).