img.c

基于嵌入式零树小波编码的C语言程序实现。

            }
         }         
      }
   }
   
   for (k=0; k<3; k++){
      mean[k] /= (img->XSize[0]*img->YSize[0]);
      M[k] = (int)(0.5 + mean[k]);
   }

   for (k = 0; k < 3; k++){
      for (l = k; l < 3; l++) {
         S[k][l] /= (img->XSize[0]*img->YSize[0]);
         S[k][l] -= mean[k]*mean[l];
         S[l][k] = S[k][l];
      }
   } 

   // S is the covariance matrix, R will be the KLT matrix
   Jacobi(mean, S, R);

#if 0
   printf("KLT matrix...\n");
   for (k=0; k<3; k++){
      for (l=0; l<3; l++){
         printf("% -11.8f  ", R[k][l]);
      }
      printf("\n");
   }
   printf("\n");
#endif

   for (k=0; k<3; k++){
      for (l=0; l<3; l++){
         T[k][l]= ROUND(R[k][l]*(1<<precision));
      }
   }

   return;
}

/*--------------------------------------------------------------------------------------------*/
void FCoImageForwardKLT(FCoImage *img, int M[3], int T[3][3], int precision)
{
   REAL R[3][3], p[3], s;
   int k, l, i, j, t;
   
   for (k=0; k<3; k++){
      for (l=0; l<3; l++){
         R[k][l] = (REAL)T[k][l]/(1<<precision);
      }
   }

   for (j=0; j<img->YSize[0]; j++){
      for (i=0; i<img->XSize[0]; i++){
         
         t = j*img->XSize[0]+i;

         for (k=0; k<3; k++){
            p[k] = img->Pixel[k][t] - M[k];
         }

         for (k=0; k<3; k++){
            s = 0;
            for (l = 0; l < 3; l++){
               s += R[k][l] * p[l];
            }
            img->Pixel[k][t] = s;
         }
      }
   }

   return;
}

/*--------------------------------------------------------------------------------------------*/
void FCoImageInverseKLT(FCoImage *img, int M[3], int T[3][3], int precision)
{
   REAL R[3][3], p[3], s;
   int k, l, i, j, t;
   
   for (k=0; k<3; k++){
      for (l=0; l<3; l++){
         R[k][l] = (REAL)T[k][l]/(1<<precision);
      }
   }

   for (j=0; j<img->YSize[0]; j++){
      for (i=0; i<img->XSize[0]; i++){
         
         t = j*img->XSize[0]+i;

         for (k=0; k<3; k++){
            p[k] = img->Pixel[k][t];
         }

         for (k=0; k<3; k++){
            s = 0;
            for (l = 0; l < 3; l++){
               s += R[l][k] * p[l];
            }
            img->Pixel[k][t] = s + M[k];
         }
      }
   }

   return;


}


/*--------------------------------------------------------------------------------------------*/
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 );
}

/*--------------------------------------------------------------------------------------------*/
// 0: PGM, 1:PPM, 2: raw/unknown
int ImageType(char *filename, int *XSize, int *YSize)
{
   FILE *fp;
   char type_p, type_c, str[256];

   if ((fp=fopen(filename, "rb"))==NULL){
		ImageWarning("Fail to read image.\n");		
      return -1;
	}
	
   fread(&type_p, sizeof(char), 1, fp);
   fread(&type_c, sizeof(char), 1, fp);
   rewind(fp);

   if ((type_p == 'P') && (type_c == '5' || type_c== '6')){
		fgets (str, 256, fp);
	   do {
         fgets (str, 256, fp);
      } while (str[0] == '#');
		
      sscanf (str, "%d%d", XSize, YSize);
		fgets (str, 256, fp);   // maximum value, always 255

      return (type_c == '5' ? 0: 1);   // PGM
	}
   else {
      *XSize = 0;
      *YSize = 0;
      return 2;
   }
}

/*--------------------------------------------------------------------------------------------*/
double ImageCompareMSE(Image *img1, Image *img2)
{
   int i;
   double Error, SumOfSquareError=0.0;

   if ((img1->XSize != img2->XSize) || (img1->YSize!=img2->YSize)){
      ImageWarning("Cannot compare images of different sizes (%dx%d and %dx%d)\n",
         img1->XSize, img1->YSize, img2->XSize, img2->YSize);
      return 0.0;
   }

   for (i=0; i< img1->XSize*img1->YSize; i++){
      Error=(int)img1->Pixel[i]-(int)img2->Pixel[i];
      SumOfSquareError+=Error*Error;
   }

   return SumOfSquareError/(double)(img1->XSize*img1->YSize);

}

/*--------------------------------------------------------------------------------------------*/
double ImageCompareMAD(Image *img1, Image *img2)
{
   int i;
   double SumOfMAD=0.0;

   if ((img1->XSize != img2->XSize) || (img1->YSize!=img2->YSize)){
      ImageWarning("Cannot compare images of different sizes (%dx%d and %dx%d)\n",
         img1->XSize, img1->YSize, img2->XSize, img2->YSize);
      return 0.0;
   }

   for (i=0; i< img1->XSize*img1->YSize; i++){
      SumOfMAD+=fabs((double)((int)img1->Pixel[i]-(int)img2->Pixel[i]));
   }

   return SumOfMAD/(double)(img1->XSize*img1->YSize);

}

/*--------------------------------------------------------------------------------------------*/
double ImageCompareSNR(Image *img1, Image *img2)
{
   int i;
   double SumOfSquareSignal=0.0, SumOfSquareError=0.0;
	double Error, RootMeanSignalError, RootMeanSquareError;

   if ((img1->XSize != img2->XSize) || (img1->YSize!=img2->YSize)){
      ImageWarning("Cannot compare images of different sizes (%dx%d and %dx%d)\n",
         img1->XSize, img1->YSize, img2->XSize, img2->YSize);
      return 0.0;
   }

   for (i=0; i< img1->XSize*img1->YSize; i++){
		SumOfSquareSignal += (img1->Pixel[i]*img1->Pixel[i]);
      Error=(int)img1->Pixel[i]-(int)img2->Pixel[i];
      SumOfSquareError+=Error*Error;
   }

   RootMeanSignalError = sqrt(SumOfSquareSignal/(double)(img1->XSize*img1->YSize));
	RootMeanSquareError = sqrt(SumOfSquareError/(double)(img1->XSize*img1->YSize));
   
   return 20*log10(RootMeanSignalError/RootMeanSquareError); /* assume maximum 255 */
}

/*--------------------------------------------------------------------------------------------*/
double ImageComparePSNR(Image *img1, Image *img2)
{
   int i;
   double Error, SumOfSquareError=0.0, RootMeanSquareError;

   if ((img1->XSize != img2->XSize) || (img1->YSize!=img2->YSize)){
      ImageWarning("Cannot compare images of different sizes (%dx%d and %dx%d)\n",
         img1->XSize, img1->YSize, img2->XSize, img2->YSize);
      return 0.0;
   }

   for (i=0; i< img1->XSize*img1->YSize; i++){
      Error=(int)img1->Pixel[i]-(int)img2->Pixel[i];
      SumOfSquareError+=Error*Error;
   }

   RootMeanSquareError = sqrt(SumOfSquareError/(double)(img1->XSize*img1->YSize));
   
   return 20*log10(255.0/RootMeanSquareError); /* assume maximum 255 */
}

/*--------------------------------------------------------------------------------------------*/
double ImageFileCompareMSE(char *PGMFileName1, char *PGMFileName2)
{
	Image *img1;
	Image *img2;
   double MSE;

	if ((img1 = ImageReadPGM(PGMFileName1))==NULL){
		ImageWarning("Fail to open PGM files.\n");		
		return 0.0;
	}

	if ((img2 = ImageReadPGM(PGMFileName2))==NULL){
		ImageDealloc(img1);
		ImageWarning("Fail to open PGM files.\n");		
		return 0.0;
	}

	if ((img1->XSize != img2->XSize) || (img1->YSize!=img2->YSize)){
      ImageWarning("Cannot compare images of different sizes (%dx%d and %dx%d)\n",
         img1->XSize, img1->YSize, img2->XSize, img2->YSize);
      return 0.0;
   }

	MSE = ImageCompareMSE(img1, img2);

	ImageDealloc(img1);
	ImageDealloc(img2);

   return MSE;
}

/*--------------------------------------------------------------------------------------------*/
/* The first image is assumed to be the original */
double ImageFileCompareSNR(char *PGMFileName1, char *PGMFileName2)
{
	Image *img1;
	Image *img2;
   double SNR;

	if ((img1 = ImageReadPGM(PGMFileName1))==NULL){
		ImageWarning("Fail to open PGM files.\n");		
		return 0.0;
	}

	if ((img2 = ImageReadPGM(PGMFileName2))==NULL){
		ImageDealloc(img1);
		ImageWarning("Fail to open PGM files.\n");		
		return 0.0;
	}

	if ((img1->XSize != img2->XSize) || (img1->YSize!=img2->YSize)){
      ImageWarning("Cannot compare images of different sizes (%dx%d and %dx%d)\n",
         img1->XSize, img1->YSize, img2->XSize, img2->YSize);
      return 0.0;
   }

	SNR = ImageCompareSNR(img1, img2);

	ImageDealloc(img1);
	ImageDealloc(img2);

   return SNR;
}

/*--------------------------------------------------------------------------------------------*/
double ImageFileComparePSNR(char *PGMFileName1, char *PGMFileName2)
{
	Image *img1;
	Image *img2;
   double PSNR;

	if ((img1 = ImageReadPGM(PGMFileName1))==NULL){
		ImageWarning("Fail to open PGM files.\n");		
		return 0.0;
	}

	if ((img2 = ImageReadPGM(PGMFileName2))==NULL){
		ImageDealloc(img1);
		ImageWarning("Fail to open PGM files.\n");		
		return 0.0;
	}

	if ((img1->XSize != img2->XSize) || (img1->YSize!=img2->YSize)){
      ImageWarning("Cannot compare images of different sizes (%dx%d and %dx%d)\n",
         img1->XSize, img1->YSize, img2->XSize, img2->YSize);
      return 0.0;
   }

	PSNR = ImageComparePSNR(img1, img2);

	ImageDealloc(img1);
	ImageDealloc(img2);

	return PSNR;
}

/*--------------------------------------------------------------------------------------------*/
double ImageFileCompareMAD(char *PGMFileName1, char *PGMFileName2)
{
	Image *img1;
	Image *img2;
   double MAD;

	if ((img1 = ImageReadPGM(PGMFileName1))==NULL){
		ImageWarning("Fail to open PGM files.\n");		
		return 0.0;
	}

	if ((img2 = ImageReadPGM(PGMFileName2))==NULL){
		ImageDealloc(img1);
		ImageWarning("Fail to open PGM files.\n");		
		return 0.0;
	}

	if ((img1->XSize != img2->XSize) || (img1->YSize!=img2->YSize)){
      ImageWarning("Cannot compare images of different sizes (%dx%d and %dx%d)\n",
         img1->XSize, img1->YSize, img2->XSize, img2->YSize);
      return 0.0;
   }

	MAD = ImageCompareMAD(img1, img2);

	ImageDealloc(img1);
	ImageDealloc(img2);

	return MAD;
}

/*--------------------------------------------------------------------------------------------*/
double FImageCompareMSE(FImage *img1, FImage *img2)
{
   int i;
   double Error, SumOfSquareError=0.0;

   if ((img1->XSize != img2->XSize) || (img1->YSize!=img2->YSize)){
      ImageWarning("Cannot compare images of different sizes (%dx%d and %dx%d)\n",
         img1->XSize, img1->YSize, img2->XSize, img2->YSize);
      return 0.0;
   }

   for (i=0; i< img1->XSize*img1->YSize; i++){
      Error=img1->Pixel[i]-img2->Pixel[i];
      SumOfSquareError+=Error*Error;
   }

   return SumOfSquareError/(double)(img1->XSize*img1->YSize);

}

/*--------------------------------------------------------------------------------------------*/
double FImageComparePSNR(FImage *img1, FImage *img2)
{
	   int i;
   double Error, SumOfSquareError=0.0, RootMeanSquareError;

   if ((img1->XSize != img2->XSize) || (img1->YSize!=img2->YSize)){
      ImageWarning("Cannot compare images of different sizes (%dx%d and %dx%d)\n",
         img1->XSize, img1->YSize, img2->XSize, img2->YSize);
      return 0.0;
   }

   for (i=0; i< img1->XSize*img1->YSize; i++){
      Error=img1->Pixel[i]-img2->Pixel[i];
      SumOfSquareError+=Error*Error;
   }

   RootMeanSquareError = sqrt(SumOfSquareError/(double)(img1->XSize*img1->YSize));
   
   return 20*log10(255.0/RootMeanSquareError); /* assume maximum 255 */
}
/*--------------------------------------------------------------------------------------------*/