如何使用scipy.space.delaunay在delaunay三角测量中找到给定点的所有邻居？

Question

我一直在寻找这个问题的答案，但找不到任何有用的答案。

我正在使用python科学计算堆栈（scipy、numpy、matplotlib），我有一组二维点，我使用scipy.spatial.Delaunay计算Delaunay训练规则（wiki）。

我需要编写一个函数，给定任何点a，它将返回所有其他点，这些点是任何单纯形（即三角形）的顶点，a也是（三角剖分中a的邻居）的顶点。然而，scipy.spatial.Delaunay（here）的文档非常糟糕，我一生都无法理解如何指定simplices，或者我将继续这样做。即使只是解释一下Delaunay输出中的neighbors、vertices和vertex_to_simplex数组是如何组织的，也足以让我开始工作。

非常感谢你的帮助。

Answer 1

上面描述的方法在所有的simplice中循环，如果有大量的点，可能需要很长的时间。更好的方法可能是使用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 2

我自己想出来的，所以这里有一个解释，任何未来的人谁是困惑这一点。

作为一个例子，让我们使用我在代码中使用的简单点格，我生成如下

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 3

以下是James Porter自己的简单单行本，使用列表理解：

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