[Scipy-svn] r4611 - branches/Interpolate1D

scipy-svn at scipy.org scipy-svn at scipy.org
Thu Aug 7 19:42:51 EDT 2008


Author: fcady
Date: 2008-08-07 18:42:44 -0500 (Thu, 07 Aug 2008)
New Revision: 4611

Modified:
   branches/Interpolate1D/TODO.txt
   branches/Interpolate1D/info.py
   branches/Interpolate1D/interpolateNd.py
Log:
added prefiltering from ndimage, which it turns out is necessary to get good data with spline order >1.  It still needs work to make it pretty and integrated, but it is functional

Modified: branches/Interpolate1D/TODO.txt
===================================================================
--- branches/Interpolate1D/TODO.txt	2008-08-07 21:25:16 UTC (rev 4610)
+++ branches/Interpolate1D/TODO.txt	2008-08-07 23:42:44 UTC (rev 4611)
@@ -36,6 +36,15 @@
     
 **include spline pre-filtering in ND interpolation,
     or understand why we don't need to do that.
+    
+    It appears that higher-order splines also smooth the data
+    increasingly.  Yes, I checked and they definitely don't
+    reproduce perfectly the input; they flatten stuff.
+    
+    The pre-filtering definitely improves the output by a whole lot;
+    I suspect it returns some sort of transform of the data.  But I
+    really wish I knew what it did; why didn't this guy document his
+    code better than he did??
 
 **README file that describes the architecture of the
     package and also includes license information.
@@ -48,9 +57,13 @@
 
 **allow newx to be in non-sorted order for 1D
     This requires rethinking the partition of newx into
-    low, high and mid
+    low, high and kind
     
-**put all extension files into their own directory
+    The main hurdle is that I can only really fill and array
+    using result[mask] = partial_result if I can guarantee
+    that everything will have the same type.  Maybe it is
+    worth compromising the generality of the function to
+    allow that.
     
 ********* LONGER TERM ************
 

Modified: branches/Interpolate1D/info.py
===================================================================
--- branches/Interpolate1D/info.py	2008-08-07 21:25:16 UTC (rev 4610)
+++ branches/Interpolate1D/info.py	2008-08-07 23:42:44 UTC (rev 4611)
@@ -27,7 +27,7 @@
     
     The following callable classes are also provided:
 
-        Interpolate1d  :   an object for interpolation of
+        Interpolate1d  :   an object for 1D interpolation of
                                 various kinds.  interp1d is a wrapper
                                 around this class.
                                 
@@ -36,15 +36,21 @@
                                 is used.  However, not all functionality
                                 of Spline is available through Interpolate1d.
                                 
-        Interpolate2d  :
+        Interpolate2d  :  an object for 2D interpolation of
+                                various kinds.  interp2d is a wrapper
+                                around this class.
         
-        Spline 2d  :  
+        Spline2d  :  an object for spline interpolation.  Interpolate1d
+                                wraps this class if spline interpolation
+                                is used.  However, not all functionality
+                                of Spline2d is available through Interpolate1d.
         
-        InterpolateNd  :
+        InterpolateNd  :  an object for interpolation of ND data.  interpNd
+                                is a wrapper around this class.
         
     These functions and classes constitute the primary api.  However, several
-    additional functions are also provided (the primary api in many cases calls
-    these functions):
+    additional less generic functions are also provided for more direct interpolation.
+    They don't have as crafted a UI as those above, but they are quick and easy:
 
         linear : linear interpolation
         logarithmic :  logarithmic interpolation

Modified: branches/Interpolate1D/interpolateNd.py
===================================================================
--- branches/Interpolate1D/interpolateNd.py	2008-08-07 21:25:16 UTC (rev 4610)
+++ branches/Interpolate1D/interpolateNd.py	2008-08-07 23:42:44 UTC (rev 4611)
@@ -132,6 +132,7 @@
         """
         
         # FIXME : include spline filtering
+        # the ndimage module says that it requires pre-filtering for 
         
         # checking format of input
         data = array(data)
@@ -191,8 +192,16 @@
         else:
             raise ValueError, "argument kind = %s not recognized" % str(kind)
                 
+        
+        # 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)
+        else:
+            self._data_array = data
+        
         # storing relevant data
-        self._data_array = data
         self.ndim = data.ndim
         self._shape = np.shape(data)
         self._spacings = spacings
@@ -295,4 +304,40 @@
         
         
     
-    
\ No newline at end of file
+    
+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.
+
+        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




More information about the Scipy-svn mailing list