6

I have my own triangulation algorithm that creates a triangulation based on both Delaunay's condition and the gradient such that the triangles align with the gradient.

This is an example output: enter image description here

The above description is not relevant to the question but is necessary for the context.

Now I want to use my triangulation with scipy.interpolate.LinearNDInterpolator to do an interpolation.

With scipy's Delaunay I would do the following

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)

This delaunay object has delaunay.points and delaunay.simplices that form the triangulation. I have the exact same information with my own triangulation, but scipy.interpolate.LinearNDInterpolator requires a scipy.spatial.Delaunay object.

I think I would need to subclass scipy.spatial.Delaunay and implement the relevant methods. However, I don't know which ones I need in order to get there.

4

1 回答 1

0

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++ 编写的，因此将代码转换为即 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
于 2021-01-24T16:44:16.383 回答