[Scipy-svn] r4620 - in branches/Interpolate1D: . tests

scipy-svn@scip... scipy-svn@scip...
Fri Aug 8 11:02:25 CDT 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