双线性插值

本文转自：http://www.cnblogs.com/yssongest/p/5303151.html

1 原理

在图像的仿射变换中，很多地方需要用到插值运算，常见的插值运算包括最邻近插值，双线性插值，双三次插值，兰索思插值等方法，其中，双线性插值由于折中的插值效果和运算速度，运用比较广泛。

越是简单的模型越适合用来举例子，假如 3x3 的256级灰度图的象素矩阵如下图所示（这个原始图把它叫做源图，Source）：
23 38 22
67 44 12
89 65 63

这个矩阵中，元素坐标(x, y)是这样确定的，x 从左到右，从 0 开始，y 从上到下，也是从 0 开始，这是图象处理中最常用的坐标系。

如果想把这副图放大为 4x4 大小的图像，那么第一步先把4x4的矩阵先画出来，如下所示，矩阵的每个像素都是未知数，等待着我们去填充（这个将要被填充的图的叫做目标图，Destination）：
? ? ? ?
? ? ? ?
? ? ? ?
? ? ? ?

然后往这个空的矩阵里面填值，要填的值从哪里来？从源图中。先填写目标图最左上角的象素，坐标为 (0,0)，那么该坐标对应源图中的坐标可以由如下公式得出：
srcX=dstX*(srcWidth/dstWidth)
srcY=dstY*(srcHeight/dstHeight)

套用公式，就可以找到对应的原图的坐标了(0x(3/4), 0x(3/4))→(0x0.75, 0x0.75)→ (0, 0)，找到了源图的对应坐标，就可以把源图中坐标为 (0, 0) 处的23象素值填进去目标图的 (0, 0) 这个位置了。

接下来，寻找目标图中坐标为(1,0)的象素对应源图中的坐标，套用公式：(1x0.75, 0x0.75)→(0.75, 0) ，结果发现，得到的坐标里有小数，计算机里的图像是数字图像，象素是最小单位，象素的坐标都是整数。这时候采用的一种策略就是采用四舍五入的方法（也可以采用直接舍掉小数位的方法），把非整数坐标转换成整数。那么按照四舍五入的方法就得到坐标 (1, 0)，把源图中坐标为 (1, 0) 处的像素值38填入目标图的坐标 (1, 0)。

依次填完每个象素，一幅放大后的图像就诞生了，像素矩阵如下所示：
23 38 22 22
67 44 12 12
89 65 63 63
89 65 63 63

这种放大图像的方法叫做最临近插值算法，这是一种最基本、最简单的图像缩放算法，效果也是最不好的，放大后的图像有很严重的马赛克，缩小后的图像有很严重的失真。

双线型内插值算法就是一种比较好的图像缩放算法，它充分的利用了源图中虚拟点四周的四个真实存在的像素值来共同决定目标图中的一个像素值，因此缩放效果比简单的最邻近插值要好很多。

双线性内插值算法描述如下:
对于一个目的像素，设置坐标通过反向变换得到的浮点坐标为(i+u, j+v)，其中i、j为浮点坐标的整数部分，u、v为浮点坐标的小数部分，是取值 [0, 1) 区间的浮点数。则这个像素的值 f(i+u, j+v) 可由原图像中坐标为 (i, j)、(i+1, j)、(i, j+1)、(i+1, j+1) 所对应的周围四个像素的值决定，即：
　　f(i+u, j+v) = (1-u)(1-v)f(i, j) + (1-u)vf(i, j+1) + u(1-v)f(i+1, j) + uvf(i+1, j+1)
其中， f(i, j) 表示源图像 (i, j) 处的的像素值，以此类推。

现在假如目标图的象素坐标为 (1, 1)，那么反推得到的对应于源图的坐标是 (0.75, 0.75)，这其实只是一个概念上的虚拟象素，实际在源图中并不存在这样一个象素，那么目标图的象素 (1, 1) 的取值不能够由这个虚拟象素来决定，而只能由源图的这四个象素共同决定：(0, 0)、(0, 1)、(1, 0)、(1, 1)，而由于 (0.75, 0.75) 离 (1, 1) 要更近一些，那么 (1, 1) 所起的决定作用更大一些，这从公式1中的系数 uv=0.75x0.75 就可以体现出来，而 (0.75, 0.75) 离 (0, 0) 最远，所以 (0,0) 所起的决定作用就要小一些，公式中系数为 (1-u)(1-v)=0.25x0.25 也体现出了这一特点。

2 计算方法

首先，在X方向上进行两次线性插值计算，然后在Y方向上进行一次插值计算。
　

在图像处理的时候，我们先根据
　　srcX = dstX* (srcWidth/dstWidth) ,
　　srcY = dstY * (srcHeight/dstHeight)
来计算目标像素在源图像中的位置，这里计算的 srcX 和 srcY 一般都是浮点数，比如f(1.2, 3.4)这个像素点是虚拟存在的，先找到与它临近的四个实际存在的像素点：
　　(1,3) (2,3)
　　(1,4) (2,4)
写成 f(i+u, j+v) 的形式，则 u=0.2, v=0.4, i=1, j=3
在沿着X方向差插值时，f(R1)=u(f(Q21)-f(Q11))+f(Q11)
沿着Y方向同理计算。

或者直接整理一步计算：
f(i+u, j+v) = (1-u)(1-v)f(i, j) + (1-u)vf(i, j+1) + u(1-v)f(i+1, j) + uvf(i+1, j+1) 。

3 加速以及优化策略

单纯按照上文实现的插值算法只能勉强完成插值的功能，速度和效果都不会理想，在具体代码实现的时候有些小技巧。参考OpenCV源码以及网上博客整理如下两点：

1. 源图像和目标图像几何中心的对齐。
2. 将浮点运算转换成整数运算。

3.1 源图像和目标图像几何中心的对齐

在计算源图像的虚拟浮点坐标的时候，
一般情况：
　　srcX = dstX* (srcWidth/dstWidth) ,
　　srcY = dstY * (srcHeight/dstHeight).
中心对齐：
　　srcX = (dstX+0.5)* (srcWidth/dstWidth) -0.5
　　srcY = (dstY+0.5) * (srcHeight/dstHeight)-0.5
　　
将公式变形，srcX=dstX* (srcWidth/dstWidth)+0.5*(srcWidth/dstWidth-1)

相当于我们在原始的浮点坐标上加上了 0.5*(srcWidth/dstWidth-1) 这样一个控制因子，这项的符号可正可负，与 srcWidth/dstWidth 的比值，也就是当前插值是扩大还是缩小图像有关。

看一个例子：假设源图像是 3x3，中心点坐标是 (1, 1)，目标图像是 9x9，中心点坐标是 (4, 4)，在进行插值映射时，尽可能希望均匀的用到源图像的像素信息，最直观的就是 (4, 4) 映射到 (1, 1) 现在直接计算 srcX=4x3/9=1.3333，也就是我们在插值的时候所利用的像素集中在图像的右下方。现在考虑中心点对齐，srcX=(4+0.5)x3/9-0.5=1，刚好满足我们的要求。

3.2 将浮点运算转换成整数运算

参考图像处理界双线性插值算法的优化。

直接进行计算的话，由于计算的srcX和srcY 都是浮点数，后续会进行大量的乘法，而图像数据量又大，速度不会理想，解决思路是：浮点运算→整数运算→"<<左右移按位运算"。

放大的主要对象是 u，v 这些浮点数，OpenCV选择的放大倍数是 2048。
如何取这个合适的放大倍数呢，要从三个方面考虑。
第一，精度问题，如果这个数取得过小，那么经过计算后可能会导致结果出现较大的误差。
第二，这个数不能太大，太大会导致计算过程超过长整形所能表达的范围。
第三，速度考虑。假如放大倍数取为 12，那么算式在最后的结果中应该需要除以 12x12=144，但是如果取为 16，则最后的除数为 16x16=256，这个数字好，我们可以用右移来实现，而右移要比普通的整除快多了。我们利用左移11位操作就可以达到放大目的。

4 代码

参考：http://www.cnblogs.com/sdxk/p/4056223.html

#include<iostream>
#include<opencv2/core/core.hpp>
#include<opencv2/highgui/highgui.hpp>
using namespace cv;
using namespace std;

void scale(Mat &srcmat, Mat &desmat, double sx, double sy)
{
	int nc = x, nl = y,srccol=0,srcrow=0;
	double alph = 0.0, beta = 0.0;
	for (int i = 0; i < nc; i++)
	{
		uchar* desdata = desmat.ptr<uchar>(i);
		for (int j = 0; j < nl; j++)
		{
			srcrow = int(i / sx);
			//下面的的几个if是判断放大后图像对应到原图像的坐标是否越界的
			if (srcrow >= srcmat.rows-1)
				srcrow = srcmat.rows - 2;
			alph = i / sx - srcrow;
			if (alph >= 1)
				alph = 1;
			srccol = int(j / sy);
			if (srccol >= srcmat.cols-1)
				srccol = srcmat.cols-2;
			beta = j / sy - srccol;
			if (beta >= 1)
				beta=1;
			for (int k = 0; k < 3; k++)
			{
				double kk = srcmat.at<Vec3b>(srcrow, srccol)[k] + beta*(srcmat.at<Vec3b>(srcrow, srccol+1)[k] - srcmat.at<Vec3b>(srcrow, srccol)[k]);
				double jj = srcmat.at<Vec3b>(srcrow+1, srccol)[k] + beta*(srcmat.at<Vec3b>(srcrow+1, srccol + 1)[k] - srcmat.at<Vec3b>(srcrow+1, srccol)[k]);
				desdata[j*3+k] = kk + alph*(jj - kk);
			}
		}
	}
}
int main(int argc ,char** argv)
{
    // Read the file
    Mat srcimg = imread("opencv.jpg", IMREAD_COLOR);
    // Check for invalid input
    if (!srcimg.data) 
    {
        cout << "Could not open or find the image" << std::endl;
        return -1;
    }
    imshow("source img", srcimg);
    
    double sx = 0, sy = 0;
    cout << "please input scale x and scale y" << endl;
    cin >> sx >> sy;
        
    int x, y;
    x = int(srcimg.rows*sx);
    y = int(srcimg.cols*sy);
    
    Mat scaimg(x, y, CV_8UC3, Scalar::all(0));
    scale(srcimg, scaimg, sx, sy);
    imshow("scaled img", scaimg);
    
    waitKey();
    return 0;
}