雙線性插值
假設源映像大小為mxn,靶心圖表像為axb。那麼兩幅映像的邊長比分別為:m/a和n/b。注意,通常這個比例不是整數,編程儲存的時候要用浮點型。靶心圖表像的第(i,j)個像素點(i行j列)可以通過邊長比對應回源映像。其對應座標為(i*m/a,j*n/b)。顯然,這個對應座標一般來說不是整數,而非整數的座標是無法在映像這種離散資料上使用的。雙線性插值通過尋找距離這個對應座標最近的四個像素點,來計算該點的值(灰階值或者RGB值)。如果你的對應座標是(2.5,4.5),那麼最近的四個像素是(2,4)、(2,5)、(3,4),(3,5)。
若映像為灰階映像,那麼(i,j)點的灰階值的數學計算模型是:
f(x,y)=b1+b2x+b3y+b4xy
其中,b1,b2,b3,b4是相關的係數。
或者:
f(i,j)=w1*p1+w2*p2+w3*p3+w4*p4;
其中,pi(i=1,2,3,4)為最近的四個像素點,wi(i=1,2,3,4)為各點相應權值。關於權值的計算,見下面的維基百科。
關於其的計算過程如下:
如圖,已知Q12,Q22,Q11,Q21,但是要插值的點為P點,這就要用雙線性插值了,首先在x軸方向上,對R1和R2兩個點進行插值,這個很簡單,然後根據R1和R2對P點進行插值,這就是所謂的雙線性插值。
附:維基百科--雙線性插值:
雙線性插值,又稱為雙線性插補。在數學上,雙線性插值是有兩個變數的插值函數的線性插值擴充,其核心思想是在兩個方向分別進行一次線性插值。
假如我們想得到未知函數 在點 的值,假設我們已知函數 在 , , , 及 四個點的值。
首先在 x 方向進行線性插值,得到
然後在 y 方向進行線性插值,得到
這樣就得到所要的結果 ,
如果選擇一個座標系統使得 的四個已知點座標分別為 (0, 0)、(0, 1)、(1, 0) 和 (1, 1),那麼插值公式就可以化簡為
或者用矩陣運算表示為
這種插值方法的結果通常不是線性,線性插值的結果與插值的順序無關。首先進行 y 方向的插值,然後進行 x 方向的插值,所得到的結果是一樣的。 opencv和Matlab中的雙線性插值
這部分的前提是,你已經明白什麼是雙線性插值並且在給定源映像和靶心圖表像尺寸的情況下,可以用筆計算出靶心圖表像某個像素點的值。當然,最好的情況是你已經用某種語言實現了網上一大堆部落格上原創或轉載的雙線性插值演算法,然後發現計算出來的結果和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
替換下面的公式:
int x=i*m/a
int y=j*n/b
利用上述公式,將得到正確的雙線性插值結果
總結:
總結一下,我得到的教訓有這麼幾條。
1.網上的一些資料有的時候並不靠譜,自己還是要多做實驗。
2.不要小瞧一些簡單的、基本的演算法,讓你寫你未必會寫,而且其中可能還藏著一些玄妙。
3.要多動手編程,多體會演算法,多看大牛寫的源碼(雖然有的時候很吃力,但是要堅持看)。
在影像處理中,雙線性插值演算法的使用頻率相當高,比如在映像的縮放中,在所有的扭曲演算法中,都可以利用該演算法改進處理的視覺效果。首先,我們看看該演算法的簡介。
在數學上,雙線性插值演算法可以看成是兩個變數間的線性插值的延伸。執行該過程的關鍵思路是先在一個方向上執行線性插值,然後再在另外一個方向上插值。下圖示意出這個過程的大概意思。
用一個簡單的數學運算式可以表示如下:
f(x,y)=f(0,0)(1-x)(1-y)+f(1,0)x(1-y)+f(0,1)(1-x)y+f(1,1)xy
合并有關項,可寫為: f(x,y)=(f(0,0)(1-x)+f(1,0)x) (1-y)+(f(0,1)(1-x)+f(1,1)x)y
由上式可以看出,這個過程存在著大量的浮點數運算,對於映像這樣大的計算使用者來說,是一個較為耗時的過程。
考慮到映像的特殊性,他的像素值的計算結果需要落在0到255之間,最多隻有256種結果,由上式可以看出,一般情況下,計算出的f(x,y)是個浮點數,我們還需要對該浮點數進行取整。因此,我們可以考慮將該過程中的所有類似於1-x、1-y的變數放大合適的倍數,得到對應的整數,最後再除以一個合適的整數作為插值的結果。
如何取這個合適的放大倍數呢,要從三個方面考慮,第一:精度問題,如果這個數取得過小,那麼經過計算後可能會導致結果出現較大的誤差。第二,這個數不能太大,太大會導致計算過程超過長整形所能表達的範圍。第三:速度考慮。假如放大倍數取為12,那麼算式在最後的結果中應該需要除以12*12=144,但是如果取為16,則最後的除數為16*16=256,這個數字好,我們可以用右移來實現,而右移要比普通的整除快多了。
綜合考慮上述三條,我們選擇2048這個數比較合適。
下面我們假定某個演算法得到了我們要取樣的座標分別PosX以及PosY,其中PosX=25.489,PosY=58.698。則這個過程的類似程式碼片段如下: 1 NewX = Int(PosX) ' 向下取整,NewX=25
2 NewY = Int(PosY) ' 向下取整,NewY=58
3 PartX = (PosX - NewX) * 2048 ' 對應運算式中的X
4 PartY = (PosY - NewY) * 2048 ' 對應運算式中的Y
5 InvX = 2048 - PartX ' 對應運算式中的1-X
6 InvY = 2048 - PartY ' 對應運算式中的1-Y
7
8 Index1 = SamStride * NewY + NewX * 3 ' 計算取樣點左上方鄰近的那個像素點的記憶體位址
9 Index2 = Index1 + SamStride ' 左下角像素點地址
10 ImageData(Speed + 2) = ((Sample(Index1 + 2) * InvX + Sample(Index1 + 5) * PartX) * InvY + (Sample(Index2 + 2) * InvX + Sample(Index2 + 5) * PartX) * PartY) \ 4194304 '處理紅色分量 11 ImageData(Speed + 1) = ((Sample(Index1 + 1) * InvX + Sample(Index1 + 4) * PartX) * InvY + (Sample(Index2 + 1) * InvX + Sample(Index2 + 4) * PartX) * PartY) \ 4194304 '處理綠色分量
12 ImageData(Speed) = ((Sample(Index1) * InvX + Sample(Index1 + 3) * PartX) * InvY + (Sample(Index2) * InvX + Sample(Index2 + 3) * PartX) * PartY) \ 4194304 '處理藍色分量
以上代碼中涉及到的變數都為整型(PosX及PosY當然為浮點型)。
代碼中Sample數組儲存了從中取樣的映像資料,SamStride為該映像的掃描行大小。
觀察上述代碼,除了有2句涉及到了浮點計算,其他都是整數之間的運算。
在Basic語言中,編譯時間如果勾選所有的進階最佳化,則\ 4194304會被編譯為>>22,即右移22位,如果使用的是C語言,則直接寫為>>22。
需要注意的是,在進行這代代碼前,需要保證PosX以及PosY在合理的範圍內,即不能超出取樣映像的寬度和高度範圍。
通過這樣的改進,速度較直接用浮點類型快至少100%以上,而處理後的效果基本沒有什麼區別。 利用Matlab實現的雙線性插值(映像放大/縮小)
引用:http://www.cnblogs.com/tiandsp/archive/2012/12/03/2800201.html
close all;
clear all;clc;m=1.8; %放大或縮小的高度n=2.3; %放大或縮小的寬度img=imread('lena.jpg');imshow(img);[h w]=size(img);imgn=zeros(h*m,w*n);rot=[m 0 0;0 n 0;0 0 1]; %變換矩陣for i=1:h*m for j=1:w*n pix=[i j 1]/rot; float_Y=pix(1)-floor(pix(1)); float_X=pix(2)-floor(pix(2)); if pix(1) < 1 %邊界處理 pix(1) = 1; end if pix(1) > h pix(1) = h; end if pix(2) < 1 pix(2) =1; end if pix(2) > w pix(2) =w; end pix_up_left=[floor(pix(1)) floor(pix(2))]; %四個相鄰的點 pix_up_right=[floor(pix(1)) ceil(pix(2))]; pix_down_left=[ceil(pix(1)) floor(pix(2))]; pix_down_right=[ceil(pix(1)) ceil(pix(2))]; value_up_left=(1-float_X)*(1-float_Y); %計算臨近四個點的權重 value_up_right=float_X*(1-float_Y); value_down_left=(1-float_X)*float_Y; value_down_right=float_X*float_Y; %按權重進行雙線性插值 imgn(i,j)=value_up_left*img(pix_up_left(1),pix_up_left(2))+ ... value_up_right*img(pix_up_right(1),pix_up_right(2))+ ... value_down_left*img(pix_down_left(1),pix_down_left(2))+ ... value_down_right*img(pix_down_right(1),pix_down_right(2)); endendfigure,imshow(uint8(imgn))
原圖 放大後