双线性内插值算法

Posted 2021-08-30 qlexcel

引出

  在图像的仿射变换中，很多地方需要用到插值运算，常见的插值运算包括最邻近插值，双线性插值，双三次插值，兰索思插值等方法，OpenCV提供了很多方法，其中，双线性插值由于折中的插值效果和运算速度，运用比较广泛。
  越是简单的模型越适合用来举例子，我们就举个简单的图像：3*3 的256级灰度图。假如图像的象素矩阵如下图所示（这个原始图把它叫做源图，Source）：

234 38  22
67  44  12
89  65  63

  这个矩阵中，元素坐标(x,y)是这样确定的，x从左到右，从0开始，y从上到下，也是从零开始，这是图象处理中最常用的坐标系。
  如果想把这副图放大为 4*4大小的图像，那么该怎么做呢？那么第一步肯定想到的是先把4*4的矩阵先画出来再说，好了矩阵画出来了，如下所示，当然，矩阵的每个像素都是未知数，等待着我们去填充（这个将要被填充的图的叫做目标图,Destination）：

? ? ? ?
? ? ? ?
? ? ? ?
? ? ? ?

  然后要往这个空的矩阵里面填值了，要填的值从哪里来来呢？是从源图中来，好，先填写目标图最左上角的象素，坐标为（0，0），那么该坐标对应源图中的坐标可以由如下公式得出srcX=dstX* (srcWidth/dstWidth) , srcY = dstY * (srcHeight/dstHeight)
  好了，套用公式，就可以找到对应的原图的坐标了(0*(3/4),0*(3/4))=>(0*0.75,0*0.75)=>(0,0)，找到了源图的对应坐标,就可以把源图中坐标为(0,0)处的234象素值填进去目标图的(0,0)这个位置了。
  接下来,如法炮制,寻找目标图中坐标为(1,0)的象素对应源图中的坐标,套用公式:
(1*0.75,0*0.75)=>(0.75,0) 结果发现,得到的坐标里面竟然有小数,这可怎么办?计算机里的图像可是数字图像,象素就是最小单位了,象素的坐标都是整数,从来没有小数坐标。这时候采用的一种策略就是采用四舍五入的方法（也可以采用直接舍掉小数位的方法），把非整数坐标转换成整数，好，那么按照四舍五入的方法就得到坐标（1，0），完整的运算过程就是这样的：(1*0.75,0*0.75)=>(0.75,0)=>(1,0) 那么就可以再填一个象素到目标矩阵中了，同样是把源图中坐标为(1,0)处的像素值38填入目标图中的坐标。
  依次填完每个象素，一幅放大后的图像就诞生了，像素矩阵如下所示：

234 38  22  22
67  44  12  12
89  65  63  63
89  65  63  63

  这种放大图像的方法叫做最临近插值算法，这是一种最基本、最简单的图像缩放算法，效果也是最不好的，放大后的图像有很严重的马赛克，缩小后的图像有很严重的失真；效果不好的根源就是其简单的最临近插值方法引入了严重的图像失真，比如，当由目标图的坐标反推得到的源图的的坐标是一个浮点数的时候，采用了四舍五入的方法，直接采用了和这个浮点数最接近的象素的值，这种方法是很不科学的，当推得坐标值为 0.75的时候，不应该就简单的取为1，既然是0.75，比1要小0.25 ，比0要大0.75 ,那么目标象素值其实应该根据这个源图中虚拟的点四周的四个真实的点来按照一定的规律计算出来的，这样才能达到更好的缩放效果。
  双线型内插值算法就是一种比较好的图像缩放算法，它充分的利用了源图中虚拟点四周的四个真实存在的像素值来共同决定目标图中的一个像素值，因此缩放效果比简单的最邻近插值要好很多。

双线性内插值算法描述

  对于一个目的像素，设置坐标通过反向变换得到的浮点坐标为(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.75×0.75就可以体现出来，而（0.75,0.75）离（0，0）最远，所以（0，0）所起的决定作用就要小一些，公式中系数为(1-u)(1-v)=0.25×0.25也体现出了这一特点。

计算方法

首先，在X方向上进行两次线性插值计算

然后在Y方向上进行一次插值计算：

在图像处理的时候，我们先根据

srcX = dstX * (srcWidth/dstWidth),
srcY = dstY * (srcHeight/dstHeight)

来计算目标像素在源图像中的位置，这里计算的srcX和srcY一般都是浮点数，比如f（1.2, 3.2）这个像素点是虚拟存在的，先找到与它临近的四个实际存在的像素点：

（1，4） （2，4）
（1，3） （2，3）

写成f(i+u,j+v)的形式，则u=0.2,v=0.2, i=1, j=3
假设他们的值为：

100 200
300 400

根据上面公式：

x1=1  x=1.2  x2=2
y1=3  y=3.2  y2=4
f(Q12)=100  f(Q22)=200
f(Q11)=300  f(Q21)=400

f(R1)=( (2-1.2)/(2-1) )*300 + ( (1.2-1)/(2-1) )*400 = 0.8*300 + 0.2*400 = 320
f(R2)=( (2-1.2)/(2-1) )*100 + ( (1.2-1)/(2-1) )*200 = 0.8*100 + 0.2*200 = 120

f(P)=( (4-3.2)/(4-3) )*f(R1) + ( (3.2-3)/(4-3) )*f(R2) = 0.8*320 + 0.2*120 = 280

也可以一步直接计算：
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)

加速以及优化策略

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

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

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

在计算源图像的虚拟浮点坐标的时候，一般情况：

srcX = dstX * (srcWidth/dstWidth) ,
srcY = dstY * (srcHeight/dstHeight)

中心对齐(OpenCV也是如此)：

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的比值也就是当前插值是扩大还是缩小图像有关，有什么作用呢？看一个例子：
  假设源图像是3*3，中心点坐标（1，1），目标图像是9*9，中心点坐标（4，4），我们在进行插值映射的时候，尽可能希望均匀的用到源图像的像素信息，最直观的就是（4,4）映射到（1,1），现在直接计算srcX=4*(3/9)=1.3333 != 1，也就是我们在插值的时候所利用的像素集中在图像的右下方，而不是均匀分布整个图像。现在考虑中心点对齐，srcX=(4+0.5)*3/9-0.5=1，刚好满足我们的要求。

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

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

  放大的主要对象是u，v这些浮点数，OpenCV选择的放大倍数是2048“如何取这个合适的放大倍数呢，要从三个方面考虑：

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

opencv代码

cv::Mat matSrc, matDst1, matDst2;  
 
matSrc = cv::imread("lena.jpg", 2 | 4);  
matDst1 = cv::Mat(cv::Size(800, 1000), matSrc.type(), cv::Scalar::all(0));  
matDst2 = cv::Mat(matDst1.size(), matSrc.type(), cv::Scalar::all(0));  
 
double scale_x = (double)matSrc.cols / matDst1.cols;  
double scale_y = (double)matSrc.rows / matDst1.rows;  
 
uchar* dataDst = matDst1.data;  
int stepDst = matDst1.step;  
uchar* dataSrc = matSrc.data;  
int stepSrc = matSrc.step;  
int iWidthSrc = matSrc.cols;  
int iHiehgtSrc = matSrc.rows;  
 
for (int j = 0; j < matDst1.rows; ++j)  
{  
    float fy = (float)((j + 0.5) * scale_y - 0.5);  
    int sy = cvFloor(fy);  
    fy -= sy;  
    sy = std::min(sy, iHiehgtSrc - 2);  
    sy = std::max(0, sy);  
 
    short cbufy[2];  
    cbufy[0] = cv::saturate_cast<short>((1.f - fy) * 2048);  
    cbufy[1] = 2048 - cbufy[0];  
 
    for (int i = 0; i < matDst1.cols; ++i)  
    {  
        float fx = (float)((i + 0.5) * scale_x - 0.5);  
        int sx = cvFloor(fx);  
        fx -= sx;  
 
        if (sx < 0) {  
            fx = 0, sx = 0;  
        }  
        if (sx >= iWidthSrc - 1) {  
            fx = 0, sx = iWidthSrc - 2;  
        }  
 
        short cbufx[2];  
        cbufx[0] = cv::saturate_cast<short>((1.f - fx) * 2048);  
        cbufx[1] = 2048 - cbufx[0];  
 
        for (int k = 0; k < matSrc.channels(); ++k)  
        {  
            *(dataDst+ j*stepDst + 3*i + k) = (*(dataSrc + sy*stepSrc + 3*sx + k) * cbufx[0] * cbufy[0] +   
                *(dataSrc + (sy+1)*stepSrc + 3*sx + k) * cbufx[0] * cbufy[1] +   
                *(dataSrc + sy*stepSrc + 3*(sx+1) + k) * cbufx[1] * cbufy[0] +   
                *(dataSrc + (sy+1)*stepSrc + 3*(sx+1) + k) * cbufx[1] * cbufy[1]) >> 22;  
        }  
    }  
}  
cv::imwrite("linear_1.jpg", matDst1);  
 
cv::resize(matSrc, matDst2, matDst1.size(), 0, 0, 1);  
cv::imwrite("linear_2.jpg", matDst2);