使用 scipy 的 RegularGridInterpolator 插值 3D 体积数据

Question

我有一个包含一些体积数据的 3D 数组，即在每个网格点我都有一些代表特定数量大小的值。我想使用 scipy 的

RegularGridInterpolator

插入该数组，但我用错了。

举一个简单的例子，我在这里使用一个我想要插值的圆。原始圆圈看起来相当粗糙，如下所示。

根据我对

RegularGridInterpolator

文档的理解，链接在这里，我认为我需要为插值函数提供新的坐标作为形状

[[x1, y1, z1],[x2, y2, z2],...]

的数组。我的想法是将新的坐标向量提升为列向量，然后将它们与

hstack

组合，但这不起作用，因为从插值器返回的数据形状错误。

这是我的代码：

from mayavi import mlab
import numpy as np
import scipy.interpolate as interp

def make_simple_3Dplot( data2plot, xVals, yVals, zVals, N_contLevels=8, fname_plot='' ):

    contLevels  = np.linspace( np.amin(data2plot), 
                               np.amax(data2plot),
                               N_contLevels)[1:].tolist()
    fig1    = mlab.figure( bgcolor=(1,1,1), fgcolor=(0,0,0),size=(800,600))

    XX, YY, ZZ = np.meshgrid(xVals, yVals, zVals, indexing='ij' )
    contPlot = mlab.contour3d( XX, YY, ZZ, 
                               data2plot, contours=contLevels,
                               transparent=True, opacity=.4,
                               figure=fig1
                             )
    mlab.xlabel('x')
    mlab.ylabel('y')
    mlab.zlabel('z')

    mlab.show()

# define original coordinates
x_min, y_min, z_min = 0, 0, 0
x_max, y_max, z_max = 10, 10, 10
Nx, Ny, Nz = 20, 30, 40
x_arr = np.linspace(x_min, x_max, Nx)
y_arr = np.linspace(y_min, y_max, Ny)
z_arr = np.linspace(z_min, z_max, Nz)

# center of circle
xc, yc, zc = 3, 5, 7

# radius of circle
rc = 2

# define original data
data_3D_original = np.zeros( (Nx, Ny, Nz) )
for ii in range(Nx):
    for jj in range(Ny):
        for kk in range(Nz):
            if np.sqrt((x_arr[ii]-xc)**2 + (y_arr[jj]-yc)**2 + (z_arr[kk]-zc)**2) < rc:
                data_3D_original[ii,jj,kk] = 1.

make_simple_3Dplot( data_3D_original, x_arr, y_arr, z_arr )

# spatial coordinates for interpolation
step_size = np.mean(np.diff(x_arr))/5.
x_interp = np.arange(x_arr[0], x_arr[-1], step_size ) 
y_interp = np.arange(y_arr[0], y_arr[-1], step_size )
z_interp = np.arange(z_arr[0], z_arr[-1], step_size )

# make interpolation function
func_interp = interp.RegularGridInterpolator( (x_arr, y_arr, z_arr), data_3D_original )

# make coordinates for interpolation, first transform vectors for coordinates
# into column vectors and then stack them together 
points = np.hstack( (x_interp[...,None], y_interp[...,None], z_interp[...,None]) )

data_3D_interp = func_interp(points)

print(data_3D_interp.shape, x_interp.shape, y_interp.shape, z_interp.shape)

除了绘图之外，输出显示为

(96,) (96,) (96,) (96,)

，而它应该是

(96,96,96) (96,) (96,) (96,)

。显然，我错过了一些东西。此外，这仅在所有坐标向量具有相同长度的情况下才有效，但在我的实际用例中并非如此。那么我在这里做错了什么？

我正在使用的相关库的版本（虽然我认为这在这里不起作用，但我也认为包含它们是一个很好的做法）：

numpy version: 1.20.3
scipy version: 1.7.2

Answer 1

离开问题一段时间后，我在洗衣服时意识到问题所在。由于我想在由

x_interp

、

y_interp

和

z_interp

跨越的完整 3D 网格上进行插值，我需要使用 numpy.meshgrid() 提供一个

full

坐标列表。

这是我需要在代码中更改的内容

XX, YY, ZZ = np.meshgrid( x_interp, y_interp, z_interp, indexing='ij' )
points = list( zip(XX.ravel(), YY.ravel(), ZZ.ravel() ) )
data_3D_interp = func_interp(points)
data_3D_interp = np.reshape(data_3D_interp, (len(x_interp), len(y_interp), len(z_interp)) )

将

step_size

设置为原始值的 1/10，会产生下图。我对结果有些失望，但我想这就是这里可能发生的事情。使用插值方法可能会改善结果，但这是我无法在我的计算机上测试的东西，因为我遇到内存错误：

numpy.core._exceptions.MemoryError: Unable to allocate ...

Answer 2

这是执行插值的最小片段：

import numpy as np
from scipy import interpolate

def model(x, y, z):
    return np.sqrt(x ** 2 + y ** 2 + z ** 2)

x = np.linspace(-1, 1, 100)
y = np.linspace(-1, 1, 101)
z = np.linspace(-1, 1, 102)

X, Y, Z = np.meshgrid(x, y, z, indexing="ij")

F = model(X, Y, Z)

interpolator = interpolate.RegularGridInterpolator((x, y, z), F)

points = np.array([
    [0, 0, 0],
    [0, 0, 1],
    [1, 0, 0],
    [0, 1, 0],
])

interpolator(points)
# array([0.01414426, 1.00005101, 1.00004901, 1.00010003])

内存问题是由于矩形 3D 网格的立方复杂性造成的。

F.nbytes / 2**20 # 7.8598 Mb 

将分辨率提高 10 倍，将消耗 1000 倍的内存，在这种情况下将消耗大约 8 Gb。