[SciPy-Dev] N-dimensional interpolation

[Pauli Virtanen]

Sun, 25 Jul 2010 21:52:24 +0000, Pauli Virtanen wrote:
[clip]
> As far as I understand, the search algorithm in scikits.delaunay walks
> the triangles by hopping to a neighbour chosen so that the target point
> is on positive side of the corresponding ridge. That generalizes to N-D
> easily.

This turned out to work quite well. A hybrid algorithm combining the two 
(first search a simplex with a positive plane distance, then continue 
with directed search) seems to work even better.

So now it should be quite robust, and I'm happy with the timings, too:

-- qhull triangulate					2.19968700409
-- qhull interpolate (meshgrid)				0.651250123978
-- qhull interpolate (random order)			22.201515913
-- scikits.delaunay triangulate				0.925380945206
-- scikits.delaunay linear_interpolate (meshgrid)	4.93817210197
-- scikits.delaunay nn_interpolate (meshgrid)		7.32182908058
-- scikits.delaunay nn_interpolate (random order)	45.4793360233
-- abs-diff-max 8.91013769433e-08
-- rel-diff-max 7.27796350818e-12

Qhull is roughly 2-3x slower in forming the triangulation in 2-D than 
scikits.delaunay, but that's acceptable. For some reason, the new 
interpolation code is faster than the linear interpolation in 
scikits.delaunay.

I'll probably try to split the code so that the Qhull geometry parts end 
up in scipy.spatial, and the interpolation routines in scipy.interpolate. 

I'll still need to test the Qhull triangulation against some pathological 
cases..

	Pauli



---------------------------------------------------------------------
import qhull
import numpy as np
import sys
import time

np.random.seed(1234)
x = np.random.rand(240*240, 2).astype(np.double)
y = np.arange(x.shape[0], dtype=np.double)

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()

# permuted order
xix = xi.reshape(nn*nn, 2)
p = np.random.permutation(nn*nn)
p2 = p.copy()
p2[p] = np.arange(nn*nn)
xix = xix[p]

# process!

print "-- qhull triangulate"
start = time.time()
ip = qhull.LinearNDInterpolator(x, y)
print time.time() - start

print "-- qhull interpolate (meshgrid)"
start = time.time()
zi = ip(xi)
print time.time() - start

print "-- qhull interpolate (random order)"
start = time.time()
zi = ip(xix)[p2].reshape(zi.shape)
print time.time() - start

from scikits.delaunay import Triangulation

print "-- scikits.delaunay triangulate"
start = time.time()
tri2 = Triangulation(x[:,0], x[:,1])
print time.time() - start

print "-- scikits.delaunay linear_interpolate (meshgrid)"
start = time.time()
ip2 = tri2.linear_interpolator(y)
zi2 = ip2[-0.1:1.1:(nn*1j),-0.1:1.1:(nn*1j)].T
print time.time() - start

print "-- scikits.delaunay nn_interpolate (meshgrid)"
start = time.time()
ip3 = tri2.nn_interpolator(y)
zi3 = ip3(xi[:,:,0].ravel(), xi[:,:,1].ravel()).reshape(zi.shape)
print time.time() - start

print "-- scikits.delaunay nn_interpolate (random order)"
start = time.time()
zi3 = ip3(xix[:,0], xix[:,1])[p2].reshape(zi.shape)
print time.time() - start

print "-- abs-diff-max", np.nanmax(abs(zi-zi2))
print "-- rel-diff-max", np.nanmax(abs(zi-zi2)/abs(zi))