Làm thế nào để tìm thấy tất cả các nước láng giềng của một điểm nhất định trong một tam giác delaunay sử dụng scipy.spatial.Delaunay?

Question

Tôi đã tìm kiếm câu trả lời cho câu hỏi này nhưng không thể tìm thấy bất kỳ điều gì hữu ích.

Tôi đang làm việc với ngăn xếp máy tính khoa học python (scipy, numpy, matplotlib) và tôi có một tập hợp các điểm 2 chiều, mà tôi tính toán tra cứu Delaunay (wiki) sử dụng scipy.spatial.Delaunay.

tôi cần phải viết một hàm, đưa ra bất kỳ điểm a, sẽ trở lại tất cả các điểm khác mà là đỉnh của bất kỳ simplex (ví dụ: hình tam giác) mà a cũng là một đỉnh của (những người hàng xóm của a trong tam giác). Tuy nhiên, các tài liệu cho scipy.spatial.Delaunay (here) là khá xấu, và tôi không thể cho cuộc sống của tôi hiểu làm thế nào các đơn giản đang được xác định hoặc tôi sẽ đi về việc này. Ngay cả chỉ là một lời giải thích về cách các mảng neighbors, vertices và vertex_to_simplex trong đầu ra Delaunay được tổ chức sẽ đủ để giúp tôi tiếp tục.

Cảm ơn rất nhiều vì đã giúp đỡ.

Answer 1

Tôi đã tự mình tìm ra, do đó, đây là giải thích cho bất kỳ ai trong tương lai bị nhầm lẫn bởi điều này.

Như một ví dụ, chúng ta hãy sử dụng mạng đơn giản của các điểm mà tôi đã làm việc với trong mã của tôi, mà tôi tạo ra như sau

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 

chi tiết ở đây không thực sự quan trọng, đủ để nói nó tạo ra một tam giác đều đặn mạng, trong đó khoảng cách giữa một điểm và bất cứ sáu người hàng xóm gần nhất là 1.

Để vẽ nó

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

Bây giờ tính toán tam giác:

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

Hãy xem xét điều này cho chúng ta.

triang.points 

đầu ra:

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

đơn giản, chỉ là một mảng của tất cả chín điểm trong mạng được minh họa ở trên. Làm thế nào chúng ta hãy xem xét:

triang.vertices 

đầu ra:

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) 

Trong mảng này, mỗi hàng đại diện cho một simplex (tam giác) trong tam giác. Ba mục trong mỗi hàng là các chỉ số của các đỉnh của đơn vị đó trong mảng điểm mà chúng ta vừa thấy. Vì vậy, ví dụ như simplex đầu tiên trong mảng này, [4, 3, 6] gồm các điểm:

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

của nó dễ dàng nhận thấy điều này bằng cách vẽ các mạng trên một mảnh giấy, ghi nhãn mỗi điểm theo chỉ số của nó, và sau đó truy tìm qua mỗi hàng trong triang.vertices.

Đây là tất cả thông tin chúng tôi cần để viết hàm tôi đã chỉ định trong câu hỏi của mình. Có vẻ như:

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

Và đó là nó! Tôi chắc chắn rằng chức năng trên có thể làm với một số tối ưu hóa, nó chỉ là những gì tôi đã đưa ra trong một vài phút. Nếu bất cứ ai có bất cứ đề nghị, cảm thấy tự do để gửi chúng. Hy vọng rằng điều này sẽ giúp ai đó trong tương lai, những người đang bối rối về điều này như tôi.

Answer 2

Đây cũng là một phiên bản đơn giản một dòng câu trả lời của riêng James Porter bằng cách sử dụng danh sách hiểu:

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

tôi cần này quá và đã xem qua các câu trả lời sau đây. Nó chỉ ra rằng nếu bạn cần những người hàng xóm cho tất cả điểm ban đầu, đó là hiệu quả hơn để tạo ra một từ điển của các nước láng giềng trong một đi (ví dụ sau đây là dành cho 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 

Những người hàng xóm của điểm v sau đó là neighbors[v]. Đối với 10.000 điểm trong này mất 370ms để chạy trên máy tính xách tay của tôi. Có lẽ những người khác có ý tưởng về việc tối ưu hóa điều này hơn nữa?

Answer 4

Các phương pháp được mô tả ở trên, chu kỳ thông qua tất cả các đơn giản, có thể mất rất nhiều thời gian, trong trường hợp có một số lượng lớn các điểm. Một cách tốt hơn có thể là sử dụng Delaunay.vertex_neighbor_vertices, đã chứa tất cả thông tin về những người hàng xóm. Thật không may, giải nén thông tin

def find_neighbors(pindex, triang): 

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

Ví dụ dưới đây như thế nào để có được các chỉ số của một số đỉnh (số 17, trong ví dụ này):

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

Dưới đây là một ellaboration về câu trả lời @astrofrog. Điều này cũng hoạt động trong hơn 2D.

Mất khoảng 300 mili giây trên bộ 2430 điểm ở chế độ 3D (khoảng 16.000 đơn giản).

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