How to find all neighbors of a given point in a delaunay triangulation using scipy.spatial.Delaunay?

Question

I have been searching for an answer to this question but cannot find anything useful.

I am working with the python scientific computing stack (scipy,numpy,matplotlib) and I have a set of 2 dimensional points, for which I compute the Delaunay traingulation (wiki) using scipy.spatial.Delaunay.

I need to write a function that, given any point a, will return all other points which are vertices of any simplex (i.e. triangle) that a is also a vertex of (the neighbors of a in the triangulation). However, the documentation for scipy.spatial.Delaunay (here) is pretty bad, and I can't for the life of me understand how the simplices are being specified or I would go about doing this. Even just an explanation of how the neighbors, vertices and vertex_to_simplex arrays in the Delaunay output are organized would be enough to get me going.

Much thanks for any help.

Omsk

Nevermind, I figured it out myself! – Omsk

On Stack Overflow we support people to answer their own questions. Could you take the effort to answer your own question and mark it as solved (by ticking the check mark at the left of your answer post)? – Superstar

tried to, apparently users w/ less than 10 reputation can't answer their own questions for 8 hours after posting :/ I'll save what I wrote in a txt file and wait until this evening to post it. – Omsk

Answer 1

I figured it out on my own, so here's an explanation for anyone future person who is confused by this.

As an example, let's use the simple lattice of points that I was working with in my code, which I generate as follows

import numpy as np
import itertools as it
from matplotlib import pyplot as plt
import scipy as sp

inputs = list(it.product([0,1,2],[0,1,2]))
i = 0
lattice = range(0,len(inputs))
for pair in inputs:
    lattice[i] = mksite(pair[0], pair[1])
    i = i +1

Details here not really important, suffice to say it generates a regular triangular lattice in which the distance between a point and any of its six nearest neighbors is 1.

To plot it

plt.plot(*np.transpose(lattice), marker = 'o', ls = '')
axes().set_aspect('equal')

Now compute the triangulation:

dela = sp.spatial.Delaunay
triang = dela(lattice)

Let's look at what this gives us.

triang.points

output:

array([[ 0.        ,  0.        ],
       [ 0.5       ,  0.8660254 ],
       [ 1.        ,  1.73205081],
       [ 1.        ,  0.        ],
       [ 1.5       ,  0.8660254 ],
       [ 2.        ,  1.73205081],
       [ 2.        ,  0.        ],
       [ 2.5       ,  0.8660254 ],
       [ 3.        ,  1.73205081]])

simple, just an array of all nine points in the lattice illustrated above. How let's look at:

triang.vertices

output:

array([[4, 3, 6],
       [5, 4, 2],
       [1, 3, 0],
       [1, 4, 2],
       [1, 4, 3],
       [7, 4, 6],
       [7, 5, 8],
       [7, 5, 4]], dtype=int32)

In this array, each row represents one simplex (triangle) in the triangulation. The three entries in each row are the indices of the vertices of that simplex in the points array we just saw. So for example the first simplex in this array, [4, 3, 6] is composed of the points:

[ 1.5       ,  0.8660254 ]
[ 1.        ,  0.        ]
[ 2.        ,  0.        ]

Its easy to see this by drawing the lattice on a piece of paper, labeling each point according to its index, and then tracing through each row in triang.vertices.

This is all the information we need to write the function I specified in my question. It looks like:

def find_neighbors(pindex, triang):
    neighbors = list()
    for simplex in triang.vertices:
        if pindex in simplex:
            neighbors.extend([simplex[i] for i in range(len(simplex)) if simplex[i] != pindex])
            '''
            this is a one liner for if a simplex contains the point we`re interested in,
            extend the neighbors list by appending all the *other* point indices in the simplex
            '''
    #now we just have to strip out all the dulicate indices and return the neighbors list:
    return list(set(neighbors))

And that's it! I'm sure the function above could do with some optimization, its just what I came up with in a few minutes. If anyone has any suggestions, feel free to post them. Hopefully this helps somebody in the future who is as confused about this as I was.

Omsk

This is a great answer and definitely one worth keeping alongside those that use vertex_neighbor_vertices. I am working with a set of modified simplices, so the other solution does not work for me. I would suggest editing out the setup for this problem though, the function at the end and a bit of explanation is sufficient. – Georgiannageorgianne

Answer 2

The methods described above cycle through all the simplices, which could take very long, in case there's a large number of points. A better way might be to use Delaunay.vertex_neighbor_vertices, which already contains all the information about the neighbors. Unfortunately, extracting the information

def find_neighbors(pindex, triang):

    return triang.vertex_neighbor_vertices[1][triang.vertex_neighbor_vertices[0][pindex]:triang.vertex_neighbor_vertices[0][pindex+1]]

The following code demonstrates how to get the indices of some vertex (number 17, in this example):

import scipy.spatial
import numpy
import pylab

x_list = numpy.random.random(200)
y_list = numpy.random.random(200)

tri = scipy.spatial.Delaunay(numpy.array([[x,y] for x,y in zip(x_list, y_list)]))

pindex = 17

neighbor_indices = find_neighbors(pindex,tri)

pylab.plot(x_list, y_list, 'b.')
pylab.plot(x_list[pindex], y_list[pindex], 'dg')
pylab.plot([x_list[i] for i in neighbor_indices],
           [y_list[i] for i in neighbor_indices], 'ro')    

pylab.show()

Erving

At least for me, this find_neighbors is returning sets of points that range from 0 to around 7 or 8 (for 2d points). I would expect 1..3 neighbors only. – Malformation

Answer 3

I know it's been a while since this question was posed. However, I just had the same problem and figured out how to solve it. Just use the (somewhat poorly documented) method vertex_neighbor_vertices of your Delaunay triangulation object (let us call it 'tri'). It will return two arrays:

    def get_neighbor_vertex_ids_from_vertex_id(vertex_id, tri):
        index_pointers, indices = tri.vertex_neighbor_vertices
        result_ids = indices[index_pointers[vertex_id]:index_pointers[vertex_id + 1]]
        return result_ids

The neighbor vertices to the point with the index vertex_id are stored somewhere in the second array that I named 'indices'. But where? This is where the first array (which I called 'index_pointers') comes in. The starting position (for the second array 'indices') is index_pointers[vertex_id], the first position past the relevant sub-array is index_pointers[vertex_id+1]. So the solution is indices[index_pointers[vertex_id]:index_pointers[vertex_id+1]]

Thine

Thank you! There is another answer which mentions using vertex_neighbour_vertices, but does not give an explanation. Your answer was just what the doctor ordered :) – Jollification

Answer 4

Here is an ellaboration on @astrofrog answer. This works also in more than 2D.

It took about 300 ms on set of 2430 points in 3D (about 16000 simplices).

from collections import defaultdict

def find_neighbors(tess):
    neighbors = defaultdict(set)

    for simplex in tess.simplices:
        for idx in simplex:
            other = set(simplex)
            other.remove(idx)
            neighbors[idx] = neighbors[idx].union(other)
    return neighbors

Keene

Answer 5

Here is also a simple one line version of James Porter's own answer using list comprehension:

find_neighbors = lambda x,triang: list(set(indx for simplex in triang.simplices if x in simplex for indx in simplex if indx !=x))

Speech

Answer 6

I needed this too and came across the following answer. It turns out that if you need the neighbors for all initial points, it's much more efficient to produce a dictionary of neighbors in one go (the following example is for 2D):

def find_neighbors(tess, points):

    neighbors = {}
    for point in range(points.shape[0]):
        neighbors[point] = []

    for simplex in tess.simplices:
        neighbors[simplex[0]] += [simplex[1],simplex[2]]
        neighbors[simplex[1]] += [simplex[2],simplex[0]]
        neighbors[simplex[2]] += [simplex[0],simplex[1]]

    return neighbors

The neighbors of point v are then neighbors[v]. For 10,000 points in this takes 370ms to run on my laptop. Maybe others have ideas on optimizing this further?

Watershed

See my answer, I've tested it for 1000 points and my solution is 6 times faster (6 ms vs 32 ms). – Amino

Answer 7

All the answers here are focused on getting the neighbors for one point (except astrofrog, but that is in 2D and this is 6x faster), however, it's equally expensive to get a mapping for all of the points → all neighbors.

You can do this with

from collections import defaultdict
from itertools import permutations
tri = Delaunay(...)
_neighbors = defaultdict(set)
for simplex in tri.vertices:
    for i, j in permutations(simplex, 2):
        _neighbors[i].add(j)

points = [tuple(p) for p in tri.points]
neighbors = {}
for k, v in _neighbors.items():
    neighbors[points[k]] = [points[i] for i in v]

This works in any dimension and this solution, finding all neighbors of all points, is faster than finding only the neighbors of one point (the excepted answer of James Porter).

Amino

Note to readers, you can also just jam in scipy sparse and use np.argwhere. Not sure of speed. – Euphorbiaceous

Answer 8

Here's mine, it takes around 30ms on a cloud of 11000 points in 2D.

It gives you a 2xP array of indices, where P is the number of pairs of neighbours that exist.

def get_delaunay_neighbour_indices(vertices: "Array['N,D', int]") -> "Array['2,P', int]":
    """
    Fine each pair of neighbouring vertices in the delaunay triangulation.
    :param vertices: The vertices of the points to perform Delaunay triangulation on
    :return: The pairs of indices of vertices
    """
    tri = Delaunay(vertices)
    spacing_indices, neighbours = tri.vertex_neighbor_vertices
    ixs = np.zeros((2, len(neighbours)), dtype=int)
    np.add.at(ixs[0], spacing_indices[1:int(np.argmax(spacing_indices))], 1)  # The argmax is unfortuantely needed when multiple final elements the same
    ixs[0, :] = np.cumsum(ixs[0, :])
    ixs[1, :] = neighbours
    assert np.max(ixs) < len(vertices)
    return ixs

Nils

The original answer had an insidious bug that showed up only rarely. I've updated the solution, which replaces ixs[0, spacing_indices[1:int(np.argmax(spacing_indices))]] = 1 with np.add.at(ixs[0], spacing_indices[1:int(np.argmax(spacing_indices))], 1), which fixes the bug. – Nils

Answer 9

We can find one simplex containing the vertex (tri.vertex_to_simplex[vertex]) and then recursively search the neighbors of this simplex (tri.neighbors) to find other simplices containing the vertex.

from scipy.spatial import Delaunay
tri = Delaunay(points)  #points is the list of input points
neighbors =[]           #array of neighbors for all vertices

for i in range(points):

    vertex = i            #vertex index
    vertexneighbors = []  #array of neighbors for vertex i
    neighbour1 = -1
    neighbour2=-1
    firstneighbour=-1
    neighbour1index = -1
    currentsimplexno= tri.vertex_to_simplex[vertex]
    for i in range(0,3):
        if (tri.simplices[currentsimplexno][i]==vertex):
            firstneighbour=tri.simplices[currentsimplexno][(i+1) % 3]
            vertexneighbors.append(firstneighbour)
            neighbour1index=(i+1) % 3
            neighbour1=tri.simplices[currentsimplexno][(i+1) % 3]
            neighbour2=tri.simplices[currentsimplexno][(i+2) % 3]
    while (neighbour2!=firstneighbour):
        vertexneighbors.append(neighbour2)
        currentsimplexno= tri.neighbors[currentsimplexno][neighbour1index]
        for i in range(0,3):
            if (tri.simplices[currentsimplexno][i]==vertex):
                neighbour1index=(i+1) % 3
                neighbour1=tri.simplices[currentsimplexno][(i+1) % 3]
                neighbour2=tri.simplices[currentsimplexno][(i+2) % 3]
    neighbors.append(vertexneighbors)
print (neighbors)

Gothurd

range(len(points)) I think? – Euphorbiaceous