用numpy和scipy插入三維體積

Question

我非常沮喪，因爲幾個小時後我似乎無法在python中做一個看似簡單的3D插值。在Matlab中所有我需要做的就是

Vi = interp3(x,y,z,V,xi,yi,zi) 

這是什麼使用SciPy的的ndimage.map_coordinate或其他numpy的方法完全等效？

由於

Answer 1

在SciPy的0.14或更高版本，有一個新的功能scipy.interpolate.RegularGridInterpolator這與interp3非常相似。

MATLAB命令Vi = interp3(x,y,z,V,xi,yi,zi)將轉化爲類似：

from numpy import array 
from scipy.interpolate import RegularGridInterpolator as rgi 
my_interpolating_function = rgi((x,y,z), V) 
Vi = my_interpolating_function(array([xi,yi,zi]).T) 

下面是一個完整的例子都證明;它會幫助你瞭解確切的差異...

MATLAB代碼：

x = linspace(1,4,11); 
y = linspace(4,7,22); 
z = linspace(7,9,33); 
V = zeros(22,11,33); 
for i=1:11 
    for j=1:22 
     for k=1:33 
      V(j,i,k) = 100*x(i) + 10*y(j) + z(k); 
     end 
    end 
end 
xq = [2,3]; 
yq = [6,5]; 
zq = [8,7]; 
Vi = interp3(x,y,z,V,xq,yq,zq); 

結果是Vi=[268 357]這的確是在這兩點​​3210和(3,5,7)值。

SciPy的CODE：

from scipy.interpolate import RegularGridInterpolator 
from numpy import linspace, zeros, array 
x = linspace(1,4,11) 
y = linspace(4,7,22) 
z = linspace(7,9,33) 
V = zeros((11,22,33)) 
for i in range(11): 
    for j in range(22): 
     for k in range(33): 
      V[i,j,k] = 100*x[i] + 10*y[j] + z[k] 
fn = RegularGridInterpolator((x,y,z), V) 
pts = array([[2,6,8],[3,5,7]]) 
print(fn(pts)) 

同樣是[268,357]。所以你會看到一些細微差別：Scipy使用x，y，z索引順序，而MATLAB使用y，x，z（奇怪）;在Scipy中，你在一個單獨的步驟中定義了一個函數，當你調用它時，座標被分組爲（x1，y1，z1），（x2，y2，z2），...而matlab使用（x1，x2，...） 。），（Y1，Y2，...），（Z1，Z2，...）。

除此之外，兩者相似並且同樣易於使用。

Answer 2

基本上，ndimage.map_coordinates作品在 「索引」 的座標（也稱爲 「體素」 或 「像素」 座標）。它的界面起初似乎有點笨拙，但它確實給你提供了很多靈活性。

如果要指定類似於matlab的interp3的插補座標，那麼您需要將輸入座標轉換爲「索引」座標。

還有一個額外的皺紋map_coordinates始終保留在輸出中的輸入數組的dtype。如果你插入一個整數數組，你會得到整數輸出，這可能是也可能不是你想要的。對於下面的代碼片段，我假設你總是需要浮點輸出。 （如果你不這樣做，這實際上更簡單。）

我會嘗試在今晚晚些時候添加更多解釋（這是相當密集的代碼）。

總之，我所擁有的interp3函數比它可能需要用於確切目的更復雜。 Howver，我記得它應該或多或少複製的interp3行爲（忽略interp3(data, zoom_factor)的「縮放」功能，該功能scipy.ndimage.zoom把手。）

import numpy as np 
from scipy.ndimage import map_coordinates 

def main(): 
    data = np.arange(5*4*3).reshape(5,4,3) 

    x = np.linspace(5, 10, data.shape[0]) 
    y = np.linspace(10, 20, data.shape[1]) 
    z = np.linspace(-100, 0, data.shape[2]) 

    # Interpolate at a single point 
    print interp3(x, y, z, data, 7.5, 13.2, -27) 

    # Interpolate a region of the x-y plane at z=-25 
    xi, yi = np.mgrid[6:8:10j, 13:18:10j] 
    print interp3(x, y, z, data, xi, yi, -25 * np.ones_like(xi)) 

def interp3(x, y, z, v, xi, yi, zi, **kwargs): 
    """Sample a 3D array "v" with pixel corner locations at "x","y","z" at the 
    points in "xi", "yi", "zi" using linear interpolation. Additional kwargs 
    are passed on to ``scipy.ndimage.map_coordinates``.""" 
    def index_coords(corner_locs, interp_locs): 
     index = np.arange(len(corner_locs)) 
     if np.all(np.diff(corner_locs) < 0): 
      corner_locs, index = corner_locs[::-1], index[::-1] 
     return np.interp(interp_locs, corner_locs, index) 

    orig_shape = np.asarray(xi).shape 
    xi, yi, zi = np.atleast_1d(xi, yi, zi) 
    for arr in [xi, yi, zi]: 
     arr.shape = -1 

    output = np.empty(xi.shape, dtype=float) 
    coords = [index_coords(*item) for item in zip([x, y, z], [xi, yi, zi])] 

    map_coordinates(v, coords, order=1, output=output, **kwargs) 

    return output.reshape(orig_shape) 

main() 

Answer 3

的確切相當於MATLAB的interp3將使用SciPy的的interpn用於一次性插補：

import numpy as np 
from scipy.interpolate import interpn 

Vi = interpn((x,y,z), V, np.array([xi,yi,zi]).T) 

兩個MATLAB和SciPy的默認方法是線性內插，並且這可以與method來改變論據。請注意，對於3維和更高維度，只有線性和最近鄰插值支持interpn，與支持三次和樣條插值的MATLAB不同。

當在同一個網格上進行多次插值調用時，最好使用插值對象RegularGridInterpolator，如在接受的答案above中那樣。 interpn內部使用RegularGridInterpolator。