雙線性插值理論與代碼執行個體

1. 雙線性插值

      假設源圖像大小為mxn，目标圖像為axb。那麼兩幅圖像的邊長比分别為：m/a和n/b。注意，通常這個比例不是整數，程式設計存儲的時候要用浮點型。目标圖像的第（i,j）個像素點（i行j列）可以通過邊長比對應回源圖像。其對應坐标為（i*m/a,j*n/b）。顯然，這個對應坐标一般來說不是整數，而非整數的坐标是無法在圖像這種離散資料上使用的。雙線性插值通過尋找距離這個對應坐标最近的四個像素點，來計算該點的值（灰階值或者RGB值）。

  若圖像為灰階圖像，那麼（i，j）點的灰階值的數學計算模型是：

f(x,y)=b1+b2x+b3y+b4xy

其中b1,b2,b3,b4是相關的系數。關于其的計算過程如下如下：

      如圖，已知Q12，Q22，Q11，Q21，但是要插值的點為P點，這就要用雙線性插值了，首先在x軸方向上，對R1和R2兩個點進行插值，這個很簡單，然後根據R1和R2對P點進行插值，這就是所謂的雙線性插值。

      雙線性插值，又稱為雙線性内插。在數學上，雙線性插值是有兩個變量的插值函數的線性插值擴充，其核心思想是在兩個方向分别進行一次線性插值。

假如我們想得到未知函數 

 在點 

 的值，假設我們已知函數 

 在 

, 

, 

, 及 

 四個點的值。

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

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

這樣就得到所要的結果 

,

如果選擇一個坐标系統使得 

 的四個已知點坐标分别為 (0, 0)、(0, 1)、(1, 0) 和 (1, 1)，那麼插值公式就可以化簡為

或者用矩陣運算表示為

與這種插值方法名稱不同的是，這種插值方法的結果通常不是線性的，它的形式是

{\displaystyle b_{1}+b_{2}x+b_{3}y+b_{4}xy.\,}

常數的數目都對應于給定的 f 的資料點數目

{\displaystyle b_{1}=f(0,0)}

{\displaystyle b_{2}=f(1,0)-f(0,0)}

{\displaystyle b_{3}=f(0,1)-f(0,0)}

{\displaystyle b_{4}=f(1,1)-f(1,0)-f(0,1)+f(0,0)}

線性插值的結果與插值的順序無關。首先進行 y 方向的插值，然後進行 x 方向的插值，所得到的結果是一樣的。

2.存在的問題

這部分的前提是，你已經明白什麼是雙線性插值并且在給定源圖像和目标圖像尺寸的情況下，可以用筆計算出目标圖像某個像素點的值。當然，最好的情況是你已經用某種語言實作了網上一大堆部落格上原創或轉載的雙線性插值算法，然後發現計算出來的結果和matlab、openCV對應的resize（）函數得到的結果完全不一樣。

那這個究竟是怎麼回事呢？

其實答案很簡單，就是坐标系的選擇問題，或者說源圖像和目标圖像之間的對應問題。

按照網上一些部落格上寫的，源圖像和目标圖像的原點（0，0）均選擇左上角，然後根據插值公式計算目标圖像每點像素，假設你需要将一幅5x5的圖像縮小成3x3，那麼源圖像和目标圖像各個像素之間的對應關系如下：

隻畫了一行，用做示意，從圖中可以很明顯的看到，如果選擇右上角為原點（0，0），那麼最右邊和最下邊的像素實際上并沒有參與計算，而且目标圖像的每個像素點計算出的灰階值也相對于源圖像偏左偏上。

那麼，讓坐标加1或者選擇右下角為原點怎麼樣呢？很不幸，還是一樣的效果，不過這次得到的圖像将偏右偏下。

最好的方法就是，兩個圖像的幾何中心重合，并且目标圖像的每個像素之間都是等間隔的，并且都和兩邊有一定的邊距，這也是matlab和openCV的做法。如下圖：

如果你不懂我上面說的什麼，沒關系，隻要在計算對應坐标的時候改為以下公式即可， 

int x=(i+0.5)*m/a-0.5

int y=(j+0.5)*n/b-0.5 

instead of  

int x=i*m/a

int y=j*n/b 

利用上述公式，将得到正确的雙線性插值結果

總結：

總結一下，我得到的教訓有這麼幾條。

1.網上的一些資料有的時候并不靠譜，自己還是要多做實驗。

2.不要小瞧一些簡單的、基本的算法，讓你寫你未必會寫，而且其中可能還藏着一些玄妙。

3.要多動手程式設計，多體會算法，多看大牛寫的源碼（雖然有的時候很吃力，但是要堅持看）。

源碼：

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);
			}
		}
	}
}//sx=1.2,sy=1.6
           

//圖像基本資訊輸出
	int image_height = image.size().height;
	int image_width = image.size().width;
	cout << "Image Info: height:" << image_height << "  width:" << image_width << endl;

	//===============================================================================================================================
	//圖像處理-----雙線性插值
	cout << "Please input the Sx and Sy:" << endl;
	float Sx, Sy;
	cin >> Sx >> Sy;
	cout << "Sx = " << Sx << "; Sy = " << Sy << ";" << endl;

	Mat final_img;
	int final_img_height, final_img_width;
	final_img_height = image.size().height * Sx;
	final_img_width = image.size().width * Sy;
	final_img.create(final_img_height, final_img_width, CV_8UC1);

	int y, x;
	int x1, x2, y1, y2;
	float temp1, temp2;
	for (y = 0; y < final_img_height; y++)
	{
		for (x = 0; x < final_img_width; x++)
		{
			x1 = (int)(x / Sx);
			x2 = x1 + 1;
			y1 = (int)(y / Sy);
			y2 = y1 + 1;

			//判斷邊界
			if (y2 >= image_height || x2 >= image_width)
			{
				final_img.at<uchar>(y, x) = image.at<uchar>(y1, x1);
				continue;
			}

			temp1 = (x2 - x / Sx) * image.at<uchar>(y1, x1) + (x / Sx - x1) * image.at<uchar>(y1, x2);
			temp2 = (x2 - x / Sx) * image.at<uchar>(y2, x1) + (x / Sx - x1) * image.at<uchar>(y2, x2);
			final_img.at<uchar>(y, x) = (uchar)((y2 - y / Sy) * temp1 + (y / Sy - y1) * temp2);
		}
	}

	//===============================================================================================================================

	//顯示處理前後的圖像
	namedWindow("original_image");
	imshow("original_image", image);
	namedWindow("final_image");
	cout << "Final image Info: height:" << final_img.size().height << "  width:" << final_img.size().width << endl;
	imshow("final_image", final_img);
           

//sx=1.2,sy=1.2