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

Question

我一直在寻找答案，但找不到有用的信息。

我正在使用Python科学计算库（scipy、numpy、matplotlib）处理一组二维点，并使用scipy.spatial.Delaunay计算Delaunay三角剖分(wiki)。

我需要编写一个函数，给定任何点a，将返回作为任何简单形式（即三角形）的顶点之一的点（即a在三角剖分的相邻点）。但是，对于scipy.spatial.Delaunay的文档(这里)很差，我无法理解简单的指定方式或者如何操作. 即使只是解释Delaunay输出中neighbors、vertices和vertex_to_simplex数组的组织方式，也足以让我开始工作。

非常感谢您的帮助。

- James Porter

---

没关系，我自己解决了！ - James Porter

在 Stack Overflow 上，我们支持人们回答自己的问题。你能否花点时间回答自己的问题并将其标记为已解决（通过在答案帖子左侧打勾）？ - Sicco

显然，声望低于10的用户在发布问题后8小时内不能回答自己的问题 :/ 我会把我写的内容保存在一个文本文件中，等到今晚再发布。 - James Porter

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

就是这样！我相信上面的函数可能需要一些优化，这只是我在几分钟内想出来的。如果有人有任何建议，请随意发布。希望这能帮助将来像我一样困惑的人。

- James Porter

---

这是一个很好的答案，绝对值得与使用vertex_neighbor_vertices的答案并列。我正在使用一组修改后的simplex，所以另一个解决方案对我不起作用。我建议编辑掉这个问题的设置，只保留最后的函数和一些解释就足够了。 - Temba

Answer 2

上述方法会循环遍历所有的简单形体，如果有大量点可能会花费很长时间。更好的方法可能是使用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()

- user2535797

---

至少对我来说，这个find_neighbors函数返回的点集范围从0到大约7或8（对于2D点）。我只期望1到3个邻居。 - user48956

Answer 3

我知道这个问题提出已经有一段时间了。然而，我最近遇到了同样的问题并找到了解决方法。只需使用您的Delaunay三角剖分对象（我们称其为“tri”）的（文档不太好的）方法vertex_neighbor_vertices。它将返回两个数组：

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

- ThommyTheCat

---

谢谢！还有另一个答案提到了使用 vertex_neighbour_vertices，但没有给出解释。你的答案正是我需要的 :) - bzm3r

Answer 4

以下是对@astrofrog答案的详细说明。这也适用于2D以上的情况。

在3D的2430个点集上（大约16000个simplex），需要约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

- glyg

Answer 5

以下是使用列表推导式简化的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))

- buergi

Answer 6

我也需要这个，然后找到了下面的答案。原来，如果你需要所有初始点的邻居，最高效的方法是一次性生成一个邻居字典（以下示例为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毫秒。也许其他人有进一步优化的想法？

- astrofrog

---

看我的答案，我已经测试了1000个点，我的解决方案比之前快了6倍（6毫秒对比32毫秒）。 - johnbaltis

Answer 7

这里所有的答案都集中在获取一个点的邻居上（除了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]

这适用于任何维度，找到所有点的所有邻居的解决方案比只找到一个点的邻居（James Porter 的期望答案）更快。

- johnbaltis

---

注意读者，你也可以直接使用Scipy稀疏矩阵并使用np.argwhere。不确定速度如何。 - mathtick

Answer 8

这是我的代码，对于一个包含11000个二维点的云，大约需要30毫秒的时间。

它会返回一个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)
    np.add.at(ixs[0], spacing_indices[1:int(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

- Peter

---

原始答案存在一个鲜有出现的阴险错误。我已经更新了解决方案，用np.add.at(ixs[0], spacing_indices[1:int(np.argmax(spacing_indices))], 1)代替了ixs[0, spacing_indices[1:int(np.argmax(spacing_indices))]] = 1，这样就解决了这个错误。 - Peter

Answer 9

我们可以找到一个包含顶点 (tri.vertex_to_simplex[vertex]) 的单形，然后递归搜索该单形的邻居 (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)

- Pragati

---

range(len(points)) 我觉得是这样的？ - mathtick