雙線性内插值算法

引出

  在圖像的仿射變換中，很多地方需要用到插值運算，常見的插值運算包括最鄰近插值，雙線性插值，雙三次插值，蘭索思插值等方法，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);