如何使用 scipy 的 Delaunay 找到两个不同集合中点的 Delaunay 邻居？

Posted 2023-03-16

【中文标题】如何使用 scipy 的 Delaunay 找到两个不同集合中点的 Delaunay 邻居？

【英文标题】：How to find the Delaunay neighbours of points in two different sets using scipy's Delaunay?

【发布时间】：2021-12-26 06:42:48

【问题描述】：

见https://docs.scipy.org/doc/scipy/reference/generated/scipy.spatial.Delaunay.html

考虑两组点。对于 X_ 中的每个点，我想在“点”中找到最近的 delaunay 邻居。我认为一种缓慢的方法是一次形成points 的Delaunay 三角剖分加上X_ 的一个点，然后以某种方式进行邻居查找。使用 scipy（或其他工具）是否有更有效的方法？

from pylab import *
from scipy.spatial import Delaunay

np.random.seed(1)

points = np.random.randn(10, 2)
X_ = np.random.randn(4, 2)

tri = Delaunay(points, incremental=True)

plt.triplot(points[:,0], points[:,1], tri.simplices, alpha=0.5)

plot(X_[:, 0], X_[:, 1], 'ro', alpha=0.5)

for i, (x, y) in enumerate(points):
    text(x, y, i)
    
for i, (x, y) in enumerate(X_):
    text(x, y, i, color='g')


tri.points, tri.points.shape

【参考方案1】：

理想情况下，最简单的解决方案可能如下所示：

Generate the Delaunay triangulation of "points"
For each point X in X_
    Add point X to the triangulation
    Get the neighbors of X in the augmented triangulation
    Remove point X from triangulation

使用scipy.spatial.Delaunay，您可以将incrementally add points 用于三角测量，但不是删除它们的机制。所以这种方法不能直接应用。

我们可以避免实际将新点添加到三角剖分中，而只需确定它将连接到哪些点而不实际更改三角剖分。基本上，我们需要识别Bowyer-Watson algorithm中的空腔，并收集空腔中三角形的顶点。如果建议的点位于其外接圆内，则 Delaunay 三角形是空腔的一部分。这可以从包含该点的三角形开始增量/局部构建，然后仅测试相邻的三角形。

这里是执行此算法的代码（针对单个测试点）。

import numpy as np
import matplotlib.pyplot as plt
from scipy.spatial import Delaunay

num_tri_points = 50

rng = np.random.default_rng(1)

# points in the Delaunay triangulation
points = rng.random((num_tri_points, 2))

# test point for identifying neighbors in the combined triangulation
X = rng.random((1, 2))*.8 + .1

tri = Delaunay(points)

# determine if point lies in the circumcircle of simplex of Delaunay triangulation tri...
# note: this code is not robust to numerical errors and could be improved with exact predicates
def in_circumcircle(tri, simplex, point):
    coords = tri.points[tri.simplices[simplex]]
    ax = coords[0][0]
    ay = coords[0][1]
    bx = coords[1][0]
    by = coords[1][1]
    cx = coords[2][0]
    cy = coords[2][1]
    d = 2 * (ax * (by - cy) + bx * (cy - ay) + cx * (ay - by))
    ux = ((ax * ax + ay * ay) * (by - cy) + (bx * bx + by * by) * (cy - ay) + (cx * cx + cy * cy) * (ay - by)) / d
    uy = ((ax * ax + ay * ay) * (cx - bx) + (bx * bx + by * by) * (ax - cx) + (cx * cx + cy * cy) * (bx - ax)) / d

    rad_sq = (ax-ux)*(ax-ux)+(ay-uy)*(ay-uy)
    point_dist_sq = (point[0]-ux)*(point[0]-ux)+(point[1]-uy)*(point[1]-uy)
    
    return point_dist_sq < rad_sq

# find the triangle containing point X
starting_tri = tri.find_simplex(X)

# remember triangles that have been tested so we don't repeat
tris_tested = set()

# collect the set of neighboring vertices in the potential combined triangulation
neighbor_vertices = set()

# queue triangles for performing the incircle check
tris_to_process = [ starting_tri[0] ]

while len(tris_to_process):        
    # get the next triangle
    tri_to_process = tris_to_process.pop()

    # remember that we have checked this triangle
    tris_tested.add(tri_to_process)

    # is the proposed point inside the circumcircle of the triangle
    if in_circumcircle(tri, tri_to_process, X[0,:]):
        # if so, add the vertices of this triangle as neighbors in the combined triangulation
        neighbor_vertices.update(tri.simplices[tri_to_process].flatten())
        # queue the neighboring triangles for processing
        for nbr in tri.neighbors[tri_to_process]:
            if nbr not in tris_tested:
                tris_to_process.append(nbr)

# plot the results    
plt.triplot(points[:,0], points[:,1], tri.simplices, alpha=0.7)
plt.plot(X[0, 0], X[0, 1], 'ro', alpha=0.5)
plt.plot(tri.points[list(neighbor_vertices)][:,0], tri.points[list(neighbor_vertices)][:,1], 'bo', alpha=0.7)
plt.show()

运行代码会给出以下图，其中包含原始三角剖分、测试点（红色）和识别出的相邻顶点（蓝色）。

【讨论】：

这很好，现在回到这个。我现在正在寻找 nD 中的实现（最快，所有技巧），即在圆周上。从一些快速阅读来看，我的猜测是快速库只是解决和硬编码（或使用 sympy）Cayley-Menger 矩阵的反转......请参阅此处的讨论。 math.stackexchange.com/questions/4056099/…

更新：这显然被称为 Tinfour 算法gwlucastrig.github.io/TinfourDocs/…

更新：见github.com/python-adaptive/adaptive/blob/master/adaptive/…

更新：我认为我们还必须区别对待船体外点的情况。可能只是用船体上的新点和两个相邻点创建 tris 并检查一些东西。我认为与向三角剖分添加一个点相同的操作。