scipy.spatial.Delaunayを使って、delaunay三角測量で与えられた点のすべての近所を見つけるには？

Question

私はこの質問に対する答えを探していますが、何か役に立つものは見つけられません。

Iは、Python科学計算スタック（scipyのダウンロード、numpyの、matplotlibの）で働いていると私はscipy.spatial.Delaunayを用いドロネーtraingulation（wiki）を計算対象の2次元の点の集合を有します。私は任意の点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 

詳細を以下のように、それは正三角形を生成言えば十分のは、私が発生、私は私のコードにして働いていた点を、単純な格子を使用してみましょう格子点間の距離と6つの最近傍のいずれか

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

それをプロットする1

あります

は今、三角測量の計算：

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) 

、各行は、三角測量の一方シンプレックス（三角形）を表します。各行の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)) 

これだけです！私は上記の関数がいくつかの最適化を行うことができると確信しています。ちょうど私が数分間で思いついたものです。誰かが何か提案があれば、投稿してください。うまくいけば、これは私のようにこれについて混乱している将来の誰かに役立ちます。ここで

Answer 2

もリストの内包表記を使用してジェームズ・ポーターの自身の答えのシンプルな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)) 

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]です。これで10,000ポイントは私のラップトップで370ミリ秒かかる。たぶん、これをさらに最適化するアイデアは他の人にありますか？

Answer 4

上記の方法は、多くの点がある場合には、非常に長くかかる可能性があるすべてのシンプルを循環します。もっと良い方法は、Delaunay.vertex_neighbor_verticesを使用することです。これにはすでに近隣に関するすべての情報が含まれています。残念ながら、（この例では数17、）次のコードは、いくつかの頂点のインデックスを取得する方法を示した情報

def find_neighbors(pindex, triang): 

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

を抽出：

ここ

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です。これは2次元以上でも機能します。

3Dの2430ポイント（約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