Difficulty with scipy.spatial.Delaunay to find all neighbors of a given point

Question

There is a very good discussion on finding the nearest neighbors of a point among points using scipy.spatial.Delaunay here: How to find all neighbors of a given point in a delaunay triangulation using scipy.spatial.Delaunay?

I followed the answers, but I have difficulty when the symmetry of the configuration is high. Here is my code:

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

##########################################################

# Gets a Delaunay triangulations and returns the 
# nearest neighbor points for each point

def find_neighbors(tri):    
    l = tri.vertex_neighbor_vertices
    neib = []
    for i in range(len(l[0])-1):
        neib.append(list(l[1][l[0][i]:l[0][i+1]]))
    return neib

# Create a square lattice
t1 = np.array([1,0])
t2 = np.array([0,1])
n = 6
lattice = np.zeros((n**2,2))
for i in range(n):
    for j in range(n):
        lattice[i*n + j] = i* t1 + j * t2


tri = Delaunay(lattice)
Neib = find_neighbors (tri)

pindex = 20
plt.title("particle index: 20")
plt.triplot(lattice[:,0], lattice[:,1], tri.simplices.copy(),"c--")
plt.plot(lattice[:,0], lattice[:,1], 'go')
plt.plot(lattice[pindex,0], lattice[pindex,1], 'bo')
plt.plot([lattice[i,0] for i in Neib[pindex]],[lattice[i,1] for i in Neib[pindex]], 'ro')
plt.show()

pindex = 21
plt.title("particle index: 21")
plt.triplot(lattice[:,0], lattice[:,1], tri.simplices.copy(),"c--")
plt.plot(lattice[:,0], lattice[:,1], 'go')
plt.plot(lattice[pindex,0], lattice[pindex,1], 'bo')
plt.plot([lattice[i,0] for i in Neib[pindex]], [lattice[i,1] for i in Neib[pindex]], 'ro')
plt.show()

from scipy.spatial import Voronoi, voronoi_plot_2d

vor = Voronoi(lattice)
fig = voronoi_plot_2d(vor)
plt.show()

and the results are below. As you see for two points with identical symmetry I get a different number of neighbors, due to the way scipy does Delaunay triangulation. Is there any code in python to obtain the neighbors through scipy.spatial.Voronoi instead? I appreciate your help.

Here is the link to the images. My reputation is less than 10 and I cannot post images right now. (alternatively, you can run the code to get the same results)

index 20

https://i.ibb.co/9bQQqmT/download.png

index 21

https://i.ibb.co/pfw5hsG/download-1.png

Voronoi diagram

https://i.ibb.co/g3NpjVL/download-2.png

score 2 · Accepted Answer · edited Jun 20 '20 at 09:12

The problem in these symmetric configuration is that the Delaunay triangulation is not uniquely defined: whenever there are four vertices on a common circle, the Voronoi diagram does not specify how it should be triangulated. So opposite points may or may not be connected in the triangulation.

Here is an exaggerated picture of the Voronoi diagram with inflated the zero length edges in the Voronoi diagram.

The red circled edge in the Voronoi diagram is actually zero length and may or may not actually appear and correspond to an edge in the Delaunay triangulation.

You probably want to just ignore neighbors associated with these degenerate Voronoi edges (so that the neighbors of all vertices are those in your image here). However, scipy .spatial.Voronoi doesn't make this really easy: you can find short Voronoi edges but not the corresponding neighboring Voronoi sites. In the Delaunay triangulation, you can look at the two triangles connected to an edge: between those two trianlges, there are four vertices. If all four vertices are cocircular (or very nearly cocircular), then the associated Voronoi edge is degenerate and that neighbor can be ignored.

Difficulty with scipy.spatial.Delaunay to find all neighbors of a given point

1 Answers1