Scipy与PyTorch中2D三次插值方法的差异探究

Scipy和PyTorch均提供多种2D图像插值方式，但二者的三次插值结果存在显著差异。经过测试：

• 在order=1/linear/bilinear模式下，scipy.interpolate.interpn、scipy.ndimage.map_coordinates、torch.nn.functional.interpolate（设置align_corners=True）、torch.nn.functional.grid_sample（设置align_corners=True）的插值结果基本一致，仅存在浮点精度级别的误差；
• 在order=3（三次插值）模式下，上述函数的插值结果均不相同。

以下是对比测试代码及order=3时的插值差异可视化说明：

# 计算模块
import numpy as np
import numpy.typing
import torch
from scipy.interpolate import interpn
from scipy.ndimage import map_coordinates
from torch.nn.functional import interpolate, grid_sample

# 可视化模块
import matplotlib.pyplot as plt
from matplotlib.colors import LogNorm


def get_interpolated_values(
    original_size: tuple[int, int],
    sample_size: tuple[int, int],
    order: int,
) -> tuple[list[np.typing.NDArray], list[str]]:
    '''返回插值后的图像列表及对应方法的标签'''

    # 初始化插值的源数据
    values = np.random.rand(*original_size)
    order2str = [
        'nearest',
        'linear',
        None,
        'cubic',
    ]

    # 初始化数据收集容器
    labels = []
    interpolated_values = []

    # Scipy interpn
    labels.append('scipy.interpn')
    interpolated_values.append(
        interpn(
            points=(
                np.mgrid[0 : original_size[0]],
                np.mgrid[0 : original_size[1]],
            ),
            values=values,
            xi=np.stack(np.mgrid[
                0 : original_size[0] - 1 : sample_size[0] * 1j,
                0 : original_size[1] - 1 : sample_size[1] * 1j,
            ], axis=-1),
            method=order2str[order],
            bounds_error=True,
        )
    )

    # Scipy map_coordinates
    labels.append('scipy.map_coordinates')
    interpolated_values.append(
        map_coordinates(
            input=values,
            coordinates=np.mgrid[
                0 : original_size[0] - 1 : sample_size[0] * 1j,
                0 : original_size[1] - 1 : sample_size[1] * 1j,
            ],
            order=order,
            mode='constant',
        )
    )

    # Torch interpolate（align_corners=True）
    value_tensor = torch.from_numpy(values)[None, None, ...]
    labels.append('torch.interpolate')
    interpolated_values.append(
        interpolate(
            input=value_tensor,
            size=sample_size,
            align_corners=True,
            mode="bi" + order2str[order],
        ).squeeze().numpy()
    )

    # Torch grid_sample（align_corners=True）
    labels.append('torch.grid_sample')
    interpolated_values.append(
        grid_sample(
            input=value_tensor,
            grid=torch.from_numpy(
                np.mgrid[
                    -1 : 1 : sample_size[0] * 1j,
                    -1 : 1 : sample_size[1] * 1j,
                ],
            ).moveaxis(0, -1).unsqueeze(0).flip(-1),
            padding_mode='zeros',
            mode="bi" + order2str[order],
            align_corners=True,
        ).squeeze().numpy()
    )

    return interpolated_values, labels


def visualize(interpolated_values, labels):

    # 可视化各插值结果
    fig, ax_row = plt.subplots(1, len(interpolated_values))
    for i_ax, ax in enumerate(ax_row):
        ax.set_title(labels[i_ax])
        im = ax.imshow(interpolated_values[i_ax], vmin=0, vmax=1)
    plt.colorbar(im)

    # 可视化各方法间的差异
    fig, ax_grid = plt.subplots(
        len(interpolated_values),
        len(interpolated_values),
    )
    fig.suptitle('插值结果差异对比')

    for i_ax, ax_row in enumerate(ax_grid):
        for j_ax, ax in enumerate(ax_row):
            if i_ax == j_ax:
                # 对角线添加方法标签
                ax.text(
                    *((n - 1) / 2 for n in interpolated_values[i_ax].shape),
                    s=labels[i_ax],
                    ha='center',
                    va='center',
                )

            # 绘制差异矩阵及颜色条
            diff = interpolated_values[i_ax] - interpolated_values[j_ax]
            im = ax.imshow(np.abs(diff))
            plt.colorbar(im)


def main():
    interpolated_values, labels = get_interpolated_values(
        original_size = (4, 4),
        sample_size = (5, 5),
        order = 3,
    )
    visualize(interpolated_values, labels)
    plt.show()


if __name__ == "__main__":
    main()

插值差异可视化说明（order=3） #

可视化结果分为两部分：

1. 第一行子图展示四种插值方法各自的输出结果；
2. 下方矩阵子图展示每两种方法插值结果的绝对差异值，对角位置标注对应方法名称。

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