您好, 欢迎来到 !    登录 | 注册 | | 设为首页 | 收藏本站

快速内插网格数据

快速内插网格数据

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

首先是scipy.ndimage.zoom。如果您只想基于对原始数据进行插值来生成更密集的规则网格,则可以采用这种方法

第二个是scipy.ndimage.map_coordinates。如果您想在数据中插入几个(或多个)任意点,但仍然利用原始数据的规则网格性质(例如,不需要四叉树),则可以采用此方法

“缩放”数组(scipy.ndimage.zoom

作为一个简单的示例(这将使用三次插值。order=1@H_301_13@用于双线性,order=0@H_301_13@最接近等):

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)
@H_301_13@

这样产生:

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

这也适用于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
@H_301_13@

这样产生:

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)
@H_301_13@

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

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

这样产生:

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

map_coordinates

首先要了解的map_coordinates@H_301_13@是它在 像素 坐标下运行(例如,就像您要为数组建立索引一样,但是值可以是浮点数)。根据您的描述,这正是您想要的,但是如果经常使人感到困惑。例如,如果您具有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]])
@H_301_13@

这样产生:

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

再一次,这是认的三次插值。使用order@H_301_13@kwarg控制插值的类型。

在这里值得注意的是,所有scipy.ndimage@H_301_13@的操作都保留了原始数组的dtype。如果要浮点结果,则需要将原始数组强制转换为浮点数:

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

您可能会注意到的另一件事是,插值坐标格式对于单点而言非常麻烦(例如,它期望使用3xN数组而不是Nx3数组)。但是,当您具有坐标序列时,可以说是更好的选择。例如,考虑沿穿过数据“多维数据集”的线进行采样的情况:

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

这样产生:

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

这也是提及如何处理边界条件的好地方。认情况下,数组外部的任何内容都设置为0。因此,序列中的最后两个值是0@H_301_13@。(即zi@H_301_13@最后两个元素> 2)。

如果我们希望将数组外的点设为-999@H_301_13@(我们不能使用,nan@H_301_13@因为这是一个整数数组。如果想要nan@H_301_13@,则需要转换为浮点数。):

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

如果我们希望它为数组外的点返回最接近的值,我们可以这样做:

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

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

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

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

认给出:

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

比较一下:

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

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

其他 2022/1/1 18:27:48 有596人围观

撰写回答


你尚未登录,登录后可以

和开发者交流问题的细节

关注并接收问题和回答的更新提醒

参与内容的编辑和改进,让解决方法与时俱进

请先登录

推荐问题


联系我
置顶