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

數學原理

假設源影象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_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;
}

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

效果展示