[Scipy-svn] r4620 - in branches/Interpolate1D: . tests
scipy-svn at scipy.org
scipy-svn at scipy.org
Fri Aug 8 12:02:25 EDT 2008
Author: fcady
Date: 2008-08-08 11:02:21 -0500 (Fri, 08 Aug 2008)
New Revision: 4620
Modified:
branches/Interpolate1D/TODO.txt
branches/Interpolate1D/interpolateNd.py
branches/Interpolate1D/tests/test_interpolate2d.py
branches/Interpolate1D/tests/test_ndimage.py
Log:
ndimage code more seemlessly incorporated into interpolateNd.py
Modified: branches/Interpolate1D/TODO.txt
===================================================================
--- branches/Interpolate1D/TODO.txt 2008-08-08 07:11:28 UTC (rev 4619)
+++ branches/Interpolate1D/TODO.txt 2008-08-08 16:02:21 UTC (rev 4620)
@@ -11,6 +11,11 @@
************ API and/or MAJOR ISSUES ***********
+**better ND interpolation
+ InterpolateNd works, but the quadratic and higher order splines are
+ sometimes really ugly. It seems they have a boundary condition that the
+ derivative must be zero.
+
**possibly allow interp2d to return a 2d array? Like in the case that
x and y are 2d arrays in meshgrid format.
Modified: branches/Interpolate1D/interpolateNd.py
===================================================================
--- branches/Interpolate1D/interpolateNd.py 2008-08-08 07:11:28 UTC (rev 4619)
+++ branches/Interpolate1D/interpolateNd.py 2008-08-08 16:02:21 UTC (rev 4620)
@@ -2,7 +2,7 @@
from numpy import array, arange, NaN
import numpy as np
-import _nd_image
+import _nd_image # extension module with interpolation code
def interpNd(data, coordinates, starting_coords=None, spacings=None, kind='linear',out=NaN):
""" A function for interpolation of 1D, real-valued data.
@@ -70,6 +70,10 @@
class InterpolateNd:
""" A callable class for interpolation of 1D, real-valued data.
+
+ Under the hood, interpolation is always done with splines. Those
+ of order 2 or higher have the boundary condition of zero derivative.
+ This is a potential shortcoming of the package; natural splines would be better.
Parameters
-----------
@@ -120,6 +124,11 @@
>>> nd.InterpolateNd(boring_data)( np.array([[2.3], [1.0], [3.9]]) )
1.0
"""
+ # Under the hood this works by interpolating entries in an array, rather than
+ # values at points in space. Known data are stored as an array; info on spacings
+ # and starting coords describes the mapping from the array to ND space. When the
+ # function is called, the coordinates are first converted into array values.
+ # The array interpolation is handled by the C extension _nd_image.
def __init__(self, data, starting_coords =None, spacings = None,
kind='linear', out=NaN):
""" data = array or list of lists
@@ -131,9 +140,6 @@
or just NaN
"""
- # FIXME : include spline filtering
- # the ndimage module says that it requires pre-filtering for
-
# checking format of input
data = array(data)
@@ -192,21 +198,21 @@
else:
raise ValueError, "argument kind = %s not recognized" % str(kind)
-
+ # Spline pre-filtering
# This step is done because it is required by the ndimage code that I'm scavenging.
# I don't fully understand why it must do this, and that's a problem. But empirically
# this step is needed in order to get good-looking data.
if self.order >1:
- self._data_array = spline_filter(data, self.order)
+ self._data_array = self._spline_filter(data, self.order)
else:
self._data_array = data
# storing relevant data
- self.ndim = data.ndim
- self._shape = np.shape(data)
+ self.ndim = self._data_array.ndim
+ self.shape_as_column_vec = np.reshape(np.shape(self._data_array),
+ (self.ndim, 1))
self._spacings = spacings
self._min_coords = starting_coords
- self._max_coords = self._min_coords + self._shape*self._spacings
self.out = out
def __call__(self, coordinates):
@@ -271,73 +277,36 @@
""" return an array of bools saying which
points are in interpolation bounds
"""
- shape_as_column_vec = np.reshape(self._shape, (self.ndim, 1))
# entry is 1 if that coordinate of a point is in its bounds
index_in_bounds = (0 <= indices) & \
- (indices <= shape_as_column_vec)
+ (indices <= self.shape_as_column_vec)
# for each point, number of coordinates that are in bounds
num_indices_in_bounds = np.sum(index_in_bounds, axis=0)
# True if each coordinate for the point is in bounds
return num_indices_in_bounds == self.ndim
-
- def _coord_in_bounds(self, coordinates):
- """ return an array of bools saying which
- points are in interpolation bounds
- """
- # entry is 1 if that coordinate of a point is in its bounds
- coord_in_bounds = (self._min_coords <= coordinates) & \
- (coordinates <= self._max_coords)
-
- # for each point, number of coordinates that are in bounds
- num_coords_in_bounds = np.sum(coord_in_bounds, axis=0)
-
- # True if each coordinate for the point is in bounds
- return num_coords_in_bounds == self.ndim
-
-
-
-
-
-
-
-
-
-import _ni_support
-def spline_filter1d(input, order = 3, axis = -1, output = np.float64,
- output_type = None):
- # takes array and everything; we can make input safe if user never touches it
- """ Calculates a one-dimensional spline filter along the given axis.
+ def _spline_filter(self, data, order = 3):
+ """ Multi-dimensional spline filter.
- The lines of the array along the given axis are filtered by a
- spline filter. The order of the spline must be >= 2 and <= 5.
- """
- if order in [0, 1]:
- output[...] = np.array(input)
- else:
- _nd_image.spline_filter1d(input, order, axis, output)
- return output
-
-def spline_filter(input, order = 3, output = np.float64,
- output_type = None):
- """ Multi-dimensional spline filter.
-
- Note: The multi-dimensional filter is implemented as a sequence of
- one-dimensional spline filters. The intermediate arrays are stored
- in the same data type as the output. Therefore, for output types
- with a limited precision, the results may be imprecise because
- intermediate results may be stored with insufficient precision.
- """
-
- output = np.zeros( np.shape(input) , dtype=np.float64 ) # place to store the data
-
- if order not in [0, 1] and input.ndim > 0:
- for axis in range(input.ndim):
- spline_filter1d(input, order, axis, output = output)
- input = output
- else:
- output[...] = input[...]
- return output
\ No newline at end of file
+ Note: The multi-dimensional filter is implemented as a sequence of
+ one-dimensional spline filters. The intermediate arrays are stored
+ in the same data type as the output. Therefore, for output types
+ with a limited precision, the results may be imprecise because
+ intermediate results may be stored with insufficient precision.
+
+ order must be from 0 to 5
+ """
+
+ # allocate memory in which to put the output
+ output = np.zeros( np.shape(data) , dtype=np.float64 )
+
+ if order in [0, 1]:
+ output = data
+ else:
+ for axis in range(data.ndim):
+ _nd_image.spline_filter1d(data, order, axis, output)
+ data = output
+ return output
\ No newline at end of file
Modified: branches/Interpolate1D/tests/test_interpolate2d.py
===================================================================
--- branches/Interpolate1D/tests/test_interpolate2d.py 2008-08-08 07:11:28 UTC (rev 4619)
+++ branches/Interpolate1D/tests/test_interpolate2d.py 2008-08-08 16:02:21 UTC (rev 4620)
@@ -9,7 +9,7 @@
from numpy import arange, meshgrid, ravel
from interpolate2d import interp2d, Interpolate2d, atleast_1d_and_contiguous
-from fitpack_wrapper2d import Spline2d
+from fitpack_wrapper import Spline2d
# unit testing
import unittest, time
Modified: branches/Interpolate1D/tests/test_ndimage.py
===================================================================
--- branches/Interpolate1D/tests/test_ndimage.py 2008-08-08 07:11:28 UTC (rev 4619)
+++ branches/Interpolate1D/tests/test_ndimage.py 2008-08-08 16:02:21 UTC (rev 4620)
@@ -13,8 +13,8 @@
class Test (unittest.TestCase):
- def assertAllclose(self, x, y):
- self.assert_(np.allclose(x, y))
+ def assertAllclose(self, x, y, err=1.0e-8):
+ self.assert_(np.allclose(x, y, atol=err))
def test_interpNd(self):
""" Make sure : the function interpNd works
@@ -67,7 +67,8 @@
X, Y = np.meshgrid(arange(10.), arange(10.))
interesting_data = X+Y
interp = nd.InterpolateNd(interesting_data, kind = 2)
- self.assertAllclose( interp(np.array([[2.3], [1.0]])) , 3.3 )
+ print "quad answer: ", interp(np.array([[2.3], [1.0]]))
+ self.assertAllclose( interp(np.array([[2.3], [1.0]])) , 3.3 , err=.1 )
def test_order0(self):
""" Make sure : block interpolation works
@@ -83,7 +84,8 @@
X, Y = np.meshgrid(arange(10.), arange(10.))
interesting_data = X+Y
interp = nd.InterpolateNd(interesting_data, kind = 3)
- self.assertAllclose( interp(np.array([[4.3], [4.1]])) , 8.4 )
+ print "cubi answer: ", interp(np.array([[4.3], [4.1]]))
+ self.assertAllclose( interp(np.array([[4.3], [4.1]])) , 8.4 , err=.1)
def test_out(self):
""" Make sure : out-of-bounds returns NaN
More information about the Scipy-svn
mailing list