双线性插值原理及Python实现

转载请注明出处：https://www.jianshu.com/p/29e5c84ea539 如果觉得有用，麻烦点个赞噢~
这里不对背景知识做讲述，只讲双线性插值的原理，以及用Python实现代码。

下面是关于双线性插值的经典说明图例：

转载请注明出处

首先，图中有5个像素点： $Q_{00}$ ， $Q_{01}$ ， $Q_{10}$ ， $Q_{11}$ ， $P$ 。其中四个红色点 $Q$ 是原图的点，绿色点 $P$ 就有意思了，理解了它的含义，对双线性插值原理就能有一个概念上的理解。
这个点 $P$ ，其实是目标图的像素点在原图上的投影！

我们先跳出来想一下：

双线性插值用来做什么？在CV领域，通常是用来改变原图的尺寸。
那么目标图的每个像素点的像素值应该取多少呢？双线性插值的思路是，对于目标图的每个像素点，找到它在原图上最相关的四个点通过插值计算得到它的像素值。
原图上哪四个点是最相关的呢？我们可以把目标图的点，投影回到原图上，投影点四周距离最近的四个原图点就是了。

回到上面的图例：

四个红色点 $Q$ 就是投影点 $P$ 的四周最近的点。通过四个红色点 $Q$ ，可以计算出投影点 $P$ 的像素值，这样目标图上的像素点的像素值也就得到了。

接下来就有两个问题。

目标图的点如何投影到原图上？

已知：

输入图的高和宽： $height_{src}$ ， $width_{src}$
目标图的高和宽： $height_{dst}$ ， $width_{dst}$

投影的公式是：
$h = h_{dst} \ast \frac{height_{src}}{height_{dst}}$
$w = w_{dst} \ast \frac{width_{src}}{width_{dst}}$

但是这个公式有一个问题，可能会导致目标图的中心跟原图的中心不对齐。例如：原图是3x3，中心点坐标(1, 1)；目标图是9x9，中心点坐标(4, 4)；通过上面的公式计算，目标图中心点在原图的投影坐标： $h = 4 * 3 / 9 = 1.3333 \neq 1$ 。
之所以会出现中心点不对齐，原因是每个像素点实际上是一个边长为1的正方形，所以对于坐标为( $h$ , $w$ )的像素点，它的中心其实是( $h+0.5$ , $w+0.5$ )。所以精确计算应该是：
$\frac{h + 0.5}{height_{src}} = \frac{h_{dst} + 0.5}{height_{dst}}$
$\frac{w + 0.5}{width_{src}} = \frac{w_{dst} + 0.5}{width_{dst}}$

转换一下得到正确的投影公式：
$h = (h_{dst} + 0.5) \ast \frac{height_{src}}{height_{dst}} - 0.5$
$w = (w_{dst} + 0.5) \ast \frac{width_{src}}{width_{dst}} - 0.5$

如何插值计算得到投影点的像素值？

已知：

四个红色点 $Q$ 的坐标值： $h_0$ ， $h_1$ ， $w_0$ ， $w_1$
四个红色点 $Q$ 的像素值： $f(Q_{00})$ ， $f(Q_{01})$ ， $f(Q_{10})$ ， $f(Q_{11})$
投影点 $P$ 的坐标值： $h$ ， $w$

思路是：每个 $Q$ 点的像素值乘以各自的权重，然后相加得到投影点 $P$ 像素值。 $Q$ 点跟 $P$ 点的距离越近，它的权重就越大。
双线性插值给出的算法很是简单粗暴：先在横轴方向上进行两次线性插值计算，然后在纵轴方向上进行一次插值计算。结合最开始那个图例看，就是先求 $R_0$ 和 $R_1$ 这两个蓝色点的像素值，然后再通过这两个值，求得 $P$ 点的像素值。
具体计算如下：
$f(R_0) \approx \frac{w_1 - w}{w_1 - w_0} f(Q_{00}) + \frac{w - w_0}{w_1 - w_0} f(Q_{01})$
$f(R_1) \approx \frac{w_1 - w}{w_1 - w_0} f(Q_{10}) + \frac{w - w_0}{w_1 - w_0} f(Q_{11})$
$f(P) \approx \frac{h_1 - h}{h_1 - h_0} f(R_0) + \frac{h - h_0}{h_1 - h_0} f(R_1)$
$\approx \frac{h_1 - h}{h_1 - h_0} \left ( \frac{w_1 - w}{w_1 - w_0} f(Q_{00}) + \frac{w - w_0}{w_1 - w_0} f(Q_{01}) \right ) + \frac{h - h_0}{h_1 - h_0} \left (\frac{w_1 - w}{w_1 - w_0} f(Q_{10}) + \frac{w - w_0}{w_1 - w_0} f(Q_{11}) \right )$
$= \frac{1}{(w_1 - w_0)(h_1 - h_0)} \left ( (h_1 - h)(w_1 - w)f(Q_{00}) + (h_1 - h)(w - w_0)f(Q_{01}) + (h - h_0)(w_1 - w)f(Q_{10}) + (h - h_0)(w - w_0)f(Q_{11}) \right )$

前面说了，四个红色点 $Q$ 是投影点 $P$ 四周最近的点，显然四个红色点彼此间的距离都是1，也即 $w_1 - w_0 = 1$ ， $h_1 - h_0 = 1$ 。上式可以写成：
$f(P) \approx (h_1 - h)(w_1 - w)f(Q_{00}) + (h_1 - h)(w - w_0)f(Q_{01}) + (h - h_0)(w_1 - w)f(Q_{10}) + (h - h_0)(w - w_0)f(Q_{11})$
再令： $u = h - h_0$ ， $v = w - w_0$ ，式子可以进一步写成：
$f(P) \approx f(Q_{00})(1 - u)(1 - v) + f(Q_{01})(1 - u)v + f(Q_{10})u(1 - v) + f(Q_{11})uv$

至此，关于双线性插值的原理就全部讲完了。下面是用Python实现的代码。
先来一个最直观的写法：

def bilinear_interpolate(src, dst_size):
    height_src, width_src, channel_src = src.shape  # (h, w, ch)
    height_dst, width_dst = dst_size  # (h, w)
    
    """
    中心对齐，投影目标图的横轴和纵轴到原图上
    """
    ws_p = np.array([(i + 0.5) / width_dst * width_src - 0.5 for i in range(width_dst)], dtype=np.float32)
    hs_p = np.array([(i + 0.5) / height_dst * height_src - 0.5 for i in range(height_dst)], dtype=np.float32)
    ws_p = np.clip(ws_p, 0, width_src-1)  # 实验发现要这样子来一下才能跟torch的输出结果一致
    hs_p = np.clip(hs_p, 0, height_src-1)
    
    """找出每个投影点在原图横轴方向的近邻点坐标对"""
    # w_0的取值范围是 0 ~ (width_src-2)，因为w_1 = w_0 + 1
    ws_0 = np.clip(np.floor(ws_p), 0, width_src-2).astype(np.int)
        
    """找出每个投影点在原图纵轴方向的近邻点坐标对"""
    # h_0的取值范围是 0 ~ (height_src-2)，因为h_1 = h_0 + 1
    hs_0 = np.clip(np.floor(hs_p), 0, height_src-2).astype(np.int)
        
    """
    计算目标图各个点的像素值
    f(h, w) = f(h_0, w_0) * (1 - u) * (1 - v)
            + f(h_0, w_1) * (1 - u) * v
            + f(h_1, w_0) * u * (1 - v)
            + f(h_1, w_1) * u * v
    """
    dst = np.zeros(shape=(height_dst, width_dst, channel_src), dtype=np.float32)
    us = hs_p - hs_0
    vs = ws_p - ws_0
    _1_us = 1 - us
    _1_vs = 1 - vs
    for h in range(height_dst):
        h_0, h_1 = hs_0[h], hs_0[h]+1  # 原图的坐标
        for w in range(width_dst):
            w_0, w_1 = ws_0[w], ws_0[w]+1 # 原图的坐标
            for c in range(channel_src):
                dst[h][w][c] = src[h_0][w_0][c] * _1_us[h] * _1_vs[w] \
                            + src[h_0][w_1][c] * _1_us[h] * vs[w] \
                            + src[h_1][w_0][c] * us[h] * _1_vs[w] \
                            + src[h_1][w_1][c] * us[h] * vs[w]
    return dst

if __name__ == '__main__':
    src = np.array([[1, 2, 3], [3, 4, 5], [6, 7, 8]])
    src = np.expand_dims(src, axis=2)
    dst = bilinear_interpolate(src, dst_size=(src.shape[0]*2, src.shape[1]*2))
    print(dst[:, :, 0])

代码里的三个for循环，其实可以改成使用numpy的矩阵运算，虽然不容易读懂代码，但速度提升是巨大的。例如，对一张(333, 600, 3)的原图，尺寸增加一倍到(666, 1200, 3)，前一份代码耗时20多秒，而下面的代码仅需182毫秒。代码如下：

"""
将for循环改成numpy的矩阵运算，提速
"""

import numpy as np
import math

def bilinear_interpolate(src, dst_size):
    height_src, width_src, channel_src = src.shape  # (h, w, ch)
    height_dst, width_dst = dst_size  # (h, w)

    """中心对齐，投影目标图的横轴和纵轴到原图上"""
    ws_p = np.array([(i + 0.5) / width_dst * width_src - 0.5 for i in range(width_dst)], dtype=np.float32)
    hs_p = np.array([(i + 0.5) / height_dst * height_src - 0.5 for i in range(height_dst)], dtype=np.float32)
    ws_p = np.clip(ws_p, 0, width_src-1)  # 实验发现要这样子来一下才能跟torch的输出结果一致
    hs_p = np.clip(hs_p, 0, height_src-1)
    ws_p = np.repeat(ws_p.reshape(1, width_dst), height_dst, axis=0)
    hs_p = np.repeat(hs_p.reshape(height_dst, 1), width_dst, axis=1)

    """找出每个投影点在原图的近邻点坐标"""
    ws_0 = np.clip(np.floor(ws_p), 0, width_src - 2).astype(np.int)
    hs_0 = np.clip(np.floor(hs_p), 0, height_src - 2).astype(np.int)
    ws_1 = ws_0 + 1
    hs_1 = hs_0 + 1

    """四个临近点的像素值"""
    f_00 = src[hs_0, ws_0, :].T
    f_01 = src[hs_0, ws_1, :].T
    f_10 = src[hs_1, ws_0, :].T
    f_11 = src[hs_1, ws_1, :].T

    """计算权重"""
    w_00 = ((hs_1 - hs_p) * (ws_1 - ws_p)).T
    w_01 = ((hs_1 - hs_p) * (ws_p - ws_0)).T
    w_10 = ((hs_p - hs_0) * (ws_1 - ws_p)).T
    w_11 = ((hs_p - hs_0) * (ws_p - ws_0)).T

    """计算目标像素值"""
    return (f_00 * w_00).T + (f_01 * w_01).T + (f_10 * w_10).T + (f_11 * w_11).T


if __name__ == '__main__':
    src = np.array([[1, 2, 3], [3, 4, 5], [6, 7, 8]])
    src = np.expand_dims(src, axis=2)
    print(src.shape)
    dst = bilinear_interpolate(src, dst_size=(src.shape[0] * 2, src.shape[1] * 2))
    print(dst.shape)
    print(dst[:, :, 0])

参考：
[1] https://www.cnblogs.com/yssongest/p/5303151.html
[2] https://github.com/ddbourgin/numpy-ml