三次卷积插值
三次卷积插值法
三次卷积插值法是一种常用的图像插值方法，用于将低分辨率的图像放大为高分辨率的图像。
它基于插值的原理，通过对原始图像的像素进行计算和重新分配，来生成更精细的图像。
具体而言，三次卷积插值法是通过对原始图像上的像素进行插值计算，生成新的像素值。
它利用相邻像素之间的灰度值关系，根据像素之间的距离和权重进行插值计算。
通常，三次卷积插值法会考虑每个像素周围的16个邻域像素的灰度值，根据这些灰度值的权重进行加权平均，来计算新的像素值。
三次卷积插值法具有较好的图像保真度和平滑性，可以有效地减少图像放大过程中的锯齿和失真。
然而，它也有一些局限性，如可能引入一定的模糊效果，对于一些细节和纹理的保留可能不如其他高级插值方法。
总的来说，三次卷积插值法是一种常用的图像放大方法，可以在一定程度上提高图像的分辨率和质量。
但在具体应用中，还需要根据实际需求和图像特点选择合适的插值方法。
'cubic' 基于三角形的三次插补法
在数学和工程领域中，三次插值是一种常见的技术，用于找到曲线或曲面上不同点之间的未知数值。
而基于三角形的三次插值法，通过构建和连接三角形的顶点来实现插值。
在本文中，我们将探讨基于三角形的三次插值法在立方体（cubic）中的应用，以及其在数学和工程领域中的重要性。
1. 立方体中的三次插值在立方体中，三次插值法可以通过连接三维空间中的顶点来实现。
具体而言，我们可以将立方体分割成许多小的三角形，在每个三角形中进行三次插值，然后将它们连接起来以获得整个立方体的插值。
这种基于三角形的方法可以帮助我们更准确地对立方体中的数据进行插值，从而在数学建模和工程设计中发挥重要作用。
2. 应用领域基于三角形的三次插值法在实际应用中有着广泛的领域，其中包括计算机图形学、地质勘探、医学成像等。
在计算机图形学中，立方体中的三次插值可以帮助我们生成更加真实和精细的图像，使得渲染和动画的效果更加逼真。
在地质勘探中，立方体中的三次插值可以帮助地质学家准确地估计地下资源的分布情况，指导资源的开发和利用。
在医学成像中，立方体中的三次插值可以帮助医生更准确地诊断病变部位，提高医疗诊断的准确性和可靠性。
3. 个人观点作为一名专业的文章写手，我对立方体中的三次插值法有着浓厚的兴趣。
这种技术不仅在学术研究中有着重要的应用，而且在工程实践中也发挥着巨大的作用。
通过对立方体中的数据进行插值，我们可以更好地理解和利用数据的潜在信息，为现实生活和工作带来更多便利和创新。
总结回顾在本文中，我们探讨了基于三角形的三次插值法在立方体中的重要应用。
立方体中的三次插值法通过连接三维空间中的顶点，实现了对立方体中数据的准确插值，为数学建模和工程设计提供了重要支持。
我们还介绍了该方法在计算机图形学、地质勘探和医学成像等领域的应用，并共享了个人对这一技术的观点和理解。
我们希望读者通过本文的阅读，能够更全面、深刻和灵活地理解立方体中的三次插值法，以及其重要的应用意义。
GIS重采样
ArcGIS中栅格数据重采样方法作者:赵立超邮箱：zhaolc@ 发布时间：2014-01-02 21:04:28 阅读(800)ArcGIS中栅格数据重采样方法处理栅格数据时，由于数据像元大小不符合要求，或者在进行栅格数据配准后，像元发生倾斜，或者对多个栅格数据进行分析时，需要使用相同的栅格分辨率，重采样到同一分辨率下，因此对栅格数据操作时经常要进行重采样操作。
如下图所示，通过重采样将栅格调整到新的分辨率下。
1.1.栅格重采样方法栅格重采样主要包括三种方法：最邻近法、双线性内插法和三次卷积插值法。
最邻近法是把原始图像中距离最近的像元值填充到新图像中；双线性内插法和三次卷积插值法都是把原始图像附近的像元值通过距离加权平均填充到新图像中。
默认情况下，采用最近邻分配重采样技术，这种方法同时适用于离散和连续值类型，而其他重采样方法只适用于连续数据。
1.1 最近邻法最邻近分配法是用于离散（分类）数据的重采样技术，因为它不会更改输入单元的值。
将输出栅格数据集中单元中心的位置定位到输入栅格后，最邻近分配法将确定输入栅格上最近的单元中心位置并将该单元的值分配给输出栅格上的单元。
最邻近分配法不会更改输入栅格数据集中单元的任何值。
输入栅格中的值2 在输出栅格中仍将为 2，决不会为 2.2 或 2.3。
由于输出单元值保持不变，因此最邻近分配法应该用于名目数据或顺序数据，其中每个值都表示一个类、一个成员或一个分类（分类数据，如土地利用、土壤或森林类型）。
考虑到根据输入栅格创建的输出栅格会在操作中旋转 45°，因此将进行重采样。
对于每个输出单元，都要从输入栅格中获取值。
在下图中，输入栅格的单元中心为灰色点。
输出单元为绿色阴影。
要处理的单元为黄色阴影。
在最邻近分配法中，将确定与要处理的单元中心（红色点）最邻近的输入栅格单元中心（橙色点)，并将其指定为要处理的单元（黄色阴影）的输出值。
对输出栅格中的每个单元都重复此过程。
三次样条曲面插值原理
三次样条曲面插值原理
三次样条曲面插值是一种用于构造二维曲面的插值方法。
其基本原理是通过已知的曲面上的若干点，计算出该曲面上的三次多项式函数，从而实现曲面的插值。
具体来说，三次样条曲面插值的原理如下：
1. 确定曲面上的插值节点：根据给定的曲面上的点的坐标，确定曲面上的插值节点。
2. 构造曲面的参数方程：利用插值节点，构造出曲面的参数方程。
三次样条曲面插值通常使用双变量的三次多项式作为参数方程。
参数方程的形式可以是Bézier曲面、B样条曲面等。
3. 确定曲面上的插值条件：根据已知的曲面上的点的坐标和曲面方程，确定曲面上的插值条件。
通常使用平滑条件（曲面上的点的正切方向相等）和代数条件（曲面上的点的坐标满足给定的条件）来确定插值条件。
4. 求解参数方程的系数：根据插值条件，求解参数方程中的系数。
可以使用线性代数的方法求解系数矩阵，得到曲面的参数化表达式。
5. 计算曲面上的点的坐标：利用参数方程和求解得到的系数，计算曲面上的点的坐标。
可以通过插值节点上的参数值，使用参数方程计算得到。
通过以上步骤，就可以构造出满足给定插值条件的三次样条曲面，从而实现曲面的插值。
高质量的快速图像缩放——二次线性插值和三次卷积插值
高质量的快速的图像缩放——二次线性插值和三次卷积插值限制条件:为了便于讨论，这里只处理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);}}。
图像缩小方法
（2.3.12）
21
⎡33 39 15 21⎤ ⎥ G=⎢ ⎢35 11 17 23⎥ ⎢ ⎣36 12 18 34⎥ ⎦
2 图像的比例放大变换
（2.3.13）
图像在缩小操作中， 是在现有的信息里如何挑选所需要的有用信息。 而在图像的放大操 作中，则需要对尺寸放大后所多出来的空格填入适当的像素值，这是信息的估计问题，所以 较图像的缩小要难一些。 由于图像相邻像素之间的相关性很强， 可以利用这个相关性来实现 图像的放大。与图像缩小相同，按比例放大不会引起图像的畸变，而不按比例放大则会产生 图像的畸变。图像放大一般采用最近邻域法、线性插值法、三次卷积法。 （1）最近邻域法 一般地，按比例将原图像放大 k 倍时，如果按照最近邻域法，则需要将一个像素值添在 新图像的 k × k 的子块中。式（2.3.14）为图像 F 的矩阵，该图像放大 3 倍得到图像 G 的矩 阵用式（2.3.15）表示。图 2.3.3 为放大 5 倍的示意图。显然，如果放大倍数太大，按照这种 方法处理会出现马赛克效应。
⎧ 1 − 2x2 + x 3 0 ≤ x ≤1 ⎪ 3 ⎪ ω ( x) = ⎨4 − 8 x + 5 x 2 − x 1 ≤ x ≤ 2 ⎪ 0 2≤ x ⎪ ⎩
三次多项式 ω ( x) 近似表示灰度内插时周围像元的灰度值对内插点灰度值的贡献大小。 可先在某一方向上内插， 如先在 X 方向上， 每 4 个值依次内插 4 次， 求出 f ( x, j − 1) ，f ( x, j ) ，
⎡ f11 ⎢f F = ⎢ 21 ⎢ f 31 ⎢ ⎣ f 41
f 12 f 22 f 32 f 42
f 13 f 23 f 33 f 43
f 14 f 24 f 34 f 44
三次样条插值算法详解
三次样条插值算法要求数据点数量较多，且在某些情况下可能存在数值不稳定性，如数据 点过多或数据点分布不均等情况。此外，该算法对于离散数据点的拟合效果可能不如其他 插值方法。
对未来研究的展望
01
02
03
改进算法稳定性
针对数值不稳定性问题， 未来研究可以探索改进算 法的数值稳定性，提高算 法的鲁棒性。
3
数据转换
对数据进行必要的转换，如标准化、归一化等， 以适应算法需求。
构建插值函数
确定插值节点
根据数据点确定插值节点，确保插值函数在节点处连续且光滑。
构造插值多项式
根据节点和数据点，构造三次多项式作为插值函数。
确定边界条件
根据实际情况确定插值函数的边界条件，如周期性、对称性等。
求解插值函数
求解线性方程组
06
结论
三次样条插值算法总结
适用性
三次样条插值算法适用于各种连续、光滑、可微的分段函数插值问题，尤其在处理具有复 杂变化趋势的数据时表现出色。
优点
该算法能够保证插值函数在分段连接处连续且具有二阶导数，从而在插值过程中保持数据 的平滑性和连续性。此外，三次样条插值算法具有简单、易实现的特点，且计算效率较高 。
根据数据点的数量和分布，合理分段，确保 拟合的精度和连续性。
求解线性方程组
使用高效的方法求解线性方程组，如高斯消 元法或迭代法。
结果输出
输出拟合得到的插值函数，以及相关的误差 分析和图表。
03
三次样条插值算法步骤
数据准备
1 2
数据收集
收集需要插值的原始数据点，确保数据准确可靠。
数据清洗
对数据进行预处理，如去除异常值、缺失值处理 等。
插值法在图像处理中的运用要点
插值方法在图像处理中的应用作者：专业姓名学号控制工程陈龙斌控制工程陈少峰控制工程殷文龙摘要本文介绍了插值方法在图像处理中的应用。
介绍了典型的最近邻插值、双线性插值、双三次插值、双信道插值、分形插值的原理。
以分形插值为重点，在图像放大领域用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.最近邻插值（Nearest-neighbor interpolation）：
最近邻插值是一种简单的插值方法，它将目标像素的值设定为原始图像中最近的像素的值。
这种方法计算速度快，但可能导致图像出现锯齿状的块状效应。
2.双线性插值（Bilinear interpolation）：
双线性插值使用了目标像素周围的四个最近邻像素的权重来计算目标像素的值。
这种方法相比最近邻插值产生了更平滑的图像，但在一些情况下可能引入一些模糊。
3.双三次插值（Bicubic interpolation）：
双三次插值在计算目标像素值时使用了周围16个最近邻像素的权重。
相较于双线性插值，双三次插值可以提供更高的图像质量，但计算复杂度也更高。
4.立方卷积插值（Cubic convolution interpolation）：
立方卷积插值是一种使用像素值的立方函数的插值方法，对图像进行平滑处理。
它相对于双三次插值在一些情况下可能产生更好的结果。
nczos 插值：
Lanczos 插值使用sinc 函数的权重来计算目标像素的值，是一种高质量的插值方法。
它通常用于需要较高图像质量的应用，如图像放大。
在选择重采样方法时，需考虑应用场景、计算资源和图像质量要求。
不同的插值方法在不同的情况下表现良好，因此最佳的选择可能取决于具体的应用需求。
常见插值算法--拉格朗日插值、三次卷积插值、三次样条插值、兰克索斯插值
常见插值算法--拉格朗⽇插值、三次卷积插值、三次样条插值、兰克索斯插值写在前⾯本⽂简单介绍了⼏种常见的插值算法并附带了相应的python代码，本⽂公式使⽤latex编写，如有错误欢迎评论指出，如果谁知道如何修改latex字号也欢迎留⾔关于⼀维、⼆维和多维插值三次卷积插值、拉格朗⽇两点插值（线性插值）、兰克索斯插值在⼆维插值时改变x和y⽅向的计算顺序不影响最终结果，这三个也是图像缩放插值时常⽤的插值算法，⽽其他插值在改变计算顺序时会产⽣明显差异，多维的情况笔者没有尝试，读者可以⾃⾏尝试或推导最近邻插值法（Nearest Neighbour Interpolation）在待求像素的四邻像素中，将距离待求像素最近的像素值赋给待求像素p_{11}p_{12}pp_{21}p_{22}python代码1def NN_interpolation(srcImg, dstH, dstW):2 scrH, scrW, _ = srcImg.shape3 dstImg = np.zeros((dstH, dstW, 3), dtype=np.uint8)4for i in range(dstH - 1):5for j in range(dstW - 1):6 scrX = round(i * (scrH / dstH))7 scrY = round(j * (scrW / dstW))8 dstImg[i, j] = srcImg[scrX, scrY]9return dstImg拉格朗⽇插值（Lagrange Interpolation）拉格朗⽇插值法需要找到k个p_i(x)函数，使得每个函数分别在在x_i处取值为1，其余点取值为0，则y_ip_i(x)可以保证在x_i处取值为y_i，在其余点取值为0，因此L_k(x)能恰好经过所有点，这样的多项式被称为拉格朗⽇插值多项式，记为L_k(x)=\sum_{i=1}^ky_ip_i(x)p_i(x)=\prod_{j \neq i}^{1 \leq j \leq k}\frac{x-x_j}{x_i-x_j}以四点即三次图像插值为例，因为横坐标间隔为1，则设四个点横坐标为-1、0、1和2，可得p_1(x)、p_2(x)、p_3(x)和p_4(x)假设y_1、y_2、y_3和y_4分别为1、2、-1、4，则可得拉格朗⽇函数如下图所⽰，待插值点横坐标范围为[0,1]在K=2时在k=2时，也被称为线性插值通⽤公式p_1=\frac{x-x_2}{x_1-x_2}p_2=\frac{x-x_1}{x_2-x_1}\begin{align} L_2x &= p_1y_1+p_2y_2 \nonumber \\ &= \frac{x-x_2}{x_1-x_2}y_1 + \frac{x-x_1}{x_2-x_1}y_2 \nonumber \end{align}图像插值像素分布如图所⽰p_{11}p_{12}pp_{21}p_{22}即当x_{i+1}=x_i+1时，设p与p_{11}的横纵坐标差分别为dx和dy\begin{align} L_2x &= \frac{x-x_2}{x_1-x_2}y_1 + \frac{x-x_1}{x_2-x_1}y_2 \nonumber \\ &= (x_2-x)y_1+(x-x_1)y_2 \nonumber \\ &= (1-dx)y_1+dxy_2 \nonumber \\ &= (y_2-y_1)dx+y_1 \nonumber \end{align}L_2'x=y_2-y_1在K=3时通⽤公式p_1=\frac{x-x_2}{x_1-x_2}\frac{x-x_3}{x_1-x_3}p_2=\frac{x-x_1}{x_2-x_1}\frac{x-x_3}{x_2-x_3}p_3=\frac{x-x_1}{x_3-x_1}\frac{x-x_2}{x_3-x_2}\begin{align} L_3x &= p_1y_1+p_2y_2+p_3y_3 \nonumber \\ &= \frac{x-x_2}{x_1-x_2}\frac{x-x_3}{x_1-x_3}y_1+\frac{x-x_1}{x_2-x_1}\frac{x-x_3}{x_2-x_3}y_2+\frac{x-x_1}{x_3-x_1}\frac{x-x_2}{x_3-x_2}y_3 \nonumber \end{align}图像插值像素分布如图所⽰p_{11}p_{12}p_{13}p_{21}p_{22}p_{23}pp_{31}p_{32}p_{33}即当x_{i+1}=x_i+1时，设p与p_{11}的横纵坐标差分别为dx和dy\begin{align} L_3x &= \frac{x-x_2}{x_1-x_2}\frac{x-x_3}{x_1-x_3}y_1 + \frac{x-x_1}{x_2-x_1}\frac{x-x_3}{x_2-x_3}y_2 + \frac{x-x_1}{x_3-x_1}\frac{x-x_2}{x_3-x_2}y_3 \nonumber \\ &= \frac{-dx(1-dx)}{(-1)\cdot(-2)}y_1 + \frac{-(1+dx)(1-dx)}{1\cdot(-1)}y_2 + \frac{(1+dx)dx}{2\cdot 1}y_3 \nonumber \\ &= (\frac{1}{2}d^2x-\frac{1}{2}dx)y_1 - (d^2x-1)y_2 + (\frac{1}{2}d^2x+\frac{1}{2}dx)y_3 \nonumber \\ &= d^2x(\frac{1}{2}y_1-y_2+\frac{1}{2}y_3)+dx(-\frac{1}{2}y_1+\frac{1}{2}y_3)+y_2 \nonumber \end{align}L_3'x=dx(y_1-2y_2+y_3)+(\frac{1}{2}y_3-\frac{1}{2}y_1)在K=4时通⽤公式p_1=\frac{x-x_2}{x_1-x_2}\frac{x-x_3}{x_1-x_3}\frac{x-x_4}{x_1-x_4}p_2=\frac{x-x_1}{x_2-x_1}\frac{x-x_3}{x_2-x_3}\frac{x-x_4}{x_2-x_4}p_3=\frac{x-x_1}{x_3-x_1}\frac{x-x_2}{x_3-x_2}\frac{x-x_4}{x_3-x_4}p_4=\frac{x-x_1}{x_4-x_1}\frac{x-x_2}{x_4-x_2}\frac{x-x_3}{x_4-x_3}\begin{align} L_4x &= p_1y_1+p_2y_2+p_3y_3+p_4y_4 \nonumber \\ &= \frac{x-x_2}{x_1-x_2}\frac{x-x_3}{x_1-x_3}\frac{x-x_4}{x_1-x_4}y_1 + \frac{x-x_1}{x_2-x_1}\frac{x-x_3} {x_2-x_3}\frac{x-x_4}{x_2-x_4}y_2 + \frac{x-x_1}{x_3-x_1}\frac{x-x_2}{x_3-x_2}\frac{x-x_4}{x_3-x_4}y_3 + \frac{x-x_1}{x_4-x_1}\frac{x-x_2}{x_4-x_2}\frac{x-x_3}{x_4-x_3}y_4\nonumber \end{align}图像插值p_{11}p_{12}p_{13}p_{14}p_{21}p_{22}p_{23}p_{24}pp_{31}p_{32}p_{33}p_{34}p_{41}p_{42}p_{43}p_{44}即当x_{i+1}=x_i+1时，设p与p_{11}的横纵坐标差分别为dx和dy\begin{align} L_4x &= \frac{x-x_2}{x_1-x_2}\frac{x-x_3}{x_1-x_3}\frac{x-x_4}{x_1-x_4}y_1 + \frac{x-x_1}{x_2-x_1}\frac{x-x_3}{x_2-x_3}\frac{x-x_4}{x_2-x_4}y_2 + \frac{x-x_1}{x_3-x_1}\frac{x-x_2}{x_3-x_2}\frac{x-x_4}{x_3-x_4}y_3 + \frac{x-x_1}{x_4-x_1}\frac{x-x_2}{x_4-x_2}\frac{x-x_3}{x_4-x_3}y_4 \nonumber \\ &= \frac{dx[-(1-dx)][-(2-dx)]}{(-1)\cdot(-2)\cdot(-3)}y_1 + \frac{(1+dx)[-(1-dx)][-(2-dx)]}{1\cdot(-1)\cdot(-2)}y_2 + \frac{(1+dx)dx[-(2-dx)]}{2\cdot 1\cdot(-1)}y_3 + \frac{(1+dx)dx[-(1-dx)]}{3\cdot 2\cdot 1}y_4 \nonumber \\ &= \frac{d^3x-3d^2x+2dx}{-6}y1 + \frac{d^3x-2d^2x-dx+2}{2}y_2 + \frac{d^3x-d^2x-2dx}{-2}y_3 + \frac{d^3x-dx}{6}y_4 \nonumber \\ &= d^3x(-\frac{1}{6}y_1+\frac{1}{2}y_2-\frac{1} {2}y_3+\frac{1}{6}y_4)+d^2x(\frac{1}{2}y_1-y_2+\frac{1}{2}y_3)+dx(-\frac{1}{3}y_1-\frac{1}{2}y_2+y_3-\frac{1}{6}y_4)+y_2 \nonumber \end{align}\begin{align} L_4'x &= d^2x(-\frac{1}{2}y_1+\frac{3}{2}y_2-\frac{3}{2}y_3+\frac{1}{2}y_4)+dx(y_1-2y_2+y_3)+(-\frac{1}{3}y_1-\frac{1}{2}y_2+y_3-\frac{1}{6}y_4) \nonumber \\ &= -[\frac{1}{2}d^2x(y_1-3y_2+3y_3-y_4)-dx(y_1-2y_2+y_3)+\frac{1}{6}(2y_1+3y_2-6y_3+y_4)] \nonumber \end{align}python代码插值核计算的时候乘法和加减法计算的顺序不同可能会导致结果存在细微的差异，读者可以⾃⾏研究⼀下1class BiLagrangeInterpolation:2 @staticmethod3def LagrangeInterpolation2(x, y1, y2):4 f1 = 1 - x5 f2 = x6 result = y1 * f1 + y2 * f27return result89 @staticmethod10def LagrangeInterpolation3(x, y1, y2, y3):11 f1 = (x ** 2 - x) / 2.012 f2 = 1 - x ** 213 f3 = (x ** 2 + x) / 2.014 result = y1 * f1 + y2 * f2 + y3 * f315return result1617 @staticmethod18def LagrangeInterpolation4(x, y1, y2, y3, y4):19 f1 = - (x ** 3 - 3 * x ** 2 + 2 * x) / 6.020 f2 = (x ** 3 - 2 * x ** 2 - x + 2) / 2.021 f3 = - (x ** 3 - x ** 2 - 2 * x) / 2.022 f4 = (x ** 3 - x) / 6.023 result = y1 * f1 + y2 * f2 + y3 * f3 + y4 * f424return result2526def biLag2_2(self, srcImg, dstH, dstW):27 dstH, dstW = int(dstH), int(dstW)28 srcH, srcW, _ = srcImg.shape29 srcImg = np.pad(srcImg, ((1, 1), (1, 1), (0, 0)), 'edge')30 dstImg = np.zeros((dstH, dstW, 3), dtype=np.uint8)31for dstY in range(dstH):32for dstX in range(dstW):33for channel in [0, 1, 2]:34# p11 p1235# p36# p21 p2237# 储存为 p(y, x)38 p = [dstY * srcH / dstH, dstX * srcW / dstW]39 p11 = [math.floor(p[0]), math.floor(p[1])]40 p12 = [p11[0], p11[1] + 1]4142 p21 = [p11[0] + 1, p11[1]]43 p22 = [p21[0], p12[1]]4445 diff_y, diff_x = p[0] - p11[0], p[1] - p11[1]46 r1 = grangeInterpolation2(diff_x, srcImg[p11[0], p11[1], channel], srcImg[p12[0], p12[1], channel])47 r2 = grangeInterpolation2(diff_x, srcImg[p21[0], p21[1], channel], srcImg[p22[0], p22[1], channel])4849 c = grangeInterpolation2(diff_y, r1, r2)5051 dstImg[dstY, dstX, channel] = np.clip(c, 0, 255)52return dstImg5354def biLag3_3(self, srcImg, dstH, dstW):55 dstH, dstW = int(dstH), int(dstW)56 srcH, srcW, _ = srcImg.shape57 srcImg = np.pad(srcImg, ((1, 1), (1, 1), (0, 0)), 'edge')58 dstImg = np.zeros((dstH, dstW, 3), dtype=np.uint8)59for dstY in range(dstH):60for dstX in range(dstW):61for channel in [0, 1, 2]:62# p11 p12 p1363#64# p21 p22 p2365# p66# p31 p32 p3367# 储存为 p(y, x)68 p = [dstY * srcH / dstH, dstX * srcW / dstW]69 p22 = [math.floor(p[0]), math.floor(p[1])]70 p21 = [p22[0], p22[1] - 1]71 p23 = [p22[0], p22[1] + 1]7273 p11 = [p21[0] - 1, p21[1]]74 p12 = [p11[0], p22[1]]75 p13 = [p11[0], p23[1]]7677 p31 = [p21[0] + 1, p21[1]]78 p32 = [p31[0], p22[1]]79 p33 = [p31[0], p23[1]]8081 diff_y, diff_x = p[0] - p22[0], p[1] - p22[1]82 r1 = grangeInterpolation3(diff_x, srcImg[p11[0], p11[1], channel], srcImg[p12[0], p12[1], channel], srcImg[p13[0], p13[1], channel])83 r2 = grangeInterpolation3(diff_x, srcImg[p21[0], p21[1], channel], srcImg[p22[0], p22[1], channel], srcImg[p23[0], p23[1], channel])84 r3 = grangeInterpolation3(diff_x, srcImg[p31[0], p31[1], channel], srcImg[p32[0], p32[1], channel], srcImg[p33[0], p33[1], channel]) 8586 c = grangeInterpolation3(diff_y, r1, r2, r3)8788 dstImg[dstY, dstX, channel] = np.clip(c, 0, 255)89return dstImg9091def biLag4_4(self, srcImg, dstH, dstW):92 dstH, dstW = int(dstH), int(dstW)93 srcH, srcW, _ = srcImg.shape94 srcImg = np.pad(srcImg, ((1, 2), (1, 2), (0, 0)), 'edge')95 dstImg = np.zeros((dstH, dstW, 3), dtype=np.uint8)96for dstY in range(dstH):97for dstX in range(dstW):98for channel in [0, 1, 2]:99# p11 p12 p13 p14100#101# p21 p22 p23 p24102# p103# p31 p32 p33 p34104#105# p41 p42 p43 p44106# 储存为 p(y, x)107 p = [dstY * srcH / dstH, dstX * srcW / dstW]108 p22 = [math.floor(p[0]), math.floor(p[1])]109 p21 = [p22[0], p22[1] - 1]110 p23 = [p22[0], p22[1] + 1]111 p24 = [p22[0], p22[1] + 2]112113 p11 = [p21[0] - 1, p21[1]]114 p12 = [p11[0], p22[1]]115 p13 = [p11[0], p23[1]]116 p14 = [p11[0], p24[1]]117118 p31 = [p21[0] + 1, p21[1]]119 p32 = [p31[0], p22[1]]120 p33 = [p31[0], p23[1]]121 p34 = [p31[0], p24[1]]122123 p41 = [p21[0] + 2, p21[1]]124 p42 = [p41[0], p22[1]]125 p43 = [p41[0], p23[1]]126 p44 = [p41[0], p24[1]]127128 diff_y, diff_x = p[0] - p22[0], p[1] - p22[1]129 r1 = grangeInterpolation4(diff_x, srcImg[p11[0], p11[1], channel], srcImg[p12[0], p12[1], channel], srcImg[p13[0], p13[1], channel], srcImg[p14[0], p14[1], channel]) 130 r2 = grangeInterpolation4(diff_x, srcImg[p21[0], p21[1], channel], srcImg[p22[0], p22[1], channel], srcImg[p23[0], p23[1], channel], srcImg[p24[0], p24[1], channel]) 131 r3 = grangeInterpolation4(diff_x, srcImg[p31[0], p31[1], channel], srcImg[p32[0], p32[1], channel], srcImg[p33[0], p33[1], channel], srcImg[p34[0], p34[1], channel]) 132 r4 = grangeInterpolation4(diff_x, srcImg[p41[0], p41[1], channel], srcImg[p42[0], p42[1], channel], srcImg[p43[0], p43[1], channel], srcImg[p44[0], p44[1], channel]) 133134 c = grangeInterpolation4(diff_y, r1, r2, r3, r4)135136 dstImg[dstY, dstX, channel] = np.clip(c, 0, 255)137return dstImg三次卷积插值法（Cubic Convolution Interpolation）使⽤上图中的卷积核进⾏加权平均计算，卷积核为u(s)，四个等距（距离为1）的采样点记为x_0、x_1、x_2和x_3，采样数值记为y_0、y_1、y_2和y_3，且保证四个点均在[-2,2]区间上，计算得到g(x)，假设y_1、y_2、y_3和y_4分别为1、2、-1、4，则可得三次卷积插值函数如下图所⽰，待插值点横坐标范围为[0,1]公式推导设u(s)=\begin{cases} A_1|s|^3+B_1|s|^2+C_1|s|+D_1, &0<|s|<1 \\ A_2|s|^3+B_2|s|^2+C_2|s|+D_2, &1<|s|<2 \\ 1, &s=0 \\ 0, &otherwise \end{cases}\because函数在s=0,1,2处连续\therefore\begin{cases} 1=u(0^+)=D_1 \\ 0=u(1^-)=A_1+B_1+C_1+D_1 \\ 0=u(1^+)=A_2+B_2+C_2+D_2 \\ 0=u(2^-)=8A_2+4B_2+2C_2+D_2 \end{cases} (1)\because函数在s=0,1,2处导函数连续\therefore\begin{cases} u'(0^-)=u'(0+) \\ u'(1^-)=u'(1+) \\ u'(2^-)=u'(2+)\end{cases} \Rightarrow \begin{cases} -C_1=C_1 \\ 3A_1+2B_1+C_1=3A_2+2B_2+C_2\\ 12A_2+4B_2+C+2=0 \end{cases} ~~~~ (2)联⽴⽅程组(1)(2)，设A_2=a，解得\begin{cases} A_1=a+2 \\ B_1=-(a+3) \\ C_1=0 \\ D_1=1 \\ A_2=a \\ B_2=-5a \\ C_2=8a \\ D_2=-4a \end{cases}\Rightarrow u(s)=\begin{cases} (a+2)|s|^3-(a+3)|s|^2+1, &0<|s|<1 \\ A_2|s|^3+B_2|s|^2+C_2|s|+D_2, &1<|s|<2\\ 1, &s=0 \\ 0, &otherwise \end{cases}\because g(x)=\sum_kC_ku(s+j-k), ~~~~k=j-1,j, j+1,j+2且0<s<1⼜\because \begin{cases}\begin{align} u(s+1)&=as^3-2as^2+as \nonumber \\ u(s)&=(a+2)s^3-(a+3)s^2+1 \nonumber \\ u(s-1)&=-(a+2)s^3+(2a+3)s^2-as \nonumber \\ u(s-2)&=-as^3+as^2 \nonumber \end{align}\end{cases}\begin{align} \therefore g(x) &= C_{j-1}u(s+1)+C_{j}u(s)+C_{j+1}u(s-1)+C_{j+2}u(s-2) \nonumber \\ &= C_{j-1}(as^3-2as^2+as)+C_j[(a+2)s^3-(a+3)s^2+1]+C_{j+1}[-(a+2)s^3+ (2a+3)s^2-as]+C_{j+2}[-a^3+as^2] \nonumber \\ &= s^3[aC_{j-1}+(a+2)C_j-(a+2)C_{j+1}-aC_{j+2}]+s^2[-2aC_{j-1}-(a+3)C_j+(2a+3)C_{j+1}+aC_{j+2}]+s[aC_{j-1}-aC_{j+1}]+C_j \nonumber \end{align} ~~(3)f在x_j处泰勒展开得到f(x)=f(x_j)+f'(x_j)(x-x_j)+\frac{1}{2}f''(x_j)(x-x_j)^2+\cdots\therefore \begin{cases} f(x_{j+1})=f(x_j)+f'(x_j)(x_{j+1}-x_j)+\frac{1}{2}f''(x_j)(x_{j+1}-x_j)^2+\cdots \\ f(x_{j+2})=f(x_j)+f'(x_j)(x_{j+2}-x_j)+\frac{1}{2}f''(x_j)(x_{j+2}-x_j)^2+\cdots \\ f(x_{j-1})=f(x_j)+f'(x_j)(x_{j-1}-x_j)+\frac{1}{2}f''(x_j)(x_{j-1}-x_j)^2+\cdots \end{cases}令x_{j+1}-x_j=h\because x_{i+1}=x_i+1\therefore x_{j+2}-x_j=2h，x_{j-1}-x_j=-h\therefore \begin{cases} f(x_{j+2})=f(x_j)+2f'(x_j)h+2f''(x_j)h^2+\cdots \\ f(x_{j+1})=f(x_j)+f'(x_j)h+\frac{1}{2}f''(x_j)h^2+\cdots \\ f(x_{j-1})=f(x_j)-f'(x_j)h+\frac{1}{2}f''(x_j)h^2+\cdots \end{cases}\therefore \begin{cases} c_{j-1}=f(x_j)-f'(x_j)h+\frac{1}{2}f''(x_j)h^2+o(h^3) \\ c_j=f(x_j) \\ c_{j+1}=f(x_j)+f'(x_j)h+\frac{1}{2}f''(x_j)h^2+o(h^3)\\ c_{j+2}=f(x_j)+2f'(x_j)h+2f''(x_j)h^2+o(h^3) \end{cases} ~~ (4)将(4)代⼊(3)，得g(x)=-(2a+1)[2hf'(x_j)+h^2f''(x_j)]s^3+[(6a+3)hf'(x_j)+\frac{4a+3}{2}h^2f''(x_j)]s^2-2ahf'(x_j)s+f(x_j)+o(h^3)\because h=1,s=x-x_J\therefore sh=x-x_j\begin{align}\therefore f(x)&= f(x_j)+f'(x_j)(x-x_j)+\frac{1}{2}f''(x_j)(x-x_j)^2+o(h^3) \nonumber \\ &= f(x_j)+f'(x_j)sh+\frac{1}{2}f''(x_j)s^2h^2+o(h^3) \nonumber \end{align}\therefore f(x)-g(x)=(2a+1)[2hf'(x_j)+h^2f''(x_j)]s^3-(2a+1)[3hf'(x_j)+h^2f''(x_j)]s^2+[(2a+1)hf'(x_j)]s+o(h^3)\because 期望f(x)-g(x)趋于0\therefore 2a+1=0 \Rightarrow a=-\frac{1}{2}\therefore u(s)=\begin{cases} \frac{3}{2}|s|^3-\frac{5}{2}|s|^2+1, &0<|s|<1 \\ -\frac{1}{2}|s|^3+\frac{5}{2}|s|^2-4|s|+2, &1<|s|<2 \\ 1, &s=0 \\ 0, &otherwise \end{cases}\therefore g(s)=s^3[-\frac{1}{2}c_{j-1}+\frac{3}{2}c_j-\frac{3}{2}c_{j+1}+\frac{1}{2}c_{j+2}]+s^2[c_{j-1}-\frac{5}{2}c_j+2c_{j+1}-\frac{1}{2}c_{j+2}]+s[-\frac{1}{2}c_{j-1}+\frac{1} {2}c_{j+1}]+c_j图像插值p_{11}p_{12}p_{13}p_{14}p_{21}p_{22}p_{23}p_{24}pp_{31}p_{32}p_{33}p_{34}p_{41}p_{42}p_{43}p_{44}python代码1class BiCubicConvInterpolation:2 @staticmethod3def CubicConvInterpolation1(p0, p1, p2, p3, s):4# ⽤g（s）公式计算，已经将四个u（s）计算完毕并整理5# as^3 + bs^2 + cs + d6 a = 0.5 * (-p0 + 3.0 * p1 - 3.0 * p2 + p3)7 b = 0.5 * (2.0 * p0 - 5.0 * p1 + 4.0 * p2 - p3)8 c = 0.5 * (-p0 + p2)9 d = p110return d + s * (c + s * (b + s * a))1112 @staticmethod13def CubicConvInterpolation2(s):14# ⽤u（s）公式计算15 s = abs(s)16if s <= 1:17return 1.5 * s ** 3 - 2.5 * s ** 2 + 118elif s <= 2:19return -0.5 * s ** 3 + 2.5 * s ** 2 - 4 * s + 220else:21return 02223def biCubic1(self, srcImg, dstH, dstW):24# p11 p12 p13 p1425#26# p21 p22 p23 p2427# p28# p31 p32 p33 p3429#30# p41 p42 p43 p4431 dstH, dstW = int(dstH), int(dstW)32 scrH, scrW, _ = srcImg.shape33 srcImg = np.pad(srcImg, ((1, 1), (1, 1), (0, 0)), 'edge')34 dstImg = np.zeros((dstH, dstW, 1), dtype=np.uint8)35for dstY in range(dstH):36for dstX in range(dstW):37for channel in [0]:38 y = dstY * scrH / dstH39 x = dstX * scrW / dstW40 y1 = math.floor(y)41 x1 = math.floor(x)4243 array = []44for i in [-1, 0, 1, 2]:45 temp = self.CubicConvInterpolation1(srcImg[y1 + i, x1 - 1, channel],46 srcImg[y1 + i, x1, channel],47 srcImg[y1 + i, x1 + 1, channel],48 srcImg[y1 + i, x1 + 2, channel],49 x - x1)50 array.append(temp)5152 temp = self.CubicConvInterpolation1(array[0], array[1], array[2], array[3], y - y1)53 dstImg[dstY, dstX, channel] = np.clip(temp, 0, 255)5455return dstImg5657def biCubic2(self, srcImg, dstH, dstW):58# p11 p12 p13 p1459#60# p21 p22 p23 p2461# p62# p31 p32 p33 p3463#64# p41 p42 p43 p4465 dstH, dstW = int(dstH), int(dstW)66 scrH, scrW, _ = srcImg.shape67 srcImg = np.pad(srcImg, ((1, 1), (1, 1), (0, 0)), 'edge')68 dstImg = np.zeros((dstH, dstW, 3), dtype=np.uint8)69for dstY in range(dstH):70for dstX in range(dstW):71for channel in [0, 1, 2]:72 y = dstY * scrH / dstH73 x = dstX * scrW / dstW74 y1 = math.floor(y)75 x1 = math.floor(x)7677 array = []78for i in [-1, 0, 1, 2]:79 temp = 080for j in [-1, 0, 1, 2]:81 temp += srcImg[y1 + i, x1 + j, channel] * self.CubicConvInterpolation2(x - (x1 + j))82 array.append(temp)8384 temp = 085for i in [-1, 0, 1, 2]:86 temp += array[i + 1] * self.CubicConvInterpolation2(y - (y1 + i))87 dstImg[dstY, dstX, channel] = np.clip(temp, 0, 255)8889return dstImg三次样条插值在n-1个区间上寻找n-1个三次曲线，使其满⾜相邻曲线在端点处值相等、⼀阶导数相等，⼆阶导数相等，在加以边界条件后可得每个曲线的⽅程，然后沿x轴依次偏移对应的距离即可得到插值结果，如仅需要特定范围内的结果，则可以⼤幅减少计算量公式推导设S_i(x)=a_i+b_i(x-x_i)+c_i(x-x_i)^2+d_i(x-x_i)^3, ~~~~i=0,1,...,n-1则 \begin{cases} S_i'(x)=b_i+2c_i(x-x_i)+3d_i(x-x_i)^2\\ S_i''(x)=2c_i+6d_i(x-x_i)\\ S_i'''(x)=6d_i\\ \end{cases} ~~~~i=0,1,...,n-1设h_i(x)=x_{i+1}-x_i，可得\begin{cases} S_i(x)=a_i+b_ih_i+c_ih_i^2+d_ih_i^3\\ S_i'(x)=b_i+2c_ih_i+3d_ih_i^2\\ S_i''(x)=2c_i+6d_ih_i\\ S_i'''(x)=6d_i\\ \end{cases} ~~~~i=0,1,...,n-1\because S_i(x)过点(x_i,y_i)\therefore S_i(x)=a_i=y+i, ~~~~i=0,1,...,n-1 ~~~~~~(1)\because S_i(x)与S_{i+1}(x)在X_{i+1}处相等\therefore S_i(x_{i+1})=S_{i+1}(x_{i+1})\Rightarrow a_i+b_ih_i+c_ih_i^2+d_ih_i^3=y_{i+1}, ~~~~i=0,1,...,n-2~~~~~~(2)\because S_i'(x)与S_{i+1}'(x)在X_{i+1}处相等\therefore S_i'(x)-S_{i+1}'(x)=0\Rightarrow b_i+2c_ih_i+3d_ih_i^2-b_{i+1}=0~~~~~~(3)\because S_i''(x)与S_{i+1}''(x)在X_{i+1}处相等\therefore S_i''(x)-S_{i+1}''(x)=0\Rightarrow 2c_i+6d_ih_i-2c_{i+1}=0, ~~~~i=0,1,...,n-2~~~~~~(4)设m_i=S_i(x_i)=2c_i，即c_i=\frac{1}{2}m_i, ~~~~i=0,1,...,n-1~~~~~~(5)将(5)代⼊(4)，得2c_i+6d_ih_i-2c_{i+1}=0\Rightarrow m_i+6h_id_i-m_{i+1}=0\Rightarrow d_i=\frac{m_{i+1}-m_i}{6h_i}, ~~~~i=0,1,...,n-2~~~~~~(6)将(1)(5)(6)代⼊(2)，得\begin{align} &a_i+b_ih_i+c_ih_i^2+d_ih_i^3=y_{i+1} \nonumber \\ \Rightarrow&y_i+b_ih_i+\frac{1}{2}m_ih_i^2+\frac{m_{i+1}-m_i}{6h_i}h_i^3=y_{i+1} \nonumber \\\Rightarrow&b_i=\frac{y_{i+1}-y_i}{h_i}-\frac{1}{2}m_ih_i-\frac{1}{6}(m_{i+1}-m_i)h_i \nonumber \\ \Rightarrow&b_i=\frac{y_{i+1}-y_i}{h_i}-\frac{1}{3}m_ih_i-\frac{1}{6}m_{i+1}h_i, ~~~~i=0,1,...,n-2~~~~~~(7) \nonumber \end{align}将(5)(6)(7)代⼊(3)，得\begin{align} &\frac{y_{i+1}-y{i}}{h_i}-\frac{1}{3}m_ih_i-\frac{1}{6}m_{i+1}h_i+2\cdot\frac{1}{2}m_ih_i+3\frac{m_{i+1}-m_i}{6h_i}h_i^2-(\frac{y_{i+2}-y_{i+1}}{h_{i+1}}-\frac{1}{3}m_{i+1}h_{i+1}-\frac{1}{6}m_{i+2}h_{i+1})=0 \nonumber \\ \Rightarrow&\frac{y_{i+1}-y{i}}{h_i}-\frac{1}{3}m_ih_i-\frac{1}{6}m_{i+1}h_i+m_ih_i+\frac{1}{2}(m_{i+1}-m_i)h_i-\frac{y_{i+2}-y_{i+1}}{h_{i+1}}+\frac{1}{3}m_{i+1}h_{i+1}+\frac{1}{6}m_{i+2}h_{i+1}=0 \nonumber \\ \Rightarrow&m_ih_i(-\frac{1}{3}+1-\frac{1}{2})+m_{i+1}h_i(-\frac{1}{6}+\frac{1} {2})+\frac{1}{3}m_{i+1}h_{i+1}+\frac{1}{6}m_{i+2}h_{i+1}=\frac{y_{i+2}-y_{i+1}}{h_{i+1}}-\frac{y_{i+1}-y_{i}}{h_{i}} \nonumber \\ \Rightarrow&\frac{1}{6}(m_ih_i+2m_{i+1}h_i+2m_{i+1}h_{i+1}+m_{i+2}h_{i+1})=\frac{y_{i+2}-y_{i+1}}{h_{i+1}}-\frac{y_{i+1}-y_{i}}{h_{i}} \nonumber \\ \Rightarrow&m_ih_i+2m_{i+1}(h_i+h_{i+1})+m_{i+2}h_{i+1}=6(\frac{y_{i+2}-y_{i+1}}{h_{i+1}}-\frac{y_{i+1}-y_{i}}{h_{i}}), ~~~~i=0,1,...,n-2~~~~~~(8) \nonumber \end{align}由(8)可知i=0,1,...,n-2，则有m_0,m_1,...,m_n，需要两个额外条件⽅程组才有解⾃然边界(Natural)m_0=0，m_n=0\begin{bmatrix} \tiny 1 & 0 & 0 & 0 & 0 & \cdots & 0\\ h_0 & 2(h_0+h_1) & h_1 & 0 & 0 & \cdots & 0\\ 0 & h_1 & 2(h_1+h_2) & h_2 & 0 & \cdots & 0\\ 0 & 0 & h_2 & 2(h_2+h_3) & h_3 & \cdots & 0\\ \vdots& & & \ddots & \ddots & \ddots & \vdots\\ 0 & \cdots & & & h_{n-2} & 2(h_{n-2}+h_{n-1}) & h_{n-1}\\ 0 & \cdots & & & 0 & 0 & 1 \end{bmatrix}\begin{bmatrix} m_0\\m_1\\m_2\\m_3\\\vdots\\m_{n-1}\\m_n \end{bmatrix}=6\begin{bmatrix} 0\\ \frac{y_2-y_1}{h_1}-\frac{y_1-y_0}{h_0}\\ \frac{y_3-y_2}{h_2}-\frac{y_2-y_1}{h_1}\\ \frac{y_4-y_3}{h_3}-\frac{y_3-y_2}{h_2}\\ \vdots\\ \frac{y_n-y_{n-1}}{h_{n-1}}-\frac{y_{n-1}-y_{n-2}}{h_{n-2}}\\ 0 \end{bmatrix}固定边界(Clamped)\begin{align} &\begin{cases} S_0'(x_0)=A\\ S_{n-1}'(x_n)=B \end{cases} \nonumber \\ \Rightarrow&\begin{cases} b_0=A\\ b_{n-1}+2c_{n-1}h_{n-1}+3d_{n-1}h_{n-1}^2=B\end{cases} \nonumber \\ \Rightarrow&\begin{cases} A=\frac{y_1-y_0}{h_0}-\frac{h_0}{2}m_0-\frac{h_0}{6}(m_1-m_0)\\ B=\frac{y_n-y_{n-1}}{h_{n-1}}-\frac{1}{3}m_{n-1}h_{n-1}+m_{n-1}h_{n-1}+\frac{1}{2}m_nh_{n-1}-\frac{1}{2}m_{n-1}h_{n-1} \end{cases} \nonumber \\ \Rightarrow&\begin{cases} 2h_0m_0+h_0m_1=6(\frac{y_1-y_0}{h_0}-A)\\ h_{n-1}m_{n-1}+2h_{n-1}m_{n}=6(B-\frac{y_n-y_{n-1}}{h_{n-1}}) \end{cases} \nonumber \\ \end{align}\begin{bmatrix} 2 & 1 & 0 & 0 & 0 & \cdots & 0\\ h_0 & 2(h_0+h_1) & h_1 & 0 & 0 & \cdots & 0\\ 0 & h_1 & 2(h_1+h_2) & h_2 & 0 & \cdots & 0\\ 0 & 0 & h_2 & 2(h_2+h_3) & h_3 & \cdots & 0\\ \vdots& & & \ddots & \ddots & \ddots & \vdots\\ 0 & \cdots & & & h_{n-2} & 2(h_{n-2}+h_{n-1}) & h_{n-1}\\ 0 & \cdots & & & 0 & 1 & 2 \end{bmatrix}\begin{bmatrix} m_0\\m_1\\m_2\\m_3\\\vdots\\m_{n-1}\\m_n \end{bmatrix}=6\begin{bmatrix} \frac{y_1-y_0}{h_0}-A\\ \frac{y_2-y_1}{h_1}-\frac{y_1-y_0}{h_0}\\ \frac{y_3-y_2}{h_2}-\frac{y_2-y_1}{h_1}\\ \frac{y_4-y_3}{h_3}-\frac{y_3-y_2}{h_2}\\ \vdots\\\frac{y_n-y_{n-1}}{h_{n-1}}-\frac{y_{n-1}-y_{n-2}}{h_{n-2}}\\ B-\frac{y_n-y_{n-1}}{h_{n-1}} \end{bmatrix}⾮节点边界(Not-A-Knot)\begin{align} &\begin{cases} S_0'''(x_1)=S_1'''(x_1)\\ S_{n-2}'''(x_{n-1})=S_{n-1}'''(x_{n-1}) \end{cases} \nonumber \\ \Rightarrow&\begin{cases} 6\cdot\frac{m_1-m_0}{6h_0}=6\cdot\frac{m_2-m_1}{6h_1}\\ 6\cdot\frac{m_{n-1}-m_{n-2}}{6h_{n-2}}=6\cdot\frac{m_n-m_{n-1}}{6h_{n-1}} \end{cases} \nonumber \\ \Rightarrow&\begin{cases} h_1(m_1-m_0)=h_0(m_2-m_1)\\ h_{n-1}(m_{n-1}-m_{n-2})=h_{n-2}(m_n-m_{n-1}) \end{cases} \nonumber \\ \Rightarrow&\begin{cases} -h_1m_0+(h_1+h_0)m_1-h_0m_2=0\\ -h_{n-1}m_{n-2}+(h_{n-1}+h_{n-2})m_{n-1}-h_{n-2}m_n=0 \end{cases} \nonumber \\ \end{align}\begin{bmatrix} -h_1 & h_0+h_1 & -h_0 & 0 & 0 & \cdots & 0\\ h_0 & 2(h_0+h_1) & h_1 & 0 & 0 & \cdots & 0\\ 0 & h_1 & 2(h_1+h_2) & h_2 & 0 & \cdots & 0\\ 0 & 0 & h_2 &2(h_2+h_3) & h_3 & \cdots & 0\\ \vdots& & & \ddots & \ddots & \ddots & \vdots\\ 0 & \cdots & & & h_{n-2} & 2(h_{n-2}+h_{n-1}) & h_{n-1}\\ 0 & \cdots & & & -h_{n-1} & h_{n-1}+h_{n-2} & -h_{n-2} \end{bmatrix}\begin{bmatrix} m_0\\m_1\\m_2\\m_3\\\vdots\\m_{n-1}\\m_n \end{bmatrix}=6\begin{bmatrix} 0\\ \frac{y_2-y_1}{h_1}-\frac{y_1-y_0}{h_0}\\ \frac{y_3-y_2}{h_2}-\frac{y_2-y_1}{h_1}\\ \frac{y_4-y_3}{h_3}-\frac{y_3-y_2}{h_2}\\ \vdots\\ \frac{y_n-y_{n-1}}{h_{n-1}}-\frac{y_{n-1}-y_{n-2}}{h_{n-2}}\\ 0 \end{bmatrix}在n=4时通⽤公式⾃然边界\begin{bmatrix} 1 & 0 & 0 & 0 \\ h_0 & 2(h_0+h_1) & h_1 & 0 \\ 0 & h_1 & 2(h_1+h_2) & h_2 \\ 0 & 0 & 0 & 1 \\ \end{bmatrix}\begin{bmatrix} m_0\\m_1\\m_2\\m_3 \end{bmatrix}=6\begin{bmatrix} 0\\ \frac{y_2-y_1}{h_1}-\frac{y_1-y_0}{h_0}\\ \frac{y_3-y_2}{h_2}-\frac{y_2-y_1}{h_1}\\ 0 \end{bmatrix}固定边界\begin{bmatrix} 2 & 1 & 0 & 0 \\ h_0 & 2(h_0+h_1) & h_1 & 0 \\ 0 & h_1 & 2(h_1+h_2) & h_2 \\ 0 & 0 & 1 & 2 \\ \end{bmatrix}\begin{bmatrix} m_0\\m_1\\m_2\\m_3 \end{bmatrix}=6\begin{bmatrix} \frac{y_1-y_0}{h_0}-A\\ \frac{y_2-y_1}{h_1}-\frac{y_1-y_0}{h_0}\\ \frac{y_3-y_2}{h_2}-\frac{y_2-y_1}{h_1}\\ B-\frac{y_3-y_2}{h_2} \end{bmatrix}⾮节点边界\begin{bmatrix} -h_1 & h_0+h_1 & -h_0 & 0 \\ h_0 & 2(h_0+h_1) & h_1 & 0 \\ 0 & h_1 & 2(h_1+h_2) & h_2 \\ 0 & -h_2 & h_1+h_2 & -h_1 \\ \end{bmatrix}\begin{bmatrix} m_0\\m_1\\m_2\\m_3 \end{bmatrix}=6\begin{bmatrix} 0\\ \frac{y_2-y_1}{h_1}-\frac{y_1-y_0}{h_0}\\ \frac{y_3-y_2}{h_2}-\frac{y_2-y_1}{h_1}\\ 0 \end{bmatrix}图像插值x_{i+1}-x_i=1 \Rightarrow h_i(x)=1在n=4时，即四个点时如下所⽰p_{11}p_{12}p_{13}p_{14}p_{21}p_{22}p_{23}p_{24}pp_{31}p_{32}p_{33}p_{34}p_{41}p_{42}p_{43}p_{44}⾃然边界（可⽤TDMA或化简计算）\begin{bmatrix} 1 & 0 & 0 & 0 \\ 1 & 4 & 1 & 0 \\ 0 & 1 & 4 & 1 \\ 0 & 0 & 0 & 1 \\ \end{bmatrix}\begin{bmatrix} m_0\\m_1\\m_2\\m_3 \end{bmatrix}=6\begin{bmatrix} 0\\ y_0+y_2-2y_1\\ y_1+y_3-2y_2\\ 0 \end{bmatrix}固定边界（只能⽤TDMA计算）\begin{bmatrix} 2 & 1 & 0 & 0 \\ 1 & 4 & 1 & 0 \\ 0 & 1 & 4 & 1 \\ 0 & 0 & 1 & 2 \\ \end{bmatrix}\begin{bmatrix} m_0\\m_1\\m_2\\m_3 \end{bmatrix}=6\begin{bmatrix} y_1-y_0-A\\ y_0+y_2-2y_1\\ y_1+y_3-2y_2\\ y_2-y_3+B \end{bmatrix}⾮节点边界（只能化简计算）\begin{bmatrix} -1 & 2 & -1 & 0 \\ 1 & 4 & 1 & 0 \\ 0 & 1 & 4 & 1 \\ 0 & -1 & 2 & -1 \\ \end{bmatrix}\begin{bmatrix} m_0\\m_1\\m_2\\m_3 \end{bmatrix}=6\begin{bmatrix} 0\\ y_0+y_2-2y_1\\ y_1+y_3-2y_2\\ 0 \end{bmatrix}python代码1class BiSplineInterpolation:2 @staticmethod3 def TDMA(a, b, c, d):4 n = len(d)56 c[0] = c[0] / b[0]7 d[0] = d[0] / b[0]89for i in range(1, n):10 coef = 1.0 / (b[i] - a[i] * c[i - 1])11 c[i] = coef * c[i]12 d[i] = coef * (d[i] - a[i] * d[i - 1])1314for i in range(n - 2, -1, -1):15 d[i] = d[i] - c[i] * d[i + 1]1617return d1819 @staticmethod20 def Simplified_Natural4(y1, y2, y3, y4):21 # 四点⾃然边界化简公式22 d1 = y1 + y3 - 2 * y223 d2 = y2 + y4 - 2 * y32425 k0 = 026 k1 = (4 * d1 - d2) * 0.427 k2 = (4 * d2 - d1) * 0.428 k3 = 02930return [k0, k1, k2, k3]3132 @staticmethod33 def Simplified_Not_A_Knot4(y1, y2, y3, y4):34 # 四点⾮节点边界化简公式35 d1 = y1 + y3 - 2 * y236 d2 = y2 + y4 - 2 * y33738 k0 = 2 * d1 - d239 k1 = d140 k2 = d241 k3 = 2 * d2 - d14243return [k0, k1, k2, k3]4445 # TDMA矩阵说明46 # a0 和 c3 没有实际意义，占位⽤47 # a0 [b0 c0 00 ] [x0] [d0]48 # [a1 b1 c1 0 ] [x1] = [d1]49 # [0 a2 b2 c2] [x2] [d2]50 # [00 a3 b3] c3 [x3] [d3]5152 def SplineInterpolationNatural4(self, x, y1, y2, y3, y4):53 # ⽤TDMA计算54 # matrix_a = [0, 1, 1, 0]55 # matrix_b = [1, 4, 4, 1]56 # matrix_c = [0, 1, 1, 0]57 # matrix_d = [0, 6 * (y1 + y3 - 2 * y2), 6 * (y2 + y4 - 2 * y3), 0]58 # matrix_x = self.TDMA(matrix_a, matrix_b, matrix_c, matrix_d)5960 # 化简计算61 matrix_x = self.Simplified_Natural4(y1, y2, y3, y4)6263 a = y264 b = y3 - y2 - matrix_x[1] / 3.0 - matrix_x[2] / 6.065 c = matrix_x[1] / 2.066 d = (matrix_x[2] - matrix_x[1]) / 6.06768 s = a + b * x + c * x * x + d * x * x * x69return s7071 def SplineInterpolationClamped4(self, x, y1, y2, y3, y4):72 # 仅有TDMA计算，⽆法化简73 A, B = 1, 17475 matrix_a = [0, 1, 1, 1]76 matrix_b = [2, 4, 4, 2]77 matrix_c = [1, 1, 1, 0]78 matrix_d = [6 * (y2 - y1 - A), 6 * (y1 + y3 - 2 * y2), 6 * (y2 + y4 - 2 * y3), 6 * (B - y4 + y3)]79 matrix_x = self.TDMA(matrix_a, matrix_b, matrix_c, matrix_d)8081 a = y282 b = y3 - y2 - matrix_x[1] / 3.0 - matrix_x[2] / 6.083 c = matrix_x[1] / 2.084 d = (matrix_x[2] - matrix_x[1]) / 6.08586 s = a + b * x + c * x * x + d * x * x * x87return s8889 def SplineInterpolationNotAKnot4(self, x, y1, y2, y3, y4):90 # ⽆法使⽤TDMA计算91 matrix_x = self.Simplified_Not_A_Knot4(y1, y2, y3, y4)9293 a = y294 b = y3 - y2 - matrix_x[1] / 3.0 - matrix_x[2] / 6.095 c = matrix_x[1] / 2.096 d = (matrix_x[2] - matrix_x[1]) / 6.09798 s = a + b * x + c * x * x + d * x * x * x99return s100101 def biSpline4(self, srcImg, dstH, dstW):102 dstH, dstW = int(dstH), int(dstW)103 srcH, srcW, _ = srcImg.shape104 srcImg = np.pad(srcImg, ((1, 2), (1, 2), (0, 0)), 'edge')105 dstImg = np.zeros((dstH, dstW, 3), dtype=np.uint8)106for dstY in range(dstH):107for dstX in range(dstW):108for channel in [0, 1, 2]:109 # p11 p12 p13 p14110 #111 # p21 p22 p23 p24112 # p113 # p31 p32 p33 p34114 #115 # p41 p42 p43 p44116 # 储存为 p(y, x)117 p = [dstY * srcH / dstH, dstX * srcW / dstW]118 p22 = [math.floor(p[0]), math.floor(p[1])]119 p21 = [p22[0], p22[1] - 1]120 p23 = [p22[0], p22[1] + 1]121 p24 = [p22[0], p22[1] + 2]122123 p11 = [p21[0] - 1, p21[1]]124 p12 = [p11[0], p22[1]]125 p13 = [p11[0], p23[1]]126 p14 = [p11[0], p24[1]]127128 p31 = [p21[0] + 1, p21[1]]129 p32 = [p31[0], p22[1]]130 p33 = [p31[0], p23[1]]131 p34 = [p31[0], p24[1]]132133 p41 = [p21[0] + 2, p21[1]]134 p42 = [p41[0], p22[1]]135 p43 = [p41[0], p23[1]]136 p44 = [p41[0], p24[1]]137138 diff_y, diff_x = p[0] - p22[0], p[1] - p22[1]139 r1 = self.SplineInterpolationNatural4(diff_x, srcImg[p11[0], p11[1], channel], srcImg[p12[0], p12[1], channel], srcImg[p13[0], p13[1], channel], srcImg[p14[0], p14[1], channel]) 140 r2 = self.SplineInterpolationNatural4(diff_x, srcImg[p21[0], p21[1], channel], srcImg[p22[0], p22[1], channel], srcImg[p23[0], p23[1], channel], srcImg[p24[0], p24[1], channel]) 141 r3 = self.SplineInterpolationNatural4(diff_x, srcImg[p31[0], p31[1], channel], srcImg[p32[0], p32[1], channel], srcImg[p33[0], p33[1], channel], srcImg[p34[0], p34[1], channel]) 142 r4 = self.SplineInterpolationNatural4(diff_x, srcImg[p41[0], p41[1], channel], srcImg[p42[0], p42[1], channel], srcImg[p43[0], p43[1], channel], srcImg[p44[0], p44[1], channel]) 143144 c = self.SplineInterpolationNatural4(diff_y, r1, r2, r3, r4)145146 dstImg[dstY, dstX, channel] = np.clip(c, 0, 255)。
三次样条插值
三次样条插值
三次样条插值是插值运算的一种，它具有计算精度高、收敛性好以及曲线拟合准确等特点，是插值运算中最常用的插值方法之
三次样条插值是以曲线为基本元素，把离散点数据连接成一个曲线，并能够在曲线上求出任意点的函数值。
它通过拟合所有离散数据点，来求出一个连续曲线，从而解决了插值法的局限性。
三次样条插值的基本原理是：在离散点的两端，曲线的曲率是零，由此可以计算出曲线的系数，从而得到曲线的表达式，这样就可以得到曲线上任意点的函数值。
三次样条插值的优点在于计算精度高、收敛性好，可以很好地拟合离散数据，并且经过插值后得到的曲线更加平滑，其结果更加可靠。
由于它的优点，三次样条插值得到了广泛的应用，如在统计分析中，用于拟合离散数据；在机械工程中，用于优化加工轨迹；在号处理中，用于滤波等。
总之，三次样条插值是插值运算的一种，它的准确性高，拟合性好，广泛应用于各种领域，是科学研究中的一种重要方法。
插值算法系列的优缺点
插值算法系列（1）三种基本插值算法及优缺点1最近邻算法基本思想：目的像素的值，通过坐标反向变换得到的浮点坐标为(i+u,j+v) (其中i、j均为浮点坐标的整数部分，u、v为浮点坐标的小数部分，用源图像中里该浮点坐标最近的像素值表示，即用最邻近的点的像素值赋予目的像素值。
特点：速度最快，缩放质量差，放大图像不连续性明显，缩小时图像失真。
2 双线性插值：双线性内插值算法描述如下:对于一个目的像素，设置坐标通过反向变换得到的浮点坐标为(i+u,j+v) (其中i、j均为浮点坐标的整数部分，u、v为浮点坐标的小数部分，是取值[0,1)区间的浮点数)，则这个像素得值f(i+u,j+v) 可由原图像中坐标为(i,j)、(i+1,j)、(i,j+1)、(i+1,j+1)所对应的周围四个像素的值决定，即：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)其中f(i,j)表示源图像(i,j)处的的像素值，以此类推。
特点：速度慢，质量好，图像平滑，但双线性插值具有低通滤波器的性质，使高频分量受损，所以会使图像轮廓变得模糊。
3 三次立方卷积插值：算法：他考虑一个浮点坐标(i+u,j+v)周围的16个邻点，目的像素值f(i+u,j+v)可由如下插值公式得到： <?xml:namespace prefix = o ns = "urn:schemas-microsoft-com:office:office" />f(i+u,j+v) = [A] * [B] * [C][A]=[ S(u + 1) S(u + 0) S(u - 1) S(u - 2) ]┏ f(i-1, j-1) f(i-1, j+0) f(i-1, j+1) f(i-1, j+2) ┓[B]=┃ f(i+0, j-1) f(i+0, j+0) f(i+0, j+1) f(i+0, j+2) ┃┃ f(i+1, j-1) f(i+1, j+0) f(i+1, j+1) f(i+1, j+2) ┃┗ f(i+2, j-1) f(i+2, j+0) f(i+2, j+1) f(i+2, j+2) ┛┏ S(v + 1) ┓[C]=┃ S(v + 0) ┃┃ S(v - 1) ┃┗ S(v - 2) ┛┏ 1-2*Abs(x)^2+Abs(x)^3 , 0 <=Abs(x) <1S(x)=｛ 4-8*Abs(x)+5*Abs(x)^2-Abs(x)^3 , 1 <=Abs(x) <2┗ 0 , Abs(x)>=2S(x)是对 Sin(x*Pi)/xPi 的逼近（Pi是圆周率——π）即对源图像进行插值扩大。
ArcGIS中栅格数据重采样方法介绍
ArcGIS中栅格数据重采样方法处理栅格数据时，由于数据像元大小不符合要求，或者在进行栅格数据配准后，像元发生倾斜，或者对多个栅格数据进行分析时，需要使用相同的栅格分辨率，重采样到同一分辨率下，因此对栅格数据操作时经常要进行重采样操作。
如下图所示，通过重采样将栅格调整到新的分辨率下。
1.栅格重采样方法栅格重采样主要包括三种方法：最邻近法、双线性内插法和三次卷积插值法。
最邻近法是把原始图像中距离最近的像元值填充到新图像中；双线性内插法和三次卷积插值法都是把原始图像附近的像元值通过距离加权平均填充到新图像中。
默认情况下，采用最近邻分配重采样技术，这种方法同时适用于离散和连续值类型，而其他重采样方法只适用于连续数据。
1.1 最近邻法最邻近分配法是用于离散（分类）数据的重采样技术，因为它不会更改输入单元的值。
将输出栅格数据集中单元中心的位置定位到输入栅格后，最邻近分配法将确定输入栅格上最近的单元中心位置并将该单元的值分配给输出栅格上的单元。
最邻近分配法不会更改输入栅格数据集中单元的任何值。
输入栅格中的值2 在输出栅格中仍将为2，决不会为2.2 或2.3。
由于输出单元值保持不变，因此最邻近分配法应该用于名目数据或顺序数据，其中每个值都表示一个类、一个成员或一个分类（分类数据，如土地利用、土壤或森林类型）。
考虑到根据输入栅格创建的输出栅格会在操作中旋转45°，因此将进行重采样。
对于每个输出单元，都要从输入栅格中获取值。
在下图中，输入栅格的单元中心为灰色点。
输出单元为绿色阴影。
要处理的单元为黄色阴影。
在最邻近分配法中，将确定与要处理的单元中心（红色点）最邻近的输入栅格单元中心（橙色点)，并将其指定为要处理的单元（黄色阴影）的输出值。
对输出栅格中的每个单元都重复此过程。
1.2 双线性插值法双线性插值法使用四个最邻近输入单元中心的值来确定输出栅格上的值。
输出单元的新值是这四个值的加权平均值，将根据它们与输出单元中心的距离进行调整。
三次立方卷积插值法
三次立方卷积插值法（Cubic Convolution Interpolation）是一种常用的数字图像插值方法。
这种方法使用周围的16个像素点来计算新的像素值，从而得到更高分辨率的图像。
以下是三次立方卷积插值法的基本步骤：
确定需要插值的点(x, y)以及周围的16个像素点。
对这16个点应用三次立方卷积核。
这个卷积核是一个分段函数，形式如下：
W(x) = { (a+2)|x|^3 - (a+3)|x|^2 + 1, |x| <= 1
{ a|x|^3 - 5a|x|^2 + 8a|x| - 4a, 1 < |x| < 2
{ 0, |x| >= 2
其中a通常取-0.5或-1。
这个函数对于|x|<1的部分是一个三次方程，对于1<|x|<2的部分是另一个三次方程，对于|x|>=2的部分是0。
3. 将卷积核应用于每个点，并将结果相加，得到新的像素值。
注意：在实际应用中，我们通常会先将图像进行预处理（如去噪、平滑等），然后再进行插值。
此外，插值后的图像可能需要进行后处理（如锐化、对比度增强等），以进一步提高图像质量。
以上步骤可以通过编程语言（如Python）和图像处理库（如OpenCV）来实现。
具体实现时，需要注意边界处理（如镜像、填充等）以及计算效率等问题。
三次卷积内插法计算例题
三次卷积内插法计算例题三次卷积内插法是一种用于数据插值的方法，它可以通过已知的数据点来估计在两个数据点之间的未知数据点的值。
下面我将以一个简单的例题来说明三次卷积内插法的计算过程。
假设我们有以下已知数据点：点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。
如果不满足这个条件，可以通过对数据进行插值或者进行其他的数据处理方法来满足这个条件。
这就是使用三次卷积内插法计算例题的详细过程。
希望这个回答能够满足你的需求。
如果你还有其他问题，欢迎继续提问。
gdal重采样算法
gdal重采样算法GDAL（Geospatial Data Abstraction Library）是一个用于处理地理空间数据的开源库，它提供了多种重采样算法，用于在数据重采样或变换过程中调整像素值。
以下是一些GDAL中常用的重采样算法：1. 最近邻插值（Nearest Neighbour Interpolation）：这是最简单的重采样方法。
对于输出图像的每个像素，该方法选择输入图像中最接近的像素值作为输出值。
这种方法速度快，但可能会产生锯齿状边缘。
2. 双线性插值（Bilinear Interpolation）：这种方法考虑了输出像素周围四个最近的输入像素。
它通过对这四个像素进行线性插值来计算输出像素的值。
这种方法在大多数情况下都能给出较好的结果，但可能不如更高级的方法那么平滑。
3. 三次卷积插值（Cubic Convolution Interpolation）：这种方法使用了一个16个最近输入像素的窗口来计算输出像素的值。
它通过对这些像素进行三次卷积插值来得到输出值。
这种方法通常能给出非常平滑的结果，但计算量较大，速度较慢。
4. 三次卷积重采样（Cubic Spline Interpolation）：这是另一种基于三次插值的方法，但它使用了不同的算法来计算输出像素的值。
这种方法通常能给出比简单三次卷积插值更好的结果，但计算量也更大。
5. 模式插值（Mode Interpolation）：这种方法计算输出像素周围输入像素的众数（最常见的值）作为输出值。
这种方法主要用于分类数据的重采样。
在GDAL中，你可以通过`gdalwarp`命令或GDAL的Python接口来指定重采样算法。
例如，使用`gdalwarp`命令时，可以使用`-r`选项来指定重采样算法，如`-r bilinear`表示使用双线性插值。
在Python中，你可以通过设置`gdal.WarpOptions`对象的`resampleAlg`属性来选择重采样算法。
双线性插值和三次卷积在图像缩放中的应用及实现
双线性插值和三次卷积在图像缩放中的应⽤及实现
双线性插值和三次卷积在图像缩放中的应⽤及实现
◆徐鹏飞
【摘要】摘要：图像缩⼩是从⼤数据量到⼩数据量的处理过程，因此只需要从多个数据中，以适当的⽅式选出所需要的信息就可以完成。
图像放⼤则是⼩数据量到⼤数据量的处理过程，因此需要对许多未知的数据进⾏估计。
好在图像的相邻像素之间的相关性很强，可以利⽤这个相关性来实现图像的放⼤。
插值算法对于缩放⽐例较⼩的情况是完全可以接受的，令⼈信服的。
⼀般缩⼩0．5倍以上或扩⼤3．0倍以下，对任何图像都是可以接受的。
本⽂采⽤双线性插值和三次卷积算法来实现。
【期刊名称】⽹络安全技术与应⽤
【年(卷),期】2017(000)012
【总页数】2
【关键词】图像缩放；双线性插值；三次卷积
1 图像的缩放
图像的缩⼩从物理意义上来说，是将描述图像的每个像素的物理尺⼨缩⼩相应的倍数就可以完成。
但如果像素的物理尺⼨不允许改变，从数码技术的⾓度来看，图像的缩⼩实际上就是通过减少像素个数来实现的。
既然图像的缩⼩是通过减少像素个数来实现，那么，就需要根据所期望缩⼩的尺⼨数据，从原图像中选择合适的像素点，使图像缩⼩之后可以尽量保持原有图像的概貌特征不丢失。
图像的缩⼩分为按⽐例缩⼩和不按⽐例缩⼩两种。
按⽐例缩⼩就是图像的长和宽按照同样的⽐例缩⼩。
不按⽐例缩⼩则是指图像缩⼩时，长和宽的缩⼩⽐例不同。
图像的放⼤，从物理含义上讲是指图像缩⼩的逆操作。
但是从信息。
三种图像重采样方法的特点和区别
图像重采样主要有三种方法，分别是最邻近法，双线性内插法和三次卷积内插法。
⑴最近邻法。
该法针对于二维图像取待采样点周围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 〕为三次多项式的值。
ArcGIS教程：栅格像元大小和重采样
ArcGIS教程：栅格像元大小和重采样不同的栅格数据集不需要使用相同的像元分辨率进行存储。
但在处理多个数据集时，最好使用相同的像元分辨率。
将多个不同分辨率的栅格数据集输入任一A rcGIS Spatia l Analyst 扩展模块工具时，会自动将一个或多个输入数据集重采样为输入数据集的最粗糙分辨率。
在默认情况下，使用最邻近分配重采样技术。
这是因为它同时适用于离散数据和连续数据，而其他重采样类型(双线性插值和三次卷积插值)只适用于连续数据。
使用重采样技术是很有必要的，因为输入像元的中心很少能与变换到所需分辨率后的像元中心对齐。
在合并不同分辨率的栅格之前，可以使用重采样工具利用双线性和三次卷积技术作为预处理步骤。
可以使用像元大小环境参数来控制默认重采样分辨率，在该参数中可指定工具是使用输入栅格的最小分辨率还是使用定义的特定像元大小。
如下图所示，在分析环境中设置的像元大小比工具中的输入栅格的像元大小要粗糙。
执行时，会首先将输入栅格重采样到较粗糙的分辨率，然后应用该工具。
执行分析时，请确定所设置的像元大小是否合适。
例如，当像元大小为5 千米时，不太可能研究鼠标移动。
五千米的像元可能更适用于研究全球变暖对地球的影响。
重采样要找到重采样输出栅格上每个像元应采用的值，必须将输出栅格中每个像元的中心映射到原始输入坐标系。
每个像元中心坐标都会反向变换，以确定该点在原始输入栅格上的位置。
确定输入位置后，将根据输入栅格中邻近的像元为输出位置分配一个值。
输出像元中心很少能与输入栅格像元中心准确对齐。
因此，多种技术已开发出来，以根据点相对于输入栅格像元中心的位置以及与这些像元相关联的值来确定输出值。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
最临近插值的的思想很简单。对于通过反向变换得到的的一个浮点坐标，对其进行简单的取整，得到一个整数型坐标，这个整数型坐标对应的像素值就是目的像素的像素值，也就是说，取浮点坐标最邻近的左上角点（对于DIB是右上角，因为它的扫描行是逆序存储的）对应的像素值。可见，最邻近插值简单且直观，但得到的图像质量不高
┗ f(i+2, j-1) f(i+2, j+0) f(i+2, j+1) f(i+2, j+2) ┛
┏ S(v + 1) ┓
[C]=┃ S(v + 0) ┃
┃ S(v - 1) ┃
┗ S(v - 2) ┛
┏ 1-2*Abs(x)^2+Abs(x)^3 , 0<=Abs(x)<1
S(x)=｛ 4-8*Abs(x)+5*Abs(x)^2-Abs(x)^3 , 1<=Abs(x)<2
┗ 0 , Abs(x)>=2
S(x)是对 Sin(x*Pi)/x 的逼近（Pi是圆周率——π）
�
双线性内插值：
对于一个目的像素，设置坐标通过反向变换得到的浮点坐标为(i+u,j+v)，其中i、j均为非负整数，u、v为[0,1)区间的浮点数，则这个像素得值 f(i+u,j+v) 可由原图像中坐标为 (i,j)、(i+1,j)、(i,j+1)、(i+1,j+1)所对应的周围四个像素的值决定，即：
三次卷积法能够克服以上两种算法的不足，计算精度高，但计算亮大，他考虑一个浮点坐标(i+u,j+v)周围的16个邻点，目的像素值f(i+u,j+v)可由如下插值公式得到：
f(i+u,j+v) = [A] * [B] * [C]
[A]=[ S(u + 1) S(u + 0) S(u - 1) S(u - 2) ]
┏ f(i-1, j-1) f(i-1, j+0) f(i-1, j+1) f(i-1, j+2) ┓
[B]=┃ f(i+0, j-1) f(i+0, j+0) f(i+0, j+1) f(i+0, j+2) ┃
┃ f(i+1, j-1) f(i+1, j+0) f(i+1, j+1) f(i+1, j+2) ┃
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)
其中f(i,j)表示源图像(i,j)处的的像素值，以此类推
这就是双线性内插值法。双线性内插值法计算量大，但缩放后低通滤波器的性质，使高频分量受损，所以可能会使图像轮廓在一定程度上变得模糊