如何使用scipy.spatial.Delaunay在delaunay三角測量中找到給定點的所有鄰居？

Question

我一直在尋找這個問題的答案，但找不到任何有用的東西。我使用python科學計算堆棧（scipy，numpy，matplotlib），我有一組2維的點，爲此我使用scipy.spatial.Delaunay計算Delaunay traingulation（wiki）。

我需要寫，因爲任何一點a功能，將返回的任何單一（即三角形）的頂點所有其他點即a也是一個頂點（的a在三角鄰居）。然而，scipy.spatial.Delaunay（here）的文檔非常糟糕，而且我不能理解單純的指定方式，或者我會這麼做。即使只是解釋Delaunay輸出中neighbors,vertices和vertex_to_simplex陣列的組織方式，也足以讓我繼續前進。

非常感謝您的幫助。

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

就是這樣！我相信上面的功能可以做一些優化，就是我在幾分鐘內想到的。如果有人有任何建議，請隨時發佈。希望這有助於未來的某些人對我和我一樣感到困惑。

Answer 2

這裏也是使用列表理解詹姆斯·波特自己的答案的一個簡單的行版本：

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

Answer 3

我需要這也和整個以下的答案來。事實證明，如果您需要所有初始點的鄰居，這是更有效的生產鄰居的字典一氣呵成（下面的例子是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]。這需要370分鐘才能在我的筆記本電腦上運行。也許其他人有進一步優化這個想法？

Answer 4

上面描述的方法循環遍歷所有的單形，這可能需要很長的時間，以防萬一有大量的點。更好的方法是使用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() 

Answer 5

這裏是@astrofrog答案的ellaboration。這也適用於2D以上。

在2430點的3D（大約16000個單體）上花了大約300毫秒。

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