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

Question

在将一些（相当大的物理）Matlab 代码转换为 Python 时，我偶然发现了这种情况。当插值相同的二维离散数据时，Python/Scipy 的

griddata()

函数给出的结果与 Matlab 的对应函数

griddata()

不同。

Matlab 示例代码：

% 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
z1 = griddata(points(:,1),points(:,2),values,xx,yy,"linear");

Python 示例代码：

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')

这是 Matlab 输出。

这是 Python 输出。

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

scipy.spatial

中的 Delaunay 三角测量工具和

scipy.interpolate.LinearNDInterpolator

函数来解决这个问题：

# # 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)

仍然没有运气。结果不一样。这会是什么原因呢？我们怎样才能从 Scipy 的

griddata()

函数得到相同的结果？

Answer 1

您的输入点形成方形网格，这使得 Delaunay 三角剖分不明确（即存在多个同样有效的三角剖分）。正方形的四个顶点之间的线性插值会给出不同的结果，具体取决于您如何将该正方形拆分为两个三角形。这可能是造成差异的最重要原因。

有效观察这种不稳定性的一种方法是向点的位置添加少量噪声几次，并绘制每次获得的 Z 值。你会看到不同的结果。位置噪声实际上不应该对插值产生如此大的影响，但它们会导致三角测量发生变化（当向点的位置添加一点点噪声时，每个正方形可能会以不同的方式分成两个三角形），并且因此您会得到显着不同的插值结果。

我将

0.001 * rand()

添加到 MATLAB 中的每个坐标，并获得以下四个图：

MATLAB 和 Python 都能产生正确的插值。但由于您的输入位于网格上，因此最好使用网格插值方法。它会快得多，并且明确（这意味着您可以用两种语言获得基本相同的结果）。