双线性插值原理及Python实现

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

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

转载请注明出处

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

我们先跳出来想一下：

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

回到上面的图例：

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

接下来就有两个问题。

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

已知：

1. 输入图的高和宽：，
2. 目标图的高和宽：，

投影的公式是：

但是这个公式有一个问题，可能会导致目标图的中心跟原图的中心不对齐。例如：原图是3x3，中心点坐标(1, 1)；目标图是9x9，中心点坐标(4, 4)；通过上面的公式计算，目标图中心点在原图的投影坐标：。
之所以会出现中心点不对齐，原因是每个像素点实际上是一个边长为1的正方形，所以对于坐标为(, )的像素点，它的中心其实是(, )。所以精确计算应该是：

转换一下得到正确的投影公式：

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

已知：

1. 四个红色点的坐标值：，，，
2. 四个红色点的像素值：，，，
3. 投影点的坐标值：，

思路是：每个点的像素值乘以各自的权重，然后相加得到投影点像素值。点跟点的距离越近，它的权重就越大。
双线性插值给出的算法很是简单粗暴：先在横轴方向上进行两次线性插值计算，然后在纵轴方向上进行一次插值计算。结合最开始那个图例看，就是先求和这两个蓝色点的像素值，然后再通过这两个值，求得点的像素值。
具体计算如下：

前面说了，四个红色点是投影点四周最近的点，显然四个红色点彼此间的距离都是1，也即，。上式可以写成：

再令：，，式子可以进一步写成：

至此，关于双线性插值的原理就全部讲完了。下面是用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