如何使用scipy.interpolate.LinearNDInterpolator自定义三角剖分？

Question

我有自己的三角剖分算法（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并实现相关方法。然而，我不知道我需要哪些方法才能实现这一点。

- johnbaltis

---

LinearNDInterpolator 还可以接受一个点数组作为其第一个参数。 - Warren Weckesser

我知道，但这将导致使用scipy.spatial.Delaunay进行新的三角剖分，而这不是我想要的。 - johnbaltis

啊，好的。https://www.youtube.com/watch?v=VrbybKWwb7c - Warren Weckesser

我有同样的问题，我认为Scipy方面没有任何开发来提供这个功能。我尝试创建一个虚拟的Delaunay对象并修改其属性内容，但是我得到了“AttributeError：无法设置属性”。好消息是，如果您正在使用2D，则Matplotlib提供此功能。使用Triangulation创建三角剖分对象，并使用CubicTriInterpolator进行插值。 - Patol75

Answer 1

我想使用triangle包提供的Delaunay三角剖分来完成同样的事情。在大约100,000个点的情况下，triangle Delaunay代码比SciPy快8倍。（我鼓励其他开发人员尝试超越它 :)）

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

相反，我尝试了@Patol75的方法（上面问题的评论），但是使用{{link1：LinearTriInterpolator}}而不是Cubic。代码可以正确运行，但比在SciPy中完成整个操作要慢。从400,000个点的云中插值400,000个点使用matplotlib代码需要大约3倍的时间，而使用scipy则更快。{{link2：Matplotlib tri code是用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

- TomNorway