网格数据的快速插值

Question

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

我想通过在粗网格中插入数据来获得更精细的网格间距.目前我正在使用scipy griddata线性插值,但它很慢(对于20x20x20阵列,约为90秒).为了我的目的,它有点过度设计,允许随机采样体积数据.有什么东西可以利用我的常规间隔数据,以及我想插入的只有一组有限的特定点吗？

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

使用MATLAB定期网格化数据的任意插值 map_coordinates

首先要强调的map_coordinates是它以像素坐标运行(例如,就像你索引数组一样,但值可以是浮点数).根据您的描述,这正是您想要的,但如果经常让人困惑.例如,如果您有x,y,z"真实世界"坐标,则需要将它们转换为基于索引的"像素"坐标.

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

如果您根据早期的RGB图像情况考虑这一点,则第一个坐标对应于"band",第二个坐标对应于"row",最后一个坐标对应于"column".什么顺序对应于完全取决于您决定如何构建数据的顺序,但我将使用这些作为"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]

再一次,这是默认的三次插值.使用orderkwarg来控制插值的类型.

值得注意的是,所有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])

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

例如,让我们沿着数组中第一个波段的第一行插入一条线,该线延伸两倍于数组的距离:

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

希望这能澄清一点!

Answer 2

Joe给出了很好的答案.基于他的建议,我创建了regulargrid包(https://pypi.python.org/pypi/regulargrid/,来源https://github.com/JohannesBuchner/regulargrid)

它通过非常快速的scipy.ndimage.map_coordinates为任意坐标刻度提供对n维笛卡尔网格的支持(根据需要).