[SciPy-Dev] N-dimensional interpolation
Pauli Virtanen
pav@iki...
Sun Jul 25 10:35:11 CDT 2010
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
More information about the SciPy-Dev
mailing list