[Yt-dev] Fwd: [SciPy-Dev] RFR: N-dimensional interpolation
Matthew Turk
matthewturk at gmail.com
Tue Aug 31 15:04:57 PDT 2010
Very cool: includes natural neighbor interpolation in N dimensions!
---------- Forwarded message ----------
From: Pauli Virtanen <pav at iki.fi>
Date: Tue, Aug 31, 2010 at 3:01 PM
Subject: [SciPy-Dev] RFR: N-dimensional interpolation
To: scipy-dev at scipy.org
Sun, 25 Jul 2010 15:35:11 +0000, Pauli Virtanen wrote:
> I took the Qhull by the horns, and wrote a straightforward `griddata`
> implementation for working in N-D:
[clip]
It's now committed to SVN (as it works-well-for-me).
Post-mortem review and testing is welcome.
http://projects.scipy.org/scipy/changeset/6653
http://projects.scipy.org/scipy/changeset/6655
http://projects.scipy.org/scipy/changeset/6657
What's in there is:
1) scipy.spatial.qhull
Delaunay decomposition and some associated low-level N-d geometry
routines.
2) scipy.interpolate.interpnd
N-dimensional interpolation:
1) Linear barycentric interpolation
2) Cubic spline interpolation (2D-only, C1 continuous,
approximately minimum-curvature).
3) scipy.interpolate.griddatand
Convenience interface to the N-d interpolation classes.
What could be added:
- More comprehensive interface to other features of Qhull
- Using qhull_restore, qhull_save to store Qhull contexts
instead of copying the relevant data?
- Optimizing the cubic interpolant
- Monotonic cubic interpolation
- Cubic interpolation in 3-d
- Natural neighbour interpolation
- etc.
***
Example:
import numpy as np
def func(x, y):
return x*(1-x)*np.cos(4*np.pi*x) * np.sin(4*np.pi*y**2)**2
grid_x, grid_y = np.mgrid[0:1:100j, 0:1:200j]
points = np.random.rand(1000, 2)
values = func(points[:,0], points[:,1])
from scipy.interpolate import griddata
grid_z0 = griddata(points, values, (grid_x, grid_y), method='nearest')
grid_z1 = griddata(points, values, (grid_x, grid_y), method='linear')
grid_z2 = griddata(points, values, (grid_x, grid_y), method='cubic')
import matplotlib.pyplot as plt
plt.subplot(221)
plt.imshow(func(grid_x, grid_y).T, extent=(0,1,0,1), origin='lower')
plt.plot(points[:,0], points[:,1], 'k.', ms=1)
plt.title('Original')
plt.subplot(222)
plt.imshow(grid_z0.T, extent=(0,1,0,1), origin='lower')
plt.title('Nearest')
plt.subplot(223)
plt.imshow(grid_z1.T, extent=(0,1,0,1), origin='lower')
plt.title('Linear')
plt.subplot(224)
plt.imshow(grid_z2.T, extent=(0,1,0,1), origin='lower')
plt.title('Cubic')
plt.gcf().set_size_inches(6, 6)
plt.show()
--
Pauli Virtanen
_______________________________________________
SciPy-Dev mailing list
SciPy-Dev at scipy.org
http://mail.scipy.org/mailman/listinfo/scipy-dev
More information about the yt-dev
mailing list