Python griddata() 和 Matlab griddata()：某些网格点的结果不同

在将一些（相当大的物理）Matlab 代码转换为 Python 时，我偶然发现了这种情况。当对相同的二维离散数据进行插值时，Python/Scipy 的 griddata() 函数给出的结果与 Matlab 的对应函数不同。 griddata() Matlab 示例代码：

Python 示例代码：

% Sample points (x,y): 7x5=35 points
points = [-3., -2.; -2.25, -2.; -1.5, -2.; -0.75, -2.; 0., -2.; 0.75, -2.; 1.5, -2.;
    -3., -1.25; -2.25, -1.25; -1.5, -1.25; -0.75, -1.25; 0., -1.25; 0.75, -1.25; 1.5, -1.25;
    -3., -0.5; -2.25, -0.5; -1.5, -0.5; -0.75, -0.5; 0., -0.5; 0.75, -0.5; 1.5, -0.5;
    -3., 0.25; -2.25, 0.25; -1.5, 0.25; -0.75, 0.25; 0., 0.25; 0.75, 0.25; 1.5, 0.25;
    -3., 1.; -2.25, 1.; -1.5, 1.; -0.75, 1.; 0., 1.; 0.75, 1.; 1.5, 1.]

% Data values at the sample points: 35 values
values = [0.12702104; -0.24534877; -0.08879096; 0.13749533; 0.2431933; -0.01159269; 0.11705169;
    0.10039714; 0.23370454; 0.0020298; -0.05512349; -0.11208613; -0.09864074; 0.45360373;
    -0.10414895; -0.18303173; 0.39306646; 0.05457107; -0.19024738; 0.06828192; 0.29422425;
    -0.18672513; -0.23610061; 0.48881179; -0.09731334; -0.05470424; 0.26214565; -0.17978073,
    -0.10337649; 0.10246465; 0.24676284; -0.05835531; 0.06804471; -0.34589734; -0.34035067];

% Grid points to interpolate at: 12x12 points
[X,Y] = meshgrid(linspace(-2,1.2,12), linspace(-1.75,0.75,12));

% Interpolation
Z = griddata(points(:,1), points(:,2), values, X, Y, "linear");

这是 Matlab 输出。

import numpy as np
import scipy

# Sample points (x, y): 7x5 = 35 points
points = np.array([
    [-3., -2.], [-2.25, -2.], [-1.5, -2.], [-0.75, -2.], [0., -2.], [0.75, -2.], [1.5, -2.],
    [-3., -1.25], [-2.25, -1.25], [-1.5, -1.25], [-0.75, -1.25], [0., -1.25], [0.75, -1.25], [1.5, -1.25],
    [-3., -0.5], [-2.25, -0.5], [-1.5, -0.5], [-0.75, -0.5], [0., -0.5], [0.75, -0.5], [1.5, -0.5],
    [-3., 0.25], [-2.25, 0.25], [-1.5, 0.25], [-0.75, 0.25], [0., 0.25], [0.75, 0.25], [1.5, 0.25],
    [-3., 1.], [-2.25, 1.], [-1.5, 1.], [-0.75, 1.], [0., 1.], [0.75, 1.], [1.5, 1.]
])

# Data values at the sample points: 35 values
values = np.array([
    0.12702104, -0.24534877, -0.08879096, 0.13749533, 0.2431933, -0.01159269, 0.11705169,
    0.10039714, 0.23370454, 0.0020298, -0.05512349, -0.11208613, -0.09864074, 0.45360373,
    -0.10414895, -0.18303173, 0.39306646, 0.05457107, -0.19024738, 0.06828192, 0.29422425,
    -0.18672513, -0.23610061, 0.48881179, -0.09731334, -0.05470424, 0.26214565, -0.17978073,
    -0.10337649, 0.10246465, 0.24676284, -0.05835531, 0.06804471, -0.34589734, -0.34035067
])

# Interpolation grid: 12x12 points 
x = np.linspace(-2, 1.2, 12)
y = np.linspace(-1.75, 0.75, 12)
X, Y = np.meshgrid(x, y)

# Interpolation
Z = scipy.interpolate.griddata(points, values, (X, Y), method='linear')

并且这是Python 输出。

输出值差异很大。然后我尝试使用

中的 Delaunay 三角测量工具和 scipy.spatial 功能来绕过它： scipy.interpolate.LinearNDInterpolator 仍然没有运气。结果不一样。这会是什么原因呢？我们如何从 Scipy 的

# # Perform Delaunay triangulation
tri = Delaunay(points, qhull_options='Qbb Qc') #Qt is always enabled
# # Interpolation function
interp = sp.LinearNDInterpolator(points, values)
ZZ = interp(X, Y)

函数中得到相同的结果？ griddata() function?

Matlab 和 Python 的 griddata 中的三角剖分方案差异

Matlab 和 Python 的 griddata 函数使用不同的默认三角剖分方案进行线性插值，这会导致在某些网格点处出现不同的结果。

Matlab

• Matlab 的 griddata 默认使用 Delaunay 三角剖分 。

Python

• Python 的 scipy.interpolate.griddata 默认使用 Qhull 库 进行三角剖分。虽然 Qhull 也能执行 Delaunay 三角剖分，但它使用的算法和参数与 Matlab 不同。

解决方法

要使 Python 的 griddata 产生与 Matlab 相似的结果，可以指定使用 scipy.spatial.Delaunay 函数进行三角剖分，并将其传递给 scipy.interpolate.LinearNDInterpolator ：

import numpy as np
from scipy.spatial import Delaunay
from scipy.interpolate import LinearNDInterpolator

# ...（的现有代码）...

# 执行 Delaunay 三角剖分
tri = Delaunay(points)  # 使用与 Matlab 相同的三角剖分方案

# 插值函数
interp = LinearNDInterpolator(tri, values)
ZZ = interp(X, Y)

其他注意事项

• 即使使用相同的三角剖分方案，由于浮点运算的细微差异，Matlab 和 Python 的结果也可能略有不同。
• 对于大型数据集， LinearNDInterpolator 可能比 griddata 更快。
• 还可以探索 scipy.interpolate 模块中的其他插值方法，例如 CloughTocher2DInterpolator （对于平滑插值）或 Rbf （对于径向基函数插值）。

通过使用 scipy.spatial.Delaunay 显式指定三角剖分方案，可以使 Python 代码的行为更接近 Matlab 的 griddata 函数，并获得更一致的插值结果。