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 aredx
,dy
anddz
. We need to map these real world coordinates to array index coordinates, so lets define three new variables:The array index input to
map_coordinates
is an array of shape(d, ...)
whered
is the number of dimensions of your original data. If you define an array such as:you can transform your real world coordinates to array index coordinates by dividing by
scaling
with a little broadcasting magic:To put it all together in an example:
Lets say we want to interpolate values along the plane
2*y - z = 0
. We take two vectors perpendicular to the planes normal vector:And get the coordinates at which we want to interpolate as:
We convert them to array index coordinates and interpoalte using
map_coordinates
:This last array is of shape
(10, 10)
and has in position[u_idx, v_idx]
the value corresponding to the coordinatecoords[:, 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.