私はこの質問に対する答えを探していましたが、役に立つものは何も見つかりませんでした。
私はPythonの科学計算スタック（scipy、numpy、matplotlib）を使用しており、2次元のポイントのセットがあり、を使用してDelaunayのtraingulation（wiki）を計算しますscipy.spatial.Delaunay。
任意の点が与えられると、（三角形分割の隣接a点）の頂点でもあるシンプレックス（つまり三角形）の頂点である他のすべての点を返す関数を作成する必要があります。ただし、（ここ）のドキュメントはかなり悪いので、シンプレックスがどのように指定されているかを一生理解することはできません。Delaunay出力の、および配列がどのように編成されているかを説明するだけでも、私は十分に理解できます。aascipy.spatial.Delaunayneighborsverticesvertex_to_simplex
助けてくれてありがとう。
私は自分でそれを理解したので、これに混乱している将来の人のためにここに説明があります。
例として、コードで使用していた単純なポイントのラティスを使用してみましょう。これは次のように生成されます。
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
ここでの詳細はそれほど重要ではありません。ある点とその6つの最近傍のいずれかとの間の距離が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]])
単純で、上に示した格子内の9つのポイントすべての配列です。どのように見てみましょう：
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)
この配列では、各行は三角形分割の1つのシンプレックス（三角形）を表します。各行の3つのエントリは、先ほど見たポイント配列内のそのシンプレックスの頂点のインデックスです。したがって、たとえば、この配列の最初のシンプレックス[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))
以上です！上記の関数は、いくつかの最適化で実行できると確信しています。これは、私が数分で思いついたものです。誰か提案があれば、遠慮なく投稿してください。うまくいけば、これは私と同じようにこれについて混乱している将来の誰かを助けるでしょう。
上記の方法は、ポイントが多数ある場合に非常に時間がかかる可能性があるすべてのシンプレックスを循環します。より良い方法は、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()
この質問が出されてからしばらく経ちました。しかし、私は同じ問題を抱えていて、それを解決する方法を見つけました。vertex_neighbor_verticesDelaunay三角形分割オブジェクトの（文書化がやや不十分な）メソッドを使用するだけです（これを「tri」と呼びます）。2つの配列を返します。
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
インデックスvertex_idを持つポイントに隣接する頂点は、「indices」という名前の2番目の配列のどこかに格納されます。しかしここで？ここで、最初の配列（「index_pointers」と呼びます）が入ります。開始位置（2番目の配列「indexes」の場合）はindex_pointers [vertex_id]で、関連するサブ配列の先の最初の位置はindex_pointers [vertex_id+1です。 ]。したがって、解決策はindexes [index_pointers [vertex_id]：index_pointers [vertex_id +1]]です。
これが@astrofrogの回答の詳細です。これは2D以上でも機能します。
3Dの2430ポイントのセットで約300ミリ秒かかりました（約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
リスト内包表記を使用したJamesPorter自身の回答の単純な1行バージョンもあります。
find_neighbors = lambda x,triang: list(set(indx for simplex in triang.simplices if x in simplex for indx in simplex if indx !=x))
私もこれが必要で、次の答えに出くわしました。すべての初期点にネイバーが必要な場合は、ネイバーのディクショナリを一度に作成する方がはるかに効率的であることがわかります（次の例は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ミリ秒かかります。たぶん他の人はこれをさらに最適化するアイデアを持っていますか？
ここでのすべての答えは、1つのポイントのネイバーを取得することに焦点を当てています（ただし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]
これはどの次元でも機能し、すべての点のすべての近傍を見つけるこのソリューションは、 1つの点の近傍のみを見つけるよりも高速です（の例外的な答えJames Porter）。
頂点（tri.vertex_to_simplex[vertex]）を含む1つのシンプレックスを見つけてから、このシンプレックス（）の近傍を再帰的に検索してtri.neighbors、頂点を含む他のシンプレックスを見つけることができます。
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)
これが私のものです。2Dで11000ポイントのクラウドで約30msかかります。
これにより、インデックスの2xP配列が得られます。ここで、Pは存在するネイバーのペアの数です。
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)
ixs[0, spacing_indices[1: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