快速插值定期采样的3D数据(在x，y和z中具有不同的间隔) [英] Fast interpolation of regularly sampled 3D data with different intervals in x,y, and z

问题描述

我有一些体积成像数据，这些数据由在x，y，z的常规网格上采样的值组成，但具有非三次体素形状(z的相邻点之间的间距大于x，y的间距).我最终希望能够在穿过体积的任意2D平面上插值，如下所示:

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:

我知道scipy.ndimage.map_coordinates，但是在我的情况下使用它不太直接，因为它隐式地假设输入数组中元素的间距在维度上相等.我可以首先根据最小的体素尺寸对输入数组进行重新采样(这样所有的体素都将成为立方体)，然后使用map_coordinates进行平面插值，但是对我的平面进行插值似乎不是一个好主意数据两次.

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.

我也知道scipy对于不规则间隔的ND数据(LinearNDInterpolator，NearestNDInterpolator等)有各种插值器，但是对于我来说，这些插值器非常慢且占用大量内存.假设我知道值 在每个维度上有规律的间隔，插值数据的最佳方法是什么?

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?

推荐答案

您可以将map_coordinates与一些代数结合使用.假设网格的间距为dx，dy和dz.我们需要将这些真实世界坐标映射到数组索引坐标，因此让我们定义三个新变量:

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

输入到map_coordinates的数组索引是形状为(d, ...)的数组，其中d是原始数据的维数.如果您定义一个数组，例如:

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])

您可以将现实世界坐标转换为数组索引坐标，方法是用一点广播魔术除以scaling:

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)

让我们说我们想沿平面2*y - z = 0插值.我们取两个垂直于平面法线向量的向量:

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, :])

我们将它们转换为数组索引坐标，并使用map_coordinates进行解释:

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)

最后一个数组的形状为(10, 10)，并且在位置[u_idx, v_idx]中具有与坐标coords[:, u_idx, v_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.

这篇关于快速插值定期采样的3D数据(在x，y和z中具有不同的间隔)的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！