Find all simplices a point is a part of in scipy.spatial.Delaunay python

Question

Is there anyway to get all the simplices/triangles a certain point is a part of in a Delaunay triangulation using scipy.spatial.Delaunay?

I know there is the find_simplex() function, that that only returns 1 triangle that a point is a part of but I would like to get all triangles that it is a part of.

So in the example, when I do find_simplex() for point 6, it only returns triangle 2, but I would like it to return the triangles 1, 2, 3, 4, 10, and 9, as point 6 is a part of all of those triangles.

Any help would be appreciated!

Answer 1

You don’t want find_simplex because it is geometric, not topological. That is, it treats a point as a location rather than as a component of the triangulation: almost all points lie in only one simplex, so that’s what it reports.

Instead, use the vertex number. The trivial answer is to use the simplices attribute:

vert=6
[i for i,s in enumerate(d.simplices) if vert in s]

With a good bit more code, it is possible to search more efficiently using the vertex_to_simplex and neighbors attributes.

Answer 2

You can get all simplices adjacent to a given vertex efficiently by

def get_simplices(self, vertex):
    "Find all simplices this `vertex` belongs to"
    visited = set()
    queue = [self.vertex_to_simplex[vertex]]
    while queue:
        simplex = queue.pop()
        for i, s in enumerate(self.neighbors[simplex]):
            if self.simplices[simplex][i] != vertex and s != -1 and s not in visited:
                queue.append(s)
        visited.add(simplex)
    return np.array(list(visited))

Example:

import scipy.spatial
import numpy as np
np.random.seed(0)
points = np.random.rand(10, 2)
tri = scipy.spatial.Delaunay(points)
vertex = 2
simplices = get_simplices(tri, vertex)
# 0, 2, 5, 9, 11
neighbors = np.unique(tri.simplices[simplices].reshape(-1)])
# 0, 1, 2, 3, 7, 8

Visualisation:

import matplotlib.pyplot as plt
plt.triplot(points[:,0], points[:,1], tri.simplices)
plt.plot(points[neighbors,0], points[neighbors,1], 'or')
plt.plot(points[vertex,0], points[vertex,1], 'ob')
plt.show()