网格数据的快速插值 [英] Fast interpolation of grid data
问题描述
我有一个大的 3d np.ndarray 数据，它表示以常规网格方式在体积上采样的物理变量(如数组 [0,0,0] 中的值表示物理坐标 (0,0,0)).
I have a large 3d np.ndarray of data that represents a physical variable sampled over a volume in a regular grid fashion (as in the value in array[0,0,0] represents the value at physical coords (0,0,0)).
我想通过在粗网格中插入数据来获得更精细的网格间距.目前我正在使用 scipy griddata 线性插值，但它非常慢(20x20x20 阵列约为 90 秒).就我的目的而言，它有点过度设计，允许对体积数据进行随机采样.有什么东西可以利用我的规则间隔数据以及我只想插入的一组有限的特定点的事实?
I would like to go to a finer grid spacing by interpolating the data in the rough grid. At the moment I'm using scipy griddata linear interpolation but it's pretty slow (~90secs for 20x20x20 array). It's a bit overengineered for my purposes, allowing random sampling of the volume data. Is there anything out there that can take advantage of my regularly spaced data and the fact that there is only a limited set of specific points I want to interpolate to?
推荐答案
好的！有两个选项可以做不同的事情，但都利用了原始数据的规则网格性质.
Sure! There are two options that do different things but both exploit the regularly-gridded nature of the original data.
第一个是scipy.ndimage.zoom.如果您只是想基于对原始数据进行插值来生成更密集的规则网格，这就是您要走的路.
The first is scipy.ndimage.zoom. If you just want to produce a denser regular grid based on interpolating the original data, this is the way to go.
第二个是scipy.ndimage.map_coordinates.如果您想在数据中插入一些(或许多)任意点，但仍要利用原始数据的规则网格性质(例如，不需要四叉树)，这是要走的路.
The second is scipy.ndimage.map_coordinates. If you'd like to interpolate a few (or many) arbitrary points in your data, but still exploit the regularly-gridded nature of the original data (e.g. no quadtree required), it's the way to go.
举一个简单的例子(这将使用三次插值.使用 order=1 表示双线性，order=0 表示最近等):
As a quick example (This will use cubic interpolation. Use order=1 for bilinear, order=0 for nearest, etc.):
import numpy as np
import scipy.ndimage as ndimage
data = np.arange(9).reshape(3,3)
print 'Original:
', data
print 'Zoomed by 2x:
', 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 倍，则会沿所有轴进行缩放.
This also works for 3D (and nD) arrays. However, be aware that if you zoom by 2x, for example, you'll zoom along all axes.
data = np.arange(27).reshape(3,3,3)
print 'Original:
', 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 图像之类的东西，您可以通过将元组序列指定为缩放因子来实现:
If you have something like a 3-band, RGB image that you'd like to zoom, you can do this by specifying a sequence of tuples as the zoom factor:
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真实世界"坐标，则需要将它们转换为基于索引的像素"坐标.
Arbitrary interpolation of regularly-gridded data using map_coordinates
The first thing to undersand about map_coordinates is that it operates in pixel coordinates (e.g. just like you'd index the array, but the values can be floats). From your description, this is exactly what you want, but if often confuses people. For example, if you have x, y, z "real-world" coordinates, you'll need to transform them to index-based "pixel" coordinates.
无论如何，假设我们想在原始数组中的位置 1.2、0.3、1.4 处插入值.
At any rate, let's say we wanted to interpolate the value in the original array at position 1.2, 0.3, 1.4.
如果您根据较早的 RGB 图像情况考虑这一点，则第一个坐标对应于波段"，第二个坐标对应于行"，最后一个对应于列".什么顺序对应什么完全取决于您决定如何构建数据，但我将使用这些作为z、y、x"坐标，因为它使与打印数组的比较更容易可视化.
If you're thinking of this in terms of the earlier RGB image case, the first coordinate corresponds to the "band", the second to the "row" and the last to the "column". What order corresponds to what depends entirely on how you decide to structure your data, but I'm going to use these as "z, y, x" coordinates, as it makes the comparison to the printed array easier to visualize.
import numpy as np
import scipy.ndimage as ndimage
data = np.arange(27).reshape(3,3,3)
print 'Original:
', 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 来控制插值的类型.
Once again, this is cubic interpolation by default. Use the order kwarg to control the type of interpolation.
这里值得注意的是，scipy.ndimage 的所有操作都保留了原始数组的 dtype.如果您想要浮点结果，则需要将原始数组转换为浮点数:
It's worth noting here that all of scipy.ndimage's operations preserve the dtype of the original array. If you want floating point results, you'll need to cast the original array as a float:
In [74]: ndimage.map_coordinates(data.astype(float), [[1.2], [0.3], [1.4]])
Out[74]: array([ 13.5965])
您可能会注意到的另一件事是内插坐标格式对于单个点来说相当麻烦(例如，它需要 3xN 数组而不是 Nx3 数组).然而，当你有坐标序列时，它可以说更好.例如，考虑沿穿过数据立方体"的一条线采样的情况:
Another thing you may notice is that the interpolated coordinates format is rather cumbersome for a single point (e.g. it expects a 3xN array instead of an Nx3 array). However, it's arguably nicer when you have sequences of coordinate. For example, consider the case of sampling along a line that passes through the "cube" of data:
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).
This is also a good place to mention how boundary conditions are handled. By default, anything outside of the array is set to 0. Thus the last two values in the sequence are 0. (i.e. zi is > 2 for the last two elements).
如果我们想要数组外的点，说 -999(我们不能使用 nan 因为这是一个整数数组.如果你想要 nan，你需要转换为浮点数.):
If we wanted the points outside the array to be, say -999 (We can't use nan as this is an integer array. If you want nan, you'll need to cast to floats.):
In [75]: ndimage.map_coordinates(data, [zi, yi, xi], cval=-999)
Out[75]: array([ 0, 1, 4, 8, 12, 17, 21, 24, -999, -999])
如果我们希望它返回数组外点的最近值，我们会这样做:
If we wanted it to return the nearest value for points outside the array, we'd do:
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" 和默认的 之外，您还可以使用 .这些都是不言自明的，但如果您感到困惑，请尝试尝试一下."reflect" 和 "wrap" 作为边界模式>常数"
You can also use "reflect" and "wrap" as boundary modes, in addition to "nearest" and the default "constant". These are fairly self-explanatory, but try experimenting a bit if you're confused.
例如，让我们沿阵列中第一个波段的第一行插入一条线，该线延伸为阵列距离的两倍:
For example, let's interpolate a line along the first row of the first band in the array that extends for twice the distance of the array:
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])
希望能澄清一点！
这篇关于网格数据的快速插值的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
