我与插值萍水相逢续(1): 双线性插值(Bilinear Interpolation)原理及使用

    上一篇博文说到我了解的插值有两个应用：一个是补偿数据的缺失，另一个是数据缩放。上一篇关于线性插值的文章中主要讲述了补偿数据缺失这一种应用，而这一篇关于双线性插值的文章，我准备用数据缩放来举例子。

1.1.关于双线性插值

    关于双线性插值，我打算从维基百科中的一幅关于双线性插值原理的图讲起。

    在这幅图中，假设有四个已知的数据点，Q11(x1,y1)、Q21(x2,y1)、Q12(x1,y2)、Q22(x2,y2)，而我们要在P(x,y)这个位置进行插值。会怎么做？

    我会做三次线性插值：第一次，在Q11和Q21间来一次，得到R1；第二次，在Q12和Q22间来一次，得到R2；第三次，在R1和R2间来一次，得到P。好的，就这么简单。这就是双线性插值的原理。

    接下来我们列数学表达式，这里我们可以参考维基百科中给出的推导：

    （https://zh.wikipedia.org/wiki/%E5%8F%8C%E7%BA%BF%E6%80%A7%E6%8F%92%E5%80%BC）

    假设(x,y)处的值等于f(x,y)，那么P处的值为f(P)，Q11、Q21、Q12、Q22处的值为f(Q11)、f(Q21)、f(Q12)、f(Q22)。那么在x方向进行的两次插值，得到R1和R2（当然也可以是y方向进行两次插值）。

    注意：第一次插值怎么做，R1怎么得到，要留心的是Q11和Q21对应的数据不是y坐标的值，而是它们各自的高度值，比如Q11处的值f(Q11)=f(x1,y1)=100。是这个样子的，所以我要把线性插值的类似的图再贴一遍，

    根据这个示意图，我们可以得到第一第二次的插值：

    通过合并同类项，我们就可以得到维基百科上的式子了：

    那么y方向上的插值也可以通过这样的方式得到：

    那么所要的结果：

    到这里，我们就可以得到需要插值的P处的值了。以上就是关于双线性插值的原理。

    在这里，有个情况下，可以简化上述公式。即假如选择一个坐标系可以使得 {\displaystyle f} 的四个已知点坐标分别为 (0, 0)、(0, 1)、(1, 0) 和 (1, 1)，那么插值公式就可以化简为

    合并同类项，可以得到：（要注意：这里的x和y不再是原来的x和y，而是在更换坐标系之后的x和y，类比到下面的算法中，是u和v）

                                                                   

    其中，

1.2.双线性插值应用于缩放的算法

    这里要介绍一篇博客http://blog.csdn.net/huang1024rui/article/details/46545329，虽然这位博主是转载的，不过很有帮助的。他做了一个真彩图像的缩放应用案例。我这里就针对一般的数据矩阵的缩放，即二维矩阵的缩放给出通用算法。如若需要进行三维数据矩阵进行缩放，可以去参考给出的url中的博文。

    下面先给出根据维基百科中最后化简的式子做的一个算法：叫算法1

-------------------------------------------------------------------------------------------------------------------

    算法准备：1）设计function输入：二维数据矩阵Data、缩放比例Rate。其中Rate>0，若0<Rate<1，则是缩小；若Rate>1，则是放大；若Rate=1，则不做缩放变换。

    2）设计function输出：返回二维矩阵Data_ret，为缩放后的数据矩阵。

    step1: 获得Data的行列维度值[m,n]，构建一个缩放过后的全零矩阵Data_ret，它的维度为[m*Rate,n*Rate]。

    step2: Data_ret中的一个位置(r_i,r_j)，映射回原矩阵Data中，在Data中的位置是(x,y)=(r_i/Rate,r_j/Rate)，对(x,y)进行向下取整，得到(i,j)。x = i +u, y = j + v，得到u和v（向下取整的时候会遇到i,j为0的情况，这种情况matlab是不允许的，需要要设置规定限值i，j的值）。

    step3:根据双线性插值的方法，满足简化的公式，可以得到Data_ret(r_i,r_j)=b1+b2*u+b3*v+b4*u*v，其中b1 = Data(i+1,j), b2 = Data(i+1,j+1) - b1, b3 = Data(i,j) - b1, b4 = Data(i,j+1) - Data(i+1,j+1) - Data(i, j) + b1（另外，整数倍率时，使用这种算法会遇到索引超出矩阵维度的情况，同样也需要进行限定，如果超出维度，就用边界两点进行线性插值，step2也是一样）。

    step4: 重复step2、3，直到Data_ret中的每个位置都插值完毕。

-------------------------------------------------------------------------------------------------------------------

    我再给出看了这一节开头提到的博客中的算法后写的算法：叫算法2。我为什么又要在这里提算法2，这是因为算法1我发现效果不怎么好，在最后的实验部分我们来看。（两个算法基于维基百科的公式，看起来是差不多的，差别在哪里呢？）

-------------------------------------------------------------------------------------------------------------------

     算法准备：1）设计function输入：二维数据矩阵Data、缩放比例Rate。其中Rate>0，若0<Rate<1，则是缩小；若Rate>1，则是放大；若Rate=1，则不做缩放变换。

    2）设计function输出：返回二维矩阵Data_ret，为缩放后的数据矩阵。

    step1: 获得Data的行列维度值[m,n]，构建一个缩放过后的全零矩阵Data_ret，它的维度为[m*Rate,n*Rate]。

    step2: Data_ret中的一个位置(r_i,r_j)，映射回原矩阵Data中，在Data中的位置是(x,y)=(r_i/Rate,r_j/Rate)，对(x,y)进行向下取整，得到(i,j)。x = i +u, y = j + v，得到u和v（向下取整的时候会遇到i,j为0的情况，这种情况matlab是不允许的，需要要设置规定限值i，j的值）。

    step3: 使用data_ret(i_r,j_r) = (1-u)*(1-v)*data(i,j) +(1-u)*v*data(i,j+1)+ u*(1-v)*data(i+1,j) +u*v*data(i+1,j+1)进行双线性插值（另外，整数倍率时，使用这种算法会遇到索引超出矩阵维度的情况，同样也需要进行限定，如果超出维度，就用边界两点进行线性插值，step2也是一样）。

    step4: 重复step2、3，直到Data_ret中的每个位置都插值完毕。

-------------------------------------------------------------------------------------------------------------------

1.3. 双线性插值应用案例

    选取MATLAB的toolbox/images/imdata中的一张图片，叫'logo.tif'，是一张Mathworks的logo，是一个二维矩阵，灰度图。我们准备将它分别放大两倍，缩小两倍。（如要使用像上一篇博文那样的插值，则使用matlab提供的interp2()函数就可以了）

    法1：使用MATLAB提供的imresize()函数

-------------------------------------------------------------------------------------------------------------------

im_gray = imread('logo.tif');
imshow(im_gray);
im_gray_bigger = imresize(im_gray,2,'bilinear');
figure;imshow(im_gray_bigger);
im_gray_smaller = imresize(im_gray,0.5,'bilinear');
figure;imshow(im_gray_smaller);

-------------------------------------------------------------------------------------------------------------------

    法2：自己编写self_imresize()函数，实现放大缩小的功能（和上一篇一样，并没有考虑太多异常处理，仅仅去实现基本的功能）

    写一个m-file，取名叫self_imresize.m，内容如下：

    先来看算法1版本：

-------------------------------------------------------------------------------------------------------------------

% 任务目标，实现一个对数据进行缩放的接口，使用双线性插值的方法
% 这种算法有一个缺陷，下一篇博文中讨论
function data_ret = self_imresize(data, rate, method)
% 判断方法是否是双线性插值
if strcmp(method,'bilinear'),
    % 进行双线性插值：1 构建一个缩放后的全零阵
    [lin_o, col_o] = size(data);
    lin_r = floor(lin_o * rate);
    col_r = floor(col_o * rate);
    data_ret = zeros(lin_r, col_r);
    % 将全零阵中的每个位置映射到原矩阵中去并且根据在原矩阵中的位置进行双线性插值
    % 打算用(x,y)记录映射到原矩阵中的位置,(i_r,j_r)为data_ret中的位置，(i,j)
    % 打算用来记录(x,y)所在的类似双线性插值的左上角的已知点位置
    for i_r = 1:lin_r,
        for j_r = 1:col_r,
           x = i_r / rate;
           y = j_r / rate;
           i = floor(x);
           j = floor(y);
           % 对于i和j这里必须要考虑一个问题，就是matlab的矩阵索引是从1开始的
           % 这里不用线性插值了，先粗糙处理一下
           % 如果i,j出现了等于0，那么就要变成1
           if i == 0,
               i = 1;
               x = x + 1;
           end
           if j == 0,
               j = 1;
               y = y + 1;
           end
           % 这里解决索引超出矩阵维度的问题，在整数倍率时候会遇到
           % 一个用下面的if条件来限定，或者在上面
           % x = i_r / rate后面减去一个0.0001就可以解决问题
           if i == lin_o,
               i = i - 1;
               x = x - 1;
           end
           if j == col_o,
               j = j - 1;
               y = y - 1;
           end
           b1 = data(i + 1, j);
           b2 = data(i + 1, j + 1) - b1;
           b3 = data(i, j) - b1;
           b4 = data(i, j + 1) - data(i + 1, j + 1) - data(i, j) + b1;
           u = x - i;
           v = y - j;
           data_ret(i_r, j_r) = b1 + b2 * u + b3 * v + b4 * u * v;
        end
    end
% 如果不是双线性插值或者没输，是不行的
else
    error('please check out the method!\n');
end
end

-------------------------------------------------------------------------------------------------------------------

    然后在命令行或者另起一个m-file：写入运行代码

-------------------------------------------------------------------------------------------------------------------

im_gray = imread('logo.tif');
figure;imshow(im_gray);
im_gray_bigger = self_imresize(im_gray,2,'bilinear');
figure;imshow(im_gray_bigger);
im_gray_smaller = self_imresize(im_gray,0.5,'bilinear');
figure;imshow(im_gray_smaller);

-------------------------------------------------------------------------------------------------------------------

    结果如下：

    我们再来看算法2版本：改写self_imresize.m文件内容如下：

-------------------------------------------------------------------------------------------------------------------

% 任务目标，实现一个对数据进行缩放的接口，使用双线性插值的方法
% 这种算法有一个缺陷，下一篇博文中讨论
function data_ret = self_imresize(data, rate, method)
% 判断方法是否是双线性插值
if strcmp(method,'bilinear'),
    % 进行双线性插值：1 构建一个缩放后的全零阵
    [lin_o, col_o] = size(data);
    lin_r = floor(lin_o * rate);
    col_r = floor(col_o * rate);
    data_ret = zeros(lin_r, col_r);
    % 将全零阵中的每个位置映射到原矩阵中去并且根据在原矩阵中的位置进行双线性插值
    % 打算用(x,y)记录映射到原矩阵中的位置,(i_r,j_r)为data_ret中的位置，(i,j)
    % 打算用来记录(x,y)所在的类似双线性插值的左上角的已知点位置
    for i_r = 1:lin_r,
        for j_r = 1:col_r,
           x = i_r / rate;
           y = j_r / rate;
           i = floor(x);
           j = floor(y);
           % 对于i和j这里必须要考虑一个问题，就是matlab的矩阵索引是从1开始的
           % 这里不用线性插值了，先粗糙处理一下
           % 如果i,j出现了等于0，那么就要变成1
           if i == 0,
               i = 1;
               x = x + 1;
           end
           if j == 0,
               j = 1;
               y = y + 1;
           end
           % 这里解决索引超出矩阵维度的问题，在整数倍率时候会遇到
           % 一个用下面的if条件来限定，或者在上面
           % x = i_r / rate后面减去一个0.0001就可以解决问题
           if i == lin_o,
               i = i - 1;
               x = x - 1;
           end
           if j == col_o,
               j = j - 1;
               y = y - 1;
           end
           u = x - i;
           v = y - j;
           data_ret(i_r,j_r) = (1-u)*(1-v)*data(i,j) +(1-u)*v*data(i,j+1)...
           + u*(1-v)*data(i+1,j) +u*v*data(i+1,j+1);      
        end
    end
% 如果不是双线性插值或者没输，是不行的
else
    error('please check out the method!\n');
end
end

-------------------------------------------------------------------------------------------------------------------

    像算法1的使用方法一样运行一样一样的运行代码，结果如下：

    发现没有，算法2的第二张图也就是放大2倍的图更清楚一点。

    当然了，光看没有用，我们来计算这两张最大的看起来清楚程度有点不一样的图的距离。使用matlab中norm()函数计算两个矩阵的距离，得到22.3060。可见，确实是不一样的。

    这时，我不禁想问，两个算法用的都是维基百科中的公式，我自己推导了也没问题，差别到底在哪里呢？

    我分析了这两个算法发现，差别就在step3做双线性插值的时候，算法1在写的时候的四个已知点我是选择的左下角的数据点作为坐标原点的，而算法2我参考了人家的（1.2节中给出url）改写的选取的坐标原点是左上角。这就是其中的区别。

    我又算了一下自己算法1的那张放大的图和matlab提供的imresize()函数做出来的图的距离，21.2307。（以上距离讲的都是欧式距离）算法2的放大2倍的图和matlab内部提供的函数做出来的图的距离，13.0983。嗯，选用左上角比选用左下角为坐标原点，要效果好。

    可以看出自己做的和matlab内部的虽然看上去差不多，但是好像有点距离。我在程序的注释中写到，此算法有缺陷，那么下一篇博文我们就学习一下这个缺陷存在什么地方，以及如何改进。届时我们再来看看和matlab内部的imresize()的距离差多少，现在我也不知道结果。

    PS: 此博文允许任何人转载