如何使用scipy.spatial.Delaunay在delaunay三角剖分中找到给定点的所有邻居？

Question

我一直在寻找这个问题的答案,但找不到任何有用的东西.

我正在使用python科学计算堆栈(scipy,numpy,matplotlib),我有一组2维点,我使用它来计算Delaunay traingulation(wiki)scipy.spatial.Delaunay.

我需要编写一个函数,给定任何一个点a,它将返回所有其他点,这些点是任何单形(即三角形)a的顶点,也是(a三角剖分中的邻居)的顶点.然而,scipy.spatial.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

我知道这个问题提出已经有一段时间了。然而，我刚刚遇到了同样的问题并想出了如何解决它。只需使用 Delaunay 三角剖分对象的（有点缺乏记录的）方法vertex_neighbor_vertices（让我们称之为“tri”）。它将返回两个数组：

    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”的第二个数组中的某个位置。但是哪里？这就是第一个数组（我称之为“index_pointers”）出现的地方。起始位置（第二个数组“indices”）是index_pointers[vertex_id]，经过相关子数组的第一个位置是index_pointers[vertex_id+1 ]。所以解决方案是indices[index_pointers[vertex_id]:index_pointers[vertex_id+1]]

Answer 3

上面描述的方法循环遍历所有单形,这可能需要很长时间,以防有大量的点.更好的方法可能是使用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()