如何使用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

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

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

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

就这样！我相信上面的函数可以做一些优化，这正是我几分钟后想到的。如果有人有任何建议，请随时张贴。希望这能帮助将来像我一样困惑的人。