基本原理：双三次插值是一种更加复杂的插值方式，它能创造出比双线性插值更平滑的图像边缘。缩放后图像中某个象素的象素值是由源图像相应像素附近的（4 x 4）个邻近象素值计算出来的，即通过一个基函数进行拟合得到一个目的像素值，具体某点v(x，y) 的像素值是使用下式计算得到：

v(x，y) =∑∑aij*x^i*y^j；其中，0≤i，j≤3；16个系数aij由16个临近像素写出的未知方程确定。

C/C++实现如下：

voidGeometryTrans::Zoom(floatratioX,floatratioY){//释放旧的输出图像缓冲区if(m_pImgDataOut!=NULL){delete[]m_pImgDataOut;m_pImgDataOut=NULL;}//输出图像的宽和高m_imgWidthOut=int(m_imgWidth*ratioX+0.5);m_imgHeightOut=int(m_imgHeight*ratioY+0.5);//输入图像每行像素字节数intlineByteIn=(m_imgWidth*m_nBitCount/8+3)/4*4;//输出图像每行像素字节数intlineByteOut=(m_imgWidthOut*m_nBitCount/8+3)/4*4;//申请缓冲区，存放输出结果m_pImgDataOut=newunsignedchar[lineByteOut*m_imgHeightOut];//每像素字节数，输入图像与输出图像相同intpixelByte=m_nBitCount/8;//输出图像在输入图像中待插值的位置坐标intcoordinateX,coordinateY;//循环变量，输出图像的坐标inti,j;//循环变量,像素的每个通道intk;//**************************************************************//对原图像进行拓展，上下分别拓展两行，左右分别拓展两列if(m_pImgData_temp!=NULL){delete[]m_pImgData_temp;m_pImgData_temp=NULL;}//拓展图像每行像素字节数intlineByteIn1=((m_imgWidth+4)*m_nBitCount/8+3)/4*4;//申请缓冲区，存放拓展后的图像m_pImgData_temp=newunsignedchar[lineByteIn1*(m_imgHeight+4)];for(k=0;k<pixelByte;k++){for(i=0;i<m_imgHeight;i++)//中间m_imgWidth*m_imgHeight的部分for(j=0;j<m_imgWidth;j++)*(m_pImgData_temp+(i+2)*lineByteIn1+(j+2)*pixelByte+k)=*(m_pImgData+i*lineByteIn+j*pixelByte+k);for(i=0;i<2;i++)//拓展上面两行for(j=0;j<m_imgWidth;j++)*(m_pImgData_temp+i*lineByteIn1+(j+2)*pixelByte+k)=*(m_pImgData+i*lineByteIn+j*pixelByte+k);for(i=0;i<2;i++)//拓展下面两行for(j=0;j<m_imgWidth;j++)*(m_pImgData_temp+(i+m_imgHeight+2)*lineByteIn1+(j+2)*pixelByte+k)=*(m_pImgData+(i+m_imgHeight-2)*lineByteIn+j*pixelByte+k);for(i=0;i<m_imgHeight;i++)//拓展左边两列for(j=0;j<2;j++)*(m_pImgData_temp+(i+2)*lineByteIn1+j*pixelByte+k)=*(m_pImgData+i*lineByteIn+j*pixelByte+k);for(i=0;i<m_imgHeight;i++)//拓展右边两列for(j=0;j<2;j++)*(m_pImgData_temp+(i+2)*lineByteIn1+(j+m_imgWidth+2)*pixelByte+k)=*(m_pImgData+i*lineByteIn+(j+m_imgWidth-2)*pixelByte+k);for(i=0;i<2;i++)//左上角部分for(j=0;j<2;j++)*(m_pImgData_temp+i*lineByteIn1+j*pixelByte+k)=*(m_pImgData+i*lineByteIn+j*pixelByte+k);for(i=0;i<2;i++)//右上角部分for(j=0;j<2;j++)*(m_pImgData_temp+i*lineByteIn1+(j+m_imgWidth+2)*pixelByte+k)=*(m_pImgData+i*lineByteIn+(j+m_imgWidth-2)*pixelByte+k);for(i=0;i<2;i++)//左下角部分for(j=0;j<2;j++)*(m_pImgData_temp+(i+m_imgHeight+2)*lineByteIn1+j*pixelByte+k)=*(m_pImgData+(i+m_imgHeight-2)*lineByteIn+j*pixelByte+k);for(i=0;i<2;i++)//右下角部分for(j=0;j<2;j++)*(m_pImgData_temp+(i+m_imgHeight+2)*lineByteIn1+(j+m_imgWidth+2)*pixelByte+k)=*(m_pImgData+(i+m_imgHeight-2)*lineByteIn+(j+m_imgWidth-2)*pixelByte+k);}//***************************************************************floatu,v;intii,jj;doubleArray[3][16];//存储临时像素值的二维数组floata[4],c[4];//系数数组//双3次多项式插值for(i=0;i<m_imgHeightOut;i++){coordinateY=i/ratioY+0.5;u=(float)(fmodf(i,ratioY)*(1/ratioY));/*计算系数矩阵*/a[0]=Sinc(1+u);a[1]=Sinc(u);a[2]=Sinc(1-u);a[3]=Sinc(2-u);for(j=0;j<m_imgWidthOut;j++){coordinateX=j/ratioX+0.5;/*计算系数矩阵*/v=(float)(fmodf(j,ratioX)*(1/ratioX));c[0]=Sinc(1+v);c[1]=Sinc(v);c[2]=Sinc(1-v);c[3]=Sinc(2-v);//输出图像坐标为(j,i)的像素映射到原图中的坐标值，即插值位置//若插值位置在输入图像范围内，则双三次多项式插值if(0<=coordinateX&&coordinateX<(m_imgWidth)&&coordinateY>=0&&coordinateY<(m_imgHeight)){for(k=0;k<pixelByte;k++)//三个通道，RGBfor(ii=0;ii<4;ii++)for(jj=0;jj<4;jj++)/*与目的像素相关的16个原始像素RGB*/{Array[k][ii*4+jj]=(double)(*(m_pImgData_temp+(coordinateY+ii)*lineByteIn1*+(coordinateX+jj)*pixelByte+k));}for(k=0;k<pixelByte;k++)*(m_pImgDataOut+i*lineByteOut+j*pixelByte+k)=(unsignedchar)ABC(a,Array[k],c);/*调用ABC求得目的像素*/}else//若不在输入图像范围内，则置255{for(k=0;k<pixelByte;k++)*(m_pImgDataOut+i*lineByteOut+j*pixelByte+k)=255;}}}}

基函数实现如下

loatGeometryTrans::Sinc(floatx)/*Sinc(x)是对Sin(x*Pi)/x的逼近（Pi是圆周率——π）*/{if(abs(x)>=0&&abs(x)<1)return(1-2*abs(x)*abs(x)+abs(x)*abs(x)*abs(x));elseif(abs(x)>=1&&abs(x)<2)return(4-8*abs(x)+5*abs(x)*abs(x)-abs(x)*abs(x)*abs(x));elsereturn0;}

矩阵内积计算函数实现如下

floatGeometryTrans::ABC(floata[],doubleb[],floatc[])/*矩阵运算函数，求得像素值，内积*/{inti,j;floatabc=0;floattmp[4];for(i=0;i<4;i++)tmp[i]=0;for(i=0;i<4;i++)for(j=0;j<4;j++)tmp[i]+=a[j]*b[j*4+i];for(i=0;i<4;i++)abc+=tmp[i]*c[i];if(abc<0)abc=0;if(abc>255)abc=255;returnabc;}