Skip to main content
 首页 » 编程设计

python之为什么 Scipy 的 ndimage.map_coordinates 对某些数组不返回值或返回错误结果

2024年11月01日30webabcd

代码返回正确值但并不总是返回值

在下面的代码中,python 为 arr_b 返回了正确的插值,但为 arr_a 返回了正确的插值。

不过,我已经研究了这个问题大约一天了,我真的不确定发生了什么。

出于某种原因,对于 arr_a,twoD_interpolate 一直返回 [0],即使我玩弄或弄乱数据和输入也是如此。

如何修复我的代码,使其实际上对 arr_a 进行插值并返回正确的结果?

import numpy as np 
from scipy.ndimage import map_coordinates 
 
 
def twoD_interpolate(arr, xmin, xmax, ymin, ymax, x1, y1): 
    """ 
    interpolate in two dimensions with "hard edges" 
    """ 
    ny, nx = arr.shape  # Note the order of ny and xy 
 
    x1 = np.atleast_1d(x1) 
    y1 = np.atleast_1d(y1) 
 
    # Mask upper and lower boundaries using @Jamies suggestion 
    np.clip(x1, xmin, xmax, out=x1) 
    np.clip(y1, ymin, ymax, out=y1) 
 
    # Change coordinates to match your array. 
    x1 = (x1 - xmin) * (xmax - xmin) / float(nx - 1) 
    y1 = (y1 - ymin) * (ymax - ymin) / float(ny - 1) 
 
    # order=1 is required to return your examples. 
    return map_coordinates(arr, np.vstack((y1, x1)), order=1) 
 
 
# test data 
arr_a = np.array([[0.7, 1.7, 2.5, 2.8, 2.9], 
                  [1.9, 2.9, 3.7, 4.0, 4.2], 
                  [1.4, 2.0, 2.5, 2.7, 3.9], 
                  [1.1, 1.3, 1.6, 1.9, 2.0], 
                  [0.6, 0.9, 1.1, 1.3, 1.4], 
                  [0.6, 0.7, 0.9, 1.1, 1.2], 
                  [0.5, 0.7, 0.9, 0.9, 1.1], 
                  [0.5, 0.6, 0.7, 0.7, 0.9], 
                  [0.5, 0.6, 0.6, 0.6, 0.7]]) 
 
 
arr_b = np.array([[6.4, 5.60, 4.8, 4.15, 3.5, 2.85, 2.2], 
                  [5.3, 4.50, 3.7, 3.05, 2.4, 1.75, 1.1], 
                  [4.7, 3.85, 3.0, 2.35, 1.7, 1.05, 0.4], 
                  [4.2, 3.40, 2.6, 1.95, 1.3, 0.65, 0.0]]) 
 
 
 
# Test the second array 
print twoD_interpolate(arr_b, 0, 6, 9, 12, 4, 11) 
 
 
# Test first area 
print twoD_interpolate( 
    arr_a, 0, 500, 0, 2000, 0, 2000) 
 
print arr_a[0] 
 
print twoD_interpolate( 
    arr_a_60, 0, 500, 0, 2000, 0, 2000)[0] 
print twoD_interpolate( 
    arr_a, 20, 100, 100, 1600, 902, 50) 
print twoD_interpolate( 
    arr_a, 100, 1600, 20, 100, 902, 50) 
print twoD_interpolate( 
    arr_a, 100, 1600, 20, 100, 50, 902) 
 
 
## Output 
[ 1.7] 
[ 0.] 
[ 0.7  1.7  2.5  2.8  2.9] 
0.0 
[ 0.] 
[ 0.] 
[ 0.] 

返回错误值的代码:

arr = np.array([[12.8, 20.0, 23.8, 26.2, 27.4, 28.6], 
                [10.0, 13.6, 15.8, 17.4, 18.2, 18.8], 
                [5.5, 7.7, 8.7, 9.5, 10.1, 10.3], 
                [3.3, 4.7, 5.1, 5.5, 5.7, 6.1]]) 
 
twoD_interpolate(arr, 0, 1, 1400, 3200, 0.5, 1684) 
# above should return 21 but is returning 3.44 

请您参考如下方法:

这实际上是我在原始问题中的错。

如果我们检查它试图插入的位置 twoD_interpolate(arr, 0, 1, 1400, 3200, 0.5, 1684) 我们得到 arr[ 170400, 0.1] 作为要查找的值,mode='nearest' 会将其裁剪为 arr[ -1 , 0.1]。请注意,我切换了 x 和 y 以获得它在数组中出现的位置。

这对应于值 arr[-1,0] = 3.3 和 arr[-1,1] = 4.7 的插值,因此插值看起来像 3.3 * .9 + 4.7 * .1 = 3.44。

问题是从容应对的。如果我们采用一个从 50 到 250 的数组:

>>> a=np.arange(50,300,50) 
>>> a 
array([ 50, 100, 150, 200, 250]) 
>>> stride=float(a.max()-a.min())/(a.shape[0]-1) 
>>> stride 
50.0 
 
>>> (75-a.min()) * stride 
1250.0   #Not what we want! 
>>> (75-a.min()) / stride 
0.5      #There we go 
>>> (175-a.min()) / stride 
2.5      #Looks good 

我们可以使用 map_coordinates 查看:

#Input array from the above. 
print map_coordinates(arr, np.array([[.5,2.5,1250]]), order=1, mode='nearest') 
[ 75 175 250] #First two are correct, last is incorrect. 

所以我们真正需要的是 (x-xmin)/stride,对于前面的示例,stride 是 1,所以这无关紧要。

代码应该是这样的:

def twoD_interpolate(arr, xmin, xmax, ymin, ymax, x1, y1): 
    """ 
    interpolate in two dimensions with "hard edges" 
    """ 
    arr = np.atleast_2d(arr) 
    ny, nx = arr.shape  # Note the order of ny and xy 
 
    x1 = np.atleast_1d(x1) 
    y1 = np.atleast_1d(y1) 
 
    # Change coordinates to match your array. 
    if nx==1: 
        x1 = np.zeros_like(x1.shape) 
    else: 
        x_stride = (xmax-xmin)/float(nx-1) 
        x1 = (x1 - xmin) / x_stride 
 
    if ny==1: 
        y1 = np.zeros_like(y1.shape) 
    else: 
        y_stride = (ymax-ymin)/float(ny-1) 
        y1 = (y1 - ymin) / y_stride 
 
    # order=1 is required to return your examples and mode=nearest prevents the need of clip. 
    return map_coordinates(arr, np.vstack((y1, x1)), order=1, mode='nearest') 

请注意,mode='nearest' 不需要 clip。

print twoD_interpolate(arr, 0, 1, 1400, 3200, 0.5, 1684) 
[ 21.024] 
 
print twoD_interpolate(arr, 0, 1, 1400, 3200, 0, 50000) 
[ 3.3] 
 
print twoD_interpolate(arr, 0, 1, 1400, 3200, .5, 50000) 
[ 5.3] 

检查一维或伪一维数组。将仅插入 x 维度,除非输入数组的形状正确:

arr = np.arange(50,300,50) 
print twoD_interpolate(arr, 50, 250, 0, 5, 75, 0) 
[75] 
 
arr = np.arange(50,300,50)[None,:] 
print twoD_interpolate(arr, 50, 250, 0, 5, 75, 0) 
[75] 
 
arr = np.arange(50,300,50) 
print twoD_interpolate(arr, 0, 5, 50, 250, 0, 75) 
[50] #Still interpolates the `x` dimension. 
 
arr = np.arange(50,300,50)[:,None] 
print twoD_interpolate(arr, 0, 5, 50, 250, 0, 75) 
[75]