我一直在寻找这个问题的答案，但找不到任何有用的东西。
我正在使用 python 科学计算堆栈（scipy、numpy、matplotlib），我有一组二维点，我使用scipy.spatial.Delaunay.
我需要编写一个函数，给定任何点a，将返回所有其他点，这些点是任何单纯形（即三角形）a的顶点，也是（a三角剖分中的邻居）的顶点。但是，scipy.spatial.Delaunay（here）的文档非常糟糕，我一生都无法理解单纯形是如何指定的，否则我会继续这样做。即使只是解释 Delaunay 输出中的neighbors,vertices和vertex_to_simplex数组是如何组织的，也足以让我继续前进。
非常感谢您的帮助。
我自己想通了，所以这里有一个解释给任何对此感到困惑的未来人。
例如，让我们使用我在代码中使用的简单点阵，生成如下
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))
就是这样！我确信上面的函数可以做一些优化，这正是我在几分钟内想出的。如果有人有任何建议，请随时发布。希望这可以帮助将来像我一样对此感到困惑的人。
上述方法循环遍历所有单纯形，如果有大量点，这可能需要很长时间。更好的方法可能是使用 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()
我知道自从提出这个问题以来已经有一段时间了。但是，我只是遇到了同样的问题并想出了如何解决它。vertex_neighbor_vertices只需使用您的 Delaunay 三角剖分对象（让我们称其为“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 的点的相邻顶点存储在我命名为“索引”的第二个数组中的某个位置。但是哪里？这是第一个数组（我称为“index_pointers”）的来源。起始位置（对于第二个数组“索引”）是 index_pointers[vertex_id]，经过相关子数组的第一个位置是 index_pointers[vertex_id+1 ]。所以解决方案是 indices[index_pointers[vertex_id]:index_pointers[vertex_id+1]]
这是@astrofrog 答案的详细说明。这也适用于二维以上。
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
这也是詹姆斯波特自己使用列表理解的答案的简单单行版本：
find_neighbors = lambda x,triang: list(set(indx for simplex in triang.simplices if x in simplex for indx in simplex if indx !=x))
我也需要这个，并遇到了以下答案。事实证明，如果您需要所有初始点的邻居，一次生成邻居字典会更有效（以下示例适用于 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 毫秒。也许其他人有进一步优化的想法？
我们可以找到一个包含顶点 ( 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)
这里的所有答案都集中在获取一个点的邻居（除了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。
这是我的，在 2D 中 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)
ixs[0, spacing_indices[1: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