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

Question

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

我正在使用 python 科学计算堆栈(scipy、numpy、matplotlib)并且我有一组二维点，为此我使用 scipy.spatial.Delaunay 计算 Delaunay traingulation ( wiki ) 。

我需要编写一个函数，给定任何点 a，将返回所有其他点，这些点是 a 也是 a 的任何单纯形(即三角形)的顶点的顶点(三角剖分中 a 的邻居)。但是，scipy.spatial.Delaunay ( here ) 的文档非常糟糕，我无法终生理解单纯形是如何指定的，否则我会着手这样做.即使只是解释 Delaunay 输出中的 neighbors、vertices 和 vertex_to_simplex 数组的组织方式也足以让我继续前进。

非常感谢您的帮助。

Answer 1

我自己弄明白了，所以这里有一个解释，供 future 对此感到困惑的人。

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

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

就是这样!我确信上面的函数可以做一些优化，这正是我在几分钟内想到的。如果有人有任何建议，请随时发布。希望这对将来像我一样对此感到困惑的人有所帮助。