二维图像双线性插值 python 快速实现

在对二维数据进行 resize / mapping / 坐标转换等操作时，经常会将原本的整数坐标变换为小数坐标，对于非整数的坐标值一种直观有效的插值方式为双线性插值。

插值简介

双线性插值，又称为双线性内插。在数学上，双线性插值是有两个变量的插值函数的线性插值扩展，其核心思想是在两个方向分别进行一次线性插值。

双线性插值作为数值分析中的一种插值算法，广泛应用在信号处理，数字图像和视频处理等方面。

假设我们出现了需要在四个相邻正方形整数点$(A,B,C,D)$坐标中间（正方形范围内）选择一个点$(a,b)$取近似值的情形。

此时我们已知的是四个点的数值 $V_A,V_B,V_C,V_D$，给定小数坐标 $E(a,b),0\le a,b\le 1$，如何插值求解E点的数值呢，解决类似问题的方法统称为插值，上图展示公式为双线性插值的计算方法。

最近邻法 (Nearest Interpolation)

一种最简便的方法为最近邻法，直接取与当前点距离最近的点的值作为插值结果：
$$V_E=V_{round(a),round(b)}$$
其中 $round$ 为四舍五入的取整操作，方法简便速度极快，但往往不够精细

双三次插值 (Bicubic interpolation)

双三次插值是用原图像中16(4*4)个点计算新图像中1个点，效果比较好，但是计算代价过大。

双线性插值 (Bilinear Interpolation)

使用一个点进行插值过于粗暴，16个点又过于繁琐，那就使用$E$​点周围4个点的数值来近似求解，这是一种平衡了计算代价和插值效果的折中方案，也是各大变换库的默认插值操作。

双线性插值

通过观察上述动图（可以动手挪一挪）可以清晰地看到，双线性插值本质就是把四个角落的数值按照正方形面积的比例线性加权后的结果。

好吧一句话已经把数学的核心部分讲完了

那么既然理解了本质，数学公式就好写了：
$$\begin{array}{lll}{V}_E& =S_{blue}{V}_A+S_{red}{V}_B+S_{yellow}{V}_C+S_{purple}{V}_D& =(X_C-X_E)(Y_C-Y_E)V_A+(X_D-X_E)(Y_E-Y_D)V_B+(X_E-X_A)(Y_E-Y_A)V_C+(X_E-X_B)(Y_B-Y_E)V_D\end{array}$$

python 实现

在实现时当然 for 循环大法可以解决一切问题，但总归是不太优雅，我们尝试使用 numpy 操作完成双线性插值

假设原始图像 image，变换后的小数坐标 X 矩阵 x_grid，Y 矩阵 y_grid，那么可以使用如下的 bilinear_by_meshgrid 函数快速双线性插值，已经处理好了边界，可以放心使用。

def bilinear_by_meshgrid(image, x_grid, y_grid):

    #               Ia, Wd                          Ic, Wb
    #           (floor_x, floor_y)              (ceil_x, floor_y)   
    #
    #                               (x, y)
    #
    #               Ib , Wc                         Id, Wa
    #           (floor_x, ceil_y)               (ceil_x, ceil_y)   
    #

    assert image.shape == x_grid.shape == y_grid.shape
    assert image.ndim == 2
    H, W = image.shape[:2]

    floor_x_grid = np.floor(x_grid).astype('int32')
    floor_y_grid = np.floor(y_grid).astype('int32')
    ceil_x_grid = floor_x_grid + 1
    ceil_y_grid = floor_y_grid + 1

    if np.max(ceil_x_grid) > W -1 or  np.max(ceil_y_grid) > H -1 or np.min(floor_x_grid) < 0 or np.min(floor_y_grid) < 0:
        print("Warning: index value out of original matrix, a crop operation will be applied.")

        floor_x_grid = np.clip(floor_x_grid, 0, W-1).astype('int32')
        ceil_x_grid = np.clip(ceil_x_grid, 0, W-1).astype('int32')
        floor_y_grid = np.clip(floor_y_grid, 0, H-1).astype('int32')
        ceil_y_grid = np.clip(ceil_y_grid, 0, H-1).astype('int32')

    Ia = image[ floor_y_grid, floor_x_grid ]
    Ib = image[ ceil_y_grid, floor_x_grid ]
    Ic = image[ floor_y_grid, ceil_x_grid ]
    Id = image[ ceil_y_grid, ceil_x_grid ]

    wa = (ceil_x_grid - x_grid) * (ceil_y_grid - y_grid)
    wb = (ceil_x_grid - x_grid) * (y_grid - floor_y_grid)
    wc = (x_grid - floor_x_grid) * (ceil_y_grid - y_grid)
    wd = (x_grid - floor_x_grid) * (y_grid - floor_y_grid)

    assert np.min(wa) >=0 and np.min(wb) >=0 and np.min(wc) >=0 and np.min(wd) >=0
    
    W = wa + wb + wc + wd
    assert np.sum(W[:, -1]) + np.sum(W[-1, :]) == 0
    
    wa[:-1, -1] = ceil_y_grid[:-1, -1] - y_grid[:-1, -1]
    wb[:-1, -1] = y_grid[:-1, -1] - floor_y_grid[:-1, -1]
    
    wb[-1, :-1] = ceil_x_grid[-1, :-1] - x_grid[-1, :-1]
    wd[-1, :-1] = x_grid[-1, :-1] - floor_x_grid[-1, :-1]
    
    wd[-1, -1] = 1
    
    W = wa + wb + wc + wd
    assert np.max(W) == np.min(W) == 1
    
    res_image = wa*Ia + wb*Ib + wc*Ic + wd*Id

    return res_image
PYTHON

该函数集成在我自己的python库 mtutils 中，可以通过：

pip install mtutils
CMAKE

直接安装，之后可以直接引用：

from mtutils import bilinear_by_meshgrid
ANGELSCRIPT

参考资料

• https://baike.baidu.com/item/双线性插值/11055945?fr=aladdin

• https://zhuanlan.zhihu.com/p/110754637

文章链接：
https://www.zywvvd.com/notes/study/image-processing/bilinear-interpolation/bilinear-interpolation/