精华内容
下载资源
问答
  • 图片缩放算法

    2020-02-18 21:12:28
    项目背景:博主之前做过一个摄像头采集数据,然后在LCD上显示视频数据的项目,假如我们摄像头采集的一帧数据的分辨率比我们的LCD的分辨率要大,那么LCD则无法显示整个图像,这时候我们就要把这么一帧图片进行缩放,...

    项目背景:博主之前做过一个摄像头采集数据,然后在LCD上显示视频数据的项目,假如我们摄像头采集的一帧数据的分辨率比我们的LCD的分辨率要大,那么LCD则无法显示整个图像,这时候我们就要把这么一帧图片进行缩放,然后再显示在LCD上。

     

    本文采用的是近邻取样插值方法缩放图片,比如原图有100个点,想缩放10倍的话,那就要取10个点,(点的编号)比如0 10 20 30 。。。100(其他点就不要了)。

    下面我们直接通过程序来分析:首先我们用  输入分辨率/压缩后的分辨率  可得到一个缩放比例值k,根据这个k值找出我们要保存原图上哪些竖线横线。横线和竖线的交点就是我们要保存的点。

    /**********************************************************************
     * 函数名称: PicZoom
     * 功能描述: 近邻取样插值方法缩放图片
     *            注意该函数会分配内存来存放缩放后的图片,用完后要用free函数释放掉
     *            "近邻取样插值"的原理请参考网友"lantianyu520"所著的"图像缩放算法"
     * 输入参数: ptOriginPic - 内含原始图片的象素数据
     *            ptBigPic    - 内含缩放后的图片的象素数据
     * 输出参数: 无
     * 返 回 值: 0 - 成功, 其他值 - 失败
     ***********************************************************************/
    int PicZoom(PT_PixelDatas ptOriginPic, PT_PixelDatas ptZoomPic)
    {
        unsigned long dwDstWidth = ptZoomPic->iWidth;//压缩后分辨率的宽度值
        unsigned long* pdwSrcXTable;
    	unsigned long x;
    	unsigned long y;
    	unsigned long dwSrcY;
    	unsigned char *pucDest;
    	unsigned char *pucSrc;
    	unsigned long dwPixelBytes = ptOriginPic->iBpp/8;
    
    	if (ptOriginPic->iBpp != ptZoomPic->iBpp)
    	{
    		return -1;
    	}
    
            //pdwSrcXTable数组用来存放要保存的是哪些“竖线”(保存上图红色竖线的位置)
        pdwSrcXTable = malloc(sizeof(unsigned long) * dwDstWidth);
        if (NULL == pdwSrcXTable)
        {
            DBG_PRINTF("malloc error!\n");
            return -1;
        }
            //pdwSrcXTable数组用来存放要保存的是哪些“竖线”(保存上图红色竖线的位置)
        for (x = 0; x < dwDstWidth; x++)//生成表 pdwSrcXTable
        {
            pdwSrcXTable[x]=(x*ptOriginPic->iWidth/ptZoomPic->iWidth);
        }
    
            //循环“压缩后图片的分辨率的高度”次,也就是开始把要留下的点放进输出buffer了
        for (y = 0; y < ptZoomPic->iHeight; y++)
        {		
            //找出每一次要保存的“横线”	
            dwSrcY = (y * ptOriginPic->iHeight / ptZoomPic->iHeight);
    
            //目的的起始地址
    		pucDest = ptZoomPic->aucPixelDatas + y*ptZoomPic->iLineBytes;
            //源的起始地址(横线的位置是每个点的起始地址,竖线的位置是偏移位置)
    		pucSrc  = ptOriginPic->aucPixelDatas + dwSrcY*ptOriginPic->iLineBytes;
    		
            //每条横线上有dwDstWidth个点要保存
            for (x = 0; x <dwDstWidth; x++)
            {
                //结合pdwSrcXTable数组找出要保存的点
    			 memcpy(pucDest+x*dwPixelBytes, pucSrc+pdwSrcXTable[x]*dwPixelBytes, dwPixelBytes);
            }
        }
    
        free(pdwSrcXTable);
    	return 0;
    }

     

     

    展开全文
  • 图片缩放算法原理

    2012-11-21 16:45:09
    图片缩放算法原理,思路清晰,容易理解,适合新手
  • 图片等比例缩放算法

    万次阅读 2017-07-07 13:56:07
    在许多语言中,都希望图片可以等比例缩小或者放大,所以在此提供一个所有语言通用的图片等比例缩小方法的算法。这里以java语言为例子  1.给两个值,设置你想要设置图片的宽和高;定义两个真正的宽和高;这里定义为...

    在许多语言中,都希望图片可以等比例缩小或者放大,但是仅仅依靠语言本身的方法,大多差强人意,所以在此提供一个所有语言通用的图片等比例缩小方法的算法。这里以java语言为例子 


    1.给两个值,设置你想要设置图片的宽和高;定义两个真正的宽和高;这里定义为double 类型 
    double setWidth,setHeight; 
    double width,height; 


    2.获取图片的宽和高 
    double imageWidth=image.getIconWidth(); 
    double imageHeight=image.getIconHeight(); 


    3.写出算法 
     
    if(setWidth/imageWidth<=setHeight/imageHeight)
    {    
    width=imageWidth*(setWidth/imageWidth);
    height=imageHeight*(setWidth/imageWidth);   
    }else{  
    width=imageWidth*(setHeight/imageHeight); 
    height=imageHeight*(setHeight/imageHeight);
    }


    (1)如果设置为int,上面的if就会变成(0<=0)了 (2)算法确实很简单,我们一般取比例值较小的 4,将上面的所得的width,height转化为int型
    展开全文
  • 图像缩放算法

    万次阅读 2016-03-30 22:23:07
    摘要:首先给出一个基本的图像缩放算法,然后一步一步的优化其速度和缩放质量; 高质量的快速的图像缩放 全文 分为:  上篇 近邻取样插值和其速度优化  中篇 二次线性插值和三次卷积插值  下篇 三次线性插值...
    转载别人的,但是这篇文章写得确实太好了,所以想分享出来,可是原创文章地址找不到了 ,很可惜。

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

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

    正文:  

      为了便于讨论,这里只处理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: 缩放原理和公式图示:

        缩放后图片             原图片
       (宽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: 缩放效果图:

                     

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

                                       
         原图        缩小图(缩小到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 

     

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

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

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


    速度测试说明:
      只测试内存数据到内存数据的缩放
      测试图片都是800*600缩放到1024*768; fps表示每秒钟的帧数,值越大表示函数越快


    A:近邻取样插值、二次线性插值、三次卷积插值 缩放效果对比

                                       
           原图         近邻取样缩放到0.6倍     近邻取样缩放到1.6倍

                                             
                    二次线性插值缩放到0.6倍   二次线性插值缩放到1.6倍

                                            
                   三次卷积插值缩放到0.6倍   三次卷积插值缩放到1.6倍

          
     原图 近邻取样缩放到8倍 二次线性插值缩放到8倍 三次卷积插值缩放到8倍 二次线性插值(近似公式)


         近邻取样插值缩放简单、速度快,但很多时候缩放出的图片质量比较差(特别是对于人物、景色等),
    图片的缩放有比较明显的锯齿;使用二次或更高次插值有利于改善缩放效果;


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

    #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];
    }


     


    二次线性插值缩放:

    C: 二次线性插值缩放原理和公式图示:

            

                 缩放后图片                 原图片
                (宽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

      聚焦看看(Sx,Sy)坐标点(Sx,Sy为浮点数)附近的情况;


             


      对于近邻取样插值的缩放算法,直接取Color0颜色作为缩放后点的颜色;
    二次线性插值需要考虑(Sx,Sy)坐标点周围的4个颜色值Color0\Color1\Color2\Color3,
    把(Sx,Sy)到A\B\C\D坐标点的距离作为系数来把4个颜色混合出缩放后点的颜色;
    ( u=Sx-floor(Sx); v=Sy-floor(Sy); 说明:floor函数的返回值为小于等于参数的最大整数 )  
      二次线性插值公式为:
     tmpColor0=Color0*(1-u) + Color2*u;
     tmpColor1=Color1*(1-u) + Color3*u;
            DstColor =tmpColor0*(1-v) + tmpColor2*v;

      展开公式为:
            pm0=(1-u)*(1-v);
            pm1=v*(1-u);
            pm2=u*(1-v);
            pm3=u*v;
      则颜色混合公式为:
            DstColor = Color0*pm0 + Color1*pm1 + Color2*pm2 + Color3*pm3;

    参数函数图示:

          

                                                  二次线性插值函数图示

    对于上面的公式,它将图片向右下各移动了半个像素,需要对此做一个修正;
      =>   Sx=(Dx+0.5)*SW/DW-0.5; Sy=(Dy+0.5)*SH/DH-0.5;
    而实际的程序,还需要考虑到边界(访问源图片可能超界)对于算法的影响,边界的处理可能有各种
    方案(不处理边界或边界回绕或边界饱和或边界映射或用背景颜色混合等;文章中默认使用边界饱和来处理超界);
    比如:边界饱和函数: 

    //访问一个点的函数,(x,y)坐标可能超出图片边界; //边界处理模式:边界饱和
    inline TARGB32 Pixels_Bound(const TPicRegion& pic,long x,long y)
    {
        
    //assert((pic.width>0)&&(pic.height>0));
        bool IsInPic=true;
        
    if (x<0) {x=0; IsInPic=false; } else if (x>=pic.width ) {x=pic.width -1; IsInPic=false; }
        
    if (y<0) {y=0; IsInPic=false; } else if (y>=pic.height) {y=pic.height-1; IsInPic=false; }
        TARGB32 result=Pixels(pic,x,y);
        
    if (!IsInPic) result.a=0;
        
    return result;
    }


     

    D: 二次线性插值缩放算法的一个参考实现:PicZoom_BilInear0
      该函数并没有做什么优化,只是一个简单的浮点实现版本;


     

        inline void Bilinear0(const TPicRegion& pic,float fx,float fy,TARGB32* result)
        {
            
    long x=(long)fx; if (x>fx) --x; //x=floor(fx);    
            long y=(long)fy; if (y>fy) --y; //y=floor(fy);
            
            TARGB32 Color0=Pixels_Bound(pic,x,y);
            TARGB32 Color2=Pixels_Bound(pic,x+1,y);
            TARGB32 Color1=Pixels_Bound(pic,x,y+1);
            TARGB32 Color3=Pixels_Bound(pic,x+1,y+1);

            
    float u=fx-x;
            
    float v=fy-y;
            
    float pm3=u*v;
            
    float pm2=u*(1-v);
            
    float pm1=v*(1-u);
            
    float pm0=(1-u)*(1-v);

            result->a=(pm0*Color0.a+pm1*Color1.a+pm2*Color2.a+pm3*Color3.a);
            result->r=(pm0*Color0.r+pm1*Color1.r+pm2*Color2.r+pm3*Color3.r);
            result->g=(pm0*Color0.g+pm1*Color1.g+pm2*Color2.g+pm3*Color3.g);
            result->b=(pm0*Color0.b+pm1*Color1.b+pm2*Color2.b+pm3*Color3.b);
        }

    void PicZoom_Bilinear0(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;
        TARGB32* pDstLine=Dst.pdata;
        
    for (unsigned long y=0;y<Dst.height;++y)
        {
            
    float srcy=(y+0.4999999)*Src.height/Dst.height-0.5;
            
    for (unsigned long x=0;x<dst_width;++x)
            {
                
    float srcx=(x+0.4999999)*Src.width/Dst.width-0.5;
                Bilinear0(Src,srcx,srcy,&pDstLine[x]);
            }
            ((TUInt8*&)pDstLine)+=Dst.byte_width;
        }
    }


    //速度测试:
    //==============================================================================
    // PicZoom_BilInear0      8.3 fps

     

    E: 把PicZoom_BilInear0的浮点计算改写为定点数实现:PicZoom_BilInear1

        inline void Bilinear1(const TPicRegion& pic,const long x_16,const long y_16,TARGB32* result)
        {
            
    long x=x_16>>16;
            
    long y=y_16>>16;
            TARGB32 Color0=Pixels_Bound(pic,x,y);
            TARGB32 Color2=Pixels_Bound(pic,x+1,y);
            TARGB32 Color1=Pixels_Bound(pic,x,y+1);
            TARGB32 Color3=Pixels_Bound(pic,x+1,y+1);

            unsigned 
    long u_8=(x_16 & 0xFFFF)>>8;
            unsigned 
    long v_8=(y_16 & 0xFFFF)>>8;
            unsigned 
    long pm3_16=(u_8*v_8);
            unsigned 
    long pm2_16=(u_8*(unsigned long)(256-v_8));
            unsigned 
    long pm1_16=(v_8*(unsigned long)(256-u_8));
            unsigned 
    long pm0_16=((256-u_8)*(256-v_8));

            result->a=((pm0_16*Color0.a+pm1_16*Color1.a+pm2_16*Color2.a+pm3_16*Color3.a)>>16);
            result->r=((pm0_16*Color0.r+pm1_16*Color1.r+pm2_16*Color2.r+pm3_16*Color3.r)>>16);
            result->g=((pm0_16*Color0.g+pm1_16*Color1.g+pm2_16*Color2.g+pm3_16*Color3.g)>>16);
            result->b=((pm0_16*Color0.b+pm1_16*Color1.b+pm2_16*Color2.b+pm3_16*Color3.b)>>16);
        }

    void PicZoom_Bilinear1(const TPicRegion& Dst,const TPicRegion& Src)
    {
        
    if (  (0==Dst.width)||(0==Dst.height)
            ||(0==Src.width)||(0==Src.height)) 
    return;

        
    long xrIntFloat_16=((Src.width)<<16)/Dst.width+1; 
        
    long yrIntFloat_16=((Src.height)<<16)/Dst.height+1;
        
    const long csDErrorX=-(1<<15)+(xrIntFloat_16>>1);
        
    const long csDErrorY=-(1<<15)+(yrIntFloat_16>>1);

        unsigned 
    long dst_width=Dst.width;

        TARGB32* pDstLine=Dst.pdata;
        
    long srcy_16=csDErrorY;
        
    long y;
        
    for (y=0;y<Dst.height;++y)
        {
            
    long srcx_16=csDErrorX;
            
    for (unsigned long x=0;x<dst_width;++x)
            {
                Bilinear1(Src,srcx_16,srcy_16,&pDstLine[x]); 
    //border
                srcx_16+=xrIntFloat_16;
            }
            srcy_16+=yrIntFloat_16;
            ((TUInt8*&)pDstLine)+=Dst.byte_width;
        }



    //速度测试:
    //==============================================================================
    // PicZoom_BilInear1     17.7 fps

    F: 二次线性插值需要考略边界访问超界的问题,我们可以将边界区域和内部区域分开处理,这样就可以优化内部的插值实现函数了:比如不需要判断访问超界、减少颜色数据复制、减少一些不必要的重复坐标计算等等


     


     

        inline void Bilinear2_Fast(TARGB32* PColor0,TARGB32* PColor1,unsigned long u_8,unsigned long v_8,TARGB32* result)
        {
            unsigned 
    long pm3_16=u_8*v_8;
            unsigned 
    long pm2_16=(u_8<<8)-pm3_16;
            unsigned 
    long pm1_16=(v_8<<8)-pm3_16;
            unsigned 
    long pm0_16=(1<<16)-pm1_16-pm2_16-pm3_16;
       
            result->a=((pm0_16*PColor0[0].a+pm2_16*PColor0[1].a+pm1_16*PColor1[0].a+pm3_16*PColor1[1].a)>>16);
            result->r=((pm0_16*PColor0[0].r+pm2_16*PColor0[1].r+pm1_16*PColor1[0].r+pm3_16*PColor1[1].r)>>16);
            result->g=((pm0_16*PColor0[0].g+pm2_16*PColor0[1].g+pm1_16*PColor1[0].g+pm3_16*PColor1[1].g)>>16);
            result->b=((pm0_16*PColor0[0].b+pm2_16*PColor0[1].b+pm1_16*PColor1[0].b+pm3_16*PColor1[1].b)>>16);
        }

        inline 
    void Bilinear2_Border(const TPicRegion& pic,const long x_16,const long y_16,TARGB32* result)
        {
            
    long x=(x_16>>16);
            
    long y=(y_16>>16);
            unsigned 
    long u_16=((unsigned short)(x_16));
            unsigned 
    long v_16=((unsigned short)(y_16));

            TARGB32 pixel[4];
            pixel[0]=Pixels_Bound(pic,x,y);
            pixel[1]=Pixels_Bound(pic,x+1,y);
            pixel[2]=Pixels_Bound(pic,x,y+1);
            pixel[3]=Pixels_Bound(pic,x+1,y+1);
            
            Bilinear2_Fast(&pixel[0],&pixel[2],u_16>>8,v_16>>8,result);
        }

    void PicZoom_Bilinear2(const TPicRegion& Dst,const TPicRegion& Src)
    {
        
    if (  (0==Dst.width)||(0==Dst.height)
            ||(0==Src.width)||(0==Src.height)) 
    return;

        
    long xrIntFloat_16=((Src.width)<<16)/Dst.width+1; 
        
    long yrIntFloat_16=((Src.height)<<16)/Dst.height+1;
        
    const long csDErrorX=-(1<<15)+(xrIntFloat_16>>1);
        
    const long csDErrorY=-(1<<15)+(yrIntFloat_16>>1);

        unsigned 
    long dst_width=Dst.width;

        
    //计算出需要特殊处理的边界
        long border_y0=-csDErrorY/yrIntFloat_16+1;              //y0+y*yr>=0; y0=csDErrorY => y>=-csDErrorY/yr
        if (border_y0>=Dst.height) border_y0=Dst.height;
        
    long border_x0=-csDErrorX/xrIntFloat_16+1;     
        
    if (border_x0>=Dst.width ) border_x0=Dst.width; 
        
    long border_y1=(((Src.height-2)<<16)-csDErrorY)/yrIntFloat_16+1; //y0+y*yr<=(height-2) => y<=(height-2-csDErrorY)/yr
        if (border_y1<border_y0) border_y1=border_y0;
        
    long border_x1=(((Src.width-2)<<16)-csDErrorX)/xrIntFloat_16+1; 
        
    if (border_x1<border_x0) border_x1=border_x0;

        TARGB32* pDstLine=Dst.pdata;
        
    long Src_byte_width=Src.byte_width;
        
    long srcy_16=csDErrorY;
        
    long y;
        
    for (y=0;y<border_y0;++y)
        {
            
    long srcx_16=csDErrorX;
            
    for (unsigned long x=0;x<dst_width;++x)
            {
                Bilinear2_Border(Src,srcx_16,srcy_16,&pDstLine[x]); 
    //border
                srcx_16+=xrIntFloat_16;
            }
            srcy_16+=yrIntFloat_16;
            ((TUInt8*&)pDstLine)+=Dst.byte_width;
        }
        
    for (y=border_y0;y<border_y1;++y)
        {
            
    long srcx_16=csDErrorX;
            
    long x;
            
    for (x=0;x<border_x0;++x)
            {
                Bilinear2_Border(Src,srcx_16,srcy_16,&pDstLine[x]);
    //border
                srcx_16+=xrIntFloat_16;
            }

            {
                unsigned 
    long v_8=(srcy_16 & 0xFFFF)>>8;
                TARGB32* PSrcLineColor= (TARGB32*)((TUInt8*)(Src.pdata)+Src_byte_width*(srcy_16>>16)) ;
                
    for (unsigned long x=border_x0;x<border_x1;++x)
                {
                    TARGB32* PColor0=&PSrcLineColor[srcx_16>>16];
                    TARGB32* PColor1=(TARGB32*)((TUInt8*)(PColor0)+Src_byte_width);
                    Bilinear2_Fast(PColor0,PColor1,(srcx_16 & 0xFFFF)>>8,v_8,&pDstLine[x]);
                    srcx_16+=xrIntFloat_16;
                }
            }

            
    for (x=border_x1;x<dst_width;++x)
            {
                Bilinear2_Border(Src,srcx_16,srcy_16,&pDstLine[x]);
    //border
                srcx_16+=xrIntFloat_16;
            }
            srcy_16+=yrIntFloat_16;
            ((TUInt8*&)pDstLine)+=Dst.byte_width;
        }
        
    for (y=border_y1;y<Dst.height;++y)
        {
            
    long srcx_16=csDErrorX;
            
    for (unsigned long x=0;x<dst_width;++x)
            {
                Bilinear2_Border(Src,srcx_16,srcy_16,&pDstLine[x]); 
    //border
                srcx_16+=xrIntFloat_16;
            }
            srcy_16+=yrIntFloat_16;
            ((TUInt8*&)pDstLine)+=Dst.byte_width;
        }
    }


    //速度测试:
    //==============================================================================
    // PicZoom_BilInear2     43.4 fps

    (F'补充:
      如果不想处理边界访问超界问题,可以考虑扩大源图片的尺寸,加一个边框 (“哨兵”优化);
    这样插值算法就不用考虑边界问题了,程序写起来也简单很多! 
      如果对缩放结果的边界像素级精度要求不是太高,我还有一个方案,一个稍微改变的缩放公式:
      Sx=Dx*(SW-1)/DW; Sy=Dy*(SH-1)/DH;  (源图片宽和高:SW>=2;SH>=2)
      证明这个公式不会造成内存访问超界: 
       要求Dx=DW-1时: sx+1=int( (dw-1)/dw*(dw-1) ) +1 <= (sw-1)
          有:  int( (sw-1)*(dw-1)/dw ) <=sw-2
                 (sw-1)*(dw-1)/dw <(sw-1)
                 (dw-1) /dw<1
                 (dw-1) <dw
      比如,按这个公式的一个简单实现: (缩放效果见前面的"二次线性插值(近似公式)"图示)

    void PicZoom_ftBilinear_Common(const TPicRegion& Dst,const TPicRegion& Src)
    {
        
    if (  (0==Dst.width)||(0==Dst.height)
            ||(2>Src.width)||(2>Src.height)) 
    return;

        
    long xrIntFloat_16=((Src.width-1)<<16)/Dst.width; 
        
    long yrIntFloat_16=((Src.height-1)<<16)/Dst.height;

        unsigned 
    long dst_width=Dst.width;
        
    long Src_byte_width=Src.byte_width;
        TARGB32* pDstLine=Dst.pdata;
        
    long srcy_16=0;
        
    for (unsigned long y=0;y<Dst.height;++y)
        {
            unsigned 
    long v_8=(srcy_16 & 0xFFFF)>>8;
            TARGB32* PSrcLineColor= (TARGB32*)((TUInt8*)(Src.pdata)+Src_byte_width*(srcy_16>>16)) ;
            
    long srcx_16=0;
            
    for (unsigned long x=0;x<dst_width;++x)
            {
                TARGB32* PColor0=&PSrcLineColor[srcx_16>>16];
                Bilinear_Fast_Common(PColor0,(TARGB32*)((TUInt8*)(PColor0)+Src_byte_width),(srcx_16 & 0xFFFF)>>8,v_8,&pDstLine[x]);
                srcx_16+=xrIntFloat_16;
            }
            srcy_16+=yrIntFloat_16;
            ((TUInt8*&)pDstLine)+=Dst.byte_width;
        }


    )

    G:利用单指令多数据处理的MMX指令一般都可以加快颜色的运算;在使用MMX改写之前,利用
    32bit寄存器(或变量)来模拟单指令多数据处理;
    数据储存原理:一个颜色数据分量只有一个字节,用2个字节来储存单个颜色分量的计算结果,
    对于很多颜色计算来说精度就够了;那么一个32bit寄存器(或变量)就可以储存2个计算出的
    临时颜色分量;从而达到了单个指令两路数据处理的目的;
    单个指令两路数据处理的计算: 
      乘法: ((0x00AA*a)<<16) | (0x00BB*a) = 0x00AA00BB * a 
        可见只要保证0x00AA*a和0x00BB*a都小于(1<<16)那么乘法可以直接使用无符号数乘法了
      加法: ((0x00AA+0x00CC)<<16) | (0x00BB+0x00DD) = 0x00AA00BB + 0x00CC00DD 
        可见只要0x00AA+0x00CC和0x00BB+0x00DD小于(1<<16)那么加法可以直接使用无符号数加法了
      (移位、减法等稍微复杂一点,因为这里没有用到就不推倒运算公式了)


     

        inline void Bilinear_Fast_Common(TARGB32* PColor0,TARGB32* PColor1,unsigned long u_8,unsigned long v_8,TARGB32* result)
        {
            unsigned 
    long pm3_8=(u_8*v_8)>>8;
            unsigned 
    long pm2_8=u_8-pm3_8;
            unsigned 
    long pm1_8=v_8-pm3_8;
            unsigned 
    long pm0_8=256-pm1_8-pm2_8-pm3_8;

            unsigned 
    long Color=*(unsigned long*)(PColor0);
            unsigned 
    long BR=(Color & 0x00FF00FF)*pm0_8;
            unsigned 
    long GA=((Color & 0xFF00FF00)>>8)*pm0_8;
                          Color=((unsigned 
    long*)(PColor0))[1];
                          GA+=((Color & 0xFF00FF00)>>8)*pm2_8;
                          BR+=(Color & 0x00FF00FF)*pm2_8;
                          Color=*(unsigned 
    long*)(PColor1);
                          GA+=((Color & 0xFF00FF00)>>8)*pm1_8;
                          BR+=(Color & 0x00FF00FF)*pm1_8;
                          Color=((unsigned 
    long*)(PColor1))[1];
                          GA+=((Color & 0xFF00FF00)>>8)*pm3_8;
                          BR+=(Color & 0x00FF00FF)*pm3_8;

            *(unsigned 
    long*)(result)=(GA & 0xFF00FF00)|((BR & 0xFF00FF00)>>8);
        }

        inline 
    void Bilinear_Border_Common(const TPicRegion& pic,const long x_16,const long y_16,TARGB32* result)
        {
            
    long x=(x_16>>16);
            
    long y=(y_16>>16);
            unsigned 
    long u_16=((unsigned short)(x_16));
            unsigned 
    long v_16=((unsigned short)(y_16));

            TARGB32 pixel[4];
            pixel[0]=Pixels_Bound(pic,x,y);
            pixel[1]=Pixels_Bound(pic,x+1,y);
            pixel[2]=Pixels_Bound(pic,x,y+1);
            pixel[3]=Pixels_Bound(pic,x+1,y+1);
            
            Bilinear_Fast_Common(&pixel[0],&pixel[2],u_16>>8,v_16>>8,result);
        }

    void PicZoom_Bilinear_Common(const TPicRegion& Dst,const TPicRegion& Src)
    {
        
    if (  (0==Dst.width)||(0==Dst.height)
            ||(0==Src.width)||(0==Src.height)) 
    return;

        
    long xrIntFloat_16=((Src.width)<<16)/Dst.width+1; 
        
    long yrIntFloat_16=((Src.height)<<16)/Dst.height+1;
        
    const long csDErrorX=-(1<<15)+(xrIntFloat_16>>1);
        
    const long csDErrorY=-(1<<15)+(yrIntFloat_16>>1);

        unsigned 
    long dst_width=Dst.width;

        
    //计算出需要特殊处理的边界
        long border_y0=-csDErrorY/yrIntFloat_16+1;              //y0+y*yr>=0; y0=csDErrorY => y>=-csDErrorY/yr
        if (border_y0>=Dst.height) border_y0=Dst.height;
        
    long border_x0=-csDErrorX/xrIntFloat_16+1;     
        
    if (border_x0>=Dst.width ) border_x0=Dst.width; 
        
    long border_y1=(((Src.height-2)<<16)-csDErrorY)/yrIntFloat_16+1; //y0+y*yr<=(height-2) => y<=(height-2-csDErrorY)/yr
        if (border_y1<border_y0) border_y1=border_y0;
        
    long border_x1=(((Src.width-2)<<16)-csDErrorX)/xrIntFloat_16+1; 
        
    if (border_x1<border_x0) border_x1=border_x0;

        TARGB32* pDstLine=Dst.pdata;
        
    long Src_byte_width=Src.byte_width;
        
    long srcy_16=csDErrorY;
        
    long y;
        
    for (y=0;y<border_y0;++y)
        {
            
    long srcx_16=csDErrorX;
            
    for (unsigned long x=0;x<dst_width;++x)
            {
                Bilinear_Border_Common(Src,srcx_16,srcy_16,&pDstLine[x]); 
    //border
                srcx_16+=xrIntFloat_16;
            }
            srcy_16+=yrIntFloat_16;
            ((TUInt8*&)pDstLine)+=Dst.byte_width;
        }
        
    for (y=border_y0;y<border_y1;++y)
        {
            
    long srcx_16=csDErrorX;
            
    long x;
            
    for (x=0;x<border_x0;++x)
            {
                Bilinear_Border_Common(Src,srcx_16,srcy_16,&pDstLine[x]);
    //border
                srcx_16+=xrIntFloat_16;
            }

            {
                unsigned 
    long v_8=(srcy_16 & 0xFFFF)>>8;
                TARGB32* PSrcLineColor= (TARGB32*)((TUInt8*)(Src.pdata)+Src_byte_width*(srcy_16>>16)) ;
                
    for (unsigned long x=border_x0;x<border_x1;++x)
                {
                    TARGB32* PColor0=&PSrcLineColor[srcx_16>>16];
                    TARGB32* PColor1=(TARGB32*)((TUInt8*)(PColor0)+Src_byte_width);
                    Bilinear_Fast_Common(PColor0,PColor1,(srcx_16 & 0xFFFF)>>8,v_8,&pDstLine[x]);
                    srcx_16+=xrIntFloat_16;
                }
            }

            
    for (x=border_x1;x<dst_width;++x)
            {
                Bilinear_Border_Common(Src,srcx_16,srcy_16,&pDstLine[x]);
    //border
                srcx_16+=xrIntFloat_16;
            }
            srcy_16+=yrIntFloat_16;
            ((TUInt8*&)pDstLine)+=Dst.byte_width;
        }
        
    for (y=border_y1;y<Dst.height;++y)
        {
            
    long srcx_16=csDErrorX;
            
    for (unsigned long x=0;x<dst_width;++x)
            {
                Bilinear_Border_Common(Src,srcx_16,srcy_16,&pDstLine[x]); 
    //border
                srcx_16+=xrIntFloat_16;
            }
            srcy_16+=yrIntFloat_16;
            ((TUInt8*&)pDstLine)+=Dst.byte_width;
        }
    }


    //速度测试:
    //==============================================================================
    // PicZoom_BilInear_Common   65.3 fps

    H:使用MMX指令改写:PicZoom_Bilinear_MMX

        inline void  Bilinear_Fast_MMX(TARGB32* PColor0,TARGB32* PColor1,unsigned long u_8,unsigned long v_8,TARGB32* result)
        {
            asm
            {    
                  MOVD      MM6,v_8
                  MOVD      MM5,u_8
                  mov       edx,PColor0
                  mov       eax,PColor1
                  PXOR      mm7,mm7

                  MOVD         MM2,dword ptr [eax]  
                  MOVD         MM0,dword ptr [eax+4]
                  PUNPCKLWD    MM5,MM5
                  PUNPCKLWD    MM6,MM6
                  MOVD         MM3,dword ptr [edx]  
                  MOVD         MM1,dword ptr [edx+4]
                  PUNPCKLDQ    MM5,MM5 
                  PUNPCKLBW    MM0,MM7
                  PUNPCKLBW    MM1,MM7
                  PUNPCKLBW    MM2,MM7
                  PUNPCKLBW    MM3,MM7
                  PSUBw        MM0,MM2
                  PSUBw        MM1,MM3
                  PSLLw        MM2,8
                  PSLLw        MM3,8
                  PMULlw       MM0,MM5
                  PMULlw       MM1,MM5
                  PUNPCKLDQ    MM6,MM6 
                  PADDw        MM0,MM2
                  PADDw        MM1,MM3

                  PSRLw        MM0,8
                  PSRLw        MM1,8
                  PSUBw        MM0,MM1
                  PSLLw        MM1,8
                  PMULlw       MM0,MM6
                  mov       eax,result
                  PADDw        MM0,MM1

                  PSRLw        MM0,8
                  PACKUSwb     MM0,MM7
                  movd      [eax],MM0 
                  
    //emms
            }
        }

        
    void Bilinear_Border_MMX(const TPicRegion& pic,const long x_16,const long y_16,TARGB32* result)
        {
            
    long x=(x_16>>16);
            
    long y=(y_16>>16);
            unsigned 
    long u_16=((unsigned short)(x_16));
            unsigned 
    long v_16=((unsigned short)(y_16));

            TARGB32 pixel[4];
            pixel[0]=Pixels_Bound(pic,x,y);
            pixel[1]=Pixels_Bound(pic,x+1,y);
            pixel[2]=Pixels_Bound(pic,x,y+1);
            pixel[3]=Pixels_Bound(pic,x+1,y+1);
            
            Bilinear_Fast_MMX(&pixel[0],&pixel[2],u_16>>8,v_16>>8,result);
        }

    void PicZoom_Bilinear_MMX(const TPicRegion& Dst,const TPicRegion& Src)
    {
        
    if (  (0==Dst.width)||(0==Dst.height)
            ||(0==Src.width)||(0==Src.height)) 
    return;

        
    long xrIntFloat_16=((Src.width)<<16)/Dst.width+1; 
        
    long yrIntFloat_16=((Src.height)<<16)/Dst.height+1;
        
    const long csDErrorX=-(1<<15)+(xrIntFloat_16>>1);
        
    const long csDErrorY=-(1<<15)+(yrIntFloat_16>>1);

        unsigned 
    long dst_width=Dst.width;

        
    //计算出需要特殊处理的边界
        long border_y0=-csDErrorY/yrIntFloat_16+1;              //y0+y*yr>=0; y0=csDErrorY => y>=-csDErrorY/yr
        if (border_y0>=Dst.height) border_y0=Dst.height;
        
    long border_x0=-csDErrorX/xrIntFloat_16+1;     
        
    if (border_x0>=Dst.width ) border_x0=Dst.width; 
        
    long border_y1=(((Src.height-2)<<16)-csDErrorY)/yrIntFloat_16+1; //y0+y*yr<=(height-2) => y<=(height-2-csDErrorY)/yr
        if (border_y1<border_y0) border_y1=border_y0;
        
    long border_x1=(((Src.width-2)<<16)-csDErrorX)/xrIntFloat_16+1; 
        
    if (border_x1<border_x0) border_x1=border_x0;

        TARGB32* pDstLine=Dst.pdata;
        
    long Src_byte_width=Src.byte_width;
        
    long srcy_16=csDErrorY;
        
    long y;
        
    for (y=0;y<border_y0;++y)
        {
            
    long srcx_16=csDErrorX;
            
    for (unsigned long x=0;x<dst_width;++x)
            {
                Bilinear_Border_MMX(Src,srcx_16,srcy_16,&pDstLine[x]); 
    //border
                srcx_16+=xrIntFloat_16;
            }
            srcy_16+=yrIntFloat_16;
            ((TUInt8*&)pDstLine)+=Dst.byte_width;
        }
        
    for (y=border_y0;y<border_y1;++y)
        {
            
    long srcx_16=csDErrorX;
            
    long x;
            
    for (x=0;x<border_x0;++x)
            {
                Bilinear_Border_MMX(Src,srcx_16,srcy_16,&pDstLine[x]);
    //border
                srcx_16+=xrIntFloat_16;
            }

            {
                unsigned 
    long v_8=(srcy_16 & 0xFFFF)>>8;
                TARGB32* PSrcLineColor= (TARGB32*)((TUInt8*)(Src.pdata)+Src_byte_width*(srcy_16>>16)) ;
                
    for (unsigned long x=border_x0;x<border_x1;++x)
                {
                    TARGB32* PColor0=&PSrcLineColor[srcx_16>>16];
                    TARGB32* PColor1=(TARGB32*)((TUInt8*)(PColor0)+Src_byte_width);
                    Bilinear_Fast_MMX(PColor0,PColor1,(srcx_16 & 0xFFFF)>>8,v_8,&pDstLine[x]);
                    srcx_16+=xrIntFloat_16;
                }
            }

            
    for (x=border_x1;x<dst_width;++x)
            {
                Bilinear_Border_MMX(Src,srcx_16,srcy_16,&pDstLine[x]);
    //border
                srcx_16+=xrIntFloat_16;
            }
            srcy_16+=yrIntFloat_16;
            ((TUInt8*&)pDstLine)+=Dst.byte_width;
        }
        
    for (y=border_y1;y<Dst.height;++y)
        {
            
    long srcx_16=csDErrorX;
            
    for (unsigned long x=0;x<dst_width;++x)
            {
                Bilinear_Border_MMX(Src,srcx_16,srcy_16,&pDstLine[x]); 
    //border
                srcx_16+=xrIntFloat_16;
            }
            srcy_16+=yrIntFloat_16;
            ((TUInt8*&)pDstLine)+=Dst.byte_width;
        }
        asm emms
    }


    //速度测试:
    //==============================================================================
    // PicZoom_BilInear_MMX 132.9 fps

     

    H' 对BilInear_MMX简单改进:PicZoom_Bilinear_MMX_Ex


     

    void PicZoom_Bilinear_MMX_Ex(const TPicRegion& Dst,const TPicRegion& Src)
    {
        
    if (  (0==Dst.width)||(0==Dst.height)
            ||(0==Src.width)||(0==Src.height)) 
    return;

        
    long xrIntFloat_16=((Src.width)<<16)/Dst.width+1; 
        
    long yrIntFloat_16=((Src.height)<<16)/Dst.height+1;
        
    const long csDErrorX=-(1<<15)+(xrIntFloat_16>>1);
        
    const long csDErrorY=-(1<<15)+(yrIntFloat_16>>1);

        unsigned 
    long dst_width=Dst.width;

        
    //计算出需要特殊处理的边界
        long border_y0=-csDErrorY/yrIntFloat_16+1;              //y0+y*yr>=0; y0=csDErrorY => y>=-csDErrorY/yr
        if (border_y0>=Dst.height) border_y0=Dst.height;
        
    long border_x0=-csDErrorX/xrIntFloat_16+1;     
        
    if (border_x0>=Dst.width ) border_x0=Dst.width; 
        
    long border_y1=(((Src.height-2)<<16)-csDErrorY)/yrIntFloat_16+1; //y0+y*yr<=(height-2) => y<=(height-2-csDErrorY)/yr
        if (border_y1<border_y0) border_y1=border_y0;
        
    long border_x1=(((Src.width-2)<<16)-csDErrorX)/xrIntFloat_16+1; 
        
    if (border_x1<border_x0) border_x1=border_x0;

        TARGB32* pDstLine=Dst.pdata;
        
    long Src_byte_width=Src.byte_width;
        
    long srcy_16=csDErrorY;
        
    long y;
        
    for (y=0;y<border_y0;++y)
        {
            
    long srcx_16=csDErrorX;
            
    for (unsigned long x=0;x<dst_width;++x)
            {
                Bilinear_Border_MMX(Src,srcx_16,srcy_16,&pDstLine[x]); 
    //border
                srcx_16+=xrIntFloat_16;
            }
            srcy_16+=yrIntFloat_16;
            ((TUInt8*&)pDstLine)+=Dst.byte_width;
        }

        
    for (y=border_y0;y<border_y1;++y)
        {
            
    long srcx_16=csDErrorX;
            
    long x;
            
    for (x=0;x<border_x0;++x)
            {
                Bilinear_Border_MMX(Src,srcx_16,srcy_16,&pDstLine[x]);
    //border
                srcx_16+=xrIntFloat_16;
            }

            {
                
    long dst_width_fast=border_x1-border_x0;
                
    if (dst_width_fast>0)
                {
                    unsigned 
    long v_8=(srcy_16 & 0xFFFF)>>8;
                    TARGB32* PSrcLineColor= (TARGB32*)((TUInt8*)(Src.pdata)+Src_byte_width*(srcy_16>>16)) ;
                    TARGB32* PSrcLineColorNext= (TARGB32*)((TUInt8*)(PSrcLineColor)+Src_byte_width) ;
                    TARGB32* pDstLine_Fast=&pDstLine[border_x0];
                    asm
                    {
                          movd         mm6,v_8
                          pxor         mm7,mm7 
    //mm7=0
                          PUNPCKLWD    MM6,MM6
                          PUNPCKLDQ    MM6,MM6
    //mm6=v_8
                        
                          mov       esi,PSrcLineColor
                          mov       ecx,PSrcLineColorNext
                          mov       edx,srcx_16
                          mov       ebx,dst_width_fast
                          mov       edi,pDstLine_Fast
                          lea       edi,[edi+ebx*4]
                          push      ebp
                          mov       ebp,xrIntFloat_16
                          neg       ebx

                    loop_start:

                              mov       eax,edx
                              shl       eax,16
                              shr       eax,24
                              
    //== movzx       eax,dh  //eax=u_8
                              MOVD      MM5,eax
                              mov       eax,edx
                              shr       eax,16     
    //srcx_16>>16

                              MOVD         MM2,dword ptr [ecx+eax*4]  
                              MOVD         MM0,dword ptr [ecx+eax*4+4]
                              PUNPCKLWD    MM5,MM5
                              MOVD         MM3,dword ptr [esi+eax*4]  
                              MOVD         MM1,dword ptr [esi+eax*4+4]
                              PUNPCKLDQ    MM5,MM5 
    //mm5=u_8
                              PUNPCKLBW    MM0,MM7
                              PUNPCKLBW    MM1,MM7
                              PUNPCKLBW    MM2,MM7
                              PUNPCKLBW    MM3,MM7
                              PSUBw        MM0,MM2
                              PSUBw        MM1,MM3
                              PSLLw        MM2,8
                              PSLLw        MM3,8
                              PMULlw       MM0,MM5
                              PMULlw       MM1,MM5
                              PADDw        MM0,MM2
                              PADDw        MM1,MM3

                              PSRLw        MM0,8
                              PSRLw        MM1,8
                              PSUBw        MM0,MM1
                              PSLLw        MM1,8
                              PMULlw       MM0,MM6
                              PADDw        MM0,MM1

                              PSRLw     MM0,8
                              PACKUSwb  MM0,MM7
                              MOVd   dword ptr    [edi+ebx*4],MM0 
    //write DstColor
                                          
                              add       edx,ebp 
    //srcx_16+=xrIntFloat_16
                              inc       ebx
                              jnz       loop_start

                          pop       ebp
                          mov       srcx_16,edx
                    }
                }
            }

            
    for (x=border_x1;x<dst_width;++x)
            {
                Bilinear_Border_MMX(Src,srcx_16,srcy_16,&pDstLine[x]);
    //border
                srcx_16+=xrIntFloat_16;
            }
            srcy_16+=yrIntFloat_16;
            ((TUInt8*&)pDstLine)+=Dst.byte_width;
        }
        
    for (y=border_y1;y<Dst.height;++y)
        {
            
    long srcx_16=csDErrorX;
            
    for (unsigned long x=0;x<dst_width;++x)
            {
                Bilinear_Border_MMX(Src,srcx_16,srcy_16,&pDstLine[x]); 
    //border
                srcx_16+=xrIntFloat_16;
            }
            srcy_16+=yrIntFloat_16;
            ((TUInt8*&)pDstLine)+=Dst.byte_width;
        }
        asm emms
    }


    //速度测试:
    //==============================================================================
    // PicZoom_Bilinear_MMX_Ex 157.0 fps

    I: 把测试成绩放在一起:


    //CPU: AMD64x2 4200+(2.37G)  zoom 800*600 to 1024*768
    //==============================================================================
    // StretchBlt                   232.7 fps   
    // PicZoom3_SSE                 711.7 fps 
    // 
    // PicZoom_BilInear0              8.3 fps
    // PicZoom_BilInear1             17.7 fps
    // PicZoom_BilInear2             43.4 fps
    // PicZoom_BilInear_Common       65.3 fps
    // PicZoom_BilInear_MMX         132.9 fps
    // PicZoom_BilInear_MMX_Ex      157.0 fps

    补充Intel Core2 4400上的测试成绩:


    //CPU: Intel Core2 4400(2.00G)  zoom 800*600 to 1024*768
    //==============================================================================
    // PicZoom3_SSE                1099.7 fps  
    // 
    // PicZoom_BilInear0             10.7 fps
    // PicZoom_BilInear1             24.2 fps
    // PicZoom_BilInear2             54.3 fps
    // PicZoom_BilInear_Common       59.8 fps
    // PicZoom_BilInear_MMX         118.4 fps
    // PicZoom_BilInear_MMX_Ex      142.9 fps

     

     

    三次卷积插值:


    J: 二次线性插值缩放出的图片很多时候让人感觉变得模糊(术语叫低通滤波),特别是在放大
    的时候;使用三次卷积插值来改善插值结果;三次卷积插值考虑映射点周围16个点(4x4)的颜色来
    计算最终的混合颜色,如图;

              
             P(0,0)所在像素为映射的点,加上它周围的15个点,按一定系数混合得到最终输出结果;

             混合公式参见PicZoom_ThreeOrder0的实现;

        插值曲线公式sin(x*PI)/(x*PI),如图:


                                 三次卷积插值曲线sin(x*PI)/(x*PI) (其中PI=3.1415926...)

     

    K:三次卷积插值缩放算法的一个参考实现:PicZoom_ThreeOrder0
      该函数并没有做过多的优化,只是一个简单的浮点实现版本;


     


            inline 
    double SinXDivX(double x) 
            {
                
    //该函数计算插值曲线sin(x*PI)/(x*PI)的值 //PI=3.1415926535897932385; 
                //下面是它的近似拟合表达式
                const float a = -1; //a还可以取 a=-2,-1,-0.75,-0.5等等,起到调节锐化或模糊程度的作用

                
    if (x<0) x=-x; //x=abs(x);
                double x2=x*x;
                
    double x3=x2*x;
                
    if (x<=1)
                  
    return (a+2)*x3 - (a+3)*x2 + 1;
                
    else if (x<=2) 
                  
    return a*x3 - (5*a)*x2 + (8*a)*x - (4*a);
                
    else
                  
    return 0;
            } 

            inline TUInt8 border_color(
    long Color)
            {
                
    if (Color<=0)
                    
    return 0;
                
    else if (Color>=255)
                    
    return 255;
                
    else
                    
    return Color;
            }
            
        
    void ThreeOrder0(const TPicRegion& pic,const float fx,const float fy,TARGB32* result)
        {
            
    long x0=(long)fx; if (x0>fx) --x0; //x0=floor(fx);    
            long y0=(long)fy; if (y0>fy) --y0; //y0=floor(fy);
            float fu=fx-x0;
            
    float fv=fy-y0;

            TARGB32 pixel[16];
            
    long i,j;

            
    for (i=0;i<4;++i)
            {
                
    for (j=0;j<4;++j)
                {
                    
    long x=x0-1+j;
                    
    long y=y0-1+i;
                    pixel[i*4+j]=Pixels_Bound(pic,x,y);
                }
            }

            
    float afu[4],afv[4];
            
    //
            afu[0]=SinXDivX(1+fu);
            afu[1]=SinXDivX(fu);
            afu[2]=SinXDivX(1-fu);
            afu[3]=SinXDivX(2-fu);
            afv[0]=SinXDivX(1+fv);
            afv[1]=SinXDivX(fv);
            afv[2]=SinXDivX(1-fv);
            afv[3]=SinXDivX(2-fv);

            
    float sR=0,sG=0,sB=0,sA=0;
            
    for (i=0;i<4;++i)
            {
                
    float aR=0,aG=0,aB=0,aA=0;
                
    for (long j=0;j<4;++j)
                {
                    aA+=afu[j]*pixel[i*4+j].a;
                    aR+=afu[j]*pixel[i*4+j].r;
                    aG+=afu[j]*pixel[i*4+j].g;
                    aB+=afu[j]*pixel[i*4+j].b;
                }
                sA+=aA*afv[i];
                sR+=aR*afv[i];
                sG+=aG*afv[i];
                sB+=aB*afv[i];
            }

            result->a=border_color((
    long)(sA+0.5));
            result->r=border_color((
    long)(sR+0.5));
            result->g=border_color((
    long)(sG+0.5));
            result->b=border_color((
    long)(sB+0.5));
        }

    void PicZoom_ThreeOrder0(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;
        TARGB32* pDstLine=Dst.pdata;
        
    for (unsigned long y=0;y<Dst.height;++y)
        {
            
    float srcy=(y+0.4999999)*Src.height/Dst.height-0.5;
            
    for (unsigned long x=0;x<dst_width;++x)
            {
                
    float srcx=(x+0.4999999)*Src.width/Dst.width-0.5;
                ThreeOrder0(Src,srcx,srcy,&pDstLine[x]);
            }
            ((TUInt8*&)pDstLine)+=Dst.byte_width;
        }
    }


    //速度测试:
    //==============================================================================
    // PicZoom_ThreeOrder0    3.6 fps

    L: 使用定点数来优化缩放函数;边界和内部分开处理;对SinXDivX做一个查找表;对border_color做一个查找表;


     

        static long SinXDivX_Table_8[(2<<8)+1];
        
    class _CAutoInti_SinXDivX_Table {
        
    private
            
    void _Inti_SinXDivX_Table()
            {
                
    for (long i=0;i<=(2<<8);++i)
                    SinXDivX_Table_8[i]=
    long(0.5+256*SinXDivX(i*(1.0/(256))))*1;
            };
        
    public:
            _CAutoInti_SinXDivX_Table() { _Inti_SinXDivX_Table(); }
        };
        
    static _CAutoInti_SinXDivX_Table __tmp_CAutoInti_SinXDivX_Table;


        
    //颜色查表
        static TUInt8 _color_table[256*3];
        
    static const TUInt8* color_table=&_color_table[256];
        
    class _CAuto_inti_color_table
        {
        
    public:
            _CAuto_inti_color_table() {
                
    for (int i=0;i<256*3;++i)
                    _color_table[i]=border_color(i-256);
            }
        };
        
    static _CAuto_inti_color_table _Auto_inti_color_table;

        
    void ThreeOrder_Fast_Common(const TPicRegion& pic,const long x_16,const long y_16,TARGB32* result)
        {
            unsigned 
    long u_8=(unsigned char)((x_16)>>8);
            unsigned 
    long v_8=(unsigned char)((y_16)>>8);
            
    const TARGB32* pixel=&Pixels(pic,(x_16>>16)-1,(y_16>>16)-1);
            
    long pic_byte_width=pic.byte_width;

            
    long au_8[4],av_8[4];
            
    //
            au_8[0]=SinXDivX_Table_8[(1<<8)+u_8];
            au_8[1]=SinXDivX_Table_8[u_8];
            au_8[2]=SinXDivX_Table_8[(1<<8)-u_8];
            au_8[3]=SinXDivX_Table_8[(2<<8)-u_8];
            av_8[0]=SinXDivX_Table_8[(1<<8)+v_8];
            av_8[1]=SinXDivX_Table_8[v_8];
            av_8[2]=SinXDivX_Table_8[(1<<8)-v_8];
            av_8[3]=SinXDivX_Table_8[(2<<8)-v_8];

            
    long sR=0,sG=0,sB=0,sA=0;
            
    for (long i=0;i<4;++i)
            {
                
    long aA=au_8[0]*pixel[0].a + au_8[1]*pixel[1].a + au_8[2]*pixel[2].a + au_8[3]*pixel[3].a;
                
    long aR=au_8[0]*pixel[0].r + au_8[1]*pixel[1].r + au_8[2]*pixel[2].r + au_8[3]*pixel[3].r;
                
    long aG=au_8[0]*pixel[0].g + au_8[1]*pixel[1].g + au_8[2]*pixel[2].g + au_8[3]*pixel[3].g;
                
    long aB=au_8[0]*pixel[0].b + au_8[1]*pixel[1].b + au_8[2]*pixel[2].b + au_8[3]*pixel[3].b;
                sA+=aA*av_8[i];
                sR+=aR*av_8[i];
                sG+=aG*av_8[i];
                sB+=aB*av_8[i];
                ((TUInt8*&)pixel)+=pic_byte_width;
            }

            result->a=color_table[sA>>16];
            result->r=color_table[sR>>16];
            result->g=color_table[sG>>16];
            result->b=color_table[sB>>16];
        }

        
    void ThreeOrder_Border_Common(const TPicRegion& pic,const long x_16,const long y_16,TARGB32* result)
        {
            
    long x0_sub1=(x_16>>16)-1;
            
    long y0_sub1=(y_16>>16)-1;
            unsigned 
    long u_16_add1=((unsigned short)(x_16))+(1<<16);
            unsigned 
    long v_16_add1=((unsigned short)(y_16))+(1<<16);

            TARGB32 pixel[16];
            
    long i;

            
    for (i=0;i<4;++i)
            {
                
    long y=y0_sub1+i;
                pixel[i*4+0]=Pixels_Bound(pic,x0_sub1+0,y);
                pixel[i*4+1]=Pixels_Bound(pic,x0_sub1+1,y);
                pixel[i*4+2]=Pixels_Bound(pic,x0_sub1+2,y);
                pixel[i*4+3]=Pixels_Bound(pic,x0_sub1+3,y);
            }
            
            TPicRegion npic;
            npic.pdata     =&pixel[0];
            npic.byte_width=4*
    sizeof(TARGB32);
            
    //npic.width     =4;
            //npic.height    =4;
            ThreeOrder_Fast_Common(npic,u_16_add1,v_16_add1,result);
        }

    void PicZoom_ThreeOrder_Common(const TPicRegion& Dst,const TPicRegion& Src)
    {
        
    if (  (0==Dst.width)||(0==Dst.height)
            ||(0==Src.width)||(0==Src.height)) 
    return;

        
    long xrIntFloat_16=((Src.width)<<16)/Dst.width+1; 
        
    long yrIntFloat_16=((Src.height)<<16)/Dst.height+1;
        
    const long csDErrorX=-(1<<15)+(xrIntFloat_16>>1);
        
    const long csDErrorY=-(1<<15)+(yrIntFloat_16>>1);

        unsigned 
    long dst_width=Dst.width;

        
    //计算出需要特殊处理的边界
        long border_y0=((1<<16)-csDErrorY)/yrIntFloat_16+1;//y0+y*yr>=1; y0=csDErrorY => y>=(1-csDErrorY)/yr
        if (border_y0>=Dst.height) border_y0=Dst.height;
        
    long border_x0=((1<<16)-csDErrorX)/xrIntFloat_16+1;
        
    if (border_x0>=Dst.width ) border_x0=Dst.width;
        
    long border_y1=(((Src.height-3)<<16)-csDErrorY)/yrIntFloat_16+1; //y0+y*yr<=(height-3) => y<=(height-3-csDErrorY)/yr
        if (border_y1<border_y0) border_y1=border_y0;
        
    long border_x1=(((Src.width-3)<<16)-csDErrorX)/xrIntFloat_16+1;; 
        
    if (border_x1<border_x0) border_x1=border_x0;

        TARGB32* pDstLine=Dst.pdata;
        
    long srcy_16=csDErrorY;
        
    long y;
        
    for (y=0;y<border_y0;++y)
        {
            
    long srcx_16=csDErrorX;
            
    for (unsigned long x=0;x<dst_width;++x)
            {
                ThreeOrder_Border_Common(Src,srcx_16,srcy_16,&pDstLine[x]); 
    //border
                srcx_16+=xrIntFloat_16;
            }
            srcy_16+=yrIntFloat_16;
            ((TUInt8*&)pDstLine)+=Dst.byte_width;
        }
        
    for (y=border_y0;y<border_y1;++y)
        {
            
    long srcx_16=csDErrorX;
            
    long x;
            
    for (x=0;x<border_x0;++x)
            {
                ThreeOrder_Border_Common(Src,srcx_16,srcy_16,&pDstLine[x]);
    //border
                srcx_16+=xrIntFloat_16;
            }
            
    for (x=border_x0;x<border_x1;++x)
            {
                ThreeOrder_Fast_Common(Src,srcx_16,srcy_16,&pDstLine[x]);
    //fast  !
                srcx_16+=xrIntFloat_16;
            }
            
    for (x=border_x1;x<dst_width;++x)
            {
                ThreeOrder_Border_Common(Src,srcx_16,srcy_16,&pDstLine[x]);
    //border
                srcx_16+=xrIntFloat_16;
            }
            srcy_16+=yrIntFloat_16;
            ((TUInt8*&)pDstLine)+=Dst.byte_width;
        }
        
    for (y=border_y1;y<Dst.height;++y)
        {
            
    long srcx_16=csDErrorX;
            
    for (unsigned long x=0;x<dst_width;++x)
            {
                ThreeOrder_Border_Common(Src,srcx_16,srcy_16,&pDstLine[x]); 
    //border
                srcx_16+=xrIntFloat_16;
            }
            srcy_16+=yrIntFloat_16;
            ((TUInt8*&)pDstLine)+=Dst.byte_width;
        }
    }


    //速度测试:
    //==============================================================================
    // PicZoom_ThreeOrder_Common    16.9 fps

     M: 用MMX来优化ThreeOrder_Common函数:ThreeOrder_MMX


     

        typedef   unsigned long TMMXData32;
        
    static TMMXData32 SinXDivX_Table_MMX[(2<<8)+1];
        
    class _CAutoInti_SinXDivX_Table_MMX {
        
    private
            
    void _Inti_SinXDivX_Table_MMX()
            {
                
    for (long i=0;i<=(2<<8);++i)
                {
                    unsigned 
    short t=long(0.5+(1<<14)*SinXDivX(i*(1.0/(256))));
                    unsigned 
    long tl=t | (((unsigned long)t)<<16);
                    SinXDivX_Table_MMX[i]=tl;
                }
            };
        
    public:
            _CAutoInti_SinXDivX_Table_MMX() { _Inti_SinXDivX_Table_MMX(); }
        };
        
    static _CAutoInti_SinXDivX_Table_MMX __tmp_CAutoInti_SinXDivX_Table_MMX;


        
    void __declspec(naked) _private_ThreeOrder_Fast_MMX()
        {
            asm
            {
                movd        mm1,dword ptr [edx]
                movd        mm2,dword ptr [edx+4]
                movd        mm3,dword ptr [edx+8]
                movd        mm4,dword ptr [edx+12]
                movd        mm5,dword ptr [(offset SinXDivX_Table_MMX)+256*4+eax*4]
                movd        mm6,dword ptr [(offset SinXDivX_Table_MMX)+eax*4]
                punpcklbw   mm1,mm7
                punpcklbw   mm2,mm7
                punpcklwd   mm5,mm5
                punpcklwd   mm6,mm6
                psllw       mm1,7
                psllw       mm2,7
                pmulhw      mm1,mm5
                pmulhw      mm2,mm6
                punpcklbw   mm3,mm7
                punpcklbw   mm4,mm7
                movd        mm5,dword ptr [(offset SinXDivX_Table_MMX)+256*4+ecx*4]
                movd        mm6,dword ptr [(offset SinXDivX_Table_MMX)+512*4+ecx*4]
                punpcklwd   mm5,mm5
                punpcklwd   mm6,mm6
                psllw       mm3,7
                psllw       mm4,7
                pmulhw      mm3,mm5
                pmulhw      mm4,mm6
                paddsw      mm1,mm2
                paddsw      mm3,mm4
                movd        mm6,dword ptr [ebx] 
    //v
                paddsw      mm1,mm3
                punpcklwd   mm6,mm6

                pmulhw      mm1,mm6
                add     edx,esi  
    //+pic.byte_width
                paddsw      mm0,mm1

                ret
            }
        }

        inline 
    void ThreeOrder_Fast_MMX(const TPicRegion& pic,const long x_16,const long y_16,TARGB32* result)
        {
            asm
            {
                mov     ecx,pic
                mov     eax,y_16
                mov     ebx,x_16
                movzx   edi,ah 
    //v_8
                mov     edx,[ecx+TPicRegion::pdata]
                shr     eax,16
                mov     esi,[ecx+TPicRegion::byte_width]
                dec     eax
                movzx   ecx,bh 
    //u_8
                shr     ebx,16
                imul    eax,esi
                lea     edx,[edx+ebx*4-4]
                add     edx,eax 
    //pixel

                mov     eax,ecx
                neg     ecx

                pxor    mm7,mm7  
    //0
                //mov     edx,pixel
                pxor    mm0,mm0  //result=0
                //lea     eax,auv_7

                lea    ebx,[(offset SinXDivX_Table_MMX)+256*4+edi*4]
                call  _private_ThreeOrder_Fast_MMX
                lea    ebx,[(offset SinXDivX_Table_MMX)+edi*4]
                call  _private_ThreeOrder_Fast_MMX
                neg    edi
                lea    ebx,[(offset SinXDivX_Table_MMX)+256*4+edi*4]
                call  _private_ThreeOrder_Fast_MMX
                lea    ebx,[(offset SinXDivX_Table_MMX)+512*4+edi*4]
                call  _private_ThreeOrder_Fast_MMX

                psraw     mm0,3
                mov       eax,result
                packuswb  mm0,mm7
                movd      [eax],mm0
                
    //emms
            }
        }

        
    void ThreeOrder_Border_MMX(const TPicRegion& pic,const long x_16,const long y_16,TARGB32* result)
        {
            unsigned 
    long x0_sub1=(x_16>>16)-1;
            unsigned 
    long y0_sub1=(y_16>>16)-1;
            
    long u_16_add1=((unsigned short)(x_16))+(1<<16);
            
    long v_16_add1=((unsigned short)(y_16))+(1<<16);

            TARGB32 pixel[16];

            
    for (long i=0;i<4;++i)
            {
                
    long y=y0_sub1+i;
                pixel[i*4+0]=Pixels_Bound(pic,x0_sub1  ,y);
                pixel[i*4+1]=Pixels_Bound(pic,x0_sub1+1,y);
                pixel[i*4+2]=Pixels_Bound(pic,x0_sub1+2,y);
                pixel[i*4+3]=Pixels_Bound(pic,x0_sub1+3,y);
            }
            
            TPicRegion npic;
            npic.pdata     =&pixel[0];
            npic.byte_width=4*
    sizeof(TARGB32);
            
    //npic.width     =4;
            //npic.height    =4;
            ThreeOrder_Fast_MMX(npic,u_16_add1,v_16_add1,result);
        }

    void PicZoom_ThreeOrder_MMX(const TPicRegion& Dst,const TPicRegion& Src)
    {
        
    if (  (0==Dst.width)||(0==Dst.height)
            ||(0==Src.width)||(0==Src.height)) 
    return;

        
    long xrIntFloat_16=((Src.width)<<16)/Dst.width+1; 
        
    long yrIntFloat_16=((Src.height)<<16)/Dst.height+1;
        
    const long csDErrorX=-(1<<15)+(xrIntFloat_16>>1);
        
    const long csDErrorY=-(1<<15)+(yrIntFloat_16>>1);

        unsigned 
    long dst_width=Dst.width;

        
    //计算出需要特殊处理的边界
        long border_y0=((1<<16)-csDErrorY)/yrIntFloat_16+1;//y0+y*yr>=1; y0=csDErrorY => y>=(1-csDErrorY)/yr
        if (border_y0>=Dst.height) border_y0=Dst.height;
        
    long border_x0=((1<<16)-csDErrorX)/xrIntFloat_16+1;
        
    if (border_x0>=Dst.width ) border_x0=Dst.width;
        
    long border_y1=(((Src.height-3)<<16)-csDErrorY)/yrIntFloat_16+1; //y0+y*yr<=(height-3) => y<=(height-3-csDErrorY)/yr
        if (border_y1<border_y0) border_y1=border_y0;
        
    long border_x1=(((Src.width-3)<<16)-csDErrorX)/xrIntFloat_16+1;; 
        
    if (border_x1<border_x0) border_x1=border_x0;

        TARGB32* pDstLine=Dst.pdata;
        
    long srcy_16=csDErrorY;
        
    long y;
        
    for (y=0;y<border_y0;++y)
        {
            
    long srcx_16=csDErrorX;
            
    for (unsigned long x=0;x<dst_width;++x)
            {
                ThreeOrder_Border_MMX(Src,srcx_16,srcy_16,&pDstLine[x]); 
    //border
                srcx_16+=xrIntFloat_16;
            }
            srcy_16+=yrIntFloat_16;
            ((TUInt8*&)pDstLine)+=Dst.byte_width;
        }
        
    for (y=border_y0;y<border_y1;++y)
        {
            
    long srcx_16=csDErrorX;
            
    long x;
            
    for (x=0;x<border_x0;++x)
            {
                ThreeOrder_Border_MMX(Src,srcx_16,srcy_16,&pDstLine[x]);
    //border
                srcx_16+=xrIntFloat_16;
            }
            
    for (x=border_x0;x<border_x1;++x)
            {
                ThreeOrder_Fast_MMX(Src,srcx_16,srcy_16,&pDstLine[x]);
    //fast MMX !
                srcx_16+=xrIntFloat_16;
            }
            
    for (x=border_x1;x<dst_width;++x)
            {
                ThreeOrder_Border_MMX(Src,srcx_16,srcy_16,&pDstLine[x]);
    //border
                srcx_16+=xrIntFloat_16;
            }
            srcy_16+=yrIntFloat_16;
            ((TUInt8*&)pDstLine)+=Dst.byte_width;
        }
        
    for (y=border_y1;y<Dst.height;++y)
        {
            
    long srcx_16=csDErrorX;
            
    for (unsigned long x=0;x<dst_width;++x)
            {
                ThreeOrder_Border_MMX(Src,srcx_16,srcy_16,&pDstLine[x]); 
    //border
                srcx_16+=xrIntFloat_16;
            }
            srcy_16+=yrIntFloat_16;
            ((TUInt8*&)pDstLine)+=Dst.byte_width;
        }
        asm emms
    }


    //速度测试:
    //==============================================================================
    // PicZoom_ThreeOrder_MMX   34.3 fps

     

     N:将测试结果放到一起:


    //CPU: AMD64x2 4200+(2.37G)  zoom 800*600 to 1024*768
    //==============================================================================
    // StretchBlt                   232.7 fps   
    // PicZoom3_SSE                 711.7 fps  
    // PicZoom_BilInear_MMX_Ex      157.0 fps
    // 
    // PicZoom_ThreeOrder0            3.6 fps
    // PicZoom_ThreeOrder_Common     16.9 fps
    // PicZoom_ThreeOrder_MMX        34.3 fps

    补充Intel Core2 4400上的测试成绩:


    //CPU: Intel Core2 4400(2.00G)  zoom 800*600 to 1024*768
    //==============================================================================
    // PicZoom3_SSE                1099.7 fps  
    // PicZoom_BilInear_MMX_Ex      142.9 fps
    // 
    // PicZoom_ThreeOrder0            4.2 fps
    // PicZoom_ThreeOrder_Common     17.6 fps
    // PicZoom_ThreeOrder_MMX        34.4 fps

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

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

    正文:

    A:对于前一篇文章中的二次线性插值、三次卷积插值算法,但它们处理缩小到0.5倍以下的
    时候效果就会越来越差;这是因为插值的时候自考虑了附近点的原因;如下图:

                          
               原图          近邻取样 缩放到0.4倍     缩放到0.2倍     缩放到0.1倍

                                                     
                         二次线性插值 缩放到0.4倍     缩放到0.2倍     缩放到0.1倍

                                                     
                         三次卷积插值 缩放到0.4倍     缩放到0.2倍     缩放到0.1倍


        可以看出:当缩小的比例很大的时候,插值算法的效果和近邻取样的效果差不多了:( ;
    一种可行的解决方案就是:缩小时考虑更多的点; 但这种解决方案有很多缺点:函数编写麻烦,
    速度也许会很慢,优化也不容易做;  还有一个方案就是预先建立一个缩放好的大小不同的图片
    列表,每一张图片都是前一张的0.5倍(这种图片列表就是MipMap链),缩放的时候根据需要缩放
    的比例从表中选择一张大小接近的图片来作为缩放的源图片; 该方案的优点:不需要编写新的
    底层缩放算法,直接使用前面优化好的插值算法; 缺点:需要预先建立MipMap链,它需要时间,
    并且它的储存需要多占用原图片的1/3空间(0.5^2+0.5^4+0.5^6+...=1/3);还有一个不太明显
    的小问题,就是在一张图片的连续的比例不同的缩放中,选择会从MipMap的一张源图片跳到另
    一张图片,视觉效果上可能会有一个小的跳跃(我在《魔兽世界》里经常看到这种效应:);一种
    改进方案就是选择MipMap图片的时候,选择出附近的两张图片作为缩放的源图片;对两张图片
    单独进行插值(和原来一致)输出两个值,然后把这两个值线性插值为最终结果;还有一个比较
    大的缺点就是当缩放比例不均匀时(比如x轴放大y轴缩小),缩放效果也不好;
    (当前很多显卡都提供了MipMap纹理和对应的插值方案,OpenGL和DirectX都提供了操作接口)

         
    ("三次线性插值和MipMap链"其实比较简单,这里只给出关键代码或算法)
     

    B: MipMap图片的生成:
         原图片缩放到0.5倍(宽和高都为原图片的1/2),在把0.5倍的图片缩放到0.25倍,....
       直到宽和高都为1个像素,如果有一个长度先到1就保持1; 缩放过程中,可以可采用前面的缩放插值算法;
       如果为了速度可以考虑这样的方案,要求原图片的宽和高必须是2的整数次方的数值,缩放时就可以直接将
       2x2的像素快速合并为一个像素(如果允许原图片宽和高为任何值,可以考虑在合并时引入Alpha通道);

    C: MipMap链图片的储存方案:

                   
                        MipMap链图片示意图         
                       
                               
             可能的一种物理储存方案(我对每张图片加了一个边框)            

    D: 定义MipMap数据结构:
       MipMap数据结构可以定义为一个TPicRegion数组和该数组的大小; 
       (MipMap图片的储存参见上面的图示)
      比如:
     

         #include <vector>
         typedef std::vector
    <TPicRegion> TMipMap;
         
    //其中,第一个元素TMipMap[0]指向原始图片,后面的依次为缩小图片;       

     
    E: MipMap的选择函数和偏好:
        在进行缩放时,根据目标图片缓冲区的大小来动态的选者MipMap中的一幅图片来作为源图片;这就需要一个
    选择函数;比如:

    long SelectBestPicIndex(const TMipMap& mip,const long dstWidth,const long dstHeight)
    {
        
    long oldS=mip[0].width*mip[0].height;
        
    long dstS=dstWidth*dstHeight;
        
    if ( (dstS>=oldS) || (mip.size()==1) )
            
    return 0;
        
    else if (dstS<=1)
            
    return mip.size()-1;
        
    else

            
    return (long)(log(oldS/dstS)*0.5+0.5);
    }

    选择函数可以增加一个偏好参数:
    mip选择偏好:0.5没有偏好,靠近0偏向选择小图片,靠近1偏向选择大图片(质量好一些)

    float public_mip_bias=0.5; //[0..1] 

    long SelectBestPicIndex(const TMipMap& mip,const long dstWidth,const long dstHeight)
    {
        
    long oldS=mip[0].width*mip[0].height;
        
    long dstS=dstWidth*dstHeight;
        
    if ( (dstS>=oldS) || (mip.size()==1) )
            
    return 0;
        
    else if (dstS<=1)
            
    return mip.size()-1;
        
    else

            
    return (long)(log(oldS/dstS)*0.5+public_mip_bias);
    }

     F:利用MipMap后的缩放效果:

                                                     
                     MipMap+近邻取样 缩放到0.4倍     缩放到0.2倍     缩放到0.1倍
                   (利用MipMap做一次近邻取样)

                                                     
                  MipMap+二次线性插值 缩放到0.4倍     缩放到0.2倍     缩放到0.1倍
                  (利用MipMap做一次二次线性插值)

                                                     
                  MipMap+三次卷积插值 缩放到0.4倍     缩放到0.2倍     缩放到0.1倍
                  (利用MipMap做一次三次卷积插值)

       
                       

    G: 在MipMap的两张图片之间插值:
      选择MipMap的时候,同时可以选择相邻的两张MipMap图片;分别进行插值算法后得到两个颜色结果;
    对两个MipMap图片产生的评价值可以作为这两个颜色的插值权重,得到最终的颜色插值结果;优点是
    缩放效果好,避免跳跃;缺点是速度慢:)  

    选择和权重函数的一个可能实现:

    struct TMipWeight {
      
    long  BigMip;
      
    long  SmallMip;
      
    float BigMipWeight;//[0..1]

    };

    TMipWeight SelectBestPicIndexEx(
    const TMipMap& mip,const long dstWidth,const long dstHeight)
    {
        
    long oldS=mip[0].width*mip[0].height;
        
    long dstS=dstWidth*dstHeight;
        TMipWeight result;
        
    if ( (dstS>=oldS) || (mip.size()==1) )
        {
            result.BigMip=0;
            result.SmallMip=0;
            result.BigMipWeight=1.0;
        }
        
    else if (dstS<=1)
        {
            result.BigMip=mip.size()-1;
            result.SmallMip=mip.size()-1;
            result.BigMipWeight=1.0;
        }
        
    else

        {
            
    float bestIndex=log(oldS/dstS)*0.5+0.5; //or + public_mip_bias
            result.BigMip=(long)bestIndex;
            
    if (bestIndex==mip.size()-1)
            {
                result.SmallMip=mip.size()-1;
                result.BigMipWeight=1.0;
            }
            
    else

            {
                result.SmallMip
    =result.BigMip+1;
                result.BigMipWeight=1.0-(bestIndex-result.BigMip);
            }
        }
        
    return
     result;
    }

     

    H:MipMap间插值效果:

                                                    
                  MipMap+两次近邻取样 缩放到0.4倍     缩放到0.2倍     缩放到0.1倍
                  (利用MipMap做两次近邻取样输出两个值,然后线性插值为最终结果)

                                                    
                         三次线性插值 缩放到0.4倍     缩放到0.2倍     缩放到0.1倍
                 (三次线性插值:利用MipMap做两次二次线性插值输出两个值,然后线性插值为最终结果)  
                       
                                                    
               MipMap+两次三次卷积插值 缩放到0.4倍    缩放到0.2倍     缩放到0.1倍
             (利用MipMap做两次三次卷积插值输出两个值,然后线性插值为最终结果)
     

     

    (图像缩放系列终于写完了,计划中写图像任意角度的高质量的快速旋转、Alpha图片混合等,尽请期待:)

     (ps: 思考中的一个图片压缩方法:利用MipMap来压缩图像数据;输入一张图片,然后生成MipMap链,保存相邻之间图片的差(数值差可能很小,很容易找好的算法压缩得很小)和最顶的一张图片(一个点);  解压的时候依次求和就得到原图片了;  该算法为无损压缩,适合于人物风景等过渡比较多的图片的压缩,不太适合线条类等相邻间颜色变化剧烈的图片;)

    展开全文
  • 图像缩放算法——差值算法

    千次阅读 2016-03-08 14:51:13
    摘要:首先给出一个基本的图像缩放算法,然后一步一步的优化其速度和缩放质量;

    转自:http://blog.chinaunix.net/space.php?uid=22915173&do=blog&id=2185545

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

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

    正文:  

      为了便于讨论,这里只处理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: 缩放原理和公式图示:

        缩放后图片             原图片
       (宽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: 缩放效果图:

                     

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

                                       
         原图        缩小图(缩小到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 

     

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

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

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


    速度测试说明:
      只测试内存数据到内存数据的缩放
      测试图片都是800*600缩放到1024*768; fps表示每秒钟的帧数,值越大表示函数越快


    A:近邻取样插值、二次线性插值、三次卷积插值 缩放效果对比

                                       
           原图         近邻取样缩放到0.6倍     近邻取样缩放到1.6倍

                                             
                    二次线性插值缩放到0.6倍   二次线性插值缩放到1.6倍

                                            
                   三次卷积插值缩放到0.6倍   三次卷积插值缩放到1.6倍

          
     原图 近邻取样缩放到8倍 二次线性插值缩放到8倍 三次卷积插值缩放到8倍 二次线性插值(近似公式)


         近邻取样插值缩放简单、速度快,但很多时候缩放出的图片质量比较差(特别是对于人物、景色等),
    图片的缩放有比较明显的锯齿;使用二次或更高次插值有利于改善缩放效果;


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

    #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];
    }


     


    二次线性插值缩放:

    C: 二次线性插值缩放原理和公式图示:

            

                 缩放后图片                 原图片
                (宽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

      聚焦看看(Sx,Sy)坐标点(Sx,Sy为浮点数)附近的情况;


             


      对于近邻取样插值的缩放算法,直接取Color0颜色作为缩放后点的颜色;
    二次线性插值需要考虑(Sx,Sy)坐标点周围的4个颜色值Color0\Color1\Color2\Color3,
    把(Sx,Sy)到A\B\C\D坐标点的距离作为系数来把4个颜色混合出缩放后点的颜色;
    ( u=Sx-floor(Sx); v=Sy-floor(Sy); 说明:floor函数的返回值为小于等于参数的最大整数 )  
      二次线性插值公式为:
     tmpColor0=Color0*(1-u) + Color2*u;
     tmpColor1=Color1*(1-u) + Color3*u;
            DstColor =tmpColor0*(1-v) + tmpColor2*v;

      展开公式为:
            pm0=(1-u)*(1-v);
            pm1=v*(1-u);
            pm2=u*(1-v);
            pm3=u*v;
      则颜色混合公式为:
            DstColor = Color0*pm0 + Color1*pm1 + Color2*pm2 + Color3*pm3;

    参数函数图示:

          

                                                  二次线性插值函数图示

    对于上面的公式,它将图片向右下各移动了半个像素,需要对此做一个修正;
      =>   Sx=(Dx+0.5)*SW/DW-0.5; Sy=(Dy+0.5)*SH/DH-0.5;
    而实际的程序,还需要考虑到边界(访问源图片可能超界)对于算法的影响,边界的处理可能有各种
    方案(不处理边界或边界回绕或边界饱和或边界映射或用背景颜色混合等;文章中默认使用边界饱和来处理超界);
    比如:边界饱和函数: 

    //访问一个点的函数,(x,y)坐标可能超出图片边界; //边界处理模式:边界饱和
    inline TARGB32 Pixels_Bound(const TPicRegion& pic,long x,long y)
    {
        
    //assert((pic.width>0)&&(pic.height>0));
        bool IsInPic=true;
        
    if (x<0) {x=0; IsInPic=false; } else if (x>=pic.width ) {x=pic.width -1; IsInPic=false; }
        
    if (y<0) {y=0; IsInPic=false; } else if (y>=pic.height) {y=pic.height-1; IsInPic=false; }
        TARGB32 result=Pixels(pic,x,y);
        
    if (!IsInPic) result.a=0;
        
    return result;
    }


     

    D: 二次线性插值缩放算法的一个参考实现:PicZoom_BilInear0
      该函数并没有做什么优化,只是一个简单的浮点实现版本;


     

        inline void Bilinear0(const TPicRegion& pic,float fx,float fy,TARGB32* result)
        {
            
    long x=(long)fx; if (x>fx) --x; //x=floor(fx);    
            long y=(long)fy; if (y>fy) --y; //y=floor(fy);
            
            TARGB32 Color0=Pixels_Bound(pic,x,y);
            TARGB32 Color2=Pixels_Bound(pic,x+1,y);
            TARGB32 Color1=Pixels_Bound(pic,x,y+1);
            TARGB32 Color3=Pixels_Bound(pic,x+1,y+1);

            
    float u=fx-x;
            
    float v=fy-y;
            
    float pm3=u*v;
            
    float pm2=u*(1-v);
            
    float pm1=v*(1-u);
            
    float pm0=(1-u)*(1-v);

            result->a=(pm0*Color0.a+pm1*Color1.a+pm2*Color2.a+pm3*Color3.a);
            result->r=(pm0*Color0.r+pm1*Color1.r+pm2*Color2.r+pm3*Color3.r);
            result->g=(pm0*Color0.g+pm1*Color1.g+pm2*Color2.g+pm3*Color3.g);
            result->b=(pm0*Color0.b+pm1*Color1.b+pm2*Color2.b+pm3*Color3.b);
        }

    void PicZoom_Bilinear0(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;
        TARGB32* pDstLine=Dst.pdata;
        
    for (unsigned long y=0;y<Dst.height;++y)
        {
            
    float srcy=(y+0.4999999)*Src.height/Dst.height-0.5;
            
    for (unsigned long x=0;x<dst_width;++x)
            {
                
    float srcx=(x+0.4999999)*Src.width/Dst.width-0.5;
                Bilinear0(Src,srcx,srcy,&pDstLine[x]);
            }
            ((TUInt8*&)pDstLine)+=Dst.byte_width;
        }
    }


    //速度测试:
    //==============================================================================
    // PicZoom_BilInear0      8.3 fps

     

    E: 把PicZoom_BilInear0的浮点计算改写为定点数实现:PicZoom_BilInear1

        inline void Bilinear1(const TPicRegion& pic,const long x_16,const long y_16,TARGB32* result)
        {
            
    long x=x_16>>16;
            
    long y=y_16>>16;
            TARGB32 Color0=Pixels_Bound(pic,x,y);
            TARGB32 Color2=Pixels_Bound(pic,x+1,y);
            TARGB32 Color1=Pixels_Bound(pic,x,y+1);
            TARGB32 Color3=Pixels_Bound(pic,x+1,y+1);

            unsigned 
    long u_8=(x_16 & 0xFFFF)>>8;
            unsigned 
    long v_8=(y_16 & 0xFFFF)>>8;
            unsigned 
    long pm3_16=(u_8*v_8);
            unsigned 
    long pm2_16=(u_8*(unsigned long)(256-v_8));
            unsigned 
    long pm1_16=(v_8*(unsigned long)(256-u_8));
            unsigned 
    long pm0_16=((256-u_8)*(256-v_8));

            result->a=((pm0_16*Color0.a+pm1_16*Color1.a+pm2_16*Color2.a+pm3_16*Color3.a)>>16);
            result->r=((pm0_16*Color0.r+pm1_16*Color1.r+pm2_16*Color2.r+pm3_16*Color3.r)>>16);
            result->g=((pm0_16*Color0.g+pm1_16*Color1.g+pm2_16*Color2.g+pm3_16*Color3.g)>>16);
            result->b=((pm0_16*Color0.b+pm1_16*Color1.b+pm2_16*Color2.b+pm3_16*Color3.b)>>16);
        }

    void PicZoom_Bilinear1(const TPicRegion& Dst,const TPicRegion& Src)
    {
        
    if (  (0==Dst.width)||(0==Dst.height)
            ||(0==Src.width)||(0==Src.height)) 
    return;

        
    long xrIntFloat_16=((Src.width)<<16)/Dst.width+1; 
        
    long yrIntFloat_16=((Src.height)<<16)/Dst.height+1;
        
    const long csDErrorX=-(1<<15)+(xrIntFloat_16>>1);
        
    const long csDErrorY=-(1<<15)+(yrIntFloat_16>>1);

        unsigned 
    long dst_width=Dst.width;

        TARGB32* pDstLine=Dst.pdata;
        
    long srcy_16=csDErrorY;
        
    long y;
        
    for (y=0;y<Dst.height;++y)
        {
            
    long srcx_16=csDErrorX;
            
    for (unsigned long x=0;x<dst_width;++x)
            {
                Bilinear1(Src,srcx_16,srcy_16,&pDstLine[x]); 
    //border
                srcx_16+=xrIntFloat_16;
            }
            srcy_16+=yrIntFloat_16;
            ((TUInt8*&)pDstLine)+=Dst.byte_width;
        }



    //速度测试:
    //==============================================================================
    // PicZoom_BilInear1     17.7 fps

    F: 二次线性插值需要考略边界访问超界的问题,我们可以将边界区域和内部区域分开处理,这样就可以优化内部的插值实现函数了:比如不需要判断访问超界、减少颜色数据复制、减少一些不必要的重复坐标计算等等


     


     

        inline void Bilinear2_Fast(TARGB32* PColor0,TARGB32* PColor1,unsigned long u_8,unsigned long v_8,TARGB32* result)
        {
            unsigned 
    long pm3_16=u_8*v_8;
            unsigned 
    long pm2_16=(u_8<<8)-pm3_16;
            unsigned 
    long pm1_16=(v_8<<8)-pm3_16;
            unsigned 
    long pm0_16=(1<<16)-pm1_16-pm2_16-pm3_16;
       
            result->a=((pm0_16*PColor0[0].a+pm2_16*PColor0[1].a+pm1_16*PColor1[0].a+pm3_16*PColor1[1].a)>>16);
            result->r=((pm0_16*PColor0[0].r+pm2_16*PColor0[1].r+pm1_16*PColor1[0].r+pm3_16*PColor1[1].r)>>16);
            result->g=((pm0_16*PColor0[0].g+pm2_16*PColor0[1].g+pm1_16*PColor1[0].g+pm3_16*PColor1[1].g)>>16);
            result->b=((pm0_16*PColor0[0].b+pm2_16*PColor0[1].b+pm1_16*PColor1[0].b+pm3_16*PColor1[1].b)>>16);
        }

        inline 
    void Bilinear2_Border(const TPicRegion& pic,const long x_16,const long y_16,TARGB32* result)
        {
            
    long x=(x_16>>16);
            
    long y=(y_16>>16);
            unsigned 
    long u_16=((unsigned short)(x_16));
            unsigned 
    long v_16=((unsigned short)(y_16));

            TARGB32 pixel[4];
            pixel[0]=Pixels_Bound(pic,x,y);
            pixel[1]=Pixels_Bound(pic,x+1,y);
            pixel[2]=Pixels_Bound(pic,x,y+1);
            pixel[3]=Pixels_Bound(pic,x+1,y+1);
            
            Bilinear2_Fast(&pixel[0],&pixel[2],u_16>>8,v_16>>8,result);
        }

    void PicZoom_Bilinear2(const TPicRegion& Dst,const TPicRegion& Src)
    {
        
    if (  (0==Dst.width)||(0==Dst.height)
            ||(0==Src.width)||(0==Src.height)) 
    return;

        
    long xrIntFloat_16=((Src.width)<<16)/Dst.width+1; 
        
    long yrIntFloat_16=((Src.height)<<16)/Dst.height+1;
        
    const long csDErrorX=-(1<<15)+(xrIntFloat_16>>1);
        
    const long csDErrorY=-(1<<15)+(yrIntFloat_16>>1);

        unsigned 
    long dst_width=Dst.width;

        
    //计算出需要特殊处理的边界
        long border_y0=-csDErrorY/yrIntFloat_16+1;              //y0+y*yr>=0; y0=csDErrorY => y>=-csDErrorY/yr
        if (border_y0>=Dst.height) border_y0=Dst.height;
        
    long border_x0=-csDErrorX/xrIntFloat_16+1;     
        
    if (border_x0>=Dst.width ) border_x0=Dst.width; 
        
    long border_y1=(((Src.height-2)<<16)-csDErrorY)/yrIntFloat_16+1; //y0+y*yr<=(height-2) => y<=(height-2-csDErrorY)/yr
        if (border_y1<border_y0) border_y1=border_y0;
        
    long border_x1=(((Src.width-2)<<16)-csDErrorX)/xrIntFloat_16+1; 
        
    if (border_x1<border_x0) border_x1=border_x0;

        TARGB32* pDstLine=Dst.pdata;
        
    long Src_byte_width=Src.byte_width;
        
    long srcy_16=csDErrorY;
        
    long y;
        
    for (y=0;y<border_y0;++y)
        {
            
    long srcx_16=csDErrorX;
            
    for (unsigned long x=0;x<dst_width;++x)
            {
                Bilinear2_Border(Src,srcx_16,srcy_16,&pDstLine[x]); 
    //border
                srcx_16+=xrIntFloat_16;
            }
            srcy_16+=yrIntFloat_16;
            ((TUInt8*&)pDstLine)+=Dst.byte_width;
        }
        
    for (y=border_y0;y<border_y1;++y)
        {
            
    long srcx_16=csDErrorX;
            
    long x;
            
    for (x=0;x<border_x0;++x)
            {
                Bilinear2_Border(Src,srcx_16,srcy_16,&pDstLine[x]);
    //border
                srcx_16+=xrIntFloat_16;
            }

            {
                unsigned 
    long v_8=(srcy_16 & 0xFFFF)>>8;
                TARGB32* PSrcLineColor= (TARGB32*)((TUInt8*)(Src.pdata)+Src_byte_width*(srcy_16>>16)) ;
                
    for (unsigned long x=border_x0;x<border_x1;++x)
                {
                    TARGB32* PColor0=&PSrcLineColor[srcx_16>>16];
                    TARGB32* PColor1=(TARGB32*)((TUInt8*)(PColor0)+Src_byte_width);
                    Bilinear2_Fast(PColor0,PColor1,(srcx_16 & 0xFFFF)>>8,v_8,&pDstLine[x]);
                    srcx_16+=xrIntFloat_16;
                }
            }

            
    for (x=border_x1;x<dst_width;++x)
            {
                Bilinear2_Border(Src,srcx_16,srcy_16,&pDstLine[x]);
    //border
                srcx_16+=xrIntFloat_16;
            }
            srcy_16+=yrIntFloat_16;
            ((TUInt8*&)pDstLine)+=Dst.byte_width;
        }
        
    for (y=border_y1;y<Dst.height;++y)
        {
            
    long srcx_16=csDErrorX;
            
    for (unsigned long x=0;x<dst_width;++x)
            {
                Bilinear2_Border(Src,srcx_16,srcy_16,&pDstLine[x]); 
    //border
                srcx_16+=xrIntFloat_16;
            }
            srcy_16+=yrIntFloat_16;
            ((TUInt8*&)pDstLine)+=Dst.byte_width;
        }
    }


    //速度测试:
    //==============================================================================
    // PicZoom_BilInear2     43.4 fps

    (F'补充:
      如果不想处理边界访问超界问题,可以考虑扩大源图片的尺寸,加一个边框 (“哨兵”优化);
    这样插值算法就不用考虑边界问题了,程序写起来也简单很多! 
      如果对缩放结果的边界像素级精度要求不是太高,我还有一个方案,一个稍微改变的缩放公式:
      Sx=Dx*(SW-1)/DW; Sy=Dy*(SH-1)/DH;  (源图片宽和高:SW>=2;SH>=2)
      证明这个公式不会造成内存访问超界: 
       要求Dx=DW-1时: sx+1=int( (dw-1)/dw*(dw-1) ) +1 <= (sw-1)
          有:  int( (sw-1)*(dw-1)/dw ) <=sw-2
                 (sw-1)*(dw-1)/dw <(sw-1)
                 (dw-1) /dw<1
                 (dw-1) <dw
      比如,按这个公式的一个简单实现: (缩放效果见前面的"二次线性插值(近似公式)"图示)

    void PicZoom_ftBilinear_Common(const TPicRegion& Dst,const TPicRegion& Src)
    {
        
    if (  (0==Dst.width)||(0==Dst.height)
            ||(2>Src.width)||(2>Src.height)) 
    return;

        
    long xrIntFloat_16=((Src.width-1)<<16)/Dst.width; 
        
    long yrIntFloat_16=((Src.height-1)<<16)/Dst.height;

        unsigned 
    long dst_width=Dst.width;
        
    long Src_byte_width=Src.byte_width;
        TARGB32* pDstLine=Dst.pdata;
        
    long srcy_16=0;
        
    for (unsigned long y=0;y<Dst.height;++y)
        {
            unsigned 
    long v_8=(srcy_16 & 0xFFFF)>>8;
            TARGB32* PSrcLineColor= (TARGB32*)((TUInt8*)(Src.pdata)+Src_byte_width*(srcy_16>>16)) ;
            
    long srcx_16=0;
            
    for (unsigned long x=0;x<dst_width;++x)
            {
                TARGB32* PColor0=&PSrcLineColor[srcx_16>>16];
                Bilinear_Fast_Common(PColor0,(TARGB32*)((TUInt8*)(PColor0)+Src_byte_width),(srcx_16 & 0xFFFF)>>8,v_8,&pDstLine[x]);
                srcx_16+=xrIntFloat_16;
            }
            srcy_16+=yrIntFloat_16;
            ((TUInt8*&)pDstLine)+=Dst.byte_width;
        }


    )

    G:利用单指令多数据处理的MMX指令一般都可以加快颜色的运算;在使用MMX改写之前,利用
    32bit寄存器(或变量)来模拟单指令多数据处理;
    数据储存原理:一个颜色数据分量只有一个字节,用2个字节来储存单个颜色分量的计算结果,
    对于很多颜色计算来说精度就够了;那么一个32bit寄存器(或变量)就可以储存2个计算出的
    临时颜色分量;从而达到了单个指令两路数据处理的目的;
    单个指令两路数据处理的计算: 
      乘法: ((0x00AA*a)<<16) | (0x00BB*a) = 0x00AA00BB * a 
        可见只要保证0x00AA*a和0x00BB*a都小于(1<<16)那么乘法可以直接使用无符号数乘法了
      加法: ((0x00AA+0x00CC)<<16) | (0x00BB+0x00DD) = 0x00AA00BB + 0x00CC00DD 
        可见只要0x00AA+0x00CC和0x00BB+0x00DD小于(1<<16)那么加法可以直接使用无符号数加法了
      (移位、减法等稍微复杂一点,因为这里没有用到就不推倒运算公式了)


     

        inline void Bilinear_Fast_Common(TARGB32* PColor0,TARGB32* PColor1,unsigned long u_8,unsigned long v_8,TARGB32* result)
        {
            unsigned 
    long pm3_8=(u_8*v_8)>>8;
            unsigned 
    long pm2_8=u_8-pm3_8;
            unsigned 
    long pm1_8=v_8-pm3_8;
            unsigned 
    long pm0_8=256-pm1_8-pm2_8-pm3_8;

            unsigned 
    long Color=*(unsigned long*)(PColor0);
            unsigned 
    long BR=(Color & 0x00FF00FF)*pm0_8;
            unsigned 
    long GA=((Color & 0xFF00FF00)>>8)*pm0_8;
                          Color=((unsigned 
    long*)(PColor0))[1];
                          GA+=((Color & 0xFF00FF00)>>8)*pm2_8;
                          BR+=(Color & 0x00FF00FF)*pm2_8;
                          Color=*(unsigned 
    long*)(PColor1);
                          GA+=((Color & 0xFF00FF00)>>8)*pm1_8;
                          BR+=(Color & 0x00FF00FF)*pm1_8;
                          Color=((unsigned 
    long*)(PColor1))[1];
                          GA+=((Color & 0xFF00FF00)>>8)*pm3_8;
                          BR+=(Color & 0x00FF00FF)*pm3_8;

            *(unsigned 
    long*)(result)=(GA & 0xFF00FF00)|((BR & 0xFF00FF00)>>8);
        }

        inline 
    void Bilinear_Border_Common(const TPicRegion& pic,const long x_16,const long y_16,TARGB32* result)
        {
            
    long x=(x_16>>16);
            
    long y=(y_16>>16);
            unsigned 
    long u_16=((unsigned short)(x_16));
            unsigned 
    long v_16=((unsigned short)(y_16));

            TARGB32 pixel[4];
            pixel[0]=Pixels_Bound(pic,x,y);
            pixel[1]=Pixels_Bound(pic,x+1,y);
            pixel[2]=Pixels_Bound(pic,x,y+1);
            pixel[3]=Pixels_Bound(pic,x+1,y+1);
            
            Bilinear_Fast_Common(&pixel[0],&pixel[2],u_16>>8,v_16>>8,result);
        }

    void PicZoom_Bilinear_Common(const TPicRegion& Dst,const TPicRegion& Src)
    {
        
    if (  (0==Dst.width)||(0==Dst.height)
            ||(0==Src.width)||(0==Src.height)) 
    return;

        
    long xrIntFloat_16=((Src.width)<<16)/Dst.width+1; 
        
    long yrIntFloat_16=((Src.height)<<16)/Dst.height+1;
        
    const long csDErrorX=-(1<<15)+(xrIntFloat_16>>1);
        
    const long csDErrorY=-(1<<15)+(yrIntFloat_16>>1);

        unsigned 
    long dst_width=Dst.width;

        
    //计算出需要特殊处理的边界
        long border_y0=-csDErrorY/yrIntFloat_16+1;              //y0+y*yr>=0; y0=csDErrorY => y>=-csDErrorY/yr
        if (border_y0>=Dst.height) border_y0=Dst.height;
        
    long border_x0=-csDErrorX/xrIntFloat_16+1;     
        
    if (border_x0>=Dst.width ) border_x0=Dst.width; 
        
    long border_y1=(((Src.height-2)<<16)-csDErrorY)/yrIntFloat_16+1; //y0+y*yr<=(height-2) => y<=(height-2-csDErrorY)/yr
        if (border_y1<border_y0) border_y1=border_y0;
        
    long border_x1=(((Src.width-2)<<16)-csDErrorX)/xrIntFloat_16+1; 
        
    if (border_x1<border_x0) border_x1=border_x0;

        TARGB32* pDstLine=Dst.pdata;
        
    long Src_byte_width=Src.byte_width;
        
    long srcy_16=csDErrorY;
        
    long y;
        
    for (y=0;y<border_y0;++y)
        {
            
    long srcx_16=csDErrorX;
            
    for (unsigned long x=0;x<dst_width;++x)
            {
                Bilinear_Border_Common(Src,srcx_16,srcy_16,&pDstLine[x]); 
    //border
                srcx_16+=xrIntFloat_16;
            }
            srcy_16+=yrIntFloat_16;
            ((TUInt8*&)pDstLine)+=Dst.byte_width;
        }
        
    for (y=border_y0;y<border_y1;++y)
        {
            
    long srcx_16=csDErrorX;
            
    long x;
            
    for (x=0;x<border_x0;++x)
            {
                Bilinear_Border_Common(Src,srcx_16,srcy_16,&pDstLine[x]);
    //border
                srcx_16+=xrIntFloat_16;
            }

            {
                unsigned 
    long v_8=(srcy_16 & 0xFFFF)>>8;
                TARGB32* PSrcLineColor= (TARGB32*)((TUInt8*)(Src.pdata)+Src_byte_width*(srcy_16>>16)) ;
                
    for (unsigned long x=border_x0;x<border_x1;++x)
                {
                    TARGB32* PColor0=&PSrcLineColor[srcx_16>>16];
                    TARGB32* PColor1=(TARGB32*)((TUInt8*)(PColor0)+Src_byte_width);
                    Bilinear_Fast_Common(PColor0,PColor1,(srcx_16 & 0xFFFF)>>8,v_8,&pDstLine[x]);
                    srcx_16+=xrIntFloat_16;
                }
            }

            
    for (x=border_x1;x<dst_width;++x)
            {
                Bilinear_Border_Common(Src,srcx_16,srcy_16,&pDstLine[x]);
    //border
                srcx_16+=xrIntFloat_16;
            }
            srcy_16+=yrIntFloat_16;
            ((TUInt8*&)pDstLine)+=Dst.byte_width;
        }
        
    for (y=border_y1;y<Dst.height;++y)
        {
            
    long srcx_16=csDErrorX;
            
    for (unsigned long x=0;x<dst_width;++x)
            {
                Bilinear_Border_Common(Src,srcx_16,srcy_16,&pDstLine[x]); 
    //border
                srcx_16+=xrIntFloat_16;
            }
            srcy_16+=yrIntFloat_16;
            ((TUInt8*&)pDstLine)+=Dst.byte_width;
        }
    }


    //速度测试:
    //==============================================================================
    // PicZoom_BilInear_Common   65.3 fps

    H:使用MMX指令改写:PicZoom_Bilinear_MMX

        inline void  Bilinear_Fast_MMX(TARGB32* PColor0,TARGB32* PColor1,unsigned long u_8,unsigned long v_8,TARGB32* result)
        {
            asm
            {    
                  MOVD      MM6,v_8
                  MOVD      MM5,u_8
                  mov       edx,PColor0
                  mov       eax,PColor1
                  PXOR      mm7,mm7

                  MOVD         MM2,dword ptr [eax]  
                  MOVD         MM0,dword ptr [eax+4]
                  PUNPCKLWD    MM5,MM5
                  PUNPCKLWD    MM6,MM6
                  MOVD         MM3,dword ptr [edx]  
                  MOVD         MM1,dword ptr [edx+4]
                  PUNPCKLDQ    MM5,MM5 
                  PUNPCKLBW    MM0,MM7
                  PUNPCKLBW    MM1,MM7
                  PUNPCKLBW    MM2,MM7
                  PUNPCKLBW    MM3,MM7
                  PSUBw        MM0,MM2
                  PSUBw        MM1,MM3
                  PSLLw        MM2,8
                  PSLLw        MM3,8
                  PMULlw       MM0,MM5
                  PMULlw       MM1,MM5
                  PUNPCKLDQ    MM6,MM6 
                  PADDw        MM0,MM2
                  PADDw        MM1,MM3

                  PSRLw        MM0,8
                  PSRLw        MM1,8
                  PSUBw        MM0,MM1
                  PSLLw        MM1,8
                  PMULlw       MM0,MM6
                  mov       eax,result
                  PADDw        MM0,MM1

                  PSRLw        MM0,8
                  PACKUSwb     MM0,MM7
                  movd      [eax],MM0 
                  
    //emms
            }
        }

        
    void Bilinear_Border_MMX(const TPicRegion& pic,const long x_16,const long y_16,TARGB32* result)
        {
            
    long x=(x_16>>16);
            
    long y=(y_16>>16);
            unsigned 
    long u_16=((unsigned short)(x_16));
            unsigned 
    long v_16=((unsigned short)(y_16));

            TARGB32 pixel[4];
            pixel[0]=Pixels_Bound(pic,x,y);
            pixel[1]=Pixels_Bound(pic,x+1,y);
            pixel[2]=Pixels_Bound(pic,x,y+1);
            pixel[3]=Pixels_Bound(pic,x+1,y+1);
            
            Bilinear_Fast_MMX(&pixel[0],&pixel[2],u_16>>8,v_16>>8,result);
        }

    void PicZoom_Bilinear_MMX(const TPicRegion& Dst,const TPicRegion& Src)
    {
        
    if (  (0==Dst.width)||(0==Dst.height)
            ||(0==Src.width)||(0==Src.height)) 
    return;

        
    long xrIntFloat_16=((Src.width)<<16)/Dst.width+1; 
        
    long yrIntFloat_16=((Src.height)<<16)/Dst.height+1;
        
    const long csDErrorX=-(1<<15)+(xrIntFloat_16>>1);
        
    const long csDErrorY=-(1<<15)+(yrIntFloat_16>>1);

        unsigned 
    long dst_width=Dst.width;

        
    //计算出需要特殊处理的边界
        long border_y0=-csDErrorY/yrIntFloat_16+1;              //y0+y*yr>=0; y0=csDErrorY => y>=-csDErrorY/yr
        if (border_y0>=Dst.height) border_y0=Dst.height;
        
    long border_x0=-csDErrorX/xrIntFloat_16+1;     
        
    if (border_x0>=Dst.width ) border_x0=Dst.width; 
        
    long border_y1=(((Src.height-2)<<16)-csDErrorY)/yrIntFloat_16+1; //y0+y*yr<=(height-2) => y<=(height-2-csDErrorY)/yr
        if (border_y1<border_y0) border_y1=border_y0;
        
    long border_x1=(((Src.width-2)<<16)-csDErrorX)/xrIntFloat_16+1; 
        
    if (border_x1<border_x0) border_x1=border_x0;

        TARGB32* pDstLine=Dst.pdata;
        
    long Src_byte_width=Src.byte_width;
        
    long srcy_16=csDErrorY;
        
    long y;
        
    for (y=0;y<border_y0;++y)
        {
            
    long srcx_16=csDErrorX;
            
    for (unsigned long x=0;x<dst_width;++x)
            {
                Bilinear_Border_MMX(Src,srcx_16,srcy_16,&pDstLine[x]); 
    //border
                srcx_16+=xrIntFloat_16;
            }
            srcy_16+=yrIntFloat_16;
            ((TUInt8*&)pDstLine)+=Dst.byte_width;
        }
        
    for (y=border_y0;y<border_y1;++y)
        {
            
    long srcx_16=csDErrorX;
            
    long x;
            
    for (x=0;x<border_x0;++x)
            {
                Bilinear_Border_MMX(Src,srcx_16,srcy_16,&pDstLine[x]);
    //border
                srcx_16+=xrIntFloat_16;
            }

            {
                unsigned 
    long v_8=(srcy_16 & 0xFFFF)>>8;
                TARGB32* PSrcLineColor= (TARGB32*)((TUInt8*)(Src.pdata)+Src_byte_width*(srcy_16>>16)) ;
                
    for (unsigned long x=border_x0;x<border_x1;++x)
                {
                    TARGB32* PColor0=&PSrcLineColor[srcx_16>>16];
                    TARGB32* PColor1=(TARGB32*)((TUInt8*)(PColor0)+Src_byte_width);
                    Bilinear_Fast_MMX(PColor0,PColor1,(srcx_16 & 0xFFFF)>>8,v_8,&pDstLine[x]);
                    srcx_16+=xrIntFloat_16;
                }
            }

            
    for (x=border_x1;x<dst_width;++x)
            {
                Bilinear_Border_MMX(Src,srcx_16,srcy_16,&pDstLine[x]);
    //border
                srcx_16+=xrIntFloat_16;
            }
            srcy_16+=yrIntFloat_16;
            ((TUInt8*&)pDstLine)+=Dst.byte_width;
        }
        
    for (y=border_y1;y<Dst.height;++y)
        {
            
    long srcx_16=csDErrorX;
            
    for (unsigned long x=0;x<dst_width;++x)
            {
                Bilinear_Border_MMX(Src,srcx_16,srcy_16,&pDstLine[x]); 
    //border
                srcx_16+=xrIntFloat_16;
            }
            srcy_16+=yrIntFloat_16;
            ((TUInt8*&)pDstLine)+=Dst.byte_width;
        }
        asm emms
    }


    //速度测试:
    //==============================================================================
    // PicZoom_BilInear_MMX 132.9 fps

     

    H' 对BilInear_MMX简单改进:PicZoom_Bilinear_MMX_Ex


     

    void PicZoom_Bilinear_MMX_Ex(const TPicRegion& Dst,const TPicRegion& Src)
    {
        
    if (  (0==Dst.width)||(0==Dst.height)
            ||(0==Src.width)||(0==Src.height)) 
    return;

        
    long xrIntFloat_16=((Src.width)<<16)/Dst.width+1; 
        
    long yrIntFloat_16=((Src.height)<<16)/Dst.height+1;
        
    const long csDErrorX=-(1<<15)+(xrIntFloat_16>>1);
        
    const long csDErrorY=-(1<<15)+(yrIntFloat_16>>1);

        unsigned 
    long dst_width=Dst.width;

        
    //计算出需要特殊处理的边界
        long border_y0=-csDErrorY/yrIntFloat_16+1;              //y0+y*yr>=0; y0=csDErrorY => y>=-csDErrorY/yr
        if (border_y0>=Dst.height) border_y0=Dst.height;
        
    long border_x0=-csDErrorX/xrIntFloat_16+1;     
        
    if (border_x0>=Dst.width ) border_x0=Dst.width; 
        
    long border_y1=(((Src.height-2)<<16)-csDErrorY)/yrIntFloat_16+1; //y0+y*yr<=(height-2) => y<=(height-2-csDErrorY)/yr
        if (border_y1<border_y0) border_y1=border_y0;
        
    long border_x1=(((Src.width-2)<<16)-csDErrorX)/xrIntFloat_16+1; 
        
    if (border_x1<border_x0) border_x1=border_x0;

        TARGB32* pDstLine=Dst.pdata;
        
    long Src_byte_width=Src.byte_width;
        
    long srcy_16=csDErrorY;
        
    long y;
        
    for (y=0;y<border_y0;++y)
        {
            
    long srcx_16=csDErrorX;
            
    for (unsigned long x=0;x<dst_width;++x)
            {
                Bilinear_Border_MMX(Src,srcx_16,srcy_16,&pDstLine[x]); 
    //border
                srcx_16+=xrIntFloat_16;
            }
            srcy_16+=yrIntFloat_16;
            ((TUInt8*&)pDstLine)+=Dst.byte_width;
        }

        
    for (y=border_y0;y<border_y1;++y)
        {
            
    long srcx_16=csDErrorX;
            
    long x;
            
    for (x=0;x<border_x0;++x)
            {
                Bilinear_Border_MMX(Src,srcx_16,srcy_16,&pDstLine[x]);
    //border
                srcx_16+=xrIntFloat_16;
            }

            {
                
    long dst_width_fast=border_x1-border_x0;
                
    if (dst_width_fast>0)
                {
                    unsigned 
    long v_8=(srcy_16 & 0xFFFF)>>8;
                    TARGB32* PSrcLineColor= (TARGB32*)((TUInt8*)(Src.pdata)+Src_byte_width*(srcy_16>>16)) ;
                    TARGB32* PSrcLineColorNext= (TARGB32*)((TUInt8*)(PSrcLineColor)+Src_byte_width) ;
                    TARGB32* pDstLine_Fast=&pDstLine[border_x0];
                    asm
                    {
                          movd         mm6,v_8
                          pxor         mm7,mm7 
    //mm7=0
                          PUNPCKLWD    MM6,MM6
                          PUNPCKLDQ    MM6,MM6
    //mm6=v_8
                        
                          mov       esi,PSrcLineColor
                          mov       ecx,PSrcLineColorNext
                          mov       edx,srcx_16
                          mov       ebx,dst_width_fast
                          mov       edi,pDstLine_Fast
                          lea       edi,[edi+ebx*4]
                          push      ebp
                          mov       ebp,xrIntFloat_16
                          neg       ebx

                    loop_start:

                              mov       eax,edx
                              shl       eax,16
                              shr       eax,24
                              
    //== movzx       eax,dh  //eax=u_8
                              MOVD      MM5,eax
                              mov       eax,edx
                              shr       eax,16     
    //srcx_16>>16

                              MOVD         MM2,dword ptr [ecx+eax*4]  
                              MOVD         MM0,dword ptr [ecx+eax*4+4]
                              PUNPCKLWD    MM5,MM5
                              MOVD         MM3,dword ptr [esi+eax*4]  
                              MOVD         MM1,dword ptr [esi+eax*4+4]
                              PUNPCKLDQ    MM5,MM5 
    //mm5=u_8
                              PUNPCKLBW    MM0,MM7
                              PUNPCKLBW    MM1,MM7
                              PUNPCKLBW    MM2,MM7
                              PUNPCKLBW    MM3,MM7
                              PSUBw        MM0,MM2
                              PSUBw        MM1,MM3
                              PSLLw        MM2,8
                              PSLLw        MM3,8
                              PMULlw       MM0,MM5
                              PMULlw       MM1,MM5
                              PADDw        MM0,MM2
                              PADDw        MM1,MM3

                              PSRLw        MM0,8
                              PSRLw        MM1,8
                              PSUBw        MM0,MM1
                              PSLLw        MM1,8
                              PMULlw       MM0,MM6
                              PADDw        MM0,MM1

                              PSRLw     MM0,8
                              PACKUSwb  MM0,MM7
                              MOVd   dword ptr    [edi+ebx*4],MM0 
    //write DstColor
                                          
                              add       edx,ebp 
    //srcx_16+=xrIntFloat_16
                              inc       ebx
                              jnz       loop_start

                          pop       ebp
                          mov       srcx_16,edx
                    }
                }
            }

            
    for (x=border_x1;x<dst_width;++x)
            {
                Bilinear_Border_MMX(Src,srcx_16,srcy_16,&pDstLine[x]);
    //border
                srcx_16+=xrIntFloat_16;
            }
            srcy_16+=yrIntFloat_16;
            ((TUInt8*&)pDstLine)+=Dst.byte_width;
        }
        
    for (y=border_y1;y<Dst.height;++y)
        {
            
    long srcx_16=csDErrorX;
            
    for (unsigned long x=0;x<dst_width;++x)
            {
                Bilinear_Border_MMX(Src,srcx_16,srcy_16,&pDstLine[x]); 
    //border
                srcx_16+=xrIntFloat_16;
            }
            srcy_16+=yrIntFloat_16;
            ((TUInt8*&)pDstLine)+=Dst.byte_width;
        }
        asm emms
    }


    //速度测试:
    //==============================================================================
    // PicZoom_Bilinear_MMX_Ex 157.0 fps

    I: 把测试成绩放在一起:


    //CPU: AMD64x2 4200+(2.37G)  zoom 800*600 to 1024*768
    //==============================================================================
    // StretchBlt                   232.7 fps   
    // PicZoom3_SSE                 711.7 fps 
    // 
    // PicZoom_BilInear0              8.3 fps
    // PicZoom_BilInear1             17.7 fps
    // PicZoom_BilInear2             43.4 fps
    // PicZoom_BilInear_Common       65.3 fps
    // PicZoom_BilInear_MMX         132.9 fps
    // PicZoom_BilInear_MMX_Ex      157.0 fps

    补充Intel Core2 4400上的测试成绩:


    //CPU: Intel Core2 4400(2.00G)  zoom 800*600 to 1024*768
    //==============================================================================
    // PicZoom3_SSE                1099.7 fps  
    // 
    // PicZoom_BilInear0             10.7 fps
    // PicZoom_BilInear1             24.2 fps
    // PicZoom_BilInear2             54.3 fps
    // PicZoom_BilInear_Common       59.8 fps
    // PicZoom_BilInear_MMX         118.4 fps
    // PicZoom_BilInear_MMX_Ex      142.9 fps

     

     

    三次卷积插值:


    J: 二次线性插值缩放出的图片很多时候让人感觉变得模糊(术语叫低通滤波),特别是在放大
    的时候;使用三次卷积插值来改善插值结果;三次卷积插值考虑映射点周围16个点(4x4)的颜色来
    计算最终的混合颜色,如图;

              
             P(0,0)所在像素为映射的点,加上它周围的15个点,按一定系数混合得到最终输出结果;

             混合公式参见PicZoom_ThreeOrder0的实现;

        插值曲线公式sin(x*PI)/(x*PI),如图:


                                 三次卷积插值曲线sin(x*PI)/(x*PI) (其中PI=3.1415926...)

     

    K:三次卷积插值缩放算法的一个参考实现:PicZoom_ThreeOrder0
      该函数并没有做过多的优化,只是一个简单的浮点实现版本;


     


            inline 
    double SinXDivX(double x) 
            {
                
    //该函数计算插值曲线sin(x*PI)/(x*PI)的值 //PI=3.1415926535897932385; 
                //下面是它的近似拟合表达式
                const float a = -1; //a还可以取 a=-2,-1,-0.75,-0.5等等,起到调节锐化或模糊程度的作用

                
    if (x<0) x=-x; //x=abs(x);
                double x2=x*x;
                
    double x3=x2*x;
                
    if (x<=1)
                  
    return (a+2)*x3 - (a+3)*x2 + 1;
                
    else if (x<=2) 
                  
    return a*x3 - (5*a)*x2 + (8*a)*x - (4*a);
                
    else
                  
    return 0;
            } 

            inline TUInt8 border_color(
    long Color)
            {
                
    if (Color<=0)
                    
    return 0;
                
    else if (Color>=255)
                    
    return 255;
                
    else
                    
    return Color;
            }
            
        
    void ThreeOrder0(const TPicRegion& pic,const float fx,const float fy,TARGB32* result)
        {
            
    long x0=(long)fx; if (x0>fx) --x0; //x0=floor(fx);    
            long y0=(long)fy; if (y0>fy) --y0; //y0=floor(fy);
            float fu=fx-x0;
            
    float fv=fy-y0;

            TARGB32 pixel[16];
            
    long i,j;

            
    for (i=0;i<4;++i)
            {
                
    for (j=0;j<4;++j)
                {
                    
    long x=x0-1+j;
                    
    long y=y0-1+i;
                    pixel[i*4+j]=Pixels_Bound(pic,x,y);
                }
            }

            
    float afu[4],afv[4];
            
    //
            afu[0]=SinXDivX(1+fu);
            afu[1]=SinXDivX(fu);
            afu[2]=SinXDivX(1-fu);
            afu[3]=SinXDivX(2-fu);
            afv[0]=SinXDivX(1+fv);
            afv[1]=SinXDivX(fv);
            afv[2]=SinXDivX(1-fv);
            afv[3]=SinXDivX(2-fv);

            
    float sR=0,sG=0,sB=0,sA=0;
            
    for (i=0;i<4;++i)
            {
                
    float aR=0,aG=0,aB=0,aA=0;
                
    for (long j=0;j<4;++j)
                {
                    aA+=afu[j]*pixel[i*4+j].a;
                    aR+=afu[j]*pixel[i*4+j].r;
                    aG+=afu[j]*pixel[i*4+j].g;
                    aB+=afu[j]*pixel[i*4+j].b;
                }
                sA+=aA*afv[i];
                sR+=aR*afv[i];
                sG+=aG*afv[i];
                sB+=aB*afv[i];
            }

            result->a=border_color((
    long)(sA+0.5));
            result->r=border_color((
    long)(sR+0.5));
            result->g=border_color((
    long)(sG+0.5));
            result->b=border_color((
    long)(sB+0.5));
        }

    void PicZoom_ThreeOrder0(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;
        TARGB32* pDstLine=Dst.pdata;
        
    for (unsigned long y=0;y<Dst.height;++y)
        {
            
    float srcy=(y+0.4999999)*Src.height/Dst.height-0.5;
            
    for (unsigned long x=0;x<dst_width;++x)
            {
                
    float srcx=(x+0.4999999)*Src.width/Dst.width-0.5;
                ThreeOrder0(Src,srcx,srcy,&pDstLine[x]);
            }
            ((TUInt8*&)pDstLine)+=Dst.byte_width;
        }
    }


    //速度测试:
    //==============================================================================
    // PicZoom_ThreeOrder0    3.6 fps

    L: 使用定点数来优化缩放函数;边界和内部分开处理;对SinXDivX做一个查找表;对border_color做一个查找表;


     

        static long SinXDivX_Table_8[(2<<8)+1];
        
    class _CAutoInti_SinXDivX_Table {
        
    private
            
    void _Inti_SinXDivX_Table()
            {
                
    for (long i=0;i<=(2<<8);++i)
                    SinXDivX_Table_8[i]=
    long(0.5+256*SinXDivX(i*(1.0/(256))))*1;
            };
        
    public:
            _CAutoInti_SinXDivX_Table() { _Inti_SinXDivX_Table(); }
        };
        
    static _CAutoInti_SinXDivX_Table __tmp_CAutoInti_SinXDivX_Table;


        
    //颜色查表
        static TUInt8 _color_table[256*3];
        
    static const TUInt8* color_table=&_color_table[256];
        
    class _CAuto_inti_color_table
        {
        
    public:
            _CAuto_inti_color_table() {
                
    for (int i=0;i<256*3;++i)
                    _color_table[i]=border_color(i-256);
            }
        };
        
    static _CAuto_inti_color_table _Auto_inti_color_table;

        
    void ThreeOrder_Fast_Common(const TPicRegion& pic,const long x_16,const long y_16,TARGB32* result)
        {
            unsigned 
    long u_8=(unsigned char)((x_16)>>8);
            unsigned 
    long v_8=(unsigned char)((y_16)>>8);
            
    const TARGB32* pixel=&Pixels(pic,(x_16>>16)-1,(y_16>>16)-1);
            
    long pic_byte_width=pic.byte_width;

            
    long au_8[4],av_8[4];
            
    //
            au_8[0]=SinXDivX_Table_8[(1<<8)+u_8];
            au_8[1]=SinXDivX_Table_8[u_8];
            au_8[2]=SinXDivX_Table_8[(1<<8)-u_8];
            au_8[3]=SinXDivX_Table_8[(2<<8)-u_8];
            av_8[0]=SinXDivX_Table_8[(1<<8)+v_8];
            av_8[1]=SinXDivX_Table_8[v_8];
            av_8[2]=SinXDivX_Table_8[(1<<8)-v_8];
            av_8[3]=SinXDivX_Table_8[(2<<8)-v_8];

            
    long sR=0,sG=0,sB=0,sA=0;
            
    for (long i=0;i<4;++i)
            {
                
    long aA=au_8[0]*pixel[0].a + au_8[1]*pixel[1].a + au_8[2]*pixel[2].a + au_8[3]*pixel[3].a;
                
    long aR=au_8[0]*pixel[0].r + au_8[1]*pixel[1].r + au_8[2]*pixel[2].r + au_8[3]*pixel[3].r;
                
    long aG=au_8[0]*pixel[0].g + au_8[1]*pixel[1].g + au_8[2]*pixel[2].g + au_8[3]*pixel[3].g;
                
    long aB=au_8[0]*pixel[0].b + au_8[1]*pixel[1].b + au_8[2]*pixel[2].b + au_8[3]*pixel[3].b;
                sA+=aA*av_8[i];
                sR+=aR*av_8[i];
                sG+=aG*av_8[i];
                sB+=aB*av_8[i];
                ((TUInt8*&)pixel)+=pic_byte_width;
            }

            result->a=color_table[sA>>16];
            result->r=color_table[sR>>16];
            result->g=color_table[sG>>16];
            result->b=color_table[sB>>16];
        }

        
    void ThreeOrder_Border_Common(const TPicRegion& pic,const long x_16,const long y_16,TARGB32* result)
        {
            
    long x0_sub1=(x_16>>16)-1;
            
    long y0_sub1=(y_16>>16)-1;
            unsigned 
    long u_16_add1=((unsigned short)(x_16))+(1<<16);
            unsigned 
    long v_16_add1=((unsigned short)(y_16))+(1<<16);

            TARGB32 pixel[16];
            
    long i;

            
    for (i=0;i<4;++i)
            {
                
    long y=y0_sub1+i;
                pixel[i*4+0]=Pixels_Bound(pic,x0_sub1+0,y);
                pixel[i*4+1]=Pixels_Bound(pic,x0_sub1+1,y);
                pixel[i*4+2]=Pixels_Bound(pic,x0_sub1+2,y);
                pixel[i*4+3]=Pixels_Bound(pic,x0_sub1+3,y);
            }
            
            TPicRegion npic;
            npic.pdata     =&pixel[0];
            npic.byte_width=4*
    sizeof(TARGB32);
            
    //npic.width     =4;
            //npic.height    =4;
            ThreeOrder_Fast_Common(npic,u_16_add1,v_16_add1,result);
        }

    void PicZoom_ThreeOrder_Common(const TPicRegion& Dst,const TPicRegion& Src)
    {
        
    if (  (0==Dst.width)||(0==Dst.height)
            ||(0==Src.width)||(0==Src.height)) 
    return;

        
    long xrIntFloat_16=((Src.width)<<16)/Dst.width+1; 
        
    long yrIntFloat_16=((Src.height)<<16)/Dst.height+1;
        
    const long csDErrorX=-(1<<15)+(xrIntFloat_16>>1);
        
    const long csDErrorY=-(1<<15)+(yrIntFloat_16>>1);

        unsigned 
    long dst_width=Dst.width;

        
    //计算出需要特殊处理的边界
        long border_y0=((1<<16)-csDErrorY)/yrIntFloat_16+1;//y0+y*yr>=1; y0=csDErrorY => y>=(1-csDErrorY)/yr
        if (border_y0>=Dst.height) border_y0=Dst.height;
        
    long border_x0=((1<<16)-csDErrorX)/xrIntFloat_16+1;
        
    if (border_x0>=Dst.width ) border_x0=Dst.width;
        
    long border_y1=(((Src.height-3)<<16)-csDErrorY)/yrIntFloat_16+1; //y0+y*yr<=(height-3) => y<=(height-3-csDErrorY)/yr
        if (border_y1<border_y0) border_y1=border_y0;
        
    long border_x1=(((Src.width-3)<<16)-csDErrorX)/xrIntFloat_16+1;; 
        
    if (border_x1<border_x0) border_x1=border_x0;

        TARGB32* pDstLine=Dst.pdata;
        
    long srcy_16=csDErrorY;
        
    long y;
        
    for (y=0;y<border_y0;++y)
        {
            
    long srcx_16=csDErrorX;
            
    for (unsigned long x=0;x<dst_width;++x)
            {
                ThreeOrder_Border_Common(Src,srcx_16,srcy_16,&pDstLine[x]); 
    //border
                srcx_16+=xrIntFloat_16;
            }
            srcy_16+=yrIntFloat_16;
            ((TUInt8*&)pDstLine)+=Dst.byte_width;
        }
        
    for (y=border_y0;y<border_y1;++y)
        {
            
    long srcx_16=csDErrorX;
            
    long x;
            
    for (x=0;x<border_x0;++x)
            {
                ThreeOrder_Border_Common(Src,srcx_16,srcy_16,&pDstLine[x]);
    //border
                srcx_16+=xrIntFloat_16;
            }
            
    for (x=border_x0;x<border_x1;++x)
            {
                ThreeOrder_Fast_Common(Src,srcx_16,srcy_16,&pDstLine[x]);
    //fast  !
                srcx_16+=xrIntFloat_16;
            }
            
    for (x=border_x1;x<dst_width;++x)
            {
                ThreeOrder_Border_Common(Src,srcx_16,srcy_16,&pDstLine[x]);
    //border
                srcx_16+=xrIntFloat_16;
            }
            srcy_16+=yrIntFloat_16;
            ((TUInt8*&)pDstLine)+=Dst.byte_width;
        }
        
    for (y=border_y1;y<Dst.height;++y)
        {
            
    long srcx_16=csDErrorX;
            
    for (unsigned long x=0;x<dst_width;++x)
            {
                ThreeOrder_Border_Common(Src,srcx_16,srcy_16,&pDstLine[x]); 
    //border
                srcx_16+=xrIntFloat_16;
            }
            srcy_16+=yrIntFloat_16;
            ((TUInt8*&)pDstLine)+=Dst.byte_width;
        }
    }


    //速度测试:
    //==============================================================================
    // PicZoom_ThreeOrder_Common    16.9 fps

     M: 用MMX来优化ThreeOrder_Common函数:ThreeOrder_MMX


     

        typedef   unsigned long TMMXData32;
        
    static TMMXData32 SinXDivX_Table_MMX[(2<<8)+1];
        
    class _CAutoInti_SinXDivX_Table_MMX {
        
    private
            
    void _Inti_SinXDivX_Table_MMX()
            {
                
    for (long i=0;i<=(2<<8);++i)
                {
                    unsigned 
    short t=long(0.5+(1<<14)*SinXDivX(i*(1.0/(256))));
                    unsigned 
    long tl=t | (((unsigned long)t)<<16);
                    SinXDivX_Table_MMX[i]=tl;
                }
            };
        
    public:
            _CAutoInti_SinXDivX_Table_MMX() { _Inti_SinXDivX_Table_MMX(); }
        };
        
    static _CAutoInti_SinXDivX_Table_MMX __tmp_CAutoInti_SinXDivX_Table_MMX;


        
    void __declspec(naked) _private_ThreeOrder_Fast_MMX()
        {
            asm
            {
                movd        mm1,dword ptr [edx]
                movd        mm2,dword ptr [edx+4]
                movd        mm3,dword ptr [edx+8]
                movd        mm4,dword ptr [edx+12]
                movd        mm5,dword ptr [(offset SinXDivX_Table_MMX)+256*4+eax*4]
                movd        mm6,dword ptr [(offset SinXDivX_Table_MMX)+eax*4]
                punpcklbw   mm1,mm7
                punpcklbw   mm2,mm7
                punpcklwd   mm5,mm5
                punpcklwd   mm6,mm6
                psllw       mm1,7
                psllw       mm2,7
                pmulhw      mm1,mm5
                pmulhw      mm2,mm6
                punpcklbw   mm3,mm7
                punpcklbw   mm4,mm7
                movd        mm5,dword ptr [(offset SinXDivX_Table_MMX)+256*4+ecx*4]
                movd        mm6,dword ptr [(offset SinXDivX_Table_MMX)+512*4+ecx*4]
                punpcklwd   mm5,mm5
                punpcklwd   mm6,mm6
                psllw       mm3,7
                psllw       mm4,7
                pmulhw      mm3,mm5
                pmulhw      mm4,mm6
                paddsw      mm1,mm2
                paddsw      mm3,mm4
                movd        mm6,dword ptr [ebx] 
    //v
                paddsw      mm1,mm3
                punpcklwd   mm6,mm6

                pmulhw      mm1,mm6
                add     edx,esi  
    //+pic.byte_width
                paddsw      mm0,mm1

                ret
            }
        }

        inline 
    void ThreeOrder_Fast_MMX(const TPicRegion& pic,const long x_16,const long y_16,TARGB32* result)
        {
            asm
            {
                mov     ecx,pic
                mov     eax,y_16
                mov     ebx,x_16
                movzx   edi,ah 
    //v_8
                mov     edx,[ecx+TPicRegion::pdata]
                shr     eax,16
                mov     esi,[ecx+TPicRegion::byte_width]
                dec     eax
                movzx   ecx,bh 
    //u_8
                shr     ebx,16
                imul    eax,esi
                lea     edx,[edx+ebx*4-4]
                add     edx,eax 
    //pixel

                mov     eax,ecx
                neg     ecx

                pxor    mm7,mm7  
    //0
                //mov     edx,pixel
                pxor    mm0,mm0  //result=0
                //lea     eax,auv_7

                lea    ebx,[(offset SinXDivX_Table_MMX)+256*4+edi*4]
                call  _private_ThreeOrder_Fast_MMX
                lea    ebx,[(offset SinXDivX_Table_MMX)+edi*4]
                call  _private_ThreeOrder_Fast_MMX
                neg    edi
                lea    ebx,[(offset SinXDivX_Table_MMX)+256*4+edi*4]
                call  _private_ThreeOrder_Fast_MMX
                lea    ebx,[(offset SinXDivX_Table_MMX)+512*4+edi*4]
                call  _private_ThreeOrder_Fast_MMX

                psraw     mm0,3
                mov       eax,result
                packuswb  mm0,mm7
                movd      [eax],mm0
                
    //emms
            }
        }

        
    void ThreeOrder_Border_MMX(const TPicRegion& pic,const long x_16,const long y_16,TARGB32* result)
        {
            unsigned 
    long x0_sub1=(x_16>>16)-1;
            unsigned 
    long y0_sub1=(y_16>>16)-1;
            
    long u_16_add1=((unsigned short)(x_16))+(1<<16);
            
    long v_16_add1=((unsigned short)(y_16))+(1<<16);

            TARGB32 pixel[16];

            
    for (long i=0;i<4;++i)
            {
                
    long y=y0_sub1+i;
                pixel[i*4+0]=Pixels_Bound(pic,x0_sub1  ,y);
                pixel[i*4+1]=Pixels_Bound(pic,x0_sub1+1,y);
                pixel[i*4+2]=Pixels_Bound(pic,x0_sub1+2,y);
                pixel[i*4+3]=Pixels_Bound(pic,x0_sub1+3,y);
            }
            
            TPicRegion npic;
            npic.pdata     =&pixel[0];
            npic.byte_width=4*
    sizeof(TARGB32);
            
    //npic.width     =4;
            //npic.height    =4;
            ThreeOrder_Fast_MMX(npic,u_16_add1,v_16_add1,result);
        }

    void PicZoom_ThreeOrder_MMX(const TPicRegion& Dst,const TPicRegion& Src)
    {
        
    if (  (0==Dst.width)||(0==Dst.height)
            ||(0==Src.width)||(0==Src.height)) 
    return;

        
    long xrIntFloat_16=((Src.width)<<16)/Dst.width+1; 
        
    long yrIntFloat_16=((Src.height)<<16)/Dst.height+1;
        
    const long csDErrorX=-(1<<15)+(xrIntFloat_16>>1);
        
    const long csDErrorY=-(1<<15)+(yrIntFloat_16>>1);

        unsigned 
    long dst_width=Dst.width;

        
    //计算出需要特殊处理的边界
        long border_y0=((1<<16)-csDErrorY)/yrIntFloat_16+1;//y0+y*yr>=1; y0=csDErrorY => y>=(1-csDErrorY)/yr
        if (border_y0>=Dst.height) border_y0=Dst.height;
        
    long border_x0=((1<<16)-csDErrorX)/xrIntFloat_16+1;
        
    if (border_x0>=Dst.width ) border_x0=Dst.width;
        
    long border_y1=(((Src.height-3)<<16)-csDErrorY)/yrIntFloat_16+1; //y0+y*yr<=(height-3) => y<=(height-3-csDErrorY)/yr
        if (border_y1<border_y0) border_y1=border_y0;
        
    long border_x1=(((Src.width-3)<<16)-csDErrorX)/xrIntFloat_16+1;; 
        
    if (border_x1<border_x0) border_x1=border_x0;

        TARGB32* pDstLine=Dst.pdata;
        
    long srcy_16=csDErrorY;
        
    long y;
        
    for (y=0;y<border_y0;++y)
        {
            
    long srcx_16=csDErrorX;
            
    for (unsigned long x=0;x<dst_width;++x)
            {
                ThreeOrder_Border_MMX(Src,srcx_16,srcy_16,&pDstLine[x]); 
    //border
                srcx_16+=xrIntFloat_16;
            }
            srcy_16+=yrIntFloat_16;
            ((TUInt8*&)pDstLine)+=Dst.byte_width;
        }
        
    for (y=border_y0;y<border_y1;++y)
        {
            
    long srcx_16=csDErrorX;
            
    long x;
            
    for (x=0;x<border_x0;++x)
            {
                ThreeOrder_Border_MMX(Src,srcx_16,srcy_16,&pDstLine[x]);
    //border
                srcx_16+=xrIntFloat_16;
            }
            
    for (x=border_x0;x<border_x1;++x)
            {
                ThreeOrder_Fast_MMX(Src,srcx_16,srcy_16,&pDstLine[x]);
    //fast MMX !
                srcx_16+=xrIntFloat_16;
            }
            
    for (x=border_x1;x<dst_width;++x)
            {
                ThreeOrder_Border_MMX(Src,srcx_16,srcy_16,&pDstLine[x]);
    //border
                srcx_16+=xrIntFloat_16;
            }
            srcy_16+=yrIntFloat_16;
            ((TUInt8*&)pDstLine)+=Dst.byte_width;
        }
        
    for (y=border_y1;y<Dst.height;++y)
        {
            
    long srcx_16=csDErrorX;
            
    for (unsigned long x=0;x<dst_width;++x)
            {
                ThreeOrder_Border_MMX(Src,srcx_16,srcy_16,&pDstLine[x]); 
    //border
                srcx_16+=xrIntFloat_16;
            }
            srcy_16+=yrIntFloat_16;
            ((TUInt8*&)pDstLine)+=Dst.byte_width;
        }
        asm emms
    }


    //速度测试:
    //==============================================================================
    // PicZoom_ThreeOrder_MMX   34.3 fps

     

     N:将测试结果放到一起:


    //CPU: AMD64x2 4200+(2.37G)  zoom 800*600 to 1024*768
    //==============================================================================
    // StretchBlt                   232.7 fps   
    // PicZoom3_SSE                 711.7 fps  
    // PicZoom_BilInear_MMX_Ex      157.0 fps
    // 
    // PicZoom_ThreeOrder0            3.6 fps
    // PicZoom_ThreeOrder_Common     16.9 fps
    // PicZoom_ThreeOrder_MMX        34.3 fps

    补充Intel Core2 4400上的测试成绩:


    //CPU: Intel Core2 4400(2.00G)  zoom 800*600 to 1024*768
    //==============================================================================
    // PicZoom3_SSE                1099.7 fps  
    // PicZoom_BilInear_MMX_Ex      142.9 fps
    // 
    // PicZoom_ThreeOrder0            4.2 fps
    // PicZoom_ThreeOrder_Common     17.6 fps
    // PicZoom_ThreeOrder_MMX        34.4 fps

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

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

    正文:

    A:对于前一篇文章中的二次线性插值、三次卷积插值算法,但它们处理缩小到0.5倍以下的
    时候效果就会越来越差;这是因为插值的时候自考虑了附近点的原因;如下图:

                          
               原图          近邻取样 缩放到0.4倍     缩放到0.2倍     缩放到0.1倍

                                                     
                         二次线性插值 缩放到0.4倍     缩放到0.2倍     缩放到0.1倍

                                                     
                         三次卷积插值 缩放到0.4倍     缩放到0.2倍     缩放到0.1倍


        可以看出:当缩小的比例很大的时候,插值算法的效果和近邻取样的效果差不多了:( ;
    一种可行的解决方案就是:缩小时考虑更多的点; 但这种解决方案有很多缺点:函数编写麻烦,
    速度也许会很慢,优化也不容易做;  还有一个方案就是预先建立一个缩放好的大小不同的图片
    列表,每一张图片都是前一张的0.5倍(这种图片列表就是MipMap链),缩放的时候根据需要缩放
    的比例从表中选择一张大小接近的图片来作为缩放的源图片; 该方案的优点:不需要编写新的
    底层缩放算法,直接使用前面优化好的插值算法; 缺点:需要预先建立MipMap链,它需要时间,
    并且它的储存需要多占用原图片的1/3空间(0.5^2+0.5^4+0.5^6+...=1/3);还有一个不太明显
    的小问题,就是在一张图片的连续的比例不同的缩放中,选择会从MipMap的一张源图片跳到另
    一张图片,视觉效果上可能会有一个小的跳跃(我在《魔兽世界》里经常看到这种效应:);一种
    改进方案就是选择MipMap图片的时候,选择出附近的两张图片作为缩放的源图片;对两张图片
    单独进行插值(和原来一致)输出两个值,然后把这两个值线性插值为最终结果;还有一个比较
    大的缺点就是当缩放比例不均匀时(比如x轴放大y轴缩小),缩放效果也不好;
    (当前很多显卡都提供了MipMap纹理和对应的插值方案,OpenGL和DirectX都提供了操作接口)

         
    ("三次线性插值和MipMap链"其实比较简单,这里只给出关键代码或算法)
     

    B: MipMap图片的生成:
         原图片缩放到0.5倍(宽和高都为原图片的1/2),在把0.5倍的图片缩放到0.25倍,....
       直到宽和高都为1个像素,如果有一个长度先到1就保持1; 缩放过程中,可以可采用前面的缩放插值算法;
       如果为了速度可以考虑这样的方案,要求原图片的宽和高必须是2的整数次方的数值,缩放时就可以直接将
       2x2的像素快速合并为一个像素(如果允许原图片宽和高为任何值,可以考虑在合并时引入Alpha通道);

    C: MipMap链图片的储存方案:

                   
                        MipMap链图片示意图         
                       
                               
             可能的一种物理储存方案(我对每张图片加了一个边框)            

    D: 定义MipMap数据结构:
       MipMap数据结构可以定义为一个TPicRegion数组和该数组的大小; 
       (MipMap图片的储存参见上面的图示)
      比如:
     

         #include <vector>
         typedef std::vector
    <TPicRegion> TMipMap;
         
    //其中,第一个元素TMipMap[0]指向原始图片,后面的依次为缩小图片;       

     
    E: MipMap的选择函数和偏好:
        在进行缩放时,根据目标图片缓冲区的大小来动态的选者MipMap中的一幅图片来作为源图片;这就需要一个
    选择函数;比如:

    long SelectBestPicIndex(const TMipMap& mip,const long dstWidth,const long dstHeight)
    {
        
    long oldS=mip[0].width*mip[0].height;
        
    long dstS=dstWidth*dstHeight;
        
    if ( (dstS>=oldS) || (mip.size()==1) )
            
    return 0;
        
    else if (dstS<=1)
            
    return mip.size()-1;
        
    else

            
    return (long)(log(oldS/dstS)*0.5+0.5);
    }

    选择函数可以增加一个偏好参数:
    mip选择偏好:0.5没有偏好,靠近0偏向选择小图片,靠近1偏向选择大图片(质量好一些)

    float public_mip_bias=0.5; //[0..1] 

    long SelectBestPicIndex(const TMipMap& mip,const long dstWidth,const long dstHeight)
    {
        
    long oldS=mip[0].width*mip[0].height;
        
    long dstS=dstWidth*dstHeight;
        
    if ( (dstS>=oldS) || (mip.size()==1) )
            
    return 0;
        
    else if (dstS<=1)
            
    return mip.size()-1;
        
    else

            
    return (long)(log(oldS/dstS)*0.5+public_mip_bias);
    }

     F:利用MipMap后的缩放效果:

                                                     
                     MipMap+近邻取样 缩放到0.4倍     缩放到0.2倍     缩放到0.1倍
                   (利用MipMap做一次近邻取样)

                                                     
                  MipMap+二次线性插值 缩放到0.4倍     缩放到0.2倍     缩放到0.1倍
                  (利用MipMap做一次二次线性插值)

                                                     
                  MipMap+三次卷积插值 缩放到0.4倍     缩放到0.2倍     缩放到0.1倍
                  (利用MipMap做一次三次卷积插值)

       
                       

    G: 在MipMap的两张图片之间插值:
      选择MipMap的时候,同时可以选择相邻的两张MipMap图片;分别进行插值算法后得到两个颜色结果;
    对两个MipMap图片产生的评价值可以作为这两个颜色的插值权重,得到最终的颜色插值结果;优点是
    缩放效果好,避免跳跃;缺点是速度慢:)  

    选择和权重函数的一个可能实现:

    struct TMipWeight {
      
    long  BigMip;
      
    long  SmallMip;
      
    float BigMipWeight;//[0..1]

    };

    TMipWeight SelectBestPicIndexEx(
    const TMipMap& mip,const long dstWidth,const long dstHeight)
    {
        
    long oldS=mip[0].width*mip[0].height;
        
    long dstS=dstWidth*dstHeight;
        TMipWeight result;
        
    if ( (dstS>=oldS) || (mip.size()==1) )
        {
            result.BigMip=0;
            result.SmallMip=0;
            result.BigMipWeight=1.0;
        }
        
    else if (dstS<=1)
        {
            result.BigMip=mip.size()-1;
            result.SmallMip=mip.size()-1;
            result.BigMipWeight=1.0;
        }
        
    else

        {
            
    float bestIndex=log(oldS/dstS)*0.5+0.5; //or + public_mip_bias
            result.BigMip=(long)bestIndex;
            
    if (bestIndex==mip.size()-1)
            {
                result.SmallMip=mip.size()-1;
                result.BigMipWeight=1.0;
            }
            
    else

            {
                result.SmallMip
    =result.BigMip+1;
                result.BigMipWeight=1.0-(bestIndex-result.BigMip);
            }
        }
        
    return
     result;
    }

     

    H:MipMap间插值效果:

                                                    
                  MipMap+两次近邻取样 缩放到0.4倍     缩放到0.2倍     缩放到0.1倍
                  (利用MipMap做两次近邻取样输出两个值,然后线性插值为最终结果)

                                                    
                         三次线性插值 缩放到0.4倍     缩放到0.2倍     缩放到0.1倍
                 (三次线性插值:利用MipMap做两次二次线性插值输出两个值,然后线性插值为最终结果)  
                       
                                                    
               MipMap+两次三次卷积插值 缩放到0.4倍    缩放到0.2倍     缩放到0.1倍
             (利用MipMap做两次三次卷积插值输出两个值,然后线性插值为最终结果)
     

     

    (图像缩放系列终于写完了,计划中写图像任意角度的高质量的快速旋转、Alpha图片混合等,尽请期待:)

     (ps: 思考中的一个图片压缩方法:利用MipMap来压缩图像数据;输入一张图片,然后生成MipMap链,保存相邻之间图片的差(数值差可能很小,很容易找好的算法压缩得很小)和最顶的一张图片(一个点);  解压的时候依次求和就得到原图片了;  该算法为无损压缩,适合于人物风景等过渡比较多的图片的压缩,不太适合线条类等相邻间颜色变化剧烈的图片;)

    展开全文
  • 其中有两种简单又常用的插值算法用来实现图像缩放,分别是最近邻插值算法和双线性插值算法。 最近邻插值算法: 最近邻插值算法的思想十分简单 设原始图像src的高为h,宽为w,图像上像素值为(x,y)。 设目标图像dst...
  • yuv图像缩放算法

    2018-07-25 11:46:39
    完整的例子,对YUV图片进行缩放,代码里有几个缩放方法
  • 阿里云OSS不支持20M以上图片缩放。 用opencv在内存里缩放,30000x30000的png图片,占用3g+内存,有点无法接受。 于是自己调研相关算法。 以300x300的png为例,正常缩放时,是将其整个读进内存,生成300x300x4的...
  • 在qt中运行,改善图片缩放效果。。在线等大神啊 毕设急求啊![图片说明](http://forum.csdn.net/PointForum/ui/scripts/csdn/Plugin/001/face/5.gif)![图片说明]...
  • 图像缩放算法小结

    万次阅读 多人点赞 2019-04-10 09:33:18
    双线性插值法效果要好于最近邻插值,只是计算量稍大一些,算法复杂些,程序运行时间也稍长些,但缩放后图像质量高,基本克服了最近邻插值灰度值不连续的特点,因为它考虑了待测采样点周围四个直接邻点对该采样点的...
  • 图像缩放算法及速度优化

    千次阅读 2015-10-14 19:01:27
    图像缩放算法及速度优化——(一)最近邻插值 图像缩放算法及速度优化——(二)双线性插值 ————————————————————以下为原文—————————————————— 第0节 简介  图像缩放...
  • 采用双线性算法实现图片缩放,今天,经过多次实现终于完成了图片缩放,来给大家分享一下
  • 图像 矩形 定点/鼠标为中心缩放算法图像 矩形 定点/鼠标为中心缩放算法条件,代码 图像 矩形 定点/鼠标为中心缩放算法 程序中 图像 或 矩形, 按照指定坐标不变缩放算法 条件,代码 图像矩形区域: Rect 缩放系数:...
  • 图像缩放算法-最近邻插值

    千次阅读 2018-09-25 19:19:27
    简单来说,用最近邻插值算法实现图像的扩大和缩小任意的尺寸,目标图的... 这种放大图像的方法叫做最临近插值算法,这是一种最基本、最简单的图像缩放算法,效果也是最不好的,放大后的图像有很严重的马赛克,缩小...
  • 在Android平台上研究了现有的图像缩放算法,分析了最近邻域算法和线性插值算法(双线性插值算法),分析其基本原理,算法流程及优劣,并决定采用线性插值算法实现基于Android平台的图像处理的图像缩放功能,通过数据...
  • 使用hammer.js实现图片手势缩放算法解释

    千次阅读 热门讨论 2019-03-15 18:00:25
    关于图片缩放,这里先讲一下算法。 transform的缩放是指长宽的缩放,不是面积的缩放。 //transform-origin的初始位置是在元素的正中间,而缩放,平移,这个位置会跟着改变。 transform的matrix值:matrix(a,b,c,d...
  • 图像的缩放算法以及python实现

    千次阅读 2020-02-29 22:09:24
    一、图像缩放算法 最邻近插值法 这种方法是图像缩放算法中最简单的一种了。优缺点也十分明显,算法速度快,但缩放效果很差,放的太大就会照成马赛克现象,放的太小会造成严重的失真现像。下面先来进行介绍 算法步骤 ...
  • FPGA视频、图像缩放算法介绍

    千次阅读 2018-01-31 11:44:50
    视频缩放算法介绍:视频缩放算法包括最邻近插值,双线性插值,双线性三次插值,Lanczos插值算法等。算法的效率是最邻近算法>双线性插值算法>双线性三次插值>Lanczos算法,而算法的效果则恰恰相反 ######双线性...
  • C# 关于图片缩放、合并的算法

    千次阅读 2011-10-09 16:23:17
    最近想做一个图片缩放的软件,能批量处理一批图片,并且能够在图片添加水印文字和水印logo。从网上参考的核心代码如下:  1.合并两个图片  ///  /// 合并两张图片,支持不透明度和透明色  //
  • 图像缩放算法 Mat-->TPicRegion

    千次阅读 2015-01-20 20:04:31
    在网上看到一篇很好的关于图像缩放算法的博客,系统的讲解了如何去对图像缩放程序进行优化,从一个基本的图像缩放算法出发,然后一步一步的优化其速度和缩放质量。但是博客中只是针对图像的数据区TPicRegion进行操作...
  • 经典图像/视频缩放算法原理及MATLAB\Simulink实现1 图像/视频缩放介绍2 图像缩放的经典算法2.1 最近邻插值法2.2 双线性插值法3 图像缩放算法的MATLAB实现3.1 最近邻算法实现3.2 双线性插值算法实现4 视频无极缩放的...

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 43,416
精华内容 17,366
关键字:

图片缩放算法