如何将scipy.interpolate.LinearNDInterpolator与自己的三角剖分一起使用

Question

我有自己的三角剖分算法，该算法基于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并实现相关方法。但是，我不知道要到达那里需要哪些。

Answer 1

triangle我想用软件包提供的 Delaunay 三角测量来做同样的事情。在大点 (~100_000) 上，三角形 Delaunay 代码比 SciPy 代码快大约八倍。（我鼓励其他开发人员尝试克服这一点:)）

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

相反，我尝试了 @Patol75 的方法（对上面的问题进行评论），但使用的LinearTriInterpolator是 Cubic 方法。代码运行正确，但比在 SciPy 中执行整个操作要慢。使用 matplotlib 代码从 400_000 个点的云中插值 400_000 个点所需的时间大约是 scipy 的 3 倍。Matplotlib tri 代码是用 C++ 编写的，因此将代码转换为 ie 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