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

Question

我想在 Python 中使用 Delaunay 三角剖分来对 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]

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

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

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

Answer 1

经过一些试验，解决方案看起来很简单(这个 post 很有帮助):

# 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 处获取内插值 p_values，我们这样做(注意:为此，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，然后您将不得不屏蔽它们(可能使用 masked arrays)。