一文彻底解决 双线性插值 的数学推导和实际应用

 

插值法是一种根据原图（source）图片信息构造目标图像（destination）的方法。常用在传统图像处理以及神经网络中的上采样。

 

例如此处我们需要将3×3的单通道图片扩展为4×4的单通道图片，此处数值应该视为对应像素的像素值，后面数学推导中的坐标，应视为像素的中点，而不是边界角点。

毕竟是有几何关系中的相似，由比例关系得：

 

                                                                                             公式 1

里面的符号在很多地方都有解释，src是源图，dst目标图，x，y分别是坐标，width、height分别是宽高。具体的对应不用记住，直接可以由两张图的相似成比例的关系推导即可。确定比例关系，很容易得到下面的公式（用不用都可以，还是上面容易记）

                                                                          公式 2

例如，上例中想求dst(0, 1)的像素值，需要知道dst(0, 1)对应原图src的哪个坐标，即用上述相似成比例即可：

                         

所以 dst(0, 1) = src(0, 0.75)

知道对应关系之后就需要确定dst(0, 1)的像素值了，因为核心想求的也是像素值。但是发现dst像素大多是与原图没有完全对应的像素坐标的（因为小数导致坐标不等嘛），此处就有不同的方法了。

如果此时直接将所求的4×4的 dst 坐标进行四舍五入，即讲像素点选择靠近理他最近的整数像素点，这个方法就是最近邻插值法。得到结果如下：

   

这个案例分辨率较小，可能看不出太多问题。但是可以想象，当较小分辨率图象插值为较大图象时，将存在严重的失真现象。因为周边局部块的像素都是一样的，并没有平滑过度，导致马赛克十分严重。

为了解决上述问题，就有了双线性插值法。

双线性插值为对两个方向进行线性插值，所以有必要了解下什么是线性插值：

                                           

此处不要将 x 和 y 看成坐标对，应该是将其看成变量和函数的值，例如此处为一条灰度随着长度均匀变化的单像素线条， x 为线条的长度坐标（自变量），y为灰度值。这样描述之后，问题变为一直x1和x2和其对应的灰度值，求解坐标为 x 的点的灰度值。同样由于线性关系，可以得到：

                                                                                                             公式 3

  简单整理有如下：

                                                                                                      公式 4

其实就理解来说，第一个公式更好理解和记忆，而第二个将灰度值y 写成 由y1和y2各带一个权重求和，更加符合我们线性插值的含义（中间点像素值由两边点像素值一定的加权组合而成），其实也是更加方便后续计算。

下面来看双线性变换法：

开始之前，注意，此时就不是单线条问题，而是二维图象，故 x，y 就是一对坐标对，而我们想进行线性插值求解的并不是y，y像素的是纵坐标，我们需要求解的是（x, y) 对应的像素值P(x，y)。双线性体现在，我们不能像上述线性插值一样找到合适的线性的关系线将代求（x,y）包含在两确定之间（此处不考虑特殊情况，如果像下图，待求点（x, y）与（x, y2) 重合，则线性就可以解决)。既然一步不能解决，就分两步，首先由水平的线性插值求到（x, y2）和（x, y1）的像素值R2、R1；再在竖直上由（x, y2）和（x, y1）求解（x, y）的像素值 P(x, y)，即完成求解。此处重申一下，我们用线性插值求解的一直都是像素值，而其各点的坐标值，其实仍然是利用 公式1 或者  公式2 求解的，此处一定不能晕乎。

                                                       

讲述完原理后，我们看下数学过程，相信理解之后很容易推导了：

       

       

      

将上述两步得到的结果一综合，就会得到哪个很长很长的式子，过程很简单。这里我就简化一下最后的结果：

   

类似与线性插值的结果，从这也能看出，这样的表示方式反映了待求点的像素值与四周的四个点的像素值相关。

到这里我们再来看下之前案例中的求解 dst(0, 1)  的像素值的双线性变化的实现：

        由公式1或者公式2，求的dst(0, 1)对应的原图为(0, 0.75)，易知其被(0，0), (0，1),  (1，0),  (1，1) 四个点包围，按照上述公式求解即可。即得到det（0，1）的像素值。这样依次求解上述 dst 中的所有点的像素值，即得到完整的插值后的图片。

此处再借用引用资料的博主的一个非常清晰图，来说明整个计算过程：

                                 

 

优化部分：

上述推导是基本的原理，下面我们看个例子，假设源图像是 3×3，中心点坐标（1，1）目标图像是9×9，中心点坐标（4，4），我们在进行插值映射的时x候，尽可能希望均匀的用到源图像的像素信息，最直观的就是（4,4）映射到（1,1）现在直接计算   srcX = 4*3/9=1.3333  !=  1 ，也就是我们在插值的时候所利用的像素集中在图像的右下方，而不是均匀分布整个图像。为此我们采用中心的策略：

如：源图像和目标图像几何中心的对齐：（用下面公式替换公式2）

                            公式5

此时，再来计算，srcX=(4+0.5)*3/9-0.5=1，符合对齐要求。

上述为什么要加0.5，然后减0.5，从计算角度，随便举几个例子很容易符合结论，对于实际含义我的理解是 这是实现了 原图中心点的对齐，下面是我的理解：

首先不要离散化看待问题，因为线性插值也是连续化的看待了像素。所以这里将dst的图象保持像素数目不变，尺寸变为和src相同的大小，这样就相当于找到了dst中像素的插值来源，即dst的像素实际是由原图src的哪个点得到。

      下面左图，为用公式2求出来的 dst 和 src 的图片的对应关系，图不重要，主要是体现对应关系，很容易知道两者的第一个像素点的中心肯定是重合的，因为符合公式2嘛。这样的话就不能对齐了，我们理想的应该是右图，图象左上角对齐，这也是公式5做的事情，在各自维度上平移一个0.5， 相当于将起点坐标对应到了样就 完整得到 dst 和 src 的对应关系了。这样做的目的就是实现中心的对应，其实也是为了方便计算，为什么这么讲，左图的时候是以像素坐标值来计算插值，而使用公式5相当于将坐标起点移到左上角，这样的计算就能够很好的满足一一对应关系了！（主要还是中心点对齐了）

                                                 

   

　

 

代码（源于参考链接）：

import cv2
import time
from math import ceil, floor
import numpy as np
 
def bilinear_interpolation(img, out_dim):
    src_h, src_w, channel = img.shape
    dst_h, dst_w = out_dim[1], out_dim[0]
    if src_h == dst_h and src_w == dst_w:
        return img.copy()
    dst_img = np.zeros((dst_h, dst_w, channel), dtype=np.uint8)  #创建一个数组，通过往里面填值，形成新的图片
    scale_x, scale_y = float(src_w) / dst_w, float(src_h) / dst_h #计算 dst与src的比例系数
    for i in range(channel):
        for dst_y in range(dst_h):
            for dst_x in range(dst_w):
 
                # 计算dst(x,y)对应回源图src的哪个坐标
                # 使用几何中心对称的方法
                # 如果不使用几何中心对成就写成： src_x = dst_x * scale_x
                src_x = (dst_x + 0.5) * scale_x - 0.5
                src_y = (dst_y + 0.5) * scale_y - 0.5
 
                # 找出用于插值的四个邻近点坐标，(x0,y0),(x0,y1),(x1,y0),(x1,y1)
                src_x0 = int(floor(src_x))
                src_x1 = min(src_x0 + 1, src_w - 1) #与边界点比，取小的一个,减1是因为是从0开始算的
                src_y0 = int(floor(src_y))
                src_y1 = min(src_y0 + 1, src_h - 1) #与边界点比，取小的一个
 
                temp0 = ((src_x1 - src_x) * img[src_y0, src_x0,i] + (src_x - src_x0) * img[src_y0, src_x1, i])
                temp1 = (src_x1 - src_x) * img[src_y1, src_x0,i] + (src_x - src_x0) * img[src_y1, src_x1, i]
                dst_img[dst_y, dst_x, i] = int((src_y1 - src_y) * temp0 + (src_y - src_y0) * temp1)
    return dst_img
 
 
if __name__ == '__main__':
    img = cv2.imread(r'C:\Users\Administrator\Desktop\pic\car.jpg')
    start = time.time()
    dst = bilinear_interpolation(img, (1000, 1000))
    print('cost {} seconds'.format(time.time() - start))
    cv2.imshow('result', dst)
    cv2.waitKey()

效果：

Reference：

https://blog.csdn.net/u014453898/article/details/105584750

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

感谢前辈的贡献！