[SciPy-Dev] N-dimensional interpolation

[Charles R Harris]

On Sun, Jul 25, 2010 at 9: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.
>
>
Do you mean scipy.interpolate?


> 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?
>
>
I'm not a computational geometry expert, but I expect the 2-D performance or
scikits.delaunay takes advantage of a special algorithm for that case (I
recall googling for that back when). Perhaps it can be adapted to the higher
dimension case or adapted for searching. Just random thoughts here, I have
no actual experience with these sorts of problems.

I certainly support getting added interpolation functionality into scipy and
barycentric interpolation has the advantage of not introducing weird
artifacts. I like it for drawing contours also.

Chuck
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/scipy-dev/attachments/20100725/2d9294e7/attachment.html>