image.c

该程序包里包括了很多跟小波变换有关的经典算法

	//分配浮点型的数据内存空间，这个是函数的输出
	pfimage=(PFIMAGE)malloc(sizeof(FIMAGE));
	if (pfimage==NULL){
		free(pimage);
		return NULL;
	}
	//打开指定路径的pgm图像
	if ((fp=fopen(PGMFileName, "rb"))==NULL){
		free(pimage);
		return NULL;
	}
	//处理文件头
	while ((ch != 'P') && (ch != '#')){ 
		ch = fgetc(fp);
	}
	
	//跳过文件附件信息
	PGMSkipComments(fp, &ch);

	//规定支持打开的pgm类型为P5
	ftype = fgetc(fp);
	//异常处理
	if (ftype != '5'){
		ImageWarning("Only binary PGM is supported.\n");
		free(pfimage);
		free(pimage);
		fclose(fp);
		return NULL;
	}
	//获取图像的宽、高数值
	pimage->xsize=(int)PGMGetVal(fp);
	pimage->ysize=(int)PGMGetVal(fp);
	//获取图像灰度值的最大值
	MaxG=(int)PGMGetVal(fp);
	
	//分配图像数据的内存空间
	if(!ImageBufferAlloc(pimage)){
		free(pimage);
		free(pfimage);
		return NULL;
	}
	
	//从文件尾开始的指针偏移，以避免读取无用的图像附加信息
	fseek(fp, -1*pimage->xsize * pimage->ysize, SEEK_END);

	//加载图像数据，如果发现数据的大小不匹配，则出错退出，释放原有的图像数据空间，关闭已打开的pgm文件
	if (fread(pimage->pixel[0], sizeof(unsigned char),
		pimage->xsize*pimage->ysize, fp) != (unsigned int)pimage->xsize*pimage->ysize){
		fclose(fp);
		ImageBufferFree(pimage);
		free(pimage);
		free(pfimage);
		return NULL;
	}
	fclose(fp);
	//将数据类型转换为浮点型
	pfimage->xsize=pimage->xsize;
	pfimage->ysize=pimage->ysize;	

	//分配内存空间
	if(!FImageBufferAlloc(pfimage)){
		free(pfimage);
		ImageFree(pimage);
		return NULL;
	}
	//将pimage中的数据转换后拷贝到pfimage中，完成浮点转换
	for (i=0; i< pfimage->xsize*pfimage->ysize; i++){
		pfimage->pixelLinear[i]=(double)pimage->pixelLinear[i];
	}
	//释放暂存用的整形数据空间
	ImageFree(pimage);

	return pfimage;
}


/*----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/	
int WritePGM(PIMAGE pimage, char *PGMFileName)
{
	FILE *fp;
	int i;
	
	if ((fp=fopen(PGMFileName, "wb+"))==NULL){
		return 0;
	}
	
	fprintf(fp, "P5\n#%s\n%d %d\n255\n", PGMFileName, pimage->xsize, pimage->ysize);

	i=fwrite(pimage->pixel[0], sizeof(unsigned char),
		(pimage->xsize)*(pimage->ysize), fp);
	
	fclose(fp);
	
	if (i!=(pimage->ysize)*(pimage->xsize)){
      return 0;
   }
   else{
      return i;
   }
}

/*----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/	
int WriteFloatToPGM(PFIMAGE pfimage, char *PGMFileName)
{
	FILE *fp;
	unsigned char **temp, *tempLinear;
	int i;
	
	temp=Byte_Alloc(pfimage->ysize, pfimage->xsize);
	tempLinear=temp[0];

	if (temp==NULL){
		return 0;
	}

	if ((fp=fopen(PGMFileName, "wb+"))==NULL){
		return 0;
	}
	
	fprintf(fp, "P5\n#%s\n%d %d\n255\n", PGMFileName, pfimage->xsize, pfimage->ysize);

	/* Copy */
	for (i=0; i<pfimage->xsize*pfimage->ysize; i++){
		tempLinear[i]=(int)((pfimage->pixelLinear[i]+0.5) < 0.0 ? 0 : 
						 ((pfimage->pixelLinear[i]+0.5) > 255.0 ? 255 : 
						  (pfimage->pixelLinear[i]+0.5)));
	}

	i=fwrite(temp[0], sizeof(unsigned char), (pfimage->xsize)*(pfimage->ysize), fp);
	fclose(fp);
	
	free(temp[0]);
	free(temp);

	if (i!=(pfimage->ysize)*(pfimage->xsize)){
      return 0;
   }
   else{
      return i;
   }
}

/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
/* Get a number from a pgm file header, skipping comments etc. */
unsigned int PGMGetVal (FILE* infile)
{
  unsigned int tmp;
  unsigned char ch;
  do { 
	  ch = fgetc(infile); 
  } while ((ch <= ' ') && (ch != '#'));

  PGMSkipComments(infile, &ch);

  ungetc(ch, infile);

  /* Do not understand this part */
  if (fscanf(infile,"%u",&tmp) != 1) {
    ImageError("%s\n","Error parsing file!");
  }

  return(tmp);
}

/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
void PGMSkipComments(FILE* infile, unsigned char *ch)
{
	while ((*ch == '#')){
		while (*ch != '\n'){ 
			*ch = fgetc(infile); 
		}
		while (*ch < ' '){
			*ch = fgetc(infile); 
		}
	}
}

/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
double ImageCompareMSE(PIMAGE pimage1, PIMAGE pimage2)
{
   int i;
   double Error, SumOfSquareError=0.0;

   if ((pimage1->xsize != pimage2->xsize) || (pimage1->ysize!=pimage2->ysize)){
      ImageWarning("Cannot compare images of different sizes (%dx%d and %dx%d)\n",
         pimage1->xsize, pimage1->ysize, pimage2->xsize, pimage2->ysize);
      return 0.0;
   }

   for (i=0; i< pimage1->xsize*pimage1->ysize; i++){
      Error=pimage1->pixelLinear[i]-pimage2->pixelLinear[i];
      SumOfSquareError+=Error*Error;
   }

   return SumOfSquareError/(double)(pimage1->xsize*pimage1->ysize);

}

/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
double ImageComparePSNR(PIMAGE pimage1, PIMAGE pimage2)
{
   int i;
   double Error, SumOfSquareError=0.0, RootMeanSquareError;

   if ((pimage1->xsize != pimage2->xsize) || (pimage1->ysize!=pimage2->ysize)){
      ImageWarning("Cannot compare images of different sizes (%dx%d and %dx%d)\n",
         pimage1->xsize, pimage1->ysize, pimage2->xsize, pimage2->ysize);
      return 0.0;
   }

   for (i=0; i< pimage1->xsize*pimage1->ysize; i++){
      Error=pimage1->pixelLinear[i]-pimage2->pixelLinear[i];
      SumOfSquareError+=Error*Error;
   }

   RootMeanSquareError = sqrt(SumOfSquareError/(double)(pimage1->xsize*pimage1->ysize));
   
   return 20*log10(255.0/RootMeanSquareError); /* assume maximum 255 */

}


/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
double FImageCompareMSE(PFIMAGE pfimage1, PFIMAGE pfimage2)
{
   int i;
   double Error, SumOfSquareError=0.0;

   if ((pfimage1->xsize != pfimage2->xsize) || (pfimage1->ysize!=pfimage2->ysize)){
      ImageWarning("Cannot compare images of different sizes (%dx%d and %dx%d)\n",
         pfimage1->xsize, pfimage1->ysize, pfimage2->xsize, pfimage2->ysize);
      return 0.0;
   }

   for (i=0; i< pfimage1->xsize*pfimage1->ysize; i++){
      Error=pfimage1->pixelLinear[i]-pfimage2->pixelLinear[i];
      SumOfSquareError+=Error*Error;
   }

   return SumOfSquareError/(double)(pfimage1->xsize*pfimage1->ysize);

}

/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
double FImageComparePSNR(PFIMAGE pfimage1, PFIMAGE pfimage2)
{
   int i;
   double Error, SumOfSquareError=0.0, RootMeanSquareError;

   if ((pfimage1->xsize != pfimage2->xsize) || (pfimage1->ysize!=pfimage2->ysize)){
      ImageWarning("Cannot compare images of different sizes (%dx%d and %dx%d)\n",
         pfimage1->xsize, pfimage1->ysize, pfimage2->xsize, pfimage2->ysize);
      return 0.0;
   }

   for (i=0; i< pfimage1->xsize*pfimage1->ysize; i++){
      Error=pfimage1->pixelLinear[i]-pfimage2->pixelLinear[i];
      SumOfSquareError+=Error*Error;
   }

   RootMeanSquareError = sqrt(SumOfSquareError/(double)(pfimage1->xsize*pfimage1->ysize));
   
   return 20*log10(255.0/RootMeanSquareError); /* assume maximum 255 */

}

/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
/* RGB to YCbCr */
void rgb2ycc(double r, double g, double b, double *y, double *cb, double *cr)
{
	(*y)  =  0.29900*r + 0.58700*g + 0.11400*b;
	(*cb) = -0.16874*r - 0.33126*g + 0.50000*b + 128.0;
	(*cr) =  0.50000*r - 0.41869*g - 0.08131*b + 128.0;

	return;
}

/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
/* YCbCr to RGB */
void ycc2rgb(double y,double cb, double cr, double *r, double *g, double *b)
{
	cb	  = cb - 128.0;
	cr	  = cr - 128.0;
	(*r) = y + 1.40200*cr;
	(*g) = y - 0.34414*cb - 0.71414*cr;
	(*b) = y + 1.77200*cb;
}


/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
/* Rescales a matrix using bicubic spline interpolation 
 * (taken from Nachtwerk audioDSP library). Coded by Makus Fick
 *
 */
PIMAGE RescaleImage(PIMAGE pimageSrc, int newWidth, int newHeight, double zeroVal)
{	
   int srcXI, srcYI;
	int xx, yy;
   int width, height;
   double xstep, ystep;
   double **matrix1, **matrix2;
   double srcX, srcXF, srcY, srcYF;
   double ip0, ip1, ip2, ip3, a0, a1, a2, a3, yd1, yd2, yd3, inv;
   PIMAGE pimageDest;

	width = pimageSrc->xsize;
   height= pimageSrc->ysize;   	
   xstep = ((double)width) / ((double)newWidth);
   ystep = ((double)height) / ((double)newHeight);
   
   pimageDest = ImageAlloc(newWidth, newHeight);

   if (pimageDest==NULL){
      return NULL;
   }
   
   if ((matrix1 = Double_Alloc(height, newWidth))==NULL){
      ImageFree(pimageDest);
      return NULL;
   }

   if ((matrix2 = Double_Alloc(newHeight, newWidth))==NULL){
      ImageFree(pimageDest);
      free(matrix1[0]);
      return NULL;
   }
   
   /* first step: decimate width */
   if( newWidth != width ){
      for (yy=0; yy < height; yy++ ){
         for (xx=0; xx<newWidth; xx++){
            srcX = xstep*xx;
            srcXI= (int)srcX;
            
            srcXF= srcX - (double)srcXI;
            
            ip0 = (( ( srcXI   < 0 ) || ( srcXI   >= width ) ) ? zeroVal : pimageSrc->pixel[yy][srcXI]);
            ip1 = (( ( srcXI+1 < 0 ) || ( srcXI+1 >= width ) ) ? zeroVal : pimageSrc->pixel[yy][srcXI+1]);
            ip2 = (( ( srcXI+2 < 0 ) || ( srcXI+2 >= width ) ) ? zeroVal : pimageSrc->pixel[yy][srcXI+2]);
				ip3 = (( ( srcXI+3 < 0 ) || ( srcXI+3 >= width ) ) ? zeroVal : pimageSrc->pixel[yy][srcXI+3]);
				
            a0  = ip0;
            a1  = ( 7.0 * ip3 - 31.5 * ip2 + 63.0 * ip1 - 38.5 * ip0) / 21.0;
            yd1 = (-3.5 * ip3 + 21.0 * ip2 - 10.5 * ip1 -  7.0 * ip0) / 21.0;
            yd2 = ip1 - ip0;
            yd3 = a1 + yd1;
            a3  = yd3 - 2.0 * yd2;
            a2  = yd2 - a3 - a1;
            
            inv = ((a3 * srcXF + a2)* srcXF + a1)* srcXF + a0;

            inv = (inv < 0.0) ? 0.0 : ((inv > 255.0) ? 255.0 : inv);
				matrix1[yy][xx] = inv;
         }
      }
   }
   else{
      for(yy=0; yy<height; yy++ ){
         for(xx=0; xx<width; xx++ ){
            matrix1[yy][xx] = pimageSrc->pixel[yy][xx];
         }
      }
   }

   /* second step: decimate height */
   if (newHeight != height){
      for( xx=0; xx<newWidth; xx++){
         for( yy=0; yy<newHeight; yy++ ){
            srcY = ystep*yy;
            srcYI	= (int)srcY;
            
            srcYF	= srcY - (double)srcYI;
            ip0 = (( ( srcYI   < 0 ) || ( srcYI   >= height ) ) ? zeroVal : matrix1[srcYI  ][xx]);
            ip1 = (( ( srcYI+1 < 0 ) || ( srcYI+1 >= height ) ) ? zeroVal : matrix1[srcYI+1][xx]);
            ip2 = (( ( srcYI+2 < 0 ) || ( srcYI+2 >= height ) ) ? zeroVal : matrix1[srcYI+2][xx]);
            ip3 = (( ( srcYI+3 < 0 ) || ( srcYI+3 >= height ) ) ? zeroVal : matrix1[srcYI+3][xx]);

            a0  = ip0;
            a1  = ( 7.0 * ip3 - 31.5 * ip2 + 63.0 * ip1 - 38.5 * ip0) / 21.0;
            yd1 = (-3.5 * ip3 + 21.0 * ip2 - 10.5 * ip1 -  7.0 * ip0) / 21.0;
				yd2 = ip1 - a0;
            yd3 = a1 + yd1;
            a3  = yd3 - 2.0 * yd2;
            a2  = yd2 - a3 - a1;

            inv = ((a3 * srcYF + a2) * srcYF + a1)* srcYF + a0;

            inv = (inv < 0.0) ? 0.0 : ((inv > 255.0) ? 255.0 : inv);
            matrix2[yy][xx] = inv;
         }
      }
   }
	else{
      for(xx=0; xx<newWidth; xx++ ){
         for(yy=0; yy<height; yy++ ){
            matrix2[yy][xx] = matrix1[yy][xx];
         }
      }
   }

   /* copy */
   for (yy=0; yy<newHeight; yy++){
      for (xx=0; xx<newWidth; xx++){
         pimageDest->pixel[yy][xx]=(unsigned char)matrix2[yy][xx];
      }
   }

   free(matrix1[0]);
   free(matrix2[0]);

   return pimageDest;
   
}

/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
void ImageError(char *fmt, ...)
{
	va_list argptr;
	
	va_start( argptr, fmt );
	fprintf(stderr, "ImageError: " );
	vprintf( fmt, argptr );
	va_end( argptr );
	exit( -1 );
}

/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
void ImageWarning(char *fmt, ...)
{
	va_list argptr;
	
	va_start( argptr, fmt );
	fprintf( stderr, "ImageWarning: " );
	vprintf( fmt, argptr );
	va_end( argptr );
}