From 7daae202d9a1ec7928603a65f5c06b4518c3596e Mon Sep 17 00:00:00 2001
From: Anton Khirnov <anton@khirnov.net>
Date: Tue, 17 Dec 2013 10:28:12 +0100
Subject: Add a convenience wrapper around scipy ndimage interpolation.

---
 cactus_utils.py | 33 +++++++++++++++++++++++++++++++++
 1 file changed, 33 insertions(+)

diff --git a/cactus_utils.py b/cactus_utils.py
index 100f31f..af5959e 100644
--- a/cactus_utils.py
+++ b/cactus_utils.py
@@ -1,6 +1,7 @@
 import h5py
 import numpy as np
 import pylab
+import scipy.ndimage.interpolation
 import scipy.interpolate
 import scipy.integrate
 import scipy.optimize
@@ -16,6 +17,38 @@ axis_data = { 'x' : 2, 'y' : 1, 'z' : 0 }
 EPSILON   = 1e-3
 
 
+def ndinterpolate(in_coords, values, out_coords):
+    """
+    A wrapper around scipy.image.interpolation.map_coordinates() providing more convenient API.
+
+    Given the values of a function on an N-dimensional regular grid, interpolate
+    its value on given coordinates.
+
+    @param in_coords An N+1-dimensional array of the coordinates of the input values.
+                     The first axis goes over the dimensions (i.e. in_coords[0] is
+                     an N-dimensional array that gives the x-coordinate for each value).
+    @param values An N-dimensional array of the values (at in_coords) of the function to
+                  be interpolated
+    @param out_coords A 2-dimensional array of the coordinates at which the output is
+                      desired. The first axis goes over the dimensions (i.e. out_coords[0]
+                      is a 1-dimensional array of length N of the x-coordinates for every
+                      desired interpolation point).
+    """
+    N = len(values.shape)
+
+    if len(in_coords.shape) != N + 1 or in_coords.shape[0] != N:
+        raise ValueError('Invalid coords/values shapes: ' + str(in_coords.shape) + str(values.shape))
+
+    out_coords = np.copy(out_coords)
+    for i in xrange(N):
+        max = np.max(in_coords[i])
+        min = np.min(in_coords[i])
+        out_coords[i] -= min
+        out_coords[i] *= values.shape[i] / (max - min)
+
+    return scipy.ndimage.interpolation.map_coordinates(values, out_coords)
+
+
 def get_default_gridfunc(datafile):
     """
     Get the name of the default grid function from a given data file
-- 
cgit v1.2.3