使用 Delaunay 三角剖分 (n-dim) 进行插值

Question

我想在 Python 中使用 Delaunay Triangulation 来插入 3D 中的点。

我拥有的是

# my array of points
points = [[1,2,3], [2,3,4], ...]
# my array of values
values = [7, 8, ...]
# an object with triangulation
tri = Delaunay(points)        
# a set of points at which I want to interpolate
p = [[1.5, 2.5, 3.5], ...]
# this gets simplexes that contain given points
s = tri.find_simplex(p)
# this gets vertices for the simplexes
v = tri.vertices[s]

我只能在这里找到一个建议使用方法进行插值的答案transform，但没有更具体。

我需要知道的是如何使用包含单纯形的顶点来获得线性插值的权重。让我们假设一个一般的 n-dim 情况，以便答案不依赖于维度。

编辑：我不想使用 LinearNDInterpolator 或类似的方法，因为我在每个点都没有一个数字作为值，而是更复杂的（数组/函数）。

Answer 1

经过一些试验，解决方案看起来很简单（这篇文章很有帮助）：

# dimension of the problem (in this example I use 3D grid,
# but the method works for any dimension n>=2)
n = 3
# my array of grid points (array of n-dimensional coordinates)
points = [[1,2,3], [2,3,4], ...]
# each point has some assigned value that will be interpolated
# (e.g. a float, but it can be a function or anything else)
values = [7, 8, ...]
# a set of points at which I want to interpolate (it must be a NumPy array)
p = np.array([[1.5, 2.5, 3.5], [1.1, 2.2, 3.3], ...])

# create an object with triangulation
tri = Delaunay(points)        
# find simplexes that contain interpolated points
s = tri.find_simplex(p)
# get the vertices for each simplex
v = tri.vertices[s]
# get transform matrices for each simplex (see explanation bellow)
m = tri.transform[s]

# for each interpolated point p, mutliply the transform matrix by 
# vector p-r, where r=m[:,n,:] is one of the simplex vertices to which 
# the matrix m is related to (again, see bellow)
b = np.einsum('ijk,ik->ij', m[:,:n,:n], p-m[:,n,:])

# get the weights for the vertices; `b` contains an n-dimensional vector
# with weights for all but the last vertices of the simplex
# (note that for n-D grid, each simplex consists of n+1 vertices);
# the remaining weight for the last vertex can be copmuted from
# the condition that sum of weights must be equal to 1
w = np.c_[b, 1-b.sum(axis=1)]

理解的关键方法是transform，它被简要记录，但是文档说明了所有需要说明的内容。对于每个单纯形，transform[:,:n,:n]包含变换矩阵，并transform[:,n,:]包含与矩阵相关的向量 r。似乎选择 r 向量作为单纯形的最后一个顶点。

另一个棘手的问题是如何获得b，因为我想做的是

for i in range(len(p)): b[i] = m[i,:n,:n].dot(p[i]-m[i,n,:])

本质上，我需要一个点积数组，同时dot给出两个数组的乘积。像上面这样的单个单纯形循环可以工作，但是可以一步完成更快，为此有numpy.einsum：

b = np.einsum('ijk,ik->ij', m[:,:n,:n], p-m[:,n,:])

现在，v包含每个单纯形的顶点索引并w保存相应的权重。为了p_values在一组点处获得插值p，我们这样做（注意：values必须是 NumPy 数组）：

values = np.array(values)
for i in range(len(p)): p_values[i] = np.inner(values[v[i]], w[i])

或者我们可以再次使用 `np.einsum' 一步完成：

p_values = np.einsum('ij,ij->i', values[v], w)

在某些内插点位于网格之外的情况下，必须小心谨慎。在这种情况下，find_simplex(p)返回-1这些点，然后您将不得不屏蔽它们（可能使用屏蔽数组）。