Interpolation/subsampling of 3D data in python without VTK

Question

What I want to do is rather simple but I havent found a straightforward approach thus far:

I have a 3D rectilinear grid with float values (therefore 3 coordinate axes -1D numpy arrays- for the centers of the grid cells and a 3D numpy array with the corresponding shape with a value for each cell center) and I want to interpolate (or you may call it subsample) this entire array to a subsampled array (e.g. size factor of 5) with linear interpolation. All the approaches I've seen this far involve 2D and then 1D interpolation or VTK tricks which Id rather not use (portability).

Could someone suggest an approach that would be the equivalent of taking 5x5x5 cells at the same time in the 3D array, averaging and returning an array 5times smaller in each direction?

Thank you in advance for any suggestions

EDIT: Here's what the data looks like, 'd' is a 3D array representing a 3D grid of cells. Each cell has a scalar float value (pressure in my case) and 'x','y' and 'z' are three 1D arrays containing the spatial coordinates of the cells of every cell (see the shapes and how the 'x' array looks like)

In [42]: x.shape
Out[42]: (181L,)

In [43]: y.shape
Out[43]: (181L,)

In [44]: z.shape
Out[44]: (421L,)

In [45]: d.shape
Out[45]: (181L, 181L, 421L)

In [46]: x
Out[46]: 
array([-0.410607  , -0.3927568 , -0.37780656, -0.36527296, -0.35475321,
       -0.34591168, -0.33846866, -0.33219107, -0.32688467, -0.3223876 ,
        ...
        0.34591168,  0.35475321,  0.36527296,  0.37780656,  0.3927568 ,
        0.410607  ])

What I want to do is create a 3D array with lets say a shape of 90x90x210 (roughly downsize by a factor of 2) by first subsampling the coordinates from the axes on arrays with the above dimensions and then 'interpolating' the 3D data to that array. Im not sure whether 'interpolating' is the right term though. Downsampling? Averaging? Here's an 2D slice of the data:

Answer 1

Here is an example of 3D interpolation on an irregular grid using scipy.interpolate.griddata.

import numpy as np
import scipy.interpolate as interpolate
import matplotlib.pyplot as plt


def func(x, y, z):
    return x ** 2 + y ** 2 + z ** 2

# Nx, Ny, Nz = 181, 181, 421
Nx, Ny, Nz = 18, 18, 42

subsample = 2
Mx, My, Mz = Nx // subsample, Ny // subsample, Nz // subsample

# Define irregularly spaced arrays
x = np.random.random(Nx)
y = np.random.random(Ny)
z = np.random.random(Nz)

# Compute the matrix D of shape (Nx, Ny, Nz).
# D could be experimental data, but here I'll define it using func
# D[i,j,k] is associated with location (x[i], y[j], z[k])
X_irregular, Y_irregular, Z_irregular = (
    x[:, None, None], y[None, :, None], z[None, None, :])
D = func(X_irregular, Y_irregular, Z_irregular)

# Create a uniformly spaced grid
xi = np.linspace(x.min(), x.max(), Mx)
yi = np.linspace(y.min(), y.max(), My)
zi = np.linspace(y.min(), y.max(), Mz)
X_uniform, Y_uniform, Z_uniform = (
    xi[:, None, None], yi[None, :, None], zi[None, None, :])

# To use griddata, I need 1D-arrays for x, y, z of length 
# len(D.ravel()) = Nx*Ny*Nz.
# To do this, I broadcast up my *_irregular arrays to each be 
# of shape (Nx, Ny, Nz)
# and then use ravel() to make them 1D-arrays
X_irregular, Y_irregular, Z_irregular = np.broadcast_arrays(
    X_irregular, Y_irregular, Z_irregular)
D_interpolated = interpolate.griddata(
    (X_irregular.ravel(), Y_irregular.ravel(), Z_irregular.ravel()),
    D.ravel(),
    (X_uniform, Y_uniform, Z_uniform),
    method='linear')

print(D_interpolated.shape)
# (90, 90, 210)

# Make plots
fig, ax = plt.subplots(2)

# Choose a z value in the uniform z-grid
# Let's take the middle value
zindex = Mz // 2
z_crosssection = zi[zindex]

# Plot a cross-section of the raw irregularly spaced data
X_irr, Y_irr = np.meshgrid(sorted(x), sorted(y))
# find the value in the irregular z-grid closest to z_crosssection
z_near_cross = z[(np.abs(z - z_crosssection)).argmin()]
ax[0].contourf(X_irr, Y_irr, func(X_irr, Y_irr, z_near_cross))
ax[0].scatter(X_irr, Y_irr, c='white', s=20)   
ax[0].set_title('Cross-section of irregular data')
ax[0].set_xlim(x.min(), x.max())
ax[0].set_ylim(y.min(), y.max())

# Plot a cross-section of the Interpolated uniformly spaced data
X_unif, Y_unif = np.meshgrid(xi, yi)
ax[1].contourf(X_unif, Y_unif, D_interpolated[:, :, zindex])
ax[1].scatter(X_unif, Y_unif, c='white', s=20)
ax[1].set_title('Cross-section of downsampled and interpolated data')
ax[1].set_xlim(x.min(), x.max())
ax[1].set_ylim(y.min(), y.max())

plt.show()

Answer 2

In short: doing interpolation in each dimension separately is the right way to go.

---

You can simply average every 5x5x5 cube and return the results. However, if your data is supposed to be continuous, you should understand that is not good subsampling practice, as it will likely induce aliasing. (Also, you can't reasonably call it "interpolation"!)

Good resampling filters need to be wider than the resampling factor in order to avoid aliasing. Since you are downsampling, you should also realize that your resampling filter needs to be scaled according to the destination resolution, not the original resolution -- in order to interpolate properly, it will likely need to be 4 or 5 times as wide as your 5x5x5 cube. This is a lot of samples -- 20*20*20 is way more than 5*5*5...

So, the reason why practical implementations of resampling typically filter each dimension separately is that it is more efficient. By taking 3 passes, you can evaluate your filter using far fewer multiply/accumulate operations per output sample.