Interpolation with Delaunay Triangulation

Question

Having a cloud point shaped like some sort of distorted paraboloid, I would like to use Delaunay Triangulation to interpolate the points. I have tried other techniques (f.ex. splines) but did not manage to enforce the desired behavior.

I was wondering if there's a quick way to use the results of scipy.spatial.Delaunay in a way where I can give the (x,y) coords and get the z-coord of the point on the simplex (triangle). From the documentation looks like I can pull out the index of the simplex but I am not sure how to take it from there.

Answer 1

You can give the Delaunay triangulation to scipy.interpolate.LinearNDInterpolator together with the set of Z-values, and it should do the job for you.

If you really want to do the interpolation yourself, you can build it up from find_simplex and transform.

Answer 2

We can compute the barycentric coordinates using the object returned by scipy.special.Delaunay() from the given sample points and interpolate the function as follows (can be done without using scipy.interpolate.LinearNDInterpolator):

from scipy.spatial import Delaunay
import matplotlib.pylab as plt

nx, ny = 160, 300

xv, yv = np.meshgrid(np.arange(nx), np.arange(ny))
xv = xv.flatten()
yv = yv.flatten()
pp = np.vstack((xv, yv)).T
    
n = 100

ix = np.random.choice(nx, n).tolist() + [0, 0, nx-1, nx-1]
iy = np.random.choice(ny, n).tolist() + [0, ny-1, 0, ny-1]
f = np.zeros((nx, ny))
for x, y in zip(ix, iy):
    f[x, y] = (x/160)**2 + (y/300)**2 + np.random.rand()
points = np.array(list(zip(ix, iy)))
tri = Delaunay(points)
ss = tri.find_simplex(pp)
ndim = tri.transform.shape[2]
print(n,len(np.unique(ss)))

out = np.zeros((nx, ny))
for i in np.unique(ss): # for all simplices (triangles)
    p = pp[ss == i] # all points in the simplex
    # compute the barycentric coordinates of the points
    b = tri.transform[i, :ndim].dot(np.transpose(p) - tri.transform[i, ndim].reshape(-1,1))
    αβγ = np.c_[np.transpose(b), 1 - b.sum(axis=0)] 
    sp = points[tri.simplices[i]]
    if len(αβγ) > 0:
        out[p[:,0], p[:,1]] = αβγ@f[sp[:,0], sp[:,1]]     
out = out / out.max()

plt.figure(figsize=(10,12))
plt.imshow(out, aspect='auto', cmap='Oranges'), plt.axis(), plt.title('interpolated output with\n barycentric coordinates', size=20)
plt.triplot(points[:,1], points[:,0], tri.simplices)
for p in points:
    plt.scatter(p[1], p[0], c=f[p[0],p[1]], s=50)
plt.title(f'Interpolation with Delaunay Triangulation', size=20)
plt.show()