How to use scipy.interpolate.LinearNDInterpolator with own triangulation

Question

I have my own triangulation algorithm that creates a triangulation based on both Delaunay's condition and the gradient such that the triangles align with the gradient.

This is an example output:

The above description is not relevant to the question but is necessary for the context.

Now I want to use my triangulation with scipy.interpolate.LinearNDInterpolator to do an interpolation.

With scipy's Delaunay I would do the following

import numpy as np
import scipy.interpolate
import scipy.spatial
points = np.random.rand(100, 2)
values = np.random.rand(100)
delaunay = scipy.spatial.Delaunay(points)
ip = scipy.interpolate.LinearNDInterpolator(delaunay, values)

This delaunay object has delaunay.points and delaunay.simplices that form the triangulation. I have the exact same information with my own triangulation, but scipy.interpolate.LinearNDInterpolator requires a scipy.spatial.Delaunay object.

I think I would need to subclass scipy.spatial.Delaunay and implement the relevant methods. However, I don't know which ones I need in order to get there.

Answer 1

I wanted to do this very same thing with the Delaunay triangulation offered by the triangle package. The triangle Delaunay code is about eight times faster than the SciPy one on large (~100_000) points. (I encourage other developers to try to beat that :) )

Unfortunately, the Scipy LinearNDInterpolator function relies heavily on specific attributes present in the SciPy Delaunay triangulation object. These are created by the _get_delaunay_info() CPython code, which is difficult to disassemble. Even knowing which attributes are needed (there seem to be many, including things like paraboloid_scale and paraboloid_shift), I'm not sure how I would extract this from a different triangulation library.

Instead, I tried @Patol75's approach (comment to the Question above), but using LinearTriInterpolator instead of the Cubic one. The code runs correctly, but is slower than doing the entire thing in SciPy. Interpolating 400_000 points from a cloud of 400_000 points takes about 3 times longer using the matplotlib code than scipy. The Matplotlib tri code is written in C++, so converting the code to i.e. CuPy is not straightforward. If we could mix the two approaches, we could reduce the total time from 3.65 sec / 10.2 sec to 1.1 seconds!

import numpy as np
np.random.seed(1)
N = 400_000
shape = (100, 100)
points = np.random.random((N, 2)) * shape # spread over 100, 100 to avoid float point errors
vals = np.random.random((N,))
interp_points1 = np.random.random((N,2)) * shape
interp_points2 = np.random.random((N,2)) * shape

triangle_input = dict(vertices=points)

### Matplotlib Tri
import triangle as tr
from matplotlib.tri import Triangulation, LinearTriInterpolator

triangle_output = tr.triangulate(triangle_input) # 280 ms
tri = tr.triangulate(triangle_input)['triangles'] # 280 ms
tri = Triangulation(*points.T, tri) # 5 ms
func = LinearTriInterpolator(tri, vals) # 9490 ms
func(*interp_points.T).data # 116 ms
# returns [0.54467719, 0.35885304, ...]
# total time 10.2 sec

### Scipy interpolate
tri = Delaunay(points) # 2720 ms
func = LinearNDInterpolator(tri, vals) # 1 ms
func(interp_points) # 925 ms

# returns [0.54467719, 0.35885304, ...]
# total time 3.65 sec