用線性插值算法實現圖像縮放
猛禽 [Mental Studio](個人專欄 )(BLOG )
http://mental.mentsu.com
在 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 ）等，有興趣的朋友可以自行參考相關書籍研究，如果你有什么研究成果，非常歡迎你為我的程序編寫插件實現。
[ Mental Studio ] 猛禽
參考文獻：
崔屹《數字圖像處理技術與應用》電子工業出版社， 1997
//上面的只能是24位的位圖，現在修改了以個Delphi版本的，支持多種位圖格式了，應該。測試，32和24位都可
procedure StretchBitmap(Dest, Src: TBitmap);
var
sw, sh, dw, dh, B, N, x, y, i, j, k, nPixelSize: DWord;
pLinePrev, pLineNext, pDest, pA, pB, pC, pD: PByte;
begin
sw := Src.Width -1;
sh := Src.Height -1;
dw := Dest.Width -1;
dh := Dest.Height -1;
//獲得顯示模式
nPixelSize := Integer(Src.PixelFormat);
if nPixelSize < 4 then
nPixelSize := 4
else if nPixelSize = 4 then
inc(nPixelSize)
else if nPixelSize > 7 then
nPixelSize := 7;
Dest.PixelFormat := TPixelFormat(nPixelSize);
nPixelSize := nPixelSize - 3;
for i := 0 to dh do
begin
pDest := Dest.ScanLine[i];
y := i * sh div dh;
N := dh - i * sh mod dh;
pLinePrev := Src.ScanLine[y];
Inc(y);
if N = dh then
pLineNext := pLinePrev
else
pLineNext := Src.ScanLine[y];
for j := 0 to dw do
begin
x := j * sw div dw * nPixelSize;
B := dw - j * sw mod dw;
pA := pLinePrev;
Inc(pA, x);
pB := pA;
Inc(pB, nPixelSize);
pC := pLineNext;
Inc(pC, x);
pD := pC;
Inc(pD, nPixelSize);
if B = dw then begin
pB := pA;
pD := pC;
end;
for k := 0 to nPixelSize -1 do
begin
pDest^ := Byte(DWord( (B * N * DWord(pA^ - pB^ - pC^ + pD^) + dw * N * pB^
+ dh * B * pC^ + (dw * dh - dh * B - dw * N)* pD^
+ dw * dh div 2) div (dw * dh) ));
Inc(pDest);
Inc(pA);
Inc(pB);
Inc(pC);
Inc(pD);
end;
end;
end;
end;