从不规则间隔点进行快速 2D 插值的实现和策略

Posted 2023-02-16

【中文标题】从不规则间隔点进行快速 2D 插值的实现和策略

【英文标题】：Implementations and strategies for fast 2D interpolation from irregularly spaced points

【发布时间】：2021-01-22 14:30:02

【问题描述】：

鉴于二维中大量（约 1000 万）不规则间隔的点，其中每个点都有一些与之相关的强度（“权重”），现有的 python 实现有哪些用于在以下位置插入值：

某个随机位置的特定点（即point = (0.5, 0.8)） 随机位置的大量点（即points = np.random.random((1_000_000, 2))） 整数位置的规则网格（即np.indices((1000, 1000)).T）

我知道 Delaunay 三角测量经常用于此目的。有没有其他方法可以这样做？

是否有任何解决方案利用多个 CPU 内核或 GPU？

例如，这里是使用scipy's LinearNDInterpolator 的方法。它似乎没有使用多个 CPU 内核。 还有other options in scipy，但是对于这个问题，我特别有兴趣了解除了 scipy 中的解决方案之外的其他解决方案。

# The %time tags are IPython magic functions that time that specific line
dimension_shape = (1000, 1000) # we spread the random [0-1] over [0-1000] to avoid floating point errors
N_points = dimension_shape[0] * dimension_shape[1]

known_points = np.random.random((N_points, 2)) * dimension_shape 
known_weights = np.random.random((N_points,))

unknown_point = (0.5, 0.8)
unknown_points = np.random.random((N_points, 2)) * dimension_shape
unknown_grid = np.indices(dimension_shape, dtype=float).T.reshape((-1, 2)) # reshape to a list of 2D points

%time tesselation = Delaunay(known_points) # create grid to know neighbours # 6 sec
%time interp_func = LinearNDInterpolator(tesselation, known_weights) # 1 ms
%time interp_func(unknown_point) # 2 sec # run it once because the scipy function needs to compile 

%time interp_func(unknown_point) # ~ns
%time interp_func(unknown_grid) # 400 ms
%time interp_func(unknown_points) # 1 min 13 sec

# Below I sort the above `unknown_points` array, and try again
%time ind = np.lexsort(np.transpose(unknown_points)[::-1]) # 306 ms
unknown_points_sorted = unknown_points[ind].copy()

%time interp_func(unknown_points_sorted) # 19 sec <- much less than 1 min!

在上面的代码中，花费大量时间的事情是构建 Delaunay 网格，以及在非常规的点网格上进行插值。请注意，首先对非常规点进行排序会显着提高速度！

不要觉得有必要从一开始就给出完整的答案。欢迎解决上述任何方面的问题。

【参考方案1】：

Scipy 非常好，我认为 Python 中没有更好的解决方案，但我可以添加一些可能对您有帮助的东西。首先，您对分数进行排序的想法非常好。所谓的“增量算法”通过一次插入一个顶点来构建 Delaunay。在现有网格中插入顶点的第一步是确定将其插入网格中的哪个三角形。为了加快速度，一些算法会从最近插入的位置开始搜索。因此，如果您的点是有序的，以便插入的每个点都相对接近前一个点，则搜索速度会快得多。如果您想了解更多详细信息，可以查找“Lawson's Walk”算法。在我自己的 Delaunay 实现中（它是用 Java 编写的，所以恐怕它对你没有帮助），我有一个基于希尔伯特空间填充曲线的排序。希尔伯特排序效果很好。但即使只是按 x/y 坐标排序也是有帮助的。

关于是否有其他方法可以在不使用 Delaunay 的情况下进行插值...您可以尝试使用反距离加权 (IDW)。 IDW 技术不需要 Delaunay，但它们确实需要某种方法来确定哪些顶点接近您希望插值的点。我已经将我的坐标空间划分为均匀间隔的 bin，将顶点存储在适当的 bin 中，然后通过查看相邻的 bin 来拉出我需要进行插值的点。它可能需要大量编码，但它会比 Delaunay 相当快并且使用更少的内存

【参考方案2】：

在 Delaunay 三角形上进行插值当然是一种可能性，但我建议对 kD-tree 中的点进行排序，使用树查询最近的邻居（在足够的半径内），然后使用 IDW 进行插值，就像以前一样建议。