三次卷积插值法

ArcGIS中栅格数据重采样方法作者:赵立超邮箱：zhaolc@ 发布时间：2014-01-02 21:04:28 阅读(800)ArcGIS中栅格数据重采样方法处理栅格数据时，由于数据像元大小不符合要求，或者在进行栅格数据配准后，像元发生倾斜，或者对多个栅格数据进行分析时，需要使用相同的栅格分辨率，重采样到同一分辨率下，因此对栅格数据操作时经常要进行重采样操作。

如下图所示，通过重采样将栅格调整到新的分辨率下。

1.1.栅格重采样方法栅格重采样主要包括三种方法：最邻近法、双线性内插法和三次卷积插值法。

最邻近法是把原始图像中距离最近的像元值填充到新图像中；双线性内插法和三次卷积插值法都是把原始图像附近的像元值通过距离加权平均填充到新图像中。

默认情况下，采用最近邻分配重采样技术，这种方法同时适用于离散和连续值类型，而其他重采样方法只适用于连续数据。

1.1 最近邻法最邻近分配法是用于离散（分类）数据的重采样技术，因为它不会更改输入单元的值。

将输出栅格数据集中单元中心的位置定位到输入栅格后，最邻近分配法将确定输入栅格上最近的单元中心位置并将该单元的值分配给输出栅格上的单元。

最邻近分配法不会更改输入栅格数据集中单元的任何值。

输入栅格中的值2 在输出栅格中仍将为 2，决不会为 2.2 或 2.3。

由于输出单元值保持不变，因此最邻近分配法应该用于名目数据或顺序数据，其中每个值都表示一个类、一个成员或一个分类（分类数据，如土地利用、土壤或森林类型）。

考虑到根据输入栅格创建的输出栅格会在操作中旋转 45°，因此将进行重采样。

对于每个输出单元，都要从输入栅格中获取值。

在下图中，输入栅格的单元中心为灰色点。

输出单元为绿色阴影。

要处理的单元为黄色阴影。

在最邻近分配法中，将确定与要处理的单元中心（红色点）最邻近的输入栅格单元中心（橙色点)，并将其指定为要处理的单元（黄色阴影）的输出值。

对输出栅格中的每个单元都重复此过程。

数字图像处理复习整理第一章1，什么是图像，模拟图像处理和数字图像处理主要区别。

图像是对客观对象的一种相似性的，生动性的描述或写真。

模拟图像：空间坐标和亮度（或色彩）都是连续性变化的图像；数字图像：空间坐标和灰度均不连续的，用离散数字（一般用整数）表示的图像。

利用光学，照相机方法对模拟图像的处理称为模拟图像处理，精度不高，稳定性差，设备笨重，操作不方便和工艺水平不高；利用计算机对数字图像进行系列操作称为数字图像处理，或计算机图像处理。

2，数字图像处理由哪些模块组成。

狭义图像处理图像分析图像理解3，数字图像处理的应用生物医学航空遥感工业应用军事公安其他第二章1，什么事图像对比度图像中最大亮度与最小亮度之比2，数字图像处理包括哪两个过程？对质量有何影响？数字图像的数据量和哪些因素有关？采样量化采样间隔越大，所得图像像素数越少，空间分辨率低，质量差，严重时出现棋盘格效应。

采样间隔越小，所得图像像素数越多，空间分辨率高，图像质量越好，但数据量大；量化等级越多，所得图像层次越丰富，灰度分辨率越高，质量越好，但数据量大。

量化等级越少，图像层次欠丰富，灰度分辨率低，质量变差。

会出现伪轮廓现象。

采样间隔量化等级3，连续图像f(x,y)与数字图像I(r,c)中各量的含义是什么？它们有何联系和区别？答：f(x,y)表示二维图像在空间（x,y）上的幅值，数字图像I(r,c)表示位于图像矩阵上第r行，第c列的元素幅值。

I(r,c)是通过对f(x,y)抽样和量化得来的。

f(x,y)各量是连续的，I(r,c)各量是离散的。

4，什么事灰度直方图？它有哪些应用？能从中获得图像的哪些信息？灰度直方图反应的是一幅图像中各灰度级像素出现的概率之间的关系。

应用：判断图像量化知否恰当；确定图像二值化阈值；计算图像中物体的面积；计算图像信息量H（熵）灰度范围，灰度级的分布，整幅图像的平均亮度。

5，熵的计算公式，灰度范围[0,L-1]6，什么是点处理？举例说明。

高质量的快速的图像缩放——二次线性插值和三次卷积插值限制条件:为了便于讨论，这里只处理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 __asmtypedef 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_BilInear1inline 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)(255-v_8));unsigned long pm1_16=(v_8*(unsigned long)(255-u_8));unsigned long pm0_16=((255-u_8)*(255-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]); //bordersrcx_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/yrif (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)/yrif (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]); //bordersrcx_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]);//bordersrcx_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]);//bordersrcx_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]); //bordersrcx_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/yrif (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)/yrif (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]); //bordersrcx_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]);//bordersrcx_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]);//bordersrcx_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]); //bordersrcx_16+=xrIntFloat_16;}srcy_16+=yrIntFloat_16;((TUInt8*&)pDstLine)+=Dst.byte_width;}}//////////////////////////////////////////////////////////////////////////////////速度测试://============================================================================== // PicZoom_BilInear_Common 65.3 fps//////////////////////////////////////////////////////////////////////////////// H: 使用MMX指令改写:PicZoom_Bilinear_MMXinline void Bilinear_Fast_MMX(TARGB32* PColor0,TARGB32* PColor1, unsigned long u_8,unsigned long v_8,TARGB32* result){asm{MOVD MM6,v_8MOVD MM5,u_8mov edx,PColor0mov eax,PColor1PXOR mm7,mm7MOVD MM2,dword ptr [eax]MOVD MM0,dword ptr [eax+4]PUNPCKLWD MM5,MM5PUNPCKLWD MM6,MM6MOVD MM3,dword ptr [edx]MOVD MM1,dword ptr [edx+4]PUNPCKLDQ MM5,MM5PUNPCKLBW MM0,MM7PUNPCKLBW MM1,MM7PUNPCKLBW MM2,MM7PUNPCKLBW MM3,MM7PSUBw MM0,MM2PSUBw MM1,MM3PSLLw MM2,8PSLLw MM3,8PMULlw MM0,MM5PMULlw MM1,MM5PUNPCKLDQ MM6,MM6PADDw MM0,MM2PADDw MM1,MM3PSRLw MM0,8PSRLw MM1,8PSUBw MM0,MM1PSLLw MM1,8PMULlw MM0,MM6mov eax,resultPADDw MM0,MM1PSRLw MM0,8PACKUSwb MM0,MM7movd [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/yrif (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)/yrif (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]);//bordersrcx_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]);//bordersrcx_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]); //bordersrcx_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_Exvoid 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/yrif (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)/yrif (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]); //bordersrcx_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]);//bordersrcx_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_8pxor mm7,mm7 //mm7=0PUNPCKLWD MM6,MM6PUNPCKLDQ MM6,MM6//mm6=v_8mov esi,PSrcLineColormov ecx,PSrcLineColorNextmov edx,srcx_16mov ebx,dst_width_fastmov edi,pDstLine_Fastlea edi,[edi+ebx*4]push ebpmov ebp,xrIntFloat_16neg ebxloop_start:mov eax,edxshl eax,16shr eax,24//== movzx eax,dh //eax=u_8MOVD MM5,eaxmov eax,edxshr eax,16 //srcx_16>>16MOVD MM2,dword ptr [ecx+eax*4]MOVD MM0,dword ptr [ecx+eax*4+4]PUNPCKLWD MM5,MM5MOVD MM3,dword ptr [esi+eax*4]MOVD MM1,dword ptr [esi+eax*4+4]PUNPCKLDQ MM5,MM5 //mm5=u_8PUNPCKLBW MM0,MM7PUNPCKLBW MM1,MM7PUNPCKLBW MM2,MM7PUNPCKLBW MM3,MM7PSUBw MM0,MM2PSUBw MM1,MM3PSLLw MM2,8PSLLw MM3,8PMULlw MM0,MM5PMULlw MM1,MM5PADDw MM0,MM2PADDw MM1,MM3PSRLw MM0,8PSRLw MM1,8PSUBw MM0,MM1PSLLw MM1,8PMULlw MM0,MM6PADDw MM0,MM1PSRLw MM0,8PACKUSwb MM0,MM7MOVd dword ptr [edi+ebx*4],MM0 //write DstColor add edx,ebp //srcx_16+=xrIntFloat_16inc ebxjnz loop_startpop ebpmov srcx_16,edx}}}for (x=border_x1;x<dst_width;++x){Bilinear_Border_MMX(Src,srcx_16,srcy_16,&pDstLine[x]);//bordersrcx_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]); //bordersrcx_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_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);elsereturn 0;}inline TUInt8 border_color(long Color){if (Color<=0)return 0;else if (Color>=255)return 255;elsereturn 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);}}。

插值方法在图像处理中的应用作者：专业姓名学号控制工程陈龙斌控制工程陈少峰控制工程殷文龙摘要本文介绍了插值方法在图像处理中的应用。

介绍了典型的最近邻插值、双线性插值、双三次插值、双信道插值、分形插值的原理。

以分形插值为重点，在图像放大领域用MATLAB进行仿真，并与其它方法的结果做了比对。

指出了各种方法的利弊，期待更进一步的研究拓展新的算法以及改进现有算法。

一、引言人类通过感觉器官从客观世界获取信息，而其中一半以上的信息都是通过视觉获得的。

图像作为人类视觉信息传递的主要媒介，具有声音、语言、文字等形式无法比拟的优势，给人以具体、直观的物体形象。

在数字化信息时代，图像处理已经成为重要的数据处理类型。

数字图像比之传统的模拟图像处理有着不可比拟的优势。

一般采用计算机处理或者硬件处理，处理的内容丰富，精度高，变通能力强，可进行非线性处理。

但是处理速度就会有所不足。

图像处理的主要内容有：几何处理、算术处理、图像增强、图像复原、图像重建、图像编码、图像识别、图像理解等。

以上这些图像处理大体上可分为图像的像质改善、图像分析和图像重建三大部分。

日常生活中，越来越多的领域需要高分辨率图像，采用图像插值技术来提高数字图像的分辨率和清晰度，从软件方面进行改进就具有十分重要的实用价值。

多媒体通信在现代网络传输中扮演重要角色，因此插值放大提高图像分辨率是一个非常重要的问题。

此外，图像变换被广泛用于遥感图像的几何校正、医学成像以及电影、电视和媒体广告等影像特技处理中。

在进行图像的一些几何变换时，通常都会出现输出像素坐标和输入栅格不重合的现象，也必须要用到图像插值。

图像插值是图像处理中图像重采样过程中的重要组成部分，而重采样过程广泛应用于改善图像质量、进行有损压缩等，因而研究图像插值具有十分重要的理论意义和实用价值。

图像插值是一个数据再生过程。

由原始图像数据再生出具有更高分辨率的图像数据。

分为图像内插值和图像间插值。

前者指将一幅较低分辨率的图像再生出一幅较高分辨率的图像。

最临近插值法原理：这种算法就是根据原图像和目标图像的尺寸，计算缩放的比例，然后根据缩放比例计算目标像素所依据的原像素，过程中自然会产生小数，这时就采用四舍五入，取与这个点最相近的点。

图解如下：如果(i+u, j+v)落在A区，即u<0.5, v<0.5，则将左上角象素的灰度值赋给待求象素，同理，落在B区则赋予右上角的象素灰度值，落在C区则赋予左下角象素的灰度值，落在D区则赋予右下角象素的灰度值。

具体Matlab源代码实现：clear all;A=imread('12.png'); %读取图像信息imshow(A); %显示原图title('原图128*128');Row = size(A,1); Col = size(A,2); %图像行数和列数nn=8; %放大倍数m = round(nn*Row); %求出变换后的坐标的最大值n = round(nn*Col);B = zeros(m,n,3); %定义变换后的图像for i = 1 : mfor j = 1 : nx = round(i/nn); y = round(j/nn); %最小临近法对图像进行插值if x==0 x = 1; endif y==0 y = 1; endif x>Row x = Row; endif y>Col y = Col;endB(i,j,:) = A(x,y,:);endendB = uint8(B); %将矩阵转换成8位无符号整数figure;imshow(B); %显示输出图片title('最邻近插值法放大8倍1024*1024');运行程序后，原图如图1所示：图1用最邻近插值法放大8倍后的图如图2所示：图2双线性内插值法：计算过程简单了解，如图，已知Q12，Q22，Q11，Q21，但是要插值的点为P 点，这就要用双线性插值了，首先在x轴方向上，对R1和R2两个点进行插值，这个很简单，然后根据R1和R2对P点进行插值，这就是所谓的双线性插值。

三次卷积内插法计算例题三次卷积内插法是一种用于数据插值的方法，它可以通过已知的数据点来估计在两个数据点之间的未知数据点的值。

下面我将以一个简单的例题来说明三次卷积内插法的计算过程。

假设我们有以下已知数据点：点1，(x1, y1) = (2, 4)。

点2，(x2, y2) = (4, 8)。

点3，(x3, y3) = (6, 6)。

现在我们要计算在点1和点2之间的一个新数据点的值，假设这个新数据点的横坐标为x。

首先，我们需要计算出x相对于已知数据点的位置。

这可以通过计算x与x1之间的距离来实现。

假设这个距离为h，则有：h = x x1。

接下来，我们需要计算三个权重系数：a0、a1、a2。

这些系数与已知数据点的纵坐标相关。

具体计算公式如下：a0 = (-y1 + 3y2 3y3 + y4) / 6。

a1 = (3y1 6y2 + 3y3) / 6。

a2 = (-3y1 + 3y3) / 6。

其中，y4为点4的纵坐标，我们需要根据已知数据点来确定它的值。

然后，我们可以计算出在点1和点2之间的新数据点的纵坐标y的估计值。

具体计算公式如下：y = a0 + a1 h + a2 h^2。

最后，我们可以得到在点1和点2之间的新数据点的坐标为(x, y)。

需要注意的是，三次卷积内插法要求已知数据点的横坐标必须是等间距的，即x2 x1 = x3 x2。

如果不满足这个条件，可以通过对数据进行插值或者进行其他的数据处理方法来满足这个条件。

这就是使用三次卷积内插法计算例题的详细过程。

希望这个回答能够满足你的需求。

如果你还有其他问题，欢迎继续提问。

图像重采样主要有三种方法，分别是最邻近法，双线性内插法和三次卷积内插法。

⑴最近邻法。

该法针对于二维图像取待采样点周围4个相邻像素点中距离最近的1个邻点的灰度值作为该点的灰度值〞如图〔1 〕。

此算法虽然计算简单，但由于仅用对该采样点影响最大的〔即最近的〕像素的灰度值作为该点的值，而没有考虑其他相邻像素的影响〔相关性〕，因此重新采样后的图像灰度值有明显的不连续性，像质损失较大。

图〔1〕图像缩放中的插值和重采样〔2〕双线性内插法作为对最近邻点法的一种改良，这种方法是“利用周围4个邻点的灰度值在两个方向上作线性内插以得到待采样点的灰度值〞。

即根据待采样点与相邻点的距离确定相应的权值计算出待采样点的灰度值。

双线性内插的示意图如图2所示，其中X、Y坐标表示像素的位置，f〔*，*〕表示像素的灰度值。

其数学表达式为：f 〔i+u，j+v〕= 〔1-u〕〔1-v〕f 〔i,j〕+〔1-u〕vf〔i,j+1〕+u〔1-v〕f〔i+1,j〕+uvf〔i+1,j+1〕〔2〕与最邻近法相比。

双线性内插法由于考虑了待采样点周围四个直接邻点对待采样点的影响，此根本克服了前者灰度不连续的缺点，但其代价是计算量有所增大。

但由于此方法仅考虑四个直接邻点灰度值的影响，而未考虑到各邻点间灰度值变化率的影响，因此具有低通滤波器的性质，使缩放后图像的高频分量受到损失，图像的轮廓变得较模糊。

用此方法缩放后的图像与原图像相比，仍然存在由于计算模型考虑不周而产生的图像质量退化与精度降低的问题。

〔3〕立方卷积法作为对双线性内插法的改良，即“不仅考虑到四个直接邻点灰度值的影响，还考虑到各邻点间灰度值变化率的影响〞，立方卷积法利用了待采样点周围更大邻域内像素的灰度值作三次插值。

此法利用了如图3所示的三次多项式S〔w〕。

S (w )的数学表达式为：I — 2 w p ui 11I xv |V 1 S<w) .4- 8 | w | + 5 | w P T w 1」I Cl w |C 2 0I u. |>2 式中，w 为自变量，S 〔w 〕为三次多项式的值。