線性插值

先講一下線性插值：已知資料 (x0, y0) 與 (x1, y1)，要計算 [x0, x1] 區間內某一位置 x 在直線上的y值（反過來也是一樣，略）：

y−y0x−x0=y1−y0x1−x0y=x1−xx1−x0y0+x−x0x1−x0y1

上面比較好理解吧，仔細看就是用x和x0，x1的距離作為一個權重，用於y0和y1的加權。雙線性插值本質上就是在兩個方向上做線性插值。

雙線性插值

在數學上，雙線性插值是有兩個變數的插值函式的線性插值擴充套件，其核心思想是在兩個方向分別進行一次線性插值[1]。見下圖：

假如我們想得到未知函式 f 在點 P = (x, y) 的值，假設我們已知函式 f 在 Q11 = (x1, y1)、Q12 = (x1, y2), Q21 = (x2, y1) 以及 Q22 = (x2, y2) 四個點的值。最常見的情況，f就是一個畫素點的畫素值

。首先在 x 方向進行線性插值，得到

然後在 y 方向進行線性插值，得到

綜合起來就是雙線性插值最後的結果：

由於影象雙線性插值只會用相鄰的4個點，因此上述公式的分母都是1。opencv中的原始碼如下，用了一些優化手段，比如用整數計算代替float（下面程式碼中的*2048就是變11位小數為整數，最後有兩個連乘，因此>>22位），以及源影象和目標影象幾何中心的對齊
SrcX=(dstX+0.5)* (srcWidth/dstWidth) -0.5
SrcY=(dstY+0.5) * (srcHeight/dstHeight)-0.5，

這個要重點說一下，源影象和目標影象的原點（0，0）均選擇左上角，然後根據插值公式計算目標影象每點畫素，假設你需要將一幅5x5的影象縮小成3x3，那麼源影象和目標影象各個畫素之間的對應關係如下。如果沒有這個中心對齊，根據基本公式去算，就會得到左邊這樣的結果；而用了對齊，就會得到右邊的結果：

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);  

好了，本篇到這裡，歡迎大家分享轉載，註明出處即可。

參考資料

=================================雙三次線性插值==================================

雙三次插值是一種更加複雜的插值方式，它能創造出比雙線性插值更平滑的影象邊緣。雙三次插值方法通常運用在一部分影象處理軟體、印表機驅動程式和數碼相機中，對原影象或原影象的某些區域進行放大。Adobe Photoshop CS 更為使用者提供了兩種不同的雙三次插值方法：雙三次插值平滑化和雙三次插值銳化。

在數值分析這個數學分支中，雙三次插值（英語：Bicubic interpolation）是二維空間中最常用的插值方法。在這種方法中，函式f在點 (x,y) 的值可以通過矩形網格中最近的十六個取樣點的加權平均得到，在這裡需要使用兩個多項式插值三次函式，每個方向使用一個。雙三次插值又叫雙立方插值，用於在影象中“插值”（Interpolating）或增加“畫素”（Pixel）數量/密度的一種方法。通常利用插值技術增加圖形資料，以便在它列印或其他形式輸出的時候，能夠增大列印面積以及（或者）解析度。目前有不同的插值技術可供選用。雙立方插值通常能產生效果最好，最精確的插補圖形，但它速度也幾乎是最慢的。“雙線性插值”（Bilinear interpolation）的速度則要快一些，但沒有前者精確。在商業性影象編輯軟體中，經常採用的是速度最快，但也是最不準確的“最近相鄰”（Nearest Neighbor）插值。其他一些插值技術通常只在高檔或單獨應用的程式中出現。顯然，無論技術多麼高階，插補過的資料肯定沒有原始資料準確。這意味著對一個圖形檔案進行插值處理後，雖然檔案長度增加了（資料量增大），但不會有原先那幅圖銳利，可能會在圖形質量上打折扣。雙三次插值通過下式進行計算：^[1]或者用一種更加緊湊的形式，計算係數的過程依賴於插值資料的特性。如果已知插值函式的導數，常用的方法就是使用四個頂點的高度以及每個頂點的三個導數。一階導數與 表示 x 與 y 方向的表面斜率，二階相互導數表示同時在 x 與 y 方向的斜率。這些值可以通過分別連續對 x 與 y 向量取微分得到。對於網格單元的每個頂點，將區域性座標（0,0, 1,0, 0,1 和 1,1) 帶入這些方程，再解這 16 個方程。

今天學習了第三種影象縮放的方法，雙三次插值法。由於理解能力比較差，看了好久的公式，還是雲裡霧裡，但是為了督促自己學習，還是把已知的部分記錄下來。

數學原理

假設源影象A大小為m*n，縮放後的目標影象B的大小為M*N。那麼根據比例我們可以得到B(X,Y)在A上的的
對應座標為A(x,y)=A(X*(m/M),Y*(n/N))。在雙線性插值法中，我們選取A(x,y)的最近四個點。而在雙立方
插值法中，我們選取的是最近的16個畫素點作為計算目標影象B(X,Y)處畫素值的引數。如圖所示：

如圖所示P點就是目標影象B在(X,Y)處對應於源影象中的位置，P的座標位置會出現小數部分，所以我們假設
P的座標為P(x+u,y+v)，其中x,y分別表示整數部分，u,v分別表示小數部分。那麼我們就可以得到如圖所示的
最近16個畫素的位置，在這裡用a(i,j)(i,j=0,1,2,3)來表示。
雙立方插值的目的就是通過找到一種關係，或者說係數，可以把這16個畫素對於P處畫素值得影響因子找出
來，從而根據這個影響因子來獲得目標影象對應點的畫素值，達到影象縮放的目的。
我在這次的學習中學習的是基於BiCubic基函式的雙三次插值法，BiCubic基函式形式如下：

我們要做的就是求出BiCubic函式中的引數x,從而獲得上面所說的16個畫素所對應的係數。在學習雙線性插
值法的時候，我們是把影象的行和列分開來理解的，那麼在這裡，我們也用這種方法描述如何求出a(i,j)對應
的係數k_ij。假設行係數為k_i,列係數為k_j。我們以a00位置為例：
首先，我們要求出當前畫素與P點的位置，比如a00距離P(x+u,y+v)的距離為(1+u,1+v)。
那麼我們可以得到：k_i_0=W(1+u),k_j_0=W(1+v).
同理我們可以得到所有行和列對應的係數：

k_i_0=W(1+u), k_i_1=W(u), k__i_2=W(1-u), k_i_3=W(2-u);
k_j_0=W(1+v), k_j_1=W(v), k_j_2=W(1-v), k_j_3=W(2-v);

這樣我們就分別得到了行和列方向上的係數。
由k_i_j=k_i*k_j我們就可以得到每個畫素a(i,j)對應的權值了。

最後通過求和公式可以得到目標圖片B(X,Y)對應的畫素值：
pixelB(X,Y)=pixelA(0,0)*k_0_0+pixelA(0,1)*k_0_1+…+pixelA(3,3)*k_3_3;
這裡其實就是個求和公式，由於不知道怎麼編輯公式，就這樣表達了。

程式實現

/**********************10-9*******************************
功能：雙三次插值縮放圖片
數學原理：假設原影象A的大小為m*n,新影象B的大小為M*N
如果我們要求B（X,Y）處的畫素值：
我們首先可以得到B（X,Y）在影象A中對應的位置（x,y）=(X*(m/M),Y*(N/n))
這個時候求得的x,y是小數值，我們可以通過這個小數值座標找到距離最近的16個畫素點，
利用所選擇的基函式，求出對應的每個畫素的權值，最終獲得pixelB（X,Y）
**********************************************************/

#include <opencv2\opencv.hpp>
#include <iostream>
#include <math.h>
using namespace std;
using namespace cv;

float a = -0.5;//BiCubic基函式
void getW_x(float w_x[4], float x);
void getW_y(float w_y[4], float y);

int main(){
    Mat image = imread("lena.jpg");//源影象

    float Row_B = image.rows*2;
    float Col_B = image.cols*2;


    Mat biggerImage(Row_B, Col_B, CV_8UC3);

    for (int i = 2; i < Row_B-4; i++){
        for (int j = 2; j < Col_B-4; j++){
            float x = i*(image.rows / Row_B);//放大後的影象的畫素位置相對於源影象的位置
            float y = j*(image.cols / Col_B);

            /*if (int(x) > 0 && int(x) < image.rows - 2 && int(y)>0 && int(y) < image.cols - 2){*/
                float w_x[4], w_y[4];//行列方向的加權係數
                getW_x(w_x, x);
                getW_y(w_y, y);

                Vec3f temp = { 0, 0, 0 };
                for (int s = 0; s <= 3; s++){
                    for (int t = 0; t <= 3; t++){
                        temp = temp + (Vec3f)(image.at<Vec3b>(int(x) + s - 1, int(y) + t - 1))*w_x[s] * w_y[t];
                    }
                }

                biggerImage.at<Vec3b>(i, j) = (Vec3b)temp;
            }
        }

    imshow("image", image);
    imshow("biggerImage", biggerImage);
    waitKey(0);

    return 0;
}
/*計算係數*/
void getW_x(float w_x[4],float x){
    int X = (int)x;//取整數部分
    float stemp_x[4];
    stemp_x[0] = 1 + (x - X);
    stemp_x[1] = x - X;
    stemp_x[2] = 1 - (x - X);
    stemp_x[3] = 2 - (x - X);

    w_x[0] = a*abs(stemp_x[0] * stemp_x[0] * stemp_x[0]) - 5 * a*stemp_x[0] * stemp_x[0] + 8 * a*abs(stemp_x[0]) - 4 * a;
    w_x[1] = (a + 2)*abs(stemp_x[1] * stemp_x[1] * stemp_x[1]) - (a + 3)*stemp_x[1] * stemp_x[1] + 1;
    w_x[2] = (a + 2)*abs(stemp_x[2] * stemp_x[2] * stemp_x[2]) - (a + 3)*stemp_x[2] * stemp_x[2] + 1;
    w_x[3] = a*abs(stemp_x[3] * stemp_x[3] * stemp_x[3]) - 5 * a*stemp_x[3] * stemp_x[3] + 8 * a*abs(stemp_x[3]) - 4 * a;
}
void getW_y(float w_y[4], float y){
    int Y = (int)y;
    float stemp_y[4];
    stemp_y[0] = 1.0 + (y - Y);
    stemp_y[1] = y - Y;
    stemp_y[2] = 1 - (y - Y);
    stemp_y[3] = 2 - (y - Y);

    w_y[0] = a*abs(stemp_y[0] * stemp_y[0] * stemp_y[0]) - 5 * a*stemp_y[0] * stemp_y[0] + 8 * a*abs(stemp_y[0]) - 4 * a;
    w_y[1] = (a + 2)*abs(stemp_y[1] * stemp_y[1] * stemp_y[1]) - (a + 3)*stemp_y[1] * stemp_y[1] + 1;
    w_y[2] = (a + 2)*abs(stemp_y[2] * stemp_y[2] * stemp_y[2]) - (a + 3)*stemp_y[2] * stemp_y[2] + 1;
    w_y[3] = a*abs(stemp_y[3] * stemp_y[3] * stemp_y[3]) - 5 * a*stemp_y[3] * stemp_y[3] + 8 * a*abs(stemp_y[3]) - 4 * a;
}

注：由於作者程式設計能力有限，希望有人能指正一下怎麼優化這裡的程式，這個程式只是實現了演算法，執行
速度慢的要死不能忍受！

效果展示