python - 在 x、y 和 z 中具有不同间隔的定期采样 3D 数据的快速插值

Question

我有一些体积成像数据，其中包含在 x、y、z 中的规则网格上采样的值，但具有非立方体素形状(z 中相邻点之间的空间大于 x、y 中的空间)。我最终希望能够在穿过体积的任意 2D 平面上插值，如下所示:

我知道 scipy.ndimage.map_coordinates，但在我的情况下使用它并不那么简单，因为它隐含地假设输入数组中元素的间距在各个维度上是相等的。我可以首先根据最小体素维度对输入数组重新采样(这样我的所有体素都将成为立方体)，然后使用 map_coordinates 在我的平面上进行插值，但它看起来不像插入我的数据两次的好主意。

我也知道 scipy 具有用于不规则间隔 ND 数据的各种插值器(LinearNDInterpolator、NearestNDInterpolator 等)，但这些就我的目的而言，它们非常缓慢且占用大量内存。鉴于我知道值 在每个维度内有规律地间隔，插入数据的最佳方法是什么？

Answer 1

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

xx = x / dx
yy = y / dy
zz = z / dz

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

scaling = np.array([dx, dy, dz])

您可以通过 scaling 使用一点广播魔法将您的真实世界坐标转换为数组索引坐标:

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 插入值。我们取两个垂直于平面法线向量的向量:

u = np.array([1, 0 ,0])
v = np.array([0, 1, 2])

并得到我们想要插值的坐标:

coords = (u[:, None, None] * np.linspace(0, 9, 10)[None, :, None] +
          v[:, None, None] * np.linspace(0, 2.5, 10)[None, None, :])

我们将它们转换为 数组索引 坐标并使用 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].

您可以基于这个想法来处理坐标不从零开始的插值，方法是在缩放之前添加一个偏移量。