Python与Matlab的griddata()函数插值结果不一致问题排查

问题背景 #

将大型物理领域Matlab代码迁移至Python时，发现对同一组二维散乱数据执行线性插值时，Matlab的griddata函数与Scipy同名函数输出结果存在显著差异。

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
Z = griddata(points(:,1), points(:,2), values, X, Y, "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的结果中，凸包外的区域直接返回NaN，且凸包内的插值数值与Matlab存在细微偏差。

尝试使用Scipy的Delaunay三角化结合LinearNDInterpolator手动实现插值，结果仍与Matlab不一致：

from scipy.spatial import Delaunay
import scipy.interpolate as sp

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

差异原因分析 #

1. Delaunay三角化实现与参数差异
   Matlab的griddata使用自研的Delaunay三角化算法，默认处理共线点、退化三角的逻辑与Scipy依赖的Qhull库不同。Scipy默认的Qhull参数更严格，会拒绝某些退化情况，导致三角剖分的拓扑结构与Matlab不一致，而线性插值的结果完全依赖三角剖分的网格。

2. 外插处理逻辑不同
   Matlab的griddata在"linear"模式下，对凸包外的点会自动执行线性外插（沿边缘三角的趋势延伸）；而Scipy的griddata和LinearNDInterpolator默认对凸包外的点返回NaN，不进行外插，这是两者视觉差异的核心原因。

3. 数值精度与计算细节
   两者在插值权重计算、浮点精度处理上存在细微差别，例如Matlab对浮点运算的舍入规则、中间结果的精度保留与Scipy略有不同，导致凸包内的数值也存在微小偏差。

解决方案：让Scipy结果与Matlab对齐 #

1. 调整Delaunay三角化参数，匹配Matlab的剖分逻辑 #

Matlab的Delaunay三角化允许共线点、保留退化三角，需要修改Scipy的Qhull参数来模拟：

from scipy.spatial import Delaunay
from scipy.interpolate import LinearNDInterpolator

# 设置Qhull参数，QJ允许共线点，Qc强制闭合凸包，QbB使用边界箱
tri = Delaunay(points, qhull_options='QJ Qc QbB')

2. 实现Matlab风格的外插逻辑 #

手动对凸包外的NaN点进行外插，模拟Matlab的线性外插行为（若要完全匹配，可选择线性外插或最近邻外插，Matlab默认是线性外插）：

# 创建基础线性插值器，不设置fill_value
interp = LinearNDInterpolator(tri, values, fill_value=None)
ZZ = interp(X, Y)

# 处理凸包外的NaN点：用最近邻插值补全（与Matlab外插效果接近）
mask = np.isnan(ZZ)
if np.any(mask):
    # 获取所有NaN点的坐标
    nan_x = X[mask]
    nan_y = Y[mask]
    from scipy.interpolate import NearestNDInterpolator
    nearest_interp = NearestNDInterpolator(points, values)
    ZZ[mask] = nearest_interp(nan_x, nan_y)

3. 统一数据精度 #

确保Python中使用与Matlab一致的float64精度，避免类型转换导致的偏差：

points = points.astype(np.float64)
values = values.astype(np.float64)
X = X.astype(np.float64)
Y = Y.astype(np.float64)

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