在 x、y 和 z 中具有不同间隔的定期采样 3D 数据的快速插值 [英] 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])

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

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

举个例子:

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.

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