I have some volumetric imaging data consisting of values sampled on a regular grid in x,y,z, but with a non-cubic voxel shape (the space between adjacent points in z is greater than in x,y). I would eventually like to be able to interpolate the values on some arbitrary 2D plane that passes through the volume, like this:
I'm aware of scipy.ndimage.map_coordinates
, but in my case using it is less straightforward because it implicitly assumes that the spacing of the elements in the input array are equal across dimensions. I could first resample my input array according to the smallest voxel dimension (so that all of my voxels would then be cubes), then use map_coordinates
to interpolate over my plane, but it doesn't seem like a great idea to interpolate my data twice.
I'm also aware that scipy
has various interpolators for irregularly-spaced ND data (LinearNDInterpolator
, NearestNDInterpolator
etc.), but these are very slow and memory-intensive for my purposes. What is the best way of interpolating my data given that I know that the values are regularly spaced within each dimension?
You can use map_coordinates
with a little bit of algebra. Lets say the spacings of your grid are dx
, dy
and dz
. We need to map these real world coordinates to array index coordinates, so lets define three new variables:
xx = x / dx
yy = y / dy
zz = z / dz
The array index input to map_coordinates
is an array of shape (d, ...)
where d
is the number of dimensions of your original data. If you define an array such as:
scaling = np.array([dx, dy, dz])
you can transform your real world coordinates to array index coordinates by dividing by scaling
with a little broadcasting magic:
idx = coords / scaling[(slice(None),) + (None,)*(coords.ndim-1)]
To put it all together in an example:
dx, dy, dz = 1, 1, 2
scaling = np.array([dx, dy, dz])
data = np.random.rand(10, 15, 5)
Lets say we want to interpolate values along the plane 2*y - z = 0
. We take two vectors perpendicular to the planes normal vector:
u = np.array([1, 0 ,0])
v = np.array([0, 1, 2])
And get the coordinates at which we want to interpolate as:
coords = (u[:, None, None] * np.linspace(0, 9, 10)[None, :, None] +
v[:, None, None] * np.linspace(0, 2.5, 10)[None, None, :])
We convert them to array index coordinates and interpoalte using map_coordinates
:
idx = coords / scaling[(slice(None),) + (None,)*(coords.ndim-1)]
new_data = ndi.map_coordinates(data, idx)
This last array is of shape (10, 10)
and has in position [u_idx, v_idx]
the value corresponding to the coordinate coords[:, u_idx, v_idx]
.
You could build on this idea to handle interpolation where your coordinates don't start at zero, by adding an offset before the scaling.