数字图像缩放之双三次插值

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

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

C/C++实现如下：

void GeometryTrans::Zoom(float ratioX, float ratioY)
{
     
    //释放旧的输出图像缓冲区
    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); 
 
    //输入图像每行像素字节数
    int lineByteIn=(m_imgWidth*m_nBitCount/8+3)/4*4;
     
    //输出图像每行像素字节数
    int lineByteOut=(m_imgWidthOut*m_nBitCount/8+3)/4*4;
 
    //申请缓冲区，存放输出结果
     
    m_pImgDataOut=new unsigned char[lineByteOut*m_imgHeightOut];
 
    //每像素字节数，输入图像与输出图像相同
    int pixelByte=m_nBitCount/8;
     
 
    //输出图像在输入图像中待插值的位置坐标
    int coordinateX, coordinateY;
     
    //循环变量，输出图像的坐标
    int i,j;
     
    //循环变量,像素的每个通道
    int k;
    //**************************************************************
    //对原图像进行拓展，上下分别拓展两行，左右分别拓展两列
    if(m_pImgData_temp!=NULL){
        delete []m_pImgData_temp;
        m_pImgData_temp=NULL;
    }
    //拓展图像每行像素字节数
    int lineByteIn1=((m_imgWidth+4)*m_nBitCount/8+3)/4*4;
    //申请缓冲区，存放拓展后的图像
     
    m_pImgData_temp=new unsigned char[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) ;
    }    
    //***************************************************************
    float u,v;
    int ii,jj;
    double Array[3][16];//存储临时像素值的二维数组
 
    float a[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++)//三个通道，RGB
                    for(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) 
                    =(unsigned char)ABC(a,Array[k],c);/*调用ABC求得目的像素*/
                     
            }
            else //若不在输入图像范围内，则置255  
            {
                for(k=0;k<pixelByte;k++)
                    *(m_pImgDataOut + i * lineByteOut + j*pixelByte+k) = 255;
            }          
        }
    }
}

基函数实现如下

loat GeometryTrans::Sinc(float x) /*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));
  else
     if(abs(x)>=1&&abs(x)<2)
      return (4-8*abs(x)+5*abs(x)*abs(x)-abs(x)*abs(x)*abs(x));
     else  
       return 0;
}

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

float GeometryTrans::ABC(float a[],double b[],float c[]) /*矩阵运算函数，求得像素值，内积*/
{
  int i,j;
  float abc=0;
  float tmp[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;
  return abc;  
}