图像处理-内插法
一、插值
1.什么是插值
插值是数值分析中的一个重要内容。由于在计算机中数据都是离散的，所以对于连续数据是没有办法直接处理的，例如：在图像的缩放以及旋转操作中，由于在现实世界数据是可以连续的所以我们队图像旋转这种操作是非常简单的。但是在计算机中执行这种操作后，新图像的某些像素点会有缺失，所以利用插值，将这些缺失的点进行填充。
我们举一个简明的例子来说明插值的原理。例如我们有以下数据
我们希望拟合出y=f(x)的函数图像，但是我的得到的数据都是离散的，那么我们应该如何确定两个点之间缺失的函数值呢？
2.Lagrang插值
插值的方法有很多，这里我们介绍Lagrange(拉格朗日)插值。假设我们拟合出的函数为p(x)首先我们应该明确我们拟合出的函数至少在我们的已知点上的数据是正确的，即:
那么我们需要构造一个多项式 其满足
经过求解可以得到一个满足上述条件的多项式(具体的求解方法这里不写了)：
不难看出当且仅当时
才等于1，其余均为0。所以我们得到一个插值多项式：
完成了数据的拟合。
举一具体的例子：
已知两点的坐标,
,利用拉格朗日插值法，拟合的函数为：
3.图像旋转中的插值
3.1 图像是如何旋转的
在计算机中图片是以矩阵或者数组的方式存放的，对图片进行任意角度的旋转看似容易但是需要有一些问题需要解决。
图像旋转的方法有大量的图片和公式，偷个懒直接从网上下了讲解的图片。
经过以上的公式运算我们就可以将图像旋转任意角度，但是运行的结果是否真的跟我们预期的一样呢？不妨写一段程序实现上述的公式，进行图像旋转：
先看一下原图
可以看到图片上是有许多像素点是没有灰度值的，那么造成这种结果的原因是什么呢？首先我们看以上的公式中含有很多三角函数，那么势必会带来许多小数，在原图的坐标向旋转图映射坐标时，通过计算会发现映射的坐标有的是小数，但是计算机中的图像矩阵坐标是没有小数的！所以需要对这些小数进行四舍五入，这样有的像素点就没有像素值了！所以我们需要进行插值。
3.2 最近邻插值法
最近邻插值法顾名思义，就是把距离像素缺失点最近的像素点的值填充进去。这种方法实现起来也非常简单，之前我们都是讲原图的坐标映射到旋转图中，现在反过来我们将旋转图的坐标映射回原图，再将原图对应点的像素复制过来，这样旋转图中的没个点都对应一个原图中的点，由于计算中会产生小数，所以有可能旋转图中的某些点都对应同一个原图中的点。我们来看一下程序运行效果。
怎么样效果是不是要提升很多呢？但是与原图相比清晰度似乎还是差了那么一点。
3.3 双线性插值法
这次需要用到我们用拉格朗日插值法的结论了。为什么我们叫线性插值？可以看到我们讲拉格朗日插值法的时候举的例子中只有两个已知点，最后拟合的函数是一个一次线性函数函数。那么双线性插值法的思想就是在图像的水平和竖直方向都做线性插值。理论的推导已经给出，所以我就偷个懒~~
假如我们想得到未知函数 在点
的值，假设我们已知函数
在
四个点的值。最常见的情况，f就是一个像素点的像素值。首先在 x 方向进行线性插值，得到
然后在 y 方向进行线性插值，得到
综合起来就是双线性插值最后的结果：
我们再来看一下程序运行结果：
对比最近邻插值，清晰度是不是又有了很大的提升呢？反正鄙人用肉眼是分辨不出与原图的差距的。
总结
其实要解释双线性插值法是没必要用拉格朗日插值法来解释的，但是图像旋转算法中还有三次插值甚至更高次的插值，如果能够理解拉格朗日插值其他高次插值也迎刃而解。例如如果我们利用一个点周围更多的像素已知的点，利用拉格朗日插值可以拟合出次数更高的函数，自然效果会更好。但是需要注意的是插值次数越高，需要利用的已知点的数量越多，程序运行越慢。
matlab代码
clear;
clc;
close all;
% read a image
img = imread('lena.jpg');
imshow(img);
title('original image');
imshow(img)
% rotation degree
degree = 60;
hh=imrotate(img,45,'bilinear');
imshow(hh)
% image size after rotation
[m, n, o] = size(img);
new_m = ceil(abs(m*cosd(degree)) + abs(n*sind(degree)));
new_n = ceil(abs(n*cosd(degree)) + abs(m*sind(degree)));
new_img_forward = zeros(new_m, new_n, o); % forward mapping
new_img_nnp = zeros(new_m, new_n, o); % reverse mapping, nearest neighbor interpolation
new_img_lp = zeros(new_m, new_n, o); % reverse mapping, bilinear interpolation
%% forward mapping
% forward mapping matrices
m1 = [1 0 0; 0 -1 0; -0.5*n 0.5*m 1];
m2 = [cosd(degree) -sind(degree) 0; sind(degree) cosd(degree) 0; 0 0 1];
m3 = [1 0 0; 0 -1 0; 0.5*new_n 0.5*new_m 1];
t=m1*m2*m3;
for i = 1:n
for j = 1:m
new_coordinate = [i j 1]*t;
col = round(new_coordinate(1));
row = round(new_coordinate(2));
new_img_forward(row, col, 1) = img(j, i, 1);
new_img_forward(row, col, 2) = img(j, i, 2);
new_img_forward(row, col, 3) = img(j, i, 3);
end
end
figure, imshow(new_img_forward/255), title('forward mapping');
%% reverse mapping
% reverse mapping matrices
rm1 = [1 0 0; 0 -1 0; -0.5*new_n 0.5*new_m 1];
rm2 = [cosd(degree) sind(degree) 0; -sind(degree) cosd(degree) 0; 0 0 1];
rm3 = [1 0 0; 0 -1 0; 0.5*n 0.5*m 1];
bbb=rm1*rm2*rm3
for i = 1:new_n
for j = 1:new_m
% rotated image's coordinates to no-rotation image's coordinates
old_coordinate = [i j 1]*bbb;
col = round(old_coordinate(1));
col_t=old_coordinate(1);
row = round(old_coordinate(2));
row_t=old_coordinate(2);
% prevent border overflow
if row < 1 || col < 1 || row > m || col > n
new_img_nnp(j, i) = 0;
new_img_lp(j, i) = 0;
else
% nearest neighbor interpolation
new_img_nnp(j, i, 1) = img(row, col, 1);
new_img_nnp(j, i, 2) = img(row, col, 2);
new_img_nnp(j, i, 3) = img(row, col, 3);
% bilinear interpolation
left = floor(col_t);
right = ceil(col_t);
top = floor(row_t);
bottom = ceil(row_t);
if bottom>m-1 || right>n-1 || left<1 ||top<1
new_img_nnp(j, i) = 0;
new_img_lp(j, i) = 0;
else
a = col_t - left;
b = row_t - top;
new_img_lp(j, i, :) = (1-a)*(1-b)*img(top, left, :) + a*(1-b)*img(top, right, :) + ...
(1-a)*(b)*img(bottom, left, :) + (a)*(b)*img(bottom, right, :);
end
end
end
end
figure, imshow(new_img_nnp/255), title('nearest neighbor interpolation');
figure, imshow(new_img_lp/255), title('bilinear interpolation');
