双线性插值法的MATLAB实现

最近做一个Fortran的数值模拟，其中有一步需要将图片进行放大或缩小，因此得自己在Fortran上实现图像伸缩的算法。一查，最常用且复杂度和效果都能接受的便是双线性插值法。然而大多数该算法的介绍，要么是照抄维基百科上的公式，要么就是说一大堆眼花缭乱的解释，公式推导写得十分繁琐，而一看代码就发现没多少干货。最后自己推了坐标变换，画了些图才确定代码该怎么写，并在MATLAB上验证后，才迁移到Fortran上。

对于图片缩放来说，首先我们具有一张原始图片，我们称为src(source)；然后我们需要缩放该图片，得到最终的图片，我们称为dst(destination)。在计算机图像中，彩色图片的数据形式是一个的矩阵。在本文中，为了简化讨论，我们只讨论灰度图，即数据形式为矩阵。

那么，图片缩放就转换为矩阵映射问题：为矩阵dst中的每个矩阵元的数值，找到其与矩阵src中矩阵元的映射关系。我们来考虑一个极端简单的例子，一张图片，将其放大为原图片的 1 倍（也就是不变啦哈哈哈）。那么此时，矩阵dst与矩阵src的矩阵元之间可以建立一个一一对应的关系，即dst(i,j)=src(i,j)。

那如果放大为原图片的 2 倍呢？如果我们也能直接通过dst(i,j)=src(i/2,j/2)来得到dst(i,j)就好了。然而矩阵元的位置只能是正整数，而i/2和j/2很可能得到一个实数，我们因此无法直接得到dst(i,j)的数值。但是，矩阵坐标(i/2,j/2)如果是非整数，则会落在四个src矩阵元之间，如图1所示。图1中的5，就是我们通过坐标缩放得到的“虚”点，这个点并不存在于矩阵src上，但我们可以通过双线性插值法，用图1中的点（1,2,3,4）来得到这个虚点5的值。

图1：双线性插值示意图【1】

设图1中的1，2，3，4这四个点的像素值分别为：。首先通过点1和点2，用线性插值得到m处的值；通过点3和点4，同样用线性插值得到n处的值。在用n处和m处的值，进行线性插值，得到“虚”点5的值。这就是双线性插值的思路，公式推导如下：

又有：，，，

故结合上面式子与式得：

以上式子中，指矩阵元的列坐标，指矩阵元的行坐标。而就是我们得到的dst(i,j)在矩阵src中的映射值了。但我们能不能直接通过i*scale和j*scale的方式来得到dst(i,j)在src中“虚”点的行列坐标呢？在参考文章【1】中说，这样子会导致原图片的最后一列数据不参与计算。因此要“几何中心对齐”，公式是这样的：。

我对此表示存疑，因为scale很小时，会产生负数的“虚”点矩阵坐标。从我的角度来说，我只希望伸缩后的图片，其中心能与原图片重合，并且算虚点坐标时不会在旁边找不出4个点。所以我的做法是，用坐标转移矩阵的方式，先将矩阵坐标转换到我们熟悉的直角坐标，如下图：

图2：矩阵的直角坐标

矩阵这样一看不就是直角坐标的右下半区嘛，我让每个矩阵元的值所在的直角坐标记为图2中小圆点的直角坐标。然后建立矩阵元直角坐标与矩阵行列坐标的转换矩阵如下。注意！图2只是为了示例坐标转换，切勿把图2当做dst矩阵。图2和图3都是src矩阵。

将dst矩阵坐标(i,j)转换为直角坐标系下坐标(x0,y0)。将这个直角坐标进行缩放，使得其缩放后的区域被最外围表示矩阵元值的圆点所组成的矩形框所框住，这样任意一个“虚”点周围都肯定有4个点可以进行双线性插值。如图3中红色框。但是缩放scale*(x0,y0)会使得缩放后的区域飘到红色框的左上角原点处。所以还需要x轴+0.5，y轴-0.5，才能把缩放后的区域挪到红色框内，所以我们得到伸缩平移矩阵如下：

图3

要强调的是，，。，以确保缩放后能被框住，且不在框上。按道理，应该取，但这么取的话，边缘会跟红色边框重合，需要在程序里判断是否超出边界，且边缘会有黑边。我为了省事，就取的。

最后我们将dst矩阵映射到src矩阵上“虚”点的直角坐标转换回矩阵坐标系就可以了：

我们就得到最终虚点在src矩阵的行列坐标(i1,j1)了。

我接下来给出用MATLAB实现的代码，并给出结果示例。

%% Lenna biger or smaller
clear;clc;
Lenna=imread('Lenna.jpg');

[r00,c00]=size(Lenna);
c00=c00/3; %因为有3个色彩通道，所以除3

r=256; %伸缩后的图片像素
c=512;

NewLenna= zeros(r,c);

scaler=(r00-1)/r;
scalec=(c00-1)/c;

M1=[0 1 -0.5;-1 0 0.5;0 0 1];
M3=[0 -1 0.5;1 0 0.5;0 0 1];
M2=[scalec 0 0.5;0 scaler -0.5;0 0 1];

for ch=1:3
    disp(ch)
    for i=1:r
    %disp(i);
    for j=1:c
        temp=M3*M2*M1*[i j 1]';
        srci=temp(1,1);
        srcj=temp(2,1);
        i0=floor(temp(1,1));
        j0=floor(temp(2,1));
        i1=i0+1;
        j1=j0+1;
        value0 = (j1 - srcj)*Lenna(i0,j0,ch) + (srcj - j0) * Lenna(i0, j1, ch);
        value1 = (j1 - srcj)*Lenna(i1,j0,ch) + (srcj - j0) * Lenna(i1, j1, ch);
        NewLenna(i, j, ch) = floor((i1 - srci) * value0 + (srci - i0) * value1);
    end
    end
end

figure()
imshow(NewLenna./256); 

根据上面的MATLAB代码，我们来试验一下著名的Lenna图。
原图如下：

Lenna316_316.jpg

Lenna256_512.jpg

Lenna256_512黑边.jpg

参考文章：
【1】：《彻底搞懂双线性插值》，胡今朝，知乎