2

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).

Community
  • 1
  • 1
mib0163
  • 393
  • 3
  • 12
  • Wow - that's a blast from the past. My final project for my degree (25 years ago) was a program to do 2D Delaunay Triangulation in 'C'. Thanks for the nostalgia trip... – SiHa May 22 '15 at 08:47

3 Answers3

7

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).

Community
  • 1
  • 1
mib0163
  • 393
  • 3
  • 12
3

You don't need to implement this from scratch, there is already built-in support in scipy for this feature:

scipy.interpolate.LinearNDInterpolator

cfh
  • 4,576
  • 1
  • 24
  • 34
  • I know, but that only works if there is a float value assigned to each point of the grid. What I interpolate is more complex. I have edited the question to make this clear. – mib0163 May 22 '15 at 14:58
  • 1
    @mib0163: I'm pretty sure LinearNDInterpolator can deal with vector-valued values, not only floats. See for instance here: http://stackoverflow.com/questions/12391747/interpolation-in-vector-valued-multi-variate-function – cfh May 22 '15 at 15:22
  • Cool, I was not aware of that. So this would be a straightforward solution. – mib0163 May 22 '15 at 16:39
0

You need an interval and a linear interpolation, i.e. the lenght of the edge and the distance of the interpolated points from the start vertex.

Micromega
  • 12,486
  • 7
  • 35
  • 72