网格数据的快速插值

Posted 2023-02-16

【中文标题】网格数据的快速插值

【英文标题】：Fast interpolation of grid data

【发布时间】：2013-06-03 17:53:27

【问题描述】：

我有一个大的 3d np.ndarray 数据，表示以规则网格方式在体积上采样的物理变量（如数组 [0,0,0] 中的值表示物理坐标（0， 0,0))。

我想通过在粗糙网格中插入数据来获得更精细的网格间距。目前我正在使用 scipy griddata 线性插值，但它非常慢（20x20x20 阵列约 90 秒）。就我的目的而言，它有点过度设计，允许对体积数据进行随机采样。有什么东西可以利用我定期间隔的数据以及我想要插值的特定点集有限的事实吗？

【问题讨论】：

fast-interpolation-of-regularly-sampled-3d-data-with-different-intervals-in-x-y 也使用 map_coordinates

【参考方案1】：

当然！有两个选项可以做不同的事情，但都利用了原始数据的规则网格特性。

第一个是scipy.ndimage.zoom。如果您只是想在对原始数据进行插值的基础上生成更密集的规则网格，那么这是可行的方法。

第二个是scipy.ndimage.map_coordinates。如果您想在数据中插入几个（或多个）任意点，但仍要利用原始数据的规则网格特性（例如，不需要四叉树），那么这是可行的方法。

---

“缩放”一个数组 (scipy.ndimage.zoom)

举个简单的例子（这将使用三次插值。使用order=1 表示双线性，order=0 表示最近等）：

import numpy as np
import scipy.ndimage as ndimage

data = np.arange(9).reshape(3,3)

print 'Original:\n', data
print 'Zoomed by 2x:\n', ndimage.zoom(data, 2)

这会产生：

Original:
[[0 1 2]
 [3 4 5]
 [6 7 8]]
Zoomed by 2x:
[[0 0 1 1 2 2]
 [1 1 1 2 2 3]
 [2 2 3 3 4 4]
 [4 4 5 5 6 6]
 [5 6 6 7 7 7]
 [6 6 7 7 8 8]]

这也适用于 3D（和 nD）数组。但是，请注意，例如，如果您缩放 2 倍，您将沿 所有 轴进行缩放。

data = np.arange(27).reshape(3,3,3)
print 'Original:\n', data
print 'Zoomed by 2x gives an array of shape:', ndimage.zoom(data, 2).shape

这会产生：

Original:
[[[ 0  1  2]
  [ 3  4  5]
  [ 6  7  8]]

 [[ 9 10 11]
  [12 13 14]
  [15 16 17]]

 [[18 19 20]
  [21 22 23]
  [24 25 26]]]
Zoomed by 2x gives an array of shape: (6, 6, 6)

如果您想要缩放 3 波段 RGB 图像，您可以通过指定元组序列作为缩放因子来实现：

print 'Zoomed by 2x along the last two axes:'
print ndimage.zoom(data, (1, 2, 2))

这会产生：

Zoomed by 2x along the last two axes:
[[[ 0  0  1  1  2  2]
  [ 1  1  1  2  2  3]
  [ 2  2  3  3  4  4]
  [ 4  4  5  5  6  6]
  [ 5  6  6  7  7  7]
  [ 6  6  7  7  8  8]]

 [[ 9  9 10 10 11 11]
  [10 10 10 11 11 12]
  [11 11 12 12 13 13]
  [13 13 14 14 15 15]
  [14 15 15 16 16 16]
  [15 15 16 16 17 17]]

 [[18 18 19 19 20 20]
  [19 19 19 20 20 21]
  [20 20 21 21 22 22]
  [22 22 23 23 24 24]
  [23 24 24 25 25 25]
  [24 24 25 25 26 26]]]

---

使用map_coordinates 对规则网格数据进行任意插值

首先要了解map_coordinates 是它在pixel 坐标中运行（例如，就像您为数组编制索引一样，但值可以是浮点数）。根据您的描述，这正是您想要的，但如果经常让人们感到困惑。例如，如果您有 x、y、z“真实世界”坐标，则需要将它们转换为基于索引的“像素”坐标。

无论如何，假设我们想要在原始数组的 1.2、0.3、1.4 位置插入值。

如果您根据早期的 RGB 图像情况来考虑这一点，则第一个坐标对应于“带”，第二个坐标对应于“行”，最后一个坐标对应于“列”。什么顺序对应什么完全取决于您决定如何构建数据，但我将使用这些作为“z、y、x”坐标，因为它使与打印数组的比较更容易可视化。

import numpy as np
import scipy.ndimage as ndimage

data = np.arange(27).reshape(3,3,3)

print 'Original:\n', data
print 'Sampled at 1.2, 0.3, 1.4:'
print ndimage.map_coordinates(data, [[1.2], [0.3], [1.4]])

这会产生：

Original:
[[[ 0  1  2]
  [ 3  4  5]
  [ 6  7  8]]

 [[ 9 10 11]
  [12 13 14]
  [15 16 17]]

 [[18 19 20]
  [21 22 23]
  [24 25 26]]]
Sampled at 1.2, 0.3, 1.4:
[14]

再一次，默认情况下这是三次插值。使用order kwarg 来控制插值的类型。

这里值得注意的是，scipy.ndimage 的所有操作都保留了原始数组的 dtype。如果您想要浮点结果，则需要将原始数组转换为浮点数：

In [74]: ndimage.map_coordinates(data.astype(float), [[1.2], [0.3], [1.4]])
Out[74]: array([ 13.5965])

您可能会注意到的另一件事是插值坐标格式对于单个点来说相当麻烦（例如，它需要一个 3xN 数组而不是 Nx3 数组）。然而，当你有坐标序列时，它可以说会更好。例如，考虑沿着穿过数据“立方体”的线进行采样的情况：

xi = np.linspace(0, 2, 10)
yi = 0.8 * xi
zi = 1.2 * xi
print ndimage.map_coordinates(data, [zi, yi, xi])

这会产生：

[ 0  1  4  8 12 17 21 24  0  0]

这也是提及如何处理边界条件的好地方。默认情况下，数组之外的任何值都设置为 0。因此序列中的最后两个值是0。 （即zi 对于最后两个元素 > 2）。

如果我们想要数组外的点，比如-999（我们不能使用nan，因为这是一个整数数组。如果你想要nan，你需要转换为浮点数。 ):

In [75]: ndimage.map_coordinates(data, [zi, yi, xi], cval=-999)
Out[75]: array([   0,    1,    4,    8,   12,   17,   21,   24, -999, -999])

如果我们希望它返回数组外点的最近值，我们会这样做：

In [76]: ndimage.map_coordinates(data, [zi, yi, xi], mode='nearest')
Out[76]: array([ 0,  1,  4,  8, 12, 17, 21, 24, 25, 25])

除了"nearest" 和默认的"constant" 之外，您还可以使用"reflect" 和"wrap" 作为边界模式。这些都是不言自明的，但是如果您感到困惑，请尝试尝试一下。

例如，让我们沿着数组中第一个波段的第一行插入一条线，该线的长度是数组距离的两倍：

xi = np.linspace(0, 5, 10)
yi, zi = np.zeros_like(xi), np.zeros_like(xi)

默认给：

In [77]: ndimage.map_coordinates(data, [zi, yi, xi])
Out[77]: array([0, 0, 1, 2, 0, 0, 0, 0, 0, 0])

比较一下：

In [78]: ndimage.map_coordinates(data, [zi, yi, xi], mode='reflect')
Out[78]: array([0, 0, 1, 2, 2, 1, 2, 1, 0, 0])

In [78]: ndimage.map_coordinates(data, [zi, yi, xi], mode='wrap')
Out[78]: array([0, 0, 1, 2, 0, 1, 1, 2, 0, 1])

希望这能澄清一点！

【讨论】：

彻底回答！感谢关于 map_coordinates 函数坐标的提示，不得不看两次才能弄清楚。

【参考方案2】：

乔的回答很好。根据他的建议，我创建了 regulargrid 包（https://pypi.python.org/pypi/regulargrid/，来源https://github.com/JohannesBuchner/regulargrid）

它通过非常快速的 scipy.ndimage.map_coordinates 为任意坐标比例提供对 n 维笛卡尔网格（根据需要）的支持。