numpy - Implementations and strategies for fast 2D interpolation from irregularly spaced points

Question

Given a large (~10 million) number of irregularly spaced points in two dimensions, where each point has some intensity ("weight") associated with it, what existing python implementations are there for interpolating the value at:

• a specific point at some random position (i.e. point = (0.5, 0.8))
• a large number of points at random positions (i.e. points = np.random.random((1_000_000, 2)))
• a regular grid at integer positions (i.e. np.indices((1000, 1000)).T)

I am aware that Delaunay triangulation is often used for this purpose. Are there alternatives to doing it this way?

Do any solutions take advantage of multiple CPU cores or GPUs?

As an example, here is an approach using scipy's LinearNDInterpolator. It does not appear to use more than one CPU core. There are also other options in scipy, but with this question I am especially interested in hearing about other solutions than the ones in scipy.

# The %time tags are IPython magic functions that time that specific line
dimension_shape = (1000, 1000) # we spread the random [0-1] over [0-1000] to avoid floating point errors
N_points = dimension_shape[0] * dimension_shape[1]

known_points = np.random.random((N_points, 2)) * dimension_shape 
known_weights = np.random.random((N_points,))

unknown_point = (0.5, 0.8)
unknown_points = np.random.random((N_points, 2)) * dimension_shape
unknown_grid = np.indices(dimension_shape, dtype=float).T.reshape((-1, 2)) # reshape to a list of 2D points

%time tesselation = Delaunay(known_points) # create grid to know neighbours # 6 sec
%time interp_func = LinearNDInterpolator(tesselation, known_weights) # 1 ms
%time interp_func(unknown_point) # 2 sec # run it once because the scipy function needs to compile 

%time interp_func(unknown_point) # ~ns
%time interp_func(unknown_grid) # 400 ms
%time interp_func(unknown_points) # 1 min 13 sec

# Below I sort the above `unknown_points` array, and try again
%time ind = np.lexsort(np.transpose(unknown_points)[::-1]) # 306 ms
unknown_points_sorted = unknown_points[ind].copy()

%time interp_func(unknown_points_sorted) # 19 sec <- much less than 1 min!

In the above code, things that take an appreciable amount of time are the construction of the Delaunay grid, and interpolation on a non-regular grid of points. Note that sorting the non-regular points first results in a significant speed improvement!

Do not feel the need to give a complete answer from the start. Tackling any aspect of the above is welcome.

question from:https://stackoverflow.com/questions/65847051/implementations-and-strategies-for-fast-2d-interpolation-from-irregularly-spaced

Answer 1

Scipy is pretty good and I don't think that there are better solutions in Python, but I can add a couple things that might be helpful to you. First off, your idea of sorting the points is a really good one. The so-called "incremental algorithms" build the Delaunay by inserting vertices one at a time. The first step in inserting a vertex in an existing mesh is to figure out which triangle in the mesh to insert it into. To speed things up, some algorithms start the search right at the point where the most recent insertion occurred. So if your points are ordered so that each point inserted is relatively close to the previous one, the search is much faster. If you want more details, you can look up the "Lawson's Walk" algorithm. In my own implementation of the Delaunay (which is in Java, so I'm afraid it won't help you), I have a sort based on the Hilbert space-filling curve. the Hilbert sort works great. But even just sorting by x/y coordinates is a help.

In terms of whether there are other ways to interpolate without using the Delaunay... You could try something using Inverse-Distance-Weighting (IDW). IDW techniques don't require the Delaunay, but they do require some way to figure out which vertices are close to the point for which you wish to interpolate. I've played with dividing my coordinate space into uniformly spaced bins, storing the vertices in the appropriate bins, and then just pulling up the points I need for an interpolation by looking at the neighboring bins. It may be a lot of coding, but it will be reasonably fast and use less memory than the Delaunay