Как найти всех соседей данной точки в триангуляции delaunay, используя scipy.spatial.Delaunay?

Ответ 1

Я сам это понял, так вот объяснение для любого будущего человека, который смущен этим.

В качестве примера позвольте использовать простую решетку точек, с которыми я работал в своем коде, который я генерирую следующим образом

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

Подробности здесь не очень важны, достаточно сказать, что он порождает правильную треугольную решетку, в которой расстояние между точкой и любым из ее шести ближайших соседей равно 1.

Чтобы построить его

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

Теперь вычислите триангуляцию:

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

Давайте посмотрим, что это дает нам.

triang.points

выход:

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]])

простой, всего лишь массив из всех девяти точек решетки, показанный выше. Как посмотреть на:

triang.vertices

выход:

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)

В этом массиве каждая строка представляет один симплекс (треугольник) в триангуляции. Три записи в каждой строке являются индексами вершин этого симплекса в массиве точек, который мы только что видели. Так, например, первый симплекс в этом массиве [4, 3, 6] состоит из точек:

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

Его легко увидеть, рисуя решетку на листе бумаги, маркируя каждую точку в соответствии с ее индексом, а затем прослеживая каждую строку в triang.vertices.

Это вся информация, которая нам нужна для написания функции, указанной в моем вопросе. Это выглядит как:

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))

И это! Я уверен, что вышеприведенная функция могла бы помочь с некоторой оптимизацией, и именно то, что я придумал за несколько минут. Если у кого есть предложения, не стесняйтесь публиковать их. Надеюсь, это поможет кому-то в будущем, кто так же запутался в этом, как и я.

Ответ 2

Описанные выше методы циклически проходят через все симплексы, которые могут занять очень много времени, если имеется большое количество точек. Лучше всего использовать Delaunay.vertex_neighbor_vertices, который уже содержит всю информацию о соседях. К сожалению, извлечение информации

def find_neighbors(pindex, triang):

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

Следующий код демонстрирует, как получить индексы некоторой вершины (номер 17 в этом примере):

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()

Ответ 3

Вот также простая однострочная версия собственного ответа Джеймса Портера с использованием понимания списка:

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

Ответ 4

Вот ответ на вопрос @astrofrog. Это работает также в более чем 2D.

Это заняло около 300 мс на множестве 2430 точек в 3D (около 16000 симплексов).

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

Ответ 5

Мне тоже было нужно это и наткнулся на следующий ответ. Оказывается, если вам нужны соседи для всех начальных точек, гораздо эффективнее производить словарь соседей за один раз (следующий пример для 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

Соседи точки v тогда neighbors[v]. Для 10 000 очков в этом требуется 370 мс для работы на моем ноутбуке. Может быть, у других есть идеи по оптимизации этого?

Ответ 6

Все ответы здесь сосредоточены на получении соседей по одной точке (за исключением astrofrog, но это в 2D, и это на 6 раз быстрее), однако одинаково дорого получить сопоставление для всех точек → всех соседей.

Вы можете сделать это с помощью

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]

Это работает в любом измерении, и это решение, находя все соседи всех точек, быстрее, чем поиск только соседей одной точки (исключенный ответ James Porter).

Ответ 7

Я знаю, что это было какое-то время, поскольку этот вопрос был задан. Однако у меня была одна и та же проблема и выяснили, как ее решить. Просто используйте (несколько плохо документированный) метод "vertex_neighbor_vertices" вашего объекта триангуляции Delaunay (назовем его "tri"). Он вернет два массива:

def get_neighbor_vertex_ids_from_vertex_id(vertex_id,tri):
    #use a less awful function name
    helper = tri.vertex_neighbor_vertices
    index_pointers = helper[0]
    indices = helper[1]
    result_ids = indices[index_pointers[vertex_id]:index_pointers[vertex_id+1]]

    return result_ids

Соседние вершины к точке с индексом vertex_id хранятся где-то во втором массиве, который я назвал "индексы". Но где? Здесь находится первый массив (который я назвал "index_pointers"). Начальная позиция (для индексов второго массива) - index_pointers [vertex_id], первая позиция за соответствующим подматрицей - index_pointers [vertex_id + 1 ]. Таким образом, решение является индексом [index_pointers [vertex_id]: index_pointers [vertex_id + 1]]