數字圖像縮放之雙三次插值

基本原理：雙三次插值是一種更加復雜的插值方式，它能創造出比雙線性插值更平滑的圖像邊緣。縮放后圖像中某個象素的象素值是由源圖像相應像素附近的（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;  
}