如何将 scipy.interpolate.LinearNDInterpolator 与自己的三角测量一起使用答案

【问题标题】：How to use scipy.interpolate.LinearNDInterpolator with own triangulation如何将 scipy.interpolate.LinearNDInterpolator 与自己的三角测量一起使用
【发布时间】：2019-02-28 06:30:39
【问题描述】：

我有 my own triangulation algorithm，它根据 Delaunay 条件和梯度创建三角剖分，使三角形与梯度对齐。

这是一个示例输出：

以上描述与问题无关，但对于上下文是必要的。

现在我想用我的三角测量与scipy.interpolate.LinearNDInterpolator 进行插值。

使用 scipy 的 Delaunay，我将执行以下操作

import numpy as np
import scipy.interpolate
import scipy.spatial
points = np.random.rand(100, 2)
values = np.random.rand(100)
delaunay = scipy.spatial.Delaunay(points)
ip = scipy.interpolate.LinearNDInterpolator(delaunay, values)

这个delaunay 对象具有构成三角剖分的delaunay.points 和delaunay.simplices。我自己的三角测量有完全相同的信息，但scipy.interpolate.LinearNDInterpolator 需要scipy.spatial.Delaunay 对象。

我想我需要继承scipy.spatial.Delaunay 并实现相关方法。但是，我不知道我需要哪些才能到达那里。

【问题讨论】：

LinearNDInterpolator 也接受点数组作为其第一个参数。
我知道，但这会导致使用 scipy.spatial.Delaunay 进行新的三角测量，这不是我想要的。
啊，对。 youtube.com/watch?v=VrbybKWwb7c
我有同样的问题，我认为 Scipy 方面没有任何开发可以提供该功能。我试图创建一个虚拟的 Delaunay 对象并修改其属性的内容，但我得到了AttributeError: can't set attribute。从好的方面来说，如果你在 2D 中工作，Matplotlib 提供了这个功能。使用Triangulation 创建三角剖分对象并使用CubicTriInterpolator 进行插值。

标签： scipy scipy-spatial

【解决方案1】：

我想用 triangle package 提供的 Delaunay 三角测量来做同样的事情。在大（~100_000）点上，三角形 Delaunay 代码比 SciPy 代码快八倍。（我鼓励其他开发者尝试打败它:)）

不幸的是，Scipy LinearNDInterpolator 函数严重依赖于 SciPy Delaunay 三角剖分对象中存在的特定属性。这些是由_get_delaunay_info() CPython 代码创建的，很难反汇编。即使知道需要哪些属性（似乎有很多，包括paraboloid_scale 和paraboloid_shift 之类的属性），我不确定如何从不同的三角测量库中提取它。

相反，我尝试了@Patol75 的方法（对上述问题的评论），但使用LinearTriInterpolator 而不是三次方。代码运行正确，但比在 SciPy 中完成整个事情要慢。使用 matplotlib 代码从 400_000 个点的云中插值 400_000 个点的时间比 scipy 长约 3 倍。 Matplotlib tri code is written in C++，因此将代码转换为即 CuPy 并不简单。如果我们可以混合使用这两种方法，我们可以将总时间从 3.65 秒/10.2 秒减少到 1.1 秒！

import numpy as np
np.random.seed(1)
N = 400_000
shape = (100, 100)
points = np.random.random((N, 2)) * shape # spread over 100, 100 to avoid float point errors
vals = np.random.random((N,))
interp_points1 = np.random.random((N,2)) * shape
interp_points2 = np.random.random((N,2)) * shape

triangle_input = dict(vertices=points)

### Matplotlib Tri
import triangle as tr
from matplotlib.tri import Triangulation, LinearTriInterpolator

triangle_output = tr.triangulate(triangle_input) # 280 ms
tri = tr.triangulate(triangle_input)['triangles'] # 280 ms
tri = Triangulation(*points.T, tri) # 5 ms
func = LinearTriInterpolator(tri, vals) # 9490 ms
func(*interp_points.T).data # 116 ms
# returns [0.54467719, 0.35885304, ...]
# total time 10.2 sec

### Scipy interpolate
tri = Delaunay(points) # 2720 ms
func = LinearNDInterpolator(tri, vals) # 1 ms
func(interp_points) # 925 ms

# returns [0.54467719, 0.35885304, ...]
# total time 3.65 sec

【讨论】：