使用Scipy RegularGridInterpolator插值3D体数据遇形状错误问题

我有一组3D体数据数组，每个网格点对应某物理量的幅值。尝试使用Scipy的RegularGridInterpolator插值时，返回的插值数据形状不符合预期（输出为(96,)，期望为(96,96,96)）；且原方式仅支持等长坐标向量，无法适配实际场景中不等长的情况。

复现代码如下：

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,)。

错误原因 #

你构造points的方式有误：np.hstack((x_interp[...,None], y_interp[...,None], z_interp[...,None]))生成的是96个点的序列（每个点是(x_i, y_i, z_i)，对应三个坐标向量的第i个元素），而非整个3D网格的所有笛卡尔积点。因此插值函数只会返回这96个点的结果，而非完整的3D体数据。

修正步骤 #

要得到完整的3D插值结果，需要：

1. 生成插值坐标的笛卡尔积网格（即所有x、y、z组合）
2. 将网格转换为RegularGridInterpolator要求的(N, 3)形状输入
3. 将插值结果重塑回目标3D形状

修正后的代码片段 #

替换原代码中points构造和插值的部分：

# 生成插值坐标的笛卡尔积网格，保持和原数据一致的indexing='ij'
XX_interp, YY_interp, ZZ_interp = np.meshgrid(x_interp, y_interp, z_interp, indexing='ij')

# 将网格转换为(N, 3)的点数组，N是总点数
points = np.stack([XX_interp.ravel(), YY_interp.ravel(), ZZ_interp.ravel()], axis=1)

# 插值后重塑为目标形状
data_3D_interp = func_interp(points).reshape(XX_interp.shape)

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

适配不等长坐标向量 #

此方法天然支持不等长的x、y、z插值坐标。例如，若x_interp长度为96，y_interp为146，z_interp为196，上述代码会自动生成(96,146,196)形状的插值结果，无需额外修改。

内容的提问来源于stack exchange，提问作者Alf