关于双线性插值

最近工作涉及到了插值算法，查了一下发现原来这东西如果抠细节的话并不简单，于是记录一下。这篇先记录的是双线性插值。

应用场景和基本原理

简单来说就是，在对图像进行缩放的时候，对于未知像素可以采用插值进行求解。这篇文章写的比较清楚，这里就不再重复了。
双线性插值的原理可以描述为，将output image的点映射回input image的坐标系（映射回去的坐标可能为小数），并找出最近的四个在input image中像素点（它们的坐标一定为整数），然后加权平均得到它的像素值。看图的话会比较清晰：

乍一看这个算法很直观，也无需专门写一篇文章来记录。但是这个直观的算法在实现的时候有一点是值得注意的，就是将output映射为input坐标系时候具体要怎么映射。查了一些文章，也有各种的实现方法，比如：

http://tech-algorithm.com/articles/bilinear-image-scaling/

http://stackoverflow.com/questions/26142288/resize-an-image-with-bilinear-interpolation-without-imresize

http://blog.sina.com.cn/s/blog_ab584cac0101h0xy.html

感觉它们的实现方法并不一致，也不知道谁的实现方法比较靠谱。所以笔者决定看matlab的imresize是怎么实现的。

在matlab中的实现

在matlab中输入：

1
2
3

% 生成一张简单的图片，并用双线性插值法放大到2倍
im = uint8([8 16;16 32]);
imresize(im, 2, 'bilinear')

返回结果会是：

ans =

    8   10   14   16
   10   13   18   20
   14   18   25   28
   16   20   28   32

现在分析这个结果是怎么得到的。所幸matlab是可以通过ctrl+D找到部分的源码实现方法的，于是笔者大致读了一下，结合网上一些资料终于把如何映射这个问题给解决了。
最主要的函数是imresize.m中的这个函数：

function [weights, indices] = contributions(in_length, out_length, ...
                                            scale, kernel, ...
                                            kernel_width, antialiasing)


if (scale < 1) && (antialiasing)
    % Use a modified kernel to simultaneously interpolate and
    % antialias.
    h = @(x) scale * kernel(scale * x);
    kernel_width = kernel_width / scale;
else
    % No antialiasing; use unmodified kernel.
    h = kernel;
end

% Output-space coordinates.
x = (1:out_length)';

% Input-space coordinates. Calculate the inverse mapping such that 0.5
% in output space maps to 0.5 in input space, and 0.5+scale in output
% space maps to 1.5 in input space.
u = x/scale + 0.5 * (1 - 1/scale);

% What is the left-most pixel that can be involved in the computation?
left = floor(u - kernel_width/2);

% What is the maximum number of pixels that can be involved in the
% computation?  Note: it's OK to use an extra pixel here; if the
% corresponding weights are all zero, it will be eliminated at the end
% of this function.
P = ceil(kernel_width) + 2;

% The indices of the input pixels involved in computing the k-th output
% pixel are in row k of the indices matrix.
indices = bsxfun(@plus, left, 0:P-1);

% The weights used to compute the k-th output pixel are in row k of the
% weights matrix.
weights = h(bsxfun(@minus, u, indices));

% Normalize the weights matrix so that each row sums to 1.
weights = bsxfun(@rdivide, weights, sum(weights, 2));

% Clamp out-of-range indices; has the effect of replicating end-points.
indices = min(max(1, indices), in_length);

% If a column in weights is all zero, get rid of it.
kill = find(~any(weights, 1));
if ~isempty(kill)
    weights(:,kill) = [];
    indices(:,kill) = [];
end

matlab的注释写得很清楚了，output的坐标用x变量表示，对应映射回input的坐标用u表示，当x=0.5时候，u也应该为0.5;当x=0.5+scale时，u为1.5。这样会得到怎么的映射效果呢，笔者画了一幅示意图：

这里的坐标系是matlab下的坐标系（左上角为原点，从1开始计数），其中红色的4个点表示input的四个点，灰色的16个点是output的点映射回到input坐标系下的坐标。图片下方有两行x和u分别是灰色点在output和input下的x方向上的坐标值。
所以根据上面的映射方法，得到output映射回input的图像特点有两个：

图片的中心是重合的（图中央的黑点为中心）
output中的相邻像素的距离为input相邻像素的距离的1/scale

另外，代码中计算weight时候用到的h函数，在bilinear中是长这个样子的（下图的左边部分，图片来自这里）：

即一个等腰三角形，代表的意思是距离某个参考点距离为1的时候，权重递减，超过1就没有权重了。

有了上面的铺垫，要计算灰色的像素值就简单了，举几个例子。
对于第一行第一个灰色点，坐标是(0.75，0.75)，(0.75±1，0.75±1)范围的点只有左上角的红点，于是它的值就取该红点的像素值。
对于第二行第一个灰色点，坐标是(0.75, 1.25)，(0.75±1, 1.25±1)范围的点有两个，即左上和左下的红点，所以根据距离权重计算为左上的点较近占3/4，左下的点距离较远占1/4，那么最后像素的值就是8x3/4+16x1/4=10
对于第二行第二个灰色点，坐标是(1.25, 1.25), (1.25±1, 1.25±1)范围的点包括了全部四个红点，于是也是根据距离来计算权值，归一化话，左上右上左下右下的点的权值分别为(27,9,9,3)/48，根据这个权重进行加权，结果为(8x27+16x9x2+32x3)/48=12.5，四舍五入取13。

对应上述matlab的结果，说明这样进行理解应该是没有问题的了。

参考页面：

在路上

关于双线性插值

应用场景和基本原理

在matlab中的实现