2

I'm trying to interpolate the value of a function at a single point from it's nearest neighbors.

I have f specified on a meshgrid, [x_grid,y_grid,z_grid] = np.meshgrid(x_range,y_range,z_range), for which I'd like to find an approximate value of a random point p_rand = (x_rand, y_rand, z_rand). What is an efficient way to find the indices of that nearest grid points and interpolate it's value? It's in 3D - nearest cube or tetrahedron of points would be fine.

anon01
  • 10,618
  • 8
  • 35
  • 58
  • 1
    `scipy` has nearest-neighbor interpolators. – hpaulj Sep 15 '15 at 19:54
  • 1
    Those interpolators create a search tree with your 'range' values, which can then can called repeatedly with the 'rand' values. Once the tree is created the interpolation (lookup) is supposed to be quite fast. – hpaulj Sep 15 '15 at 20:11
  • like @hpaulj said, `scipy`'s implementation of nearest neighbour interpolator uses a native C backend of spatial tree, with the option of parallelisation. The performance actually is pretty good. – Cong Ma Sep 15 '15 at 20:37
  • @CongMa So I just tried calling `interpolate.griddata(original_points,original_data, random_point)` twice in a row (as I thought hpaulj was suggesting) but it's not much faster. How do I reuse the interpolated grid so it's just evaluating a handful of times around the random point? – anon01 Sep 15 '15 at 20:41
  • Have you tried [`scipy.interpolate.NearestNDInterpolator`](http://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.NearestNDInterpolator.html#scipy.interpolate.NearestNDInterpolator)? – Cong Ma Sep 15 '15 at 20:43

1 Answers1

4

To expand on the comments of @hpaulj and mine... The class you're looking for is scipy.interpolate.NearestNDInterpolator

This class is based on scipy's own scipy.spatial.cKDTree. The cKDTree class implements the k-dimensional space-partition tree, or "k-d tree" data structure, which trades construction time and space for fast search.

To use scipy.interpolate.NearestNDInterpolator, you initialise an instance like

from scipy.interpolate import NearestNDInterpolator
interpolator = NearestNDInterpolator(your_points, your_values_on_the_points)

After creating interpolator, use it to evaluate the interpolant value at random_point by

interpolant = interpolator(random_point)

Once created, interpolator can be reused for different input points, which is a Good Thing (tm). You can also evaluate the interpolant value for multiple points by passing all of them into the call. [1]

If you look at the source, you'll discover how the interpolator is implemented using cKDTree.

[1] Actually there's a potential optimisation: if you need "vectorised" evaluation for many points, the underlying cKDTree's query() method supports parallelisation, done in native C code running in threads. Although scipy's own implementation of NearestNDInterpolator doesn't use this feature, probably catering to greatest common divisor, you can override this by making your own subclass that uses parallelisation with a suitable choice of n_jobs parameter.

Note: the really good thing about using a k-d tree based interpolator is that its application can be extended to a "grid" in arbitrary shape (not necessarily rectangular).

EDIT:

Oh, so you meant to use the nearest neighbors for linear interpolation?

Then very sorry, I misread your question!

But then, you have two choices.

  1. If your grid is sufficiently regular, and its construction (starting value / ending value / step known to you in all dimensions), it is not hard to write a function findneighbor() that locates the neighbors given the query point's coordinates. Then you do the vanilla linear interpolation.

  2. If your "grid" is not very regular, and you just have a lot of grind point coordinates (which may not lie on rectangular lattices), you can still use scipy.spatial.cKDTree to locate N nearest neighbors (possibly N = 1 + (dimension of the grid)). After that, you interpolate on that N points.

Cong Ma
  • 10,692
  • 3
  • 31
  • 47
  • well... I'm a bit resistant to doing nearest interpolation, but it is waaaay faster than linear interpolation. Thanks! – anon01 Sep 15 '15 at 23:39
  • @anon0909 I guess the doc assumes familiarity with k-nearest neighbor in general. If you feel this can be improved, you can file an issue to `scipy`. – Cong Ma Sep 21 '15 at 07:42