從`scipy.ndimage.map_coordinates`輸出向量/數組

Question

基本上，是否有可能讓scipy.ndimage.map_coordinates返回一個多值結構，而不僅僅是一個標量？我希望能夠插入一次以在一個點上檢索5個值，而不是必須插入5次。

這是我在MWE上的嘗試，以證明問題。我將從一個標量的3D插值開始。現在我不會在點之間，因爲那不是重點。

import numpy as np 
from scipy import ndimage 
coords = np.array([[1.,1.,1.]]) 
a = np.arange(3*3*3).reshape(3,3,3) 
ndimage.map_coordinates(a,coords.T) # array([13.]) 

現在，假設我要a有對價值的，不只是一個。我的直覺是

a = np.arange(3*3*3*2).reshape(3,3,3,2) 
a[1,1,1] # array([26.,27.]) 
ndimage.map_coordinates(a[:,:,:],coords.T) # I'd like array([26.,27.]) 

而是期望的輸出，我得到如下：

RuntimeError        Traceback (most recent call last) 
(...)/<ipython-input-84-77334fb7469f> in <module>() 
----> 1 ndimage.map_coordinates(a[:,:,:],np.array([[1.,1.,1.]]).T) 

/usr/lib/python2.7/dist-packages/scipy/ndimage/interpolation.pyc in map_coordinates(input, coordinates, output, order, mode, cval, prefilter) 
    287   raise RuntimeError('input and output rank must be > 0') 
    288  if coordinates.shape[0] != input.ndim: 
--> 289   raise RuntimeError('invalid shape for coordinate array') 
    290  mode = _extend_mode_to_code(mode) 
    291  if prefilter and order > 1: 

RuntimeError: invalid shape for coordinate array 

我找不到任何結構（a，coords等的形狀的排列），這給了我正在尋找的答案。另外，如果有比使用map_coordinates更好的方法，請繼續。我認爲scipy.interpolate.interp1d可能是要走的路，但我找不到任何文檔或者它可能會做什麼的暗示...

Answer 1

這是不可能的，我想。

但張量積插值並不難：

import numpy as np 
from scipy.interpolate import interp1d 

def interpn(*args, **kw): 
    """Interpolation on N-D. 

    ai = interpn(x, y, z, ..., a, xi, yi, zi, ...) 
    where the arrays x, y, z, ... define a rectangular grid 
    and a.shape == (len(x), len(y), len(z), ...) 
    """ 
    method = kw.pop('method', 'cubic') 
    if kw: 
     raise ValueError("Unknown arguments: " % kw.keys()) 
    nd = (len(args)-1)//2 
    if len(args) != 2*nd+1: 
     raise ValueError("Wrong number of arguments") 
    q = args[:nd] 
    qi = args[nd+1:] 
    a = args[nd] 
    for j in range(nd): 
     a = interp1d(q[j], a, axis=j, kind=method)(qi[j]) 
    return a 

import matplotlib.pyplot as plt 

x = np.linspace(0, 1, 6) 
y = np.linspace(0, 1, 7) 
k = np.array([0, 1]) 
z = np.cos(2*x[:,None,None] + k[None,None,:]) * np.sin(3*y[None,:,None]) 

xi = np.linspace(0, 1, 60) 
yi = np.linspace(0, 1, 70) 
zi = interpn(x, y, z, xi, yi, method='linear') 

plt.subplot(221) 
plt.imshow(z[:,:,0].T, interpolation='nearest') 

plt.subplot(222) 
plt.imshow(zi[:,:,0].T, interpolation='nearest') 

plt.subplot(223) 
plt.imshow(z[:,:,1].T, interpolation='nearest') 

plt.subplot(224) 
plt.imshow(zi[:,:,1].T, interpolation='nearest') 

plt.show()