I would like to use Delaunay Triangulation in Python to interpolate the points in 3D.
What I have is
# 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]
I was only able to find one answer here that suggest to use transform method for the interpolation, but without being any more specific.
What I need to know is how to use the vertices of the containing simplex to get the weights for the linear interpolation. Let's assume a general n-dim case so that the answer does not depend on the dimension.
EDIT: I do not want to use LinearNDInterpolator or similar approach because I do not have a number at each point as a value but something more complex (array/function).
After some experimenting, the solution looks simple (this post was quite helpful):
# 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)]
The key method to understand is transform, which is briefly documented, however the documentation says all it needs to be said. For each simplex, transform[:,:n,:n] contains the transformation matrix, and transform[:,n,:] contains the vector r to which the matrix is related to. It seems that r vector is chosen as the last vertex of the simplex.
Another tricky point is how to get b, because what I want to do is something like
for i in range(len(p)): b[i] = m[i,:n,:n].dot(p[i]-m[i,n,:])
Essentially, I need an array of dot products, while dot gives a product of two arrays. The loop over the individual simplexes like above would work, but a it can be done faster in one step, for which there is numpy.einsum:
b = np.einsum('ijk,ik->ij', m[:,:n,:n], p-m[:,n,:])
Now, v contains indices of vertex points for each simplex and w holds corresponding weights. To get the interpolated values p_values at set of points p, we do (note: values must be NumPy array for this):
values = np.array(values)
for i in range(len(p)): p_values[i] = np.inner(values[v[i]], w[i])
Or we may do this in a single step using `np.einsum' again:
p_values = np.einsum('ij,ij->i', values[v], w)
Some care must be taken in situations, when some of the interpolated points lie outside the grid. In such case, find_simplex(p) returns -1 for those points and then you will have to mask out them (using masked arrays perhaps).