SciPy Delaunay三角剖分改变了单纯形的多个点，以实现参数的微小变化

Posted 2023-03-16

【中文标题】SciPy Delaunay三角剖分改变了单纯形的多个点，以实现参数的微小变化

【英文标题】：SciPy Delaunay triangulation changes multiple points of simplex for tiny change in parameters

【发布时间】：2021-11-04 17:05:11

【问题描述】：

我正在使用 Delaunay 三角剖分法在一些函数值中插值，这些函数值是在常规 4 维网格上的一组参数处评估的。有时，当参数值发生少量变化并将其带到新的单纯形时，单纯形中不止一个点会发生变化。我希望随着我不断改变参数之一，我会通过一次仅更改单纯形中的一个点来从单纯形转移到单纯形（我的代码中通常也是这种情况）。相反，请考虑以下脚本：

import numpy as np
from scipy.spatial import Delaunay

# hideous construction to get the desired 4d grid of points
# with points at [-1, -0.5, 0, 0.5, 1] along each axis
X = np.vstack(list(map(np.ravel, np.meshgrid(*[np.linspace(-1, 1, 5) for i in range(4)])))).T

tri = Delaunay(X)

delta = 1e-10

print(np.sort(tri.vertices[tri.find_simplex([-0.25, -0.05, 0.5+delta, 0.1])]))
print(np.sort(tri.vertices[tri.find_simplex([-0.25, -0.05, 0.5-delta, 0.1])]))

产生

[192 292 317 318 322]
[167 292 293 313 317]

请注意，这两个单纯形相差 3 个点，我期望有一个，而且我还没有设计一个 2-D 或 3-D 示例，其中不止一个顶点会发生变化。

我 99% 确定这是因为我的分数在常规网格上，但我找不到关于为什么或如何避免问题的详细答案。我知道三角测量不是独特的，但这从根本上不是问题。各种技巧似乎改变了我遇到此问题的位置，但我还没有找到防止问题出现在任何地方的“修复”。

编辑

我已经设法找到了一个 3D 示例，这让我可以将问题可视化。

import numpy as np
from scipy.spatial import Delaunay

X = np.vstack(list(map(np.ravel, np.meshgrid(*[np.linspace(-1, 1, 5) for i in range(3)])))).T
tri = Delaunay(X)

delta = 1e-6
x = np.array([-0.25, 0, 0.07])

fig = pl.figure()
ax = fig.add_subplot(projection='3d')

ax.scatter(*x, c='k')

x[1] = delta
s = X[tri.vertices[tri.find_simplex(x)]]

for i, si in enumerate(s):
    for j, sj in enumerate(s[i:]):
        ax.plot3D(*np.vstack([si, sj]).T, c='C0')

x[1] = -delta
s = X[tri.vertices[tri.find_simplex(x)]]
for i, si in enumerate(s):
    for j, sj in enumerate(s[i:]):
        ax.plot3D(*np.vstack([si, sj]).T, c='C1')

ax.set_xlabel('$x$')
ax.set_ylabel('$y$')
ax.set_zlabel('$z$')

这里是输出的两个方面。

蓝色和橙色的单纯形包含在它穿过 y=0（从正到负）之前和之后的插值点。我的假设是两个单纯形在 y=0 平面上具有相同的三角形面，但显然这是不正确的。在这种退化的情况下，这是 Delaunay 三角剖分的基础还是与实施有关？有没有办法避免它？ QHull 选项 Qx（在 SciPy 的 D>4 的默认选项中）似乎对这个示例有所帮助，但我不确定全局。

【参考方案1】：

您的问题并不是关于 scipy.spatial 中三角测量的实现。它更多地是关于 Delaunay triangulation 作为数学对象的数学运算。

D 维度中的 Delaunay 三角剖分非常明确，......当点是“一般位置”时。这意味着没有来自输入点的D+2 点位于公共球面上。如果发生这种情况，有人会说 Delaunay 三角剖分是“退化的”。当三角剖分退化时，Delaunay三角剖分定义不明确，在保留Delaunay性质的情况下，存在多种对点的凸包进行三角剖分的方法。

您观察到的是什么：您的点位于规则网格上，这是一个非常退化的点集（对于 Delaunay 属性）。坐标的任何细微修改都可以触发多个单纯形的翻转，从而恢复 Delaunay 属性。

也许您可以通过查看 Delaunay 三角剖分的对偶对象来理解这种行为：它的 Voronoi diagram。对于接近规则网格的点集，该图是退化的：它的 Voronoi 边的长度为零，或长度接近于零。点坐标的任何微小修改都可以改变 Voronoi 图的拓扑结构（因此也改变了 Delaunay 三角剖分）。

【讨论】：

我知道三角剖分是退化的，但是一旦计算出来，它就不应该改变（对吧？）。我的问题是，当我更改我插值的点时，我假设当我的插值点穿过单纯形的一个（超）平面时，包含的单纯形一次更改一个点。

啊，我刚刚（终于）设法在 3D 中找到了一个示例，现在我看到包含的单纯形可以一次更改一个点，我将在问题中添加这一点。