interpolate 3D volume with numpy and or scipy
Solution 1
In scipy 0.14 or later, there is a new function scipy.interpolate.RegularGridInterpolator which closely resembles interp3.
The MATLAB command Vi = interp3(x,y,z,V,xi,yi,zi) would translate to something like:
from numpy import array
from scipy.interpolate import RegularGridInterpolator as rgi
my_interpolating_function = rgi((x,y,z), V)
Vi = my_interpolating_function(array([xi,yi,zi]).T)
Here is a full example demonstrating both; it will help you understand the exact differences...
MATLAB CODE:
x = linspace(1,4,11);
y = linspace(4,7,22);
z = linspace(7,9,33);
V = zeros(22,11,33);
for i=1:11
for j=1:22
for k=1:33
V(j,i,k) = 100*x(i) + 10*y(j) + z(k);
end
end
end
xq = [2,3];
yq = [6,5];
zq = [8,7];
Vi = interp3(x,y,z,V,xq,yq,zq);
The result is Vi=[268 357] which is indeed the value at those two points (2,6,8) and (3,5,7).
SCIPY CODE:
from scipy.interpolate import RegularGridInterpolator
from numpy import linspace, zeros, array
x = linspace(1,4,11)
y = linspace(4,7,22)
z = linspace(7,9,33)
V = zeros((11,22,33))
for i in range(11):
for j in range(22):
for k in range(33):
V[i,j,k] = 100*x[i] + 10*y[j] + z[k]
fn = RegularGridInterpolator((x,y,z), V)
pts = array([[2,6,8],[3,5,7]])
print(fn(pts))
Again it's [268,357]. So you see some slight differences: Scipy uses x,y,z index order while MATLAB uses y,x,z (strangely); In Scipy you define a function in a separate step and when you call it, the coordinates are grouped like (x1,y1,z1),(x2,y2,z2),... while matlab uses (x1,x2,...),(y1,y2,...),(z1,z2,...).
Other than that, the two are similar and equally easy to use.
Solution 2
The exact equivalent to MATLAB's interp3 would be using scipy's interpn for one-off interpolation:
import numpy as np
from scipy.interpolate import interpn
Vi = interpn((x,y,z), V, np.array([xi,yi,zi]).T)
The default method for both MATLAB and scipy is linear interpolation, and this can be changed with the method argument. Note that only linear and nearest-neighbor interpolation is supported by interpn for 3 dimensions and above, unlike MATLAB which supports cubic and spline interpolation as well.
When making multiple interpolation calls on the same grid it is preferable to use the interpolation object RegularGridInterpolator, as in the accepted answer above. interpn uses RegularGridInterpolator internally.
Solution 3
Basically, ndimage.map_coordinates works in "index" coordinates (a.k.a. "voxel" or "pixel" coordinates). The interface to it seems a bit clunky at first, but it does give you a lot of flexibility.
If you want to specify the interpolated coordinates similar to matlab's interp3, then you'll need to convert your intput coordinates into "index" coordinates.
There's also the additional wrinkle that map_coordinates always preserves the dtype of the input array in the output. If you interpolate an integer array, you'll get integer output, which may or may not be what you want. For the code snippet below, I'll assume that you always want floating point output. (If you don't, it's actually simpler.)
I'll try to add more explanation later tonight (this is rather dense code).
All in all, the interp3 function I have is more complex than it may need to be for your exact purposes. Howver, it should more or less replicate the behavior of interp3 as I remember it (ignoring the "zooming" functionality of interp3(data, zoom_factor), which scipy.ndimage.zoom handles.)
import numpy as np
from scipy.ndimage import map_coordinates
def main():
data = np.arange(5*4*3).reshape(5,4,3)
x = np.linspace(5, 10, data.shape[0])
y = np.linspace(10, 20, data.shape[1])
z = np.linspace(-100, 0, data.shape[2])
# Interpolate at a single point
print interp3(x, y, z, data, 7.5, 13.2, -27)
# Interpolate a region of the x-y plane at z=-25
xi, yi = np.mgrid[6:8:10j, 13:18:10j]
print interp3(x, y, z, data, xi, yi, -25 * np.ones_like(xi))
def interp3(x, y, z, v, xi, yi, zi, **kwargs):
"""Sample a 3D array "v" with pixel corner locations at "x","y","z" at the
points in "xi", "yi", "zi" using linear interpolation. Additional kwargs
are passed on to ``scipy.ndimage.map_coordinates``."""
def index_coords(corner_locs, interp_locs):
index = np.arange(len(corner_locs))
if np.all(np.diff(corner_locs) < 0):
corner_locs, index = corner_locs[::-1], index[::-1]
return np.interp(interp_locs, corner_locs, index)
orig_shape = np.asarray(xi).shape
xi, yi, zi = np.atleast_1d(xi, yi, zi)
for arr in [xi, yi, zi]:
arr.shape = -1
output = np.empty(xi.shape, dtype=float)
coords = [index_coords(*item) for item in zip([x, y, z], [xi, yi, zi])]
map_coordinates(v, coords, order=1, output=output, **kwargs)
return output.reshape(orig_shape)
main()
Solution 4
The question is old, but I think it needs some clarification, as nobody pointed out that the operation requested (trilinear interpolation) can be easily implemented from scratch with consistent saving of computational time (around 10 times faster) w.r.t. scipy.interpolate's RegularGridInterpolator.
Code
import numpy as np
from itertools import product
def trilinear_interpolation(x_volume, y_volume, z_volume, volume, x_needed, y_needed, z_needed):
"""
Trilinear interpolation (from Wikipedia)
:param x_volume: x points of the volume grid
:type crack_type: list or numpy.ndarray
:param y_volume: y points of the volume grid
:type crack_type: list or numpy.ndarray
:param x_volume: z points of the volume grid
:type crack_type: list or numpy.ndarray
:param volume: volume
:type crack_type: list or numpy.ndarray
:param x_needed: desired x coordinate of volume
:type crack_type: float
:param y_needed: desired y coordinate of volume
:type crack_type: float
:param z_needed: desired z coordinate of volume
:type crack_type: float
:return volume_needed: desired value of the volume, i.e. volume(x_needed, y_needed, z_needed)
:type volume_needed: float
"""
# dimensinoal check
if np.shape(volume) != (len(x_volume), len(y_volume), len(z_volume)):
raise ValueError(f'dimension mismatch, volume must be a ({len(x_volume)}, {len(y_volume)}, {len(z_volume)}) list or numpy.ndarray')
# check of the indices needed for the correct control volume definition
i = searchsorted(x_volume, x_needed)
j = searchsorted(y_volume, y_needed)
k = searchsorted(z_volume, z_needed)
# control volume definition
control_volume_coordinates = np.array(
[[x_volume[i - 1], y_volume[j - 1], z_volume[k - 1]], [x_volume[i], y_volume[j], z_volume[k]]])
xd = (np.array([x_needed, y_needed, z_needed]) - control_volume_coordinates[0]) / (control_volume_coordinates[1] - control_volume_coordinates[0])
# interpolation along x
c2 = [[0, 0], [0, 0]]
for m, n in product([0, 1], [0, 1]):
c2[m][n] = volume[i - 1][j - 1 + m][k - 1 + n] * (1 - xd[0]) + volume[i][j - 1 + m][k - 1 + n] * xd[0]
# interpolation along y
c1 = [0, 0]
c1[0] = c2[0][0] * (1 - xd[1]) + c2[1][0] * xd[1]
c1[1] = c2[0][1] * (1 - xd[1]) + c2[1][1] * xd[1]
# interpolation along z
volume_needed = c1[0] * (1 - xd[2]) + c1[1] * xd[2]
return volume_needed
def searchsorted(l, x):
for i in l:
if i >= x: break
return l.index(i)
from scipy.interpolate import RegularGridInterpolator
def trilin_interp_regular_grid(x_volume, y_volume, z_volume, volume, x_needed, y_needed, z_needed):
# dimensinoal check
if np.shape(volume) != (len(x_volume), len(y_volume), len(z_volume)):
raise ValueError(f'dimension mismatch, volume must be a ({len(x_volume)}, {len(y_volume)}, {len(z_volume)}) list or numpy.ndarray')
# trilinear interpolation on a regular grid
fn = RegularGridInterpolator((x_volume,y_volume,z_volume), volume)
volume_needed = fn(np.array([x_needed, y_needed, z_needed]))
return volume_needed
Edit
For the sake of speed, you can compile the code provided above using @jit:
import numpy as np
from numba import njit
@njit(fastmath=True)
def trilinear_interpolation_jit(
x_volume,
y_volume,
z_volume,
volume,
x_needed,
y_needed,
z_needed
):
"""
Trilinear interpolation (from Wikipedia)
:param x_volume: x points of the volume grid
:type crack_type: list or numpy.ndarray
:param y_volume: y points of the volume grid
:type crack_type: list or numpy.ndarray
:param x_volume: z points of the volume grid
:type crack_type: list or numpy.ndarray
:param volume: volume
:type crack_type: list or numpy.ndarray
:param x_needed: desired x coordinate of volume
:type crack_type: float
:param y_needed: desired y coordinate of volume
:type crack_type: float
:param z_needed: desired z coordinate of volume
:type crack_type: float
:return volume_needed: desired value of the volume, i.e. volume(x_needed, y_needed, z_needed)
:type volume_needed: float
"""
# dimensinoal check
assert np.shape(volume) == (
len(x_volume), len(y_volume), len(z_volume)
), "Incompatible lengths"
# check of the indices needed for the correct control volume definition
i = np.searchsorted(x_volume, x_needed)
j = np.searchsorted(y_volume, y_needed)
k = np.searchsorted(z_volume, z_needed)
# control volume definition
control_volume_coordinates = np.array(
[
[
x_volume[i - 1],
y_volume[j - 1],
z_volume[k - 1]
],
[
x_volume[i],
y_volume[j],
z_volume[k]
]
]
)
xd = (
np.array([x_needed, y_needed, z_needed]) - control_volume_coordinates[0]
) / (
control_volume_coordinates[1] - control_volume_coordinates[0]
)
# interpolation along x
c2 = [[0., 0.], [0., 0.]]
for m, n in [(0, 0), (0, 1), (1, 0), (1, 1)]:
c2[m][n] = volume[i - 1][j - 1 + m][k - 1 + n] \
* (1. - xd[0]) + volume[i][j - 1 + m][k - 1 + n] * xd[0]
# interpolation along y
c1 = [0., 0.]
c1[0] = c2[0][0] * (1. - xd[1]) + c2[1][0] * xd[1]
c1[1] = c2[0][1] * (1. - xd[1]) + c2[1][1] * xd[1]
# interpolation along z
volume_needed = c1[0] * (1. - xd[2]) + c1[1] * xd[2]
return volume_needed
@njit(fastmath=True)
def trilint_jit(
x_volume,
y_volume,
z_volume,
volume,
x_needed,
y_needed,
z_needed
):
trilint_size = x_needed.size * y_needed.size * z_needed.size
jitted_trilint = np.zeros(trilint_size)
m = 0
for x in range(0, len(x_needed)):
for y in range(0, len(y_needed)):
for z in range(0, len(z_needed)):
jitted_trilint[m]=trilinear_interpolation_jit(
x_volume,
y_volume,
z_volume,
the_volume,
x_needed[x],
y_needed[y],
z_needed[z]
)
m = m + 1
return jitted_trilint
Example
import numpy as np
import time
x_volume = np.array([100., 1000.])
y_volume = np.array([0.2, 0.4, 0.6, 0.8, 1])
z_volume = np.array([0, 0.2, 0.5, 0.8, 1.])
the_volume = np.array(
[[[0.902, 0.985, 1.12, 1.267, 1.366],
[0.822, 0.871, 0.959, 1.064, 1.141],
[0.744, 0.77, 0.824, 0.897, 0.954],
[0.669, 0.682, 0.715, 0.765, 0.806],
[0.597, 0.607, 0.631, 0.667, 0.695]],
[[1.059, 1.09, 1.384, 1.682, 1.881],
[0.948, 0.951, 1.079, 1.188, 1.251],
[0.792, 0.832, 0.888, 0.940, 0.971],
[0.726, 0.733, 0.754, 0.777, 0.792],
[0.642, 0.656, 0.675, 0.691, 0.700]]])
x_needed = np.linspace(100, 1000, 10)
y_needed = np.linspace(0.3, 1, 60)
z_needed = np.linspace(0, 1, 7)
start = time.time()
jitted_trilint = trilint_jit(
x_volume, y_volume, z_volume, the_volume, x_needed, y_needed, z_needed
)
end = time.time()
print('---')
print(f"NUMBA: {end - start}")
print('---')
start = time.time()
manual_trilint = []
for x in x_needed:
for y in y_needed:
for z in z_needed:
manual_trilint.append(
trilinear_interpolation(
x_volume, y_volume, z_volume, the_volume, x, y, z
)
)
end = time.time()
print('---')
print(f"Manual: {end - start}")
print('---')
start = time.time()
auto_trilint = []
for x in x_needed:
for y in y_needed:
for z in z_needed:
auto_trilint.append(
trilin_interp_regular_grid(
x_volume, y_volume, z_volume, the_volume, x, y, z
)
)
end = time.time()
print('---')
print(f"Auto: {end - start}")
print('---')
Which on my laptop outputs:
---
NUMBA: 0.007552146911621094
---
---
Manual: 0.06470489501953125
---
---
Auto: 1.1136455535888672
---
user1301295
Updated on October 30, 2020Comments
-
user1301295 over 3 years
I am extremely frustrated because after several hours I can't seem to be able to do a seemingly easy 3D interpolation in python. In Matlab all I had to do was
Vi = interp3(x,y,z,V,xi,yi,zi)What is the exact equivalent of this using scipy's ndimage.map_coordinate or other numpy methods?
Thanks
-
BayesianMonk over 3 yearsThis is actually not very time efficient when it comes to interpolate a large number of points. This could benefit from broadcasting to gain a large speed up -
Pietro D'Antuono about 2 years@BayesianMonk you are right. The other alternative for speed-up is Numba, for which I added an edit to my answer.