如何利用ndimage.map_coordinates建立插补问题？

Question

根据scipy.ndimage.map_coordinates的文档，

坐标数组用于为输出中的每个点查找输入中的对应坐标。输入在这些坐标下的值是通过要求的顺序的样条插值来确定的。

输出的形状是由坐标阵列的形状通过放下第一轴得到的。数组沿第一轴的值是输入数组中的坐标，在输入数组中找到输出值。

我有一个在三维网格(t，x，y)上定义的离散的三维函数；在这个三维网格的每一点上，这个函数都有一个唯一的值，除非它的值是零。

我有另一组数组，格式是熊猫数据，它有三列: t_new、x_new和y_new。

我想使用scipy.ndimage.map_coordinates插值函数，以便在上述数据集中显示的新数据集上计算其值。

由于我收到了以下错误消息，所以我确信我没有正确设置map_coordinates：

File "D:\Users\username\Anaconda3\lib\site-packages\scipy\ndimage\interpolation.py", line 437, in map_coordinates

raise RuntimeError('invalid shape for coordinate array')

下面是我对插值函数的定义：

from scipy.ndimage import map_coordinates



def interpolator_3d(df, func_values):

    # The coordinates at which input is evaluated

    coordinates = df[['t_new', 'x_new', 'y_new']].values.T    # (3, 1273)

    # list of input array [[t0, x0, y0, value0], [t1, x1, y1, value1], ...]

    input_arr = func_values                                 # (1780020000, 4)

    return map_coordinates(input_arr, coordinates)

Answer 1

关于如何使用map_coordinates，至少有两个问题。请记住，此功能是为图像重采样而设计的。

1. ，如果你有一个3D函数，数组input_arr应该是三维的.map_coordinates将使用索引作为t、x和y坐标。函数的值v必须存储在每个相应的位置。如果您的原始函数有另一个基本网格，那么您必须将所有内容都相应地规范化到数组索引之前和之后。这需要一个等距的网格作为输入，

，

1. ，coordinates，coordinates必须是一个数组，例如，表单[[t_new_0, t_new_1, ...], [x_new_0, x_new_1 ...], [y_new_0, y_new_1, ...]]。结果将是一个插值样本的列表[[v_new_0, v_new_1, ...]]。通常，如果input_array是n-dimensional，则coordinates必须是包含相同形状S的n数组的列表。结果将是形状S.

数组的列表。

以n=3维和5个样本为例，将其插值为一维形状：

import numpy as np
from scipy import ndimage

a = np.arange(64.).reshape((4, 4, 4))
print(a)

out = ndimage.map_coordinates(a, [
  [0.5, 1.0, 1.5, 2.0, 2.5], [0.1, 0.2, 0.3, 0.4, 0.5], [2.0, 1.9, 1.8, 1.7, 1.6]
])
print(out)

输出：

[[[ 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. 27.]
  [28. 29. 30. 31.]]
 [[32. 33. 34. 35.]
  [36. 37. 38. 39.]
  [40. 41. 42. 43.]
  [44. 45. 46. 47.]]
 [[48. 49. 50. 51.]
  [52. 53. 54. 55.]
  [56. 57. 58. 59.]
  [60. 61. 62. 63.]]]

[ 7.6688, 18.148 , 26.3424, 34.6304, 45.3904]

更新：

这意味着，如果您的input_array具有具有长度1780020000 = 19778 * 500 * 180的表单[[t0, x0, y0, value0], [t1, x1, y1, value1], ...]，则必须相应地将其转换为形状(19778, 500, 180)数组。

t_max, x_max, y_max, _ = np.max(func_values, axis=0).astype(int) + 1  # 19778, 500, 180
input_arr = np.zeros((t_max, x_max, y_max), dtype=float)

for t, x, y, v in func_values:
  input_arr[int(t), int(x), int(y)] = v