快速内插网格数据 [英] Fast interpolation of grid data

问题描述

我有一个大的3d np.nd数据数组，它代表以规则网格方式在某个体积上采样的物理变量(如array [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.

第一个是

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.

第二个是

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:\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倍，则会沿所有轴进行缩放.

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:\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图像之类的东西，可以通过指定一个元组序列作为缩放因子来实现:

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

---

使用

---

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图像情况，则第一个坐标对应于"band"，第二个坐标对应于"row"，最后一个坐标对应"column".什么顺序对应什么完全取决于您决定如何构造数据的方式，但是我将使用这些顺序作为"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:\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控制插值的类型.

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"和默认的"constant"外，您还可以使用"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])

将此与以下内容进行比较:

Compare this to:

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

希望这可以澄清一些事情！

Hopefully that clarifies things a bit!

这篇关于快速内插网格数据的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！