[轉帖]線性插值演算法實現映像縮放、旋轉詳解

  線性插值演算法實現映像縮放、旋轉詳解
 
      這是一篇關於圖形處理線性插值演算法細節的文章，轉載此文的目的在於能給那些對影像處理演算法感興趣的網友一些啟示，對於大量的入門網友來說，這樣的文章或許有些讓人眼暈，但我相信哪怕只理解一些表皮的影像處理演算法知識，以後在使用軟體處理圖片時便能做到“心裡有數”，還是有所助宜的。
　　在Windows中做過映像方面程式的人應該都知道Windows的GDI有一個API函數：StretchBlt，對應在VCL中是TCanvas類的StretchDraw方法。它可以很簡單地實現映像的縮放操作。但問題是它是用了速度最快，最簡單但效果也是最差的“最近鄰域法”，雖然在大多數情況下，它也夠用了，但對於要求較高的情況就不行了。

　　不久前做了一個小玩意兒，用於管理我用DC拍的一堆照片，其中有一個外掛程式提供了縮放功能，目前的版本就是用了StretchDraw，有時效果不能令人滿意，我一直想加入兩個更好的：線性插值法和三次樣條法。經過研究發現三次樣條法的計算量實在太大，不太實用，所以決定就只做線性插值法的版本了。

　　從數位影像處理的基本理論，我們可以知道：映像的變形變換就是源映像到靶心圖表像的座標變換。簡單的想法就是把源映像的每個點座標通過變形運算轉為靶心圖表像的相應點的新座標，但是這樣會導致一個問題就是目標點的座標通常不會是整數，而且像放大操作會導致靶心圖表像中沒有被源映像的點映射到，這是所謂“向前映射”方法的缺點。所以一般都是採用“逆向映射”法。

　　但是逆向映射法同樣會出現映射到源映像座標時不是整數的問題。這裡就需要“重採樣濾波器”。這個術語看起來很專業，其實不過是因為它借用了電子訊號處理中的慣用說法（在大多數情況下，它的功能類似於電子訊號處理中的帶通濾波器），理解起來也不複雜，就是如何確定這個非整數座標處的點應該是什麼顏色的問題。前面說到的三種方法：最近鄰域法，線性插值法和三次樣條法都是所謂的“重採樣濾波器”。

　　所謂“最近鄰域法”就是把這個非整數座標作一個四捨五入，取最近的整數點座標處的點的顏色。而“線性插值法”就是根據周圍最接近的幾個點（對於平面映像來說，共有四點）的顏色作線性插值計算（對於平面映像來說就是二維線性插值）來估計這點的顏色，在大多數情況下，它的準確度要高於最近鄰域法，當然效果也要好得多，最明顯的就是在放大時，映像邊緣的鋸齒比最近鄰域法小非常多。當然它同時還帶業個問題：就是映像會顯得比較柔和。這個濾波器用專業術語來說（呵呵，賣弄一下偶的專業^_^）叫做：帶阻效能好，但有帶通損失，通帶曲線的矩形係數不高。至於三次樣條法我就不說了，複雜了一點，可自行參考數位影像處理方面的專業書籍，如本文的參考文獻。

　　再來討論一下座標變換的演算法。簡單的空間變換可以用一個變換矩陣來表示：

[x’,y’,w’]=[u,v,w]*T

其中：x’,y’為靶心圖表像座標，u,v為源映像座標，w,w’稱為齊次座標，通常設為1，T為一個3X3的變換矩陣。

這種表示方法雖然很數學化，但是用這種形式可以很方便地表示多種不同的變換，如平移，旋轉，縮放等。對於縮放來說，相當於：

[Su 0 0 ]

[x, y, 1] = [u, v, 1] * | 0 Sv 0 |

[0 0 1 ]

其中Su,Sv分別是X軸方向和Y軸方向上的縮放率，大於1時放大，大於0小於1時縮小，小於0時反轉。

矩陣是不是看上去比較暈？其實把上式按矩陣乘法展開就是：

{ x = u * Su

{ y = v * Sv

就這麼簡單。^_^

　　有了上面三個方面的準備，就可以開始編寫代碼實現了。思路很簡單：首先用兩重迴圈遍曆靶心圖表像的每個點座標，通過上面的變換式（注意：因為是用逆向映射，相應的變換式應該是：u = x / Su 和v = y / Sv）取得源座標。因為源座標不是整數座標，需要進行二維線性插值運算：

P = n*b*PA + n * ( 1 – b )*PB + ( 1 – n ) * b * PC + ( 1 – n ) * ( 1 – b ) * PD

　　其中：n為v（映射後相應點在源映像中的Y軸座標，一般不是整數）下面最接近的行的Y軸座標與v的差；同樣b也類似，不過它是X軸座標。PA-PD分別是(u,v)點周圍最接近的四個（左上，右上，左下，右下）源映像點的顏色（用TCanvas的Pixels屬性）。P為(u,v)點的插值顏色，即(x,y)點的近似顏色。

　　這段代碼我就不寫了，因為它的效率實在太低：要對靶心圖表像的每一個點的RGB進行上面那一串複雜的浮點運算。所以一定要進行最佳化。對於VCL應用來說，有個比較簡單的最佳化方法就是用TBitmap的ScanLine屬性，按行進行處理，可以避免Pixels的像素級操作，對效能可以有很大的改善。這已經是算是用VCL進行影像處理的基本最佳化常識了。不過這個方法並不總是管用的，比如作映像旋轉的時候，這時需要更多的技巧。

　　無論如何，浮點運算的開銷都是比整數大很多的，這個也是一定要最佳化掉的。從上面可以看出，浮點數是在變換時引入的，而變換參數Su,Sv通常就是浮點數，所以就從它下手最佳化。一般來說，Su,Sv可以表示成分數的形式：

Su = ( double )Dw / Sw; Sv = ( double )Dh / Sh

　　其中Dw, Dh為靶心圖表像的寬度和高度，Sw, Sh為源映像的寬度和高度（因為都是整數，為求得浮點結果，需要進行類型轉換）。

　　將新的Su, Sv代入前面的變換公式和插值公式，可以匯出新的插值公式：

因為：

b = 1 – x * Sw % Dw / ( double )Dw; n = 1 – y * Sh % Dh / ( double )Dh

設：

B = Dw – x * Sw % Dw; N = Dh – y * Sh % Dh

則：

b = B / ( double )Dw; n = N / ( double )Dh

用整數的B，N代替浮點的b, n，轉換插值公式：

P = ( B * N * ( PA – PB – PC + PD ) + Dw * N * PB + DH * B * PC + ( Dw * Dh – Dh * B – Dw * N ) * PD ) / ( double )( Dw * Dh )

這裡最終結果P是浮點數，對其四捨五入即可得到結果。為完全消除浮點數，可以用這樣的方法進行四捨五入：

P = ( B * N … * PD + Dw * Dh / 2 ) / ( Dw * Dh )

這樣，P就直接是四捨五入後的整數值，全部的計算都是整數運算了。

簡單最佳化後的代碼如下：

int __fastcall TResizeDlg::Stretch_Linear(Graphics::TBitmap * aDest, Graphics::TBitmap * aSrc)

{

int sw = aSrc->Width - 1, sh = aSrc->Height - 1, dw = aDest->Width - 1, dh = aDest->Height - 1;

int B, N, x, y;

int nPixelSize = GetPixelSize( aDest->PixelFormat );

BYTE * pLinePrev, *pLineNext;

BYTE * pDest;

BYTE * pA, *pB, *pC, *pD;

for ( int i = 0; i <= dh; ++i )

{

pDest = ( BYTE * )aDest->ScanLine[i];

y = i * sh / dh;

N = dh - i * sh % dh;

pLinePrev = ( BYTE * )aSrc->ScanLine[y++];

pLineNext = ( N == dh ) ? pLinePrev : ( BYTE * )aSrc->ScanLine[y];

for ( int j = 0; j <= dw; ++j )

{

x = j * sw / dw * nPixelSize;

B = dw - j * sw % dw;

pA = pLinePrev + x;

pB = pA + nPixelSize;

pC = pLineNext + x;

pD = pC + nPixelSize;

if ( B == dw )

{

pB = pA;

pD = pC;

}

for ( int k = 0; k < nPixelSize; ++k )

*pDest++ = ( BYTE )( int )(

( B * N * ( *pA++ - *pB - *pC + *pD ) + dw * N * *pB++

+ dh * B * *pC++ + ( dw * dh - dh * B - dw * N ) * *pD++

+ dw * dh / 2 ) / ( dw * dh )

);

}

}

return 0;

}

　　應該說還是比較簡潔的。因為寬度高度都是從0開始算，所以要減一，GetPixelSize是根據PixelFormat屬性來判斷每個像素有多少位元組，此代碼只支援24或32位色的情況（對於15或16位色需要按位拆開—因為不拆開的話會在計算中出現不期望的進位或借位，導致映像顏色混亂—處理較麻煩；對於8位及8位以下索引色需要查調色盤，並且需要重索引，也很麻煩，所以都不支援；但8位灰階映像可以支援）。另外代碼中加入一些在映像邊緣時防止訪問越界的代碼。

　　通過比較，在PIII-733的機器上，靶心圖表像小於1024x768的情況下，基本感覺不出速度比StretchDraw有明顯的慢（用浮點時感覺比較明顯）。效果也相當令人滿意，不論是縮小還是放大，映像品質比StretchDraw方法有明顯提高。

　　不過由於採用了整數運算，有一個問題必須加以重視，那就是溢出的問題：由於式中的分母是dw * dh，而結果應該是一個Byte即8位位元，有符號整數最大可表示31位位元，所以dw * dh的值不能超過23位位元，即按2:1的寬高比計算靶心圖表像解析度不能超過4096*2048。當然這個也是可以通過用無符號數（可以增加一位）及降低計算精度等方法來實現擴充的，有興趣的朋友可以自己試試。

　　當然這段代碼還遠沒有最佳化到極致，而且還有很多問題沒有深入研究，比如抗混疊（anti-aliasing）等，有興趣的朋友可以自行參考相關書籍研究。（轉自：迪派影像）