Sure! There are two options that do different things, but both use a regular source data structure.
The first is scipy.ndimage.zoom . If you just want to create a denser regular grid based on the interpolation of the source data, this is the way to go.
Second scipy.ndimage.map_coordinates . If you want to interpolate several (or many) arbitrary points in your data, but still use the regular nature of the source data (for example, you don't need a quadrant), this is the way to go.
As a quick example (this will use cubic interpolation. Use order=1 for bilinear, order=0 for nearest, etc.):
import numpy as np import scipy.ndimage as ndimage data = np.arange(9).reshape(3,3) print 'Original:\n', data print 'Zoomed by 2x:\n', ndimage.zoom(data, 2)
This gives:
Original: [[0 1 2] [3 4 5] [6 7 8]] Zoomed by 2x: [[0 0 1 1 2 2] [1 1 1 2 2 3] [2 2 3 3 4 4] [4 4 5 5 6 6] [5 6 6 7 7 7] [6 6 7 7 8 8]]
This also works for 3D (and nD) arrays. However, keep in mind that if you zoom in 2x, for example, you will scale along all axes.
data = np.arange(27).reshape(3,3,3) print 'Original:\n', data print 'Zoomed by 2x gives an array of shape:', ndimage.zoom(data, 2).shape
This gives:
Original: [[[ 0 1 2] [ 3 4 5] [ 6 7 8]] [[ 9 10 11] [12 13 14] [15 16 17]] [[18 19 20] [21 22 23] [24 25 26]]] Zoomed by 2x gives an array of shape: (6, 6, 6)
If you have something like a 3-way RGB image that you want to enlarge, you can do this by specifying a sequence of tuples as a zoom factor:
print 'Zoomed by 2x along the last two axes:' print ndimage.zoom(data, (1, 2, 2))
This gives:
Zoomed by 2x along the last two axes: [[[ 0 0 1 1 2 2] [ 1 1 1 2 2 3] [ 2 2 3 3 4 4] [ 4 4 5 5 6 6] [ 5 6 6 7 7 7] [ 6 6 7 7 8 8]] [[ 9 9 10 10 11 11] [10 10 10 11 11 12] [11 11 12 12 13 13] [13 13 14 14 15 15] [14 15 15 16 16 16] [15 15 16 16 17 17]] [[18 18 19 19 20 20] [19 19 19 20 20 21] [20 20 21 21 22 22] [22 22 23 23 24 24] [23 24 24 25 25 25] [24 24 25 25 26 26]]]
Arbitrary interpolation of data with a regular grid using map_coordinates
The first thing about map_coordinates is that it works in pixel coordinates (for example, just like you indexed an array, but values can be floats). From your description, this is exactly what you want, but often confuses people. For example, if you have the "real-world" x, y, z coordinates, you will need to convert them to the "pixel" coordinates based on the index.
In any case, let's say we wanted to interpolate the value in the original array at positions 1.2, 0.3, 1.4.
If you think about this in the earlier case of the RGB image, the first coordinate corresponds to the "strip", the second to the "row", and the last to the "column". Which order corresponds to the one that completely depends on how you decide to structure your data, but I will use them as the coordinates "z, y, x", as this simplifies the visualization of the comparison with the print array.
import numpy as np import scipy.ndimage as ndimage data = np.arange(27).reshape(3,3,3) print 'Original:\n', data print 'Sampled at 1.2, 0.3, 1.4:' print ndimage.map_coordinates(data, [[1.2], [0.3], [1.4]])
This gives:
Original: [[[ 0 1 2] [ 3 4 5] [ 6 7 8]] [[ 9 10 11] [12 13 14] [15 16 17]] [[18 19 20] [21 22 23] [24 25 26]]] Sampled at 1.2, 0.3, 1.4: [14]
Again, this is the default cubic interpolation. Use order kwarg to control the type of interpolation.
It is worth noting here that all scipy.ndimage operations save the dtype of the original array. If you want to get floating point results, you will need to distinguish the original array as a float:
In [74]: ndimage.map_coordinates(data.astype(float), [[1.2], [0.3], [1.4]]) Out[74]: array([ 13.5965])
Another thing you may notice is that the format of the interpolated coordinates is rather cumbersome for a single point (for example, it expects a 3xN array instead of an Nx3 array). However, this is possibly better when you have a sequence of coordinates. For example, consider the case of sampling along a line that passes through a “cube” of data:
xi = np.linspace(0, 2, 10) yi = 0.8 * xi zi = 1.2 * xi print ndimage.map_coordinates(data, [zi, yi, xi])
This gives:
[ 0 1 4 8 12 17 21 24 0 0]
It is also a good place to mention how boundary conditions are handled. By default, everything outside the array is set to 0. Thus, the last two values in the sequence are 0 . (i.e., zi > 2 for the last two elements).
If we wanted the points outside the array to be, say -999 (we cannot use nan , since it is an integer array. If you want nan , you will need to reset to float.):
In [75]: ndimage.map_coordinates(data, [zi, yi, xi], cval=-999) Out[75]: array([ 0, 1, 4, 8, 12, 17, 21, 24, -999, -999])
If we wanted to return the closest value to points outside the array, we would do:
In [76]: ndimage.map_coordinates(data, [zi, yi, xi], mode='nearest') Out[76]: array([ 0, 1, 4, 8, 12, 17, 21, 24, 25, 25])
You can also use "reflect" and "wrap" as border modes in addition to the "nearest" and the default is "constant" . They are pretty clear, but try experimenting a bit if you're confused.
For example, let interpolate a line along the first line of the first strip in the array, which extends twice the distance of the array:
xi = np.linspace(0, 5, 10) yi, zi = np.zeros_like(xi), np.zeros_like(xi)
Default value:
In [77]: ndimage.map_coordinates(data, [zi, yi, xi]) Out[77]: array([0, 0, 1, 2, 0, 0, 0, 0, 0, 0])
Compare this to:
In [78]: ndimage.map_coordinates(data, [zi, yi, xi], mode='reflect') Out[78]: array([0, 0, 1, 2, 2, 1, 2, 1, 0, 0]) In [78]: ndimage.map_coordinates(data, [zi, yi, xi], mode='wrap') Out[78]: array([0, 0, 1, 2, 0, 1, 1, 2, 0, 1])
Hope this clears things up a bit!