用线性插值算法实现图像缩放
猛禽 [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;