图形图像处理－之－高质量的快速的图像缩放上篇近邻取样插值和其速度优化

ゝ一世哀愁。 2022-06-12 11:47 142阅读 0赞

转自 [http://blog.csdn.net/housisong/article/details/1404896][http_blog.csdn.net_housisong_article_details_1404896]

图形图像处理－之－高质量的快速的图像缩放 上篇 近邻取样插值和其速度优化  
                   [HouSisong@GMail.com][HouSisong_GMail.com]  2006.11.22

(2009.03.07  可以到这里下载缩放算法的完整的可以编译的项目源代码:  [http://blog.csdn.net/housisong/archive/2009/03/07/3967270.aspx][http_blog.csdn.net_housisong_archive_2009_03_07_3967270.aspx]  )

( 2007.06.06 更新测试数据,编译器由vc6改为vc2005,CPU由赛扬2G改为AMD64x2 4200+(2.1G) )

(2007.01.02更新)

tag:图像缩放,速度优化,定点数优化,近邻取样插值,二次线性插值,三次线性插值,  
   MipMap链,三次卷积插值,MMX,SSE,SSE2,CPU缓存优化

摘要:首先给出一个基本的图像缩放算法，然后一步一步的优化其速度和缩放质量；

高质量的快速的图像缩放 全文 分为:  
     上篇 近邻取样插值和其速度优化  
     中篇 二次线性插值和三次卷积插值  
     下篇 三次线性插值和MipMap链  
     补充 使用SSE2优化

正文：

为了便于讨论，这里只处理32bit的ARGB颜色；  
  代码使用C++;涉及到汇编优化的时候假定为x86平台;使用的编译器为vc2005;  
  为了代码的可读性,没有加入异常处理代码;  
  测试使用的CPU为AMD64x2 4200+(2.37G)  和 Intel Core2 4400(2.00G);

速度测试说明:  
  只测试内存数据到内存数据的缩放  
  测试图片都是800\*600缩放到1024\*768; fps表示每秒钟的帧数,值越大表示函数越快  
  
  
//Windows GDI相关函数参考速度:  
//==============================================================================  
// BitBlt             544.7 fps  //is copy 800\*600 to 800\*600  
// BitBlt             331.6 fps  //is copy 1024\*1024 to 1024\*1024  
// StretchBlt         232.7 fps  //is zoom 800\*600 to 1024\*1024

A: 首先定义图像数据结构:

\#define asm \_\_asm  
  
typedef unsigned  char TUInt8;  //  \[0..255\]   
 struct TARGB32       // 32 bit color   
\{  
    TUInt8  B,G,R,A;           //  A is alpha   
\};  
  
struct TPicRegion   // 一块颜色数据区的描述，便于参数传递   
\{  
    TARGB32\*    pdata;          // 颜色数据首地址   
     long        byte\_width;     // 一行数据的物理宽度(字节宽度)；  
                 // abs(byte\_width)有可能大于等于width\*sizeof(TARGB32);   
     long        width;          // 像素宽度   
     long        height;         // 像素高度   
\};  
  
// 那么访问一个点的函数可以写为：   
inline TARGB32& Pixels( const TPicRegion& pic, const  long x, const  long y)  
\{  
     return ( (TARGB32\*)((TUInt8\*)pic.pdata+pic.byte\_width\*y) )\[x\];  
\}

B: 缩放原理和公式图示:

![zoom_map.GIF][]

缩放后图片             原图片  
   (宽DW,高DH)          (宽SW,高SH)

(Sx-0)/(SW-0)=(Dx-0)/(DW-0)   (Sy-0)/(SH-0)=(Dy-0)/(DH-0)  
 =>   Sx=Dx\*SW/DW                    Sy=Dy\*SH/DH

C: 缩放算法的一个参考实现

//给出一个最简单的缩放函数(插值方式为近邻取样,而且我“尽力”把它写得慢一些了:D)  
//Src.PColorData指向源数据区,Dst.PColorData指向目的数据区  
//函数将大小为Src.Width\*Src.Height的图片缩放到Dst.Width\*Dst.Height的区域中

void PicZoom0( const TPicRegion& Dst, const TPicRegion& Src)  
\{  
     if (  (0==Dst.width)||(0==Dst.height)  
        ||(0==Src.width)||(0==Src.height))  return;  
     for ( long x=0;x<Dst.width;++x)  
    \{  
         for ( long y=0;y<Dst.height;++y)  
        \{  
             long srcx=(x\*Src.width/Dst.width);  
             long srcy=(y\*Src.height/Dst.height);  
            Pixels(Dst,x,y)=Pixels(Src,srcx,srcy);  
        \}  
    \}  
\}

//速度测试:  
//==============================================================================  
// PicZoom0            19.4 fps

D: 优化PicZoom0函数

a.PicZoom0函数并没有按照颜色数据在内存中的排列顺序读写(内部循环递增y行  
索引)，将造成CPU缓存预读失败和内存颠簸导致巨大的性能损失,(很多硬件都有这种特性,  
包括缓存、内存、显存、硬盘等,优化顺序访问，随机访问时会造成巨大的性能损失)  
所以先交换x,y循环的顺序:

void PicZoom1( const TPicRegion& Dst, const TPicRegion& Src)  
\{  
     if (  (0==Dst.width)||(0==Dst.height)  
        ||(0==Src.width)||(0==Src.height))  return;  
     for ( long y=0;y<Dst.height;++y)  
    \{  
         for ( long x=0;x<Dst.width;++x)  
        \{  
             long srcx=(x\*Src.width/Dst.width);  
             long srcy=(y\*Src.height/Dst.height);  
            Pixels(Dst,x,y)=Pixels(Src,srcx,srcy);  
        \}  
    \}  
\}

//速度测试:  
//==============================================================================  
// PicZoom1            30.1 fps

b.“(x\*Src.Width/Dst.Width)”表达式中有一个除法运算，它属于很慢的操作(比一般  
的加减运算慢几十倍!),使用定点数的方法来优化它；

void PicZoom2( const TPicRegion& Dst, const TPicRegion& Src)  
\{  
     if (  (0==Dst.width)||(0==Dst.height)  
        ||(0==Src.width)||(0==Src.height))  return;  
    //函数能够处理的最大图片尺寸65536\*65536  
    unsigned long xrIntFloat\_16=(Src.width<<16)/Dst.width+1;  // 16.16格式定点数  
    unsigned long yrIntFloat\_16=(Src.height<<16)/Dst.height+1;  // 16.16格式定点数

//可证明: (Dst.width-1)\*xrIntFloat\_16<Src.width成立

for (unsigned long y=0;y<Dst.height;++y)  
    \{  
        for (unsigned long x=0;x<Dst.width;++x)  
        \{  
            unsigned long srcx=(x\*xrIntFloat\_16)>>16;  
            unsigned long srcy=(y\*yrIntFloat\_16)>>16;  
            Pixels(Dst,x,y)=Pixels(Src,srcx,srcy);  
        \}  
    \}

//速度测试:  
//==============================================================================  
// PicZoom2           185.8 fps

c.  在x的循环中y一直不变，那么可以提前计算与y相关的值; 1.可以发现srcy的值和x变量无关，可以提前到x轴循环之前；2.展开Pixels函数，优化与y相关的指针计算；

void PicZoom3(const TPicRegion& Dst,const TPicRegion& Src)  
\{  
     if (  (0==Dst.width)||(0==Dst.height)  
        ||(0==Src.width)||(0==Src.height))  return;  
    unsigned long xrIntFloat\_16=(Src.width<<16)/Dst.width+1;   
    unsigned long yrIntFloat\_16=(Src.height<<16)/Dst.height+1;

unsigned long dst\_width=Dst.width;  
    TARGB32\* pDstLine=Dst.pdata;  
    unsigned long srcy\_16=0;  
    for (unsigned long y=0;y<Dst.height;++y)  
    \{  
        TARGB32\* pSrcLine=((TARGB32\*)((TUInt8\*)Src.pdata+Src.byte\_width\*(srcy\_16>>16)));  
        unsigned long srcx\_16=0;  
        for (unsigned long x=0;x<dst\_width;++x)  
        \{  
            pDstLine\[x\]=pSrcLine\[srcx\_16>>16\];  
            srcx\_16+=xrIntFloat\_16;  
        \}  
        srcy\_16+=yrIntFloat\_16;  
        ((TUInt8\*&)pDstLine)+=Dst.byte\_width;  
    \}  
\}

//速度测试:  
//==============================================================================  
// PicZoom3           414.4 fps

d.定点数优化使函数能够处理的最大图片尺寸和缩放结果(肉眼不可察觉的误差)受到了一  
定的影响,这里给出一个使用浮点运算的版本,可以在有这种需求的场合使用:

void PicZoom3\_float( const TPicRegion& Dst, const TPicRegion& Src)  
\{  
     // 注意: 该函数需要FPU支持   
     if (  (0==Dst.width)||(0==Dst.height)  
        ||(0==Src.width)||(0==Src.height))  return;  
    double xrFloat=1.000000001/((double)Dst.width/Src.width);  
    double yrFloat=1.000000001/((double)Dst.height/Src.height);  
  
    unsigned  short RC\_Old;  
    unsigned  short RC\_Edit;  
    asm   // 设置FPU的取整方式  为了直接使用fist浮点指令   
    \{  
        FNSTCW  RC\_Old              //  保存协处理器控制字,用来恢复   
        FNSTCW  RC\_Edit             //  保存协处理器控制字,用来修改   
        FWAIT  
        OR      RC\_Edit, 0x0F00     //  改为 RC=11  使FPU向零取整        
        FLDCW   RC\_Edit             //  载入协处理器控制字,RC场已经修改   
    \}  
  
    unsigned  long dst\_width=Dst.width;  
    TARGB32\* pDstLine=Dst.pdata;  
     double srcy=0;  
     for (unsigned  long y=0;y<Dst.height;++y)  
    \{  
        TARGB32\* pSrcLine=((TARGB32\*)((TUInt8\*)Src.pdata+Src.byte\_width\*(( long)srcy)));  
         /\*\*//\*   
        double srcx=0;  
        for (unsigned long x=0;x<dst\_width;++x)  
        \{  
            pDstLine\[x\]=pSrcLine\[(unsigned long)srcx\];//因为默认的浮点取整是一个很慢  
                                     //的操作! 所以才使用了直接操作FPU的内联汇编代码。  
            srcx+=xrFloat;  
        \} \*/  
        asm fld       xrFloat             // st0==xrFloat   
        asm fldz                          // st0==0   st1==xrFloat   
        unsigned  long srcx=0;  
         for ( long x=0;x<dst\_width;++x)  
        \{  
            asm fist dword ptr srcx       // srcx=(long)st0   
            pDstLine\[x\]=pSrcLine\[srcx\];  
            asm fadd  st,st(1)            // st0+=st1   st1==xrFloat   
        \}  
        asm fstp      st  
        asm fstp      st  
  
        srcy+=yrFloat;  
        ((TUInt8\*&)pDstLine)+=Dst.byte\_width;  
    \}  
  
    asm   // 恢复FPU的取整方式   
    \{  
        FWAIT  
        FLDCW   RC\_Old   
    \}  
\}

//速度测试:  
//==============================================================================  
// PicZoom3\_float     286.2 fps

e.注意到这样一个事实:每一行的缩放比例是固定的;那么可以预先建立一个缩放映射表格  
  来处理缩放映射算法(PicZoom3\_Table和PicZoom3\_float的实现等价);

void PicZoom3\_Table( const TPicRegion& Dst, const TPicRegion& Src)  
\{  
     if (  (0==Dst.width)||(0==Dst.height)  
        ||(0==Src.width)||(0==Src.height))  return;  
    unsigned  long dst\_width=Dst.width;  
    unsigned  long\* SrcX\_Table =  new unsigned  long\[dst\_width\];  
     for (unsigned  long x=0;x<dst\_width;++x) // 生成表 SrcX\_Table   
    \{  
        SrcX\_Table\[x\]=(x\*Src.width/Dst.width);  
    \}  
  
    TARGB32\* pDstLine=Dst.pdata;  
     for (unsigned  long y=0;y<Dst.height;++y)  
    \{  
        unsigned  long srcy=(y\*Src.height/Dst.height);  
        TARGB32\* pSrcLine=((TARGB32\*)((TUInt8\*)Src.pdata+Src.byte\_width\*srcy));  
         for (unsigned  long x=0;x<dst\_width;++x)  
            pDstLine\[x\]=pSrcLine\[SrcX\_Table\[x\]\];  
        ((TUInt8\*&)pDstLine)+=Dst.byte\_width;  
    \}  
  
    delete \[\] SrcX\_Table;  
\}

//速度测试:  
//==============================================================================  
// PicZoom3\_Table     390.1 fps

f.为了加快缩放，可以采用根据缩放比例动态生成函数的方式来得到更快的缩放函数；这  
  有点像编译器的工作原理；要实现它需要的工作量比较大(或比较晦涩)就不再实现了；  
  (动态生成是一种不错的思路，但个人觉得对于缩放，实现它的必要性不大)

g.现代CPU中，在读取数据和写入数据时，都有自动的缓存机制；很容易知道，算法中生  
  成的数据不会很快再次使用，所以不需要写入缓存的帮助；在SSE指令集中增加了movntq  
  等指令来完成这个功能；  
  (尝试过利用CPU显式prefetcht0、prefetchnta预读指令或直接的mov读取指令等速度反  
   而略有下降:(   但预读在copy算法中速度优化效果很明显 )

void PicZoom3\_SSE( const TPicRegion& Dst, const TPicRegion& Src)  
\{  
     // 警告: 函数需要CPU支持MMX和movntq指令   
     if (  (0==Dst.width)||(0==Dst.height)  
        ||(0==Src.width)||(0==Src.height))  return;  
    unsigned long xrIntFloat\_16=(Src.width<<16)/Dst.width+1;   
    unsigned long yrIntFloat\_16=(Src.height<<16)/Dst.height+1;  
  
    unsigned  long dst\_width=Dst.width;  
    TARGB32\* pDstLine=Dst.pdata;  
    unsigned  long srcy\_16=0;  
     for (unsigned  long y=0;y<Dst.height;++y)  
    \{  
        TARGB32\* pSrcLine=((TARGB32\*)((TUInt8\*)Src.pdata+Src.byte\_width\*(srcy\_16>>16)));  
  
        asm  
        \{  
            push      ebp  
            mov       esi,pSrcLine  
            mov       edi,pDstLine  
            mov       edx,xrIntFloat\_16  
            mov       ecx,dst\_width  
            xor       ebp,ebp            // srcx\_16=0   
  
            and    ecx, (not 3)     //循环 4次展开   
            TEST   ECX,ECX    // nop   
            jle    EndWriteLoop  
  
            lea       edi,\[edi+ecx\*4\]  
            neg       ecx  
  
               // todo: 预读   
  
                WriteLoop:  
                        mov       eax,ebp  
                        shr       eax,16             // srcx\_16>>16   
                        lea       ebx,\[ebp+edx\]  
                        movd      mm0,\[esi+eax\*4\]  
                        shr       ebx,16             // srcx\_16>>16   
                        PUNPCKlDQ mm0,\[esi+ebx\*4\]  
                        lea       ebp,\[ebp+edx\*2\]  
                         
                         //  movntq qword ptr \[edi+ecx\*4\], mm0   // 不使用缓存的写入指令   
                        asm \_emit 0x0F asm \_emit 0xE7 asm \_emit 0x04 asm \_emit 0x8F    
  
                        mov       eax,ebp  
                        shr       eax,16             // srcx\_16>>16   
                        lea       ebx,\[ebp+edx\]  
                        movd      mm1,\[esi+eax\*4\]  
                        shr       ebx,16             // srcx\_16>>16   
                        PUNPCKlDQ mm1,\[esi+ebx\*4\]  
                        lea       ebp,\[ebp+edx\*2\]  
                          
                         //  movntq qword ptr \[edi+ecx\*4+8\], mm1  // 不使用缓存的写入指令   
                        asm \_emit 0x0F asm \_emit 0xE7 asm \_emit 0x4C asm \_emit 0x8F asm \_emit 0x08  
  
                        add ecx, 4  
                        jnz WriteLoop  
  
                         // sfence  // 刷新写入  
                        asm \_emit 0x0F asm \_emit 0xAE asm \_emit 0xF8    
                        emms  
                EndWriteLoop:  
  
            mov    ebx,ebp  
            pop    ebp  
  
             // 处理边界  循环次数为0,1,2,3；(这个循环可以展开,做一个跳转表,略)   
            mov    ecx,dst\_width  
            and    ecx,3  
            TEST   ECX,ECX  
            jle    EndLineZoom  
  
            lea       edi,\[edi+ecx\*4\]  
            neg       ecx  
      StartBorder:  
            mov       eax,ebx  
            shr       eax,16             // srcx\_16>>16   
            mov       eax,\[esi+eax\*4\]  
            mov       \[edi+ecx\*4\],eax  
            add       ebx,edx  
  
            inc       ECX  
            JNZ       StartBorder  
      EndLineZoom:  
        \}  
  
         //  
        srcy\_16+=yrIntFloat\_16;  
        ((TUInt8\*&)pDstLine)+=Dst.byte\_width;  
    \}  
\}

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

//鉴于有读者反映汇编代码阅读困难，这里给出一个使用intel提供的函数调用方式的实现，  
//读者可以相互对照来阅读代码  
//要编译PicZoom3\_SSE\_mmh,需要\#include <mmintrin.h> \#include <xmmintrin.h>  
//并且需要编译器支持  
//函数PicZoom3\_SSE\_mmh速度为 593.7 fps

void PicZoom3\_SSE\_mmh(const TPicRegion& Dst,const TPicRegion& Src)  
\{  
    //警告: 函数需要CPU支持MMX和movntq指令

if (  (0==Dst.width)||(0==Dst.height)  
        ||(0==Src.width)||(0==Src.height)) return;

unsigned long xrIntFloat\_16=(Src.width<<16)/Dst.width+1;   
    unsigned long yrIntFloat\_16=(Src.height<<16)/Dst.height+1;

unsigned long dst\_width=Dst.width;  
    TARGB32\* pDstLine=Dst.pdata;  
    unsigned long srcy\_16=0;  
    unsigned long for4count=dst\_width/4\*4;  
    for (unsigned long y=0;y<Dst.height;++y)  
    \{  
        TARGB32\* pSrcLine=((TARGB32\*)((TUInt8\*)Src.pdata+Src.byte\_width\*(srcy\_16>>16)));  
        unsigned long srcx\_16=0;  
        unsigned long x;  
        for (x=0;x<for4count;x+=4)//循环4次展开  
        \{  
            \_\_m64 m0=\_m\_from\_int(\*(int\*)(&pSrcLine\[srcx\_16>>16\]));  
            srcx\_16+=xrIntFloat\_16;  
            m0=\_m\_punpckldq(m0, \_m\_from\_int(\*(int\*)(&pSrcLine\[srcx\_16>>16\])) );  
            srcx\_16+=xrIntFloat\_16;  
            \_\_m64 m1=\_m\_from\_int(\*(int\*)(&pSrcLine\[srcx\_16>>16\]));  
            srcx\_16+=xrIntFloat\_16;  
            m1=\_m\_punpckldq(m1, \_m\_from\_int(\*(int\*)(&pSrcLine\[srcx\_16>>16\])) );  
            srcx\_16+=xrIntFloat\_16;  
            \_mm\_stream\_pi((\_\_m64 \*)&pDstLine\[x\],m0); //不使用缓存的写入指令  
            \_mm\_stream\_pi((\_\_m64 \*)&pDstLine\[x+2\],m1); //不使用缓存的写入指令  
        \}  
        for (x=for4count;x<dst\_width;++x)//处理边界  
        \{  
            pDstLine\[x\]=pSrcLine\[srcx\_16>>16\];  
            srcx\_16+=xrIntFloat\_16;  
        \}  
        srcy\_16+=yrIntFloat\_16;  
        ((TUInt8\*&)pDstLine)+=Dst.byte\_width;  
    \}  
    \_m\_empty();  
\}

//速度测试:  
//==============================================================================  
// PicZoom3\_SSE     711.7 fps

E: 缩放效果图：

![test0.GIF][]            ![test1.GIF][]

原图       放大图(x轴放大8倍,y轴放大12倍)

![zaka0.JPG][]             ![zaka1.JPG][]                  ![zaka2.JPG][]  
     原图        缩小图(缩小到0.66倍)      放大图(放大到1.6倍)

F: 把测试成绩放在一起：  
  
  
//CPU: AMD64x2 4200+(2.1G)  zoom 800\*600 to 1024\*768  
//==============================================================================  
// BitBlt             544.7 fps  //is copy 800\*600 to 800\*600  
// BitBlt             331.6 fps  //is copy 1024\*1024 to 1024\*1024  
// StretchBlt         232.7 fps  //is zoom 800\*600 to 1024\*1024  
//   
// PicZoom0            19.4 fps  
// PicZoom1            30.1 fps  
// PicZoom2           185.8 fps  
// PicZoom3           414.4 fps  
// PicZoom3\_float     286.2 fps  
// PicZoom3\_Table     390.1 fps  
// PicZoom3\_SSE       711.7 fps

补充Intel Core2 4400上的测试成绩:  
  
//CPU: Intel Core2 4400(2.00G)  zoom 800\*600 to 1024\*768  
//==============================================================================  
// PicZoom0            15.0 fps  
// PicZoom1            63.9 fps  
// PicZoom2           231.2 fps  
// PicZoom3           460.5 fps  
// PicZoom3\_float     422.5 fps  
// PicZoom3\_Table     457.6 fps  
// PicZoom3\_SSE      1099.7 fps

[http_blog.csdn.net_housisong_article_details_1404896]: http://blog.csdn.net/housisong/article/details/1404896
[HouSisong_GMail.com]: mailto:HouSisong@263.net
[http_blog.csdn.net_housisong_archive_2009_03_07_3967270.aspx]: http://blog.csdn.net/housisong/archive/2009/03/07/3967270.aspx
[zoom_map.GIF]: /images/20220612/fcd75e5375a442c2a14aa3201c2f2b95.png
[test0.GIF]: https://p-blog.csdn.net/images/p_blog_csdn_net/housisong/test0.GIF
[test1.GIF]: /images/20220612/f5b9bd622f344a3d849e61d35f1ba385.png
[zaka0.JPG]: https://p-blog.csdn.net/images/p_blog_csdn_net/housisong/zaka0.JPG
[zaka1.JPG]: https://p-blog.csdn.net/images/p_blog_csdn_net/housisong/zaka1.JPG
[zaka2.JPG]: /images/20220612/437814dec70d440b83d73479c248c3f2.png