[SciPy-Dev] N-dimensional interpolation

[Travis Oliphant]

Fantastic.   Excellent work.   This is something we have been missing for a long time.

Travis 

--
(mobile phone of)
Travis Oliphant
Enthought, Inc.
1-512-536-1057
http://www.enthought.com

On Jul 25, 2010, at 10:35 AM, Pauli Virtanen <pav at iki.fi> wrote:

> Hi all,
> 
> I took the Qhull by the horns, and wrote a straightforward `griddata` 
> implementation for working in N-D:
> 
>    http://github.com/pv/scipy-work/tree/qhull
>    http://github.com/pv/scipy-work/blob/qhull/scipy/interpolate/qhull.pyx
> 
> Sample (4-D):
> -------------------------------------------------------------
> import numpy as np
> import matplotlib.pyplot as plt
> from scipy.interpolate import griddata
> 
> def f(x, y, z, u):
>    return np.cos(x + np.cos(y + np.cos(z + np.cos(u))))
> 
> points = np.random.rand(2000, 4)
> values = f(*points.T)
> 
> p = np.linspace(0, 1, 1500)
> z = griddata(points, values, np.c_[p, p, p, p], method='linear')
> 
> plt.plot(p, z, '-', p, f(p, p, p, p), '-')
> plt.show()
> -------------------------------------------------------------
> 
> It performs the N-D Delaunay tesselation with Qhull, and after that does 
> linear barycentric interpolation on the simplex containing each point. 
> (The nearest-neighbor "interpolation" is, on the other hand, implemented 
> using scipy.spatial.KDTree.)
> 
> I ended up writing some custom simplex-walking code for locating the 
> simplex containing the points -- this is much faster than a brute-force 
> search. However, it is slightly different from what `qh_findbestfacet` 
> does. [If you're a computational geometry expert, see below...]
> 
> Speed-wise, Qhull appears to lose in 2-D to `scikits.delaunay` by a 
> factor of about 10 in triangulation in my tests; this, however, includes 
> the time taken by computing the barycentric coordinate transforms. 
> Interpolation, however, seems to be faster.
> 
>    ***
> 
> I think this would be a nice addition to `scipy.optimize`. I'd like to 
> get it in for 0.9.
> 
> Later on, we can specialize the `griddata` function to perform better in 
> low dimensions (ndim=2, 3) and add more interpolation methods there; for 
> example natural neighbors, interpolating splines, etc. In ndim > 3, the 
> `griddata` is however already feature-equivalent to Octave and MATLAB.
> 
> Comments?
> 
>    Pauli
> 
> 
> PS. A question for computational geometry experts:
> 
> `qh_findbestfacet` finds the facet on the lower convex hull whose 
> oriented hyperplane distance to the point lifted on the paraboloid is 
> maximal --- however, this does not always give the simplex containing the 
> point in the projected Delaunay tesselation. There's a counterexample in 
> qhull.pyx:_find_simplex docstring.
> 
> On the other hand, Qhull documentation claims that this is the way to 
> locate the Delaunay simplex containing a point. What gives?
> 
> It's clear, however, that if a simplex contains the point, then the 
> hyperplane distance is positive. So currently, I just walk around 
> positive-distance simplices checking the inside condition for each and 
> hoping for the best. If nothing is found, I just fall back to brute force 
> search.
> 
> This seems to work well in practice. However, if a point is outside the 
> triangulation (but inside the bounding box), I have to fall back to a 
> brute-force search. Is there a better option here?
> 
>    ***
> 
> Another sample (2-D, object-oriented interface):
> -------------------------------------------------------------
> import time
> import numpy as np
> import matplotlib.pyplot as plt
> from scipy.interpolate import LinearNDInterpolator
> 
> # some random data
> x = np.array([(0,0), (0,1), (1, 1), (1, 0)], dtype=np.double)
> np.random.seed(0)
> xr = np.random.rand(200*200, 2).astype(np.double)
> x = np.r_[xr, x]
> y = np.arange(x.shape[0], dtype=np.double)
> 
> # evaluate on a grid
> nn = 500
> xx = np.linspace(-0.1, 1.1, nn)
> yy = np.linspace(-0.1, 1.1, nn)
> xx, yy = np.broadcast_arrays(xx[:,None], yy[None,:])
> 
> 
> xi = np.array([xx, yy]).transpose(1,2,0).copy()
> start = time.time()
> ip = LinearNDInterpolator(x, y)
> print "Triangulation", time.time() - start
> start = time.time()
> zi = ip(xi)
> print "Interpolation", time.time() - start
> 
> from scikits.delaunay import Triangulation
> 
> start = time.time()
> tri2 = Triangulation(x[:,0], x[:,1])
> print "scikits.delaunay triangulation", time.time() - start
> 
> start = time.time()
> ip2 = tri2.linear_interpolator(y)
> zi2 = ip2[-0.1:1.1:500j,-0.1:1.1:500j].T
> print "scikits.delaunay interpolation", time.time() - start
> 
> print "rel-difference", np.nanmax(abs(zi - zi2))/np.nanmax(np.abs(zi))
> 
> plt.figure()
> plt.imshow(zi)
> plt.clim(y.min(), y.max())
> plt.colorbar()
> plt.figure()
> plt.imshow(zi2)
> plt.clim(y.min(), y.max())
> plt.colorbar()
> plt.show()
> -------------------------------------------------------------
> 
> -- 
> Pauli Virtanen
> 
> _______________________________________________
> SciPy-Dev mailing list
> SciPy-Dev at scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-dev