img.c

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

}


//--------------------------------------------------------------------------------------------
// Read PPM image
//--------------------------------------------------------------------------------------------
CoImage *CoImageReadPPM(char *PPMFileName)
{
	FILE *fp;
	CoImage *img;
	int v, i, c;	
	char type_p, type_c, str[256];
	int XSize[3], YSize[3];

	if ((fp=fopen(PPMFileName, "rb"))==NULL){
		ImageWarning("Fail to open file [%s].\n", PPMFileName);
		return NULL;
	}
	
   fread(&type_p, sizeof(char), 1, fp);
   fread(&type_c, sizeof(char), 1, fp);
   rewind(fp);
	
   if (type_p!='P'){
      ImageWarning("Not PGM/PPM file.\n");
		fclose(fp);
		return NULL;
   }

   if (type_c != '5' && type_c != '6'){
		ImageWarning("Only binary PPM/PGM is supported.\n");
		fclose(fp);
		return NULL;
	}
   
   fgets (str, 256, fp);
   do {     // skip the comments
      fgets (str, 256, fp);
   } while (str[0] == '#');

   sscanf (str, "%d%d", &XSize[0], &YSize[0]);   
   fgets (str, 256, fp);   // maximum value, always 255
   
	XSize[1] = XSize[2] = XSize[0];
   YSize[1] = YSize[2] = YSize[0];
	
   if (type_c == '5'){  // PGM		
      img = CoImageAlloc(XSize, YSize, 1);
      fseek(fp, -1*img->XSize[0]*img->YSize[0]*img->color, SEEK_END);
      if ((int)fread(img->Pixel[0], sizeof(unsigned char),	
         img->XSize[0]*img->YSize[0], fp) != img->XSize[0]*img->YSize[0]){	
         
         ImageWarning("Fail to read PPM image data.\n");
         fclose(fp);
         CoImageDealloc(img);
         return NULL;
      }
	}
	else {   // PPM
		img = CoImageAlloc(XSize, YSize, 3);
      fseek(fp, -1*img->XSize[0]*img->YSize[0]*img->color, SEEK_END);
      for (i=0; i<img->XSize[0]*img->YSize[0]; i++){
         for (c=0; c<3; c++){
            if ((v = fgetc(fp)) == EOF){
               ImageWarning("Fail to read PPM image data\n");
               fclose(fp);
               CoImageDealloc(img);
               return NULL;
            }
            img->Pixel[c][i] = v;
         }
      }
	}
		
	fclose(fp);
	return img;
}

//--------------------------------------------------------------------------------------------
LCoImage *LCoImageReadPPM(char *PPMFileName)
{
	LCoImage *img;
   CoImage *temp = CoImageReadPPM(PPMFileName);
   
   if (temp!=NULL){
      img = LCoImageAlloc(temp->XSize, temp->YSize, temp->color);      
      CoImageToLCoImage(img, temp);      
      CoImageDealloc(temp);      
      return img;
   }
   else{
      return NULL;
   }
}

//--------------------------------------------------------------------------------------------
FCoImage *FCoImageReadPPM(char *PPMFileName)
{
	FCoImage *img;
   CoImage *temp = CoImageReadPPM(PPMFileName);
   
   if (temp!=NULL){
      img = FCoImageAlloc(temp->XSize, temp->YSize, temp->color);      
      CoImageToFCoImage(img, temp);      
      CoImageDealloc(temp);      
      return img;
   }
   else{
      return NULL;
   }
}

//--------------------------------------------------------------------------------------------
// Write PPM image
//--------------------------------------------------------------------------------------------
int CoImageWritePPM(CoImage *img, char *PPMFileName)
{
	FILE *fp;
	int i, c;

	if ((fp = fopen(PPMFileName,"wb")) == NULL){
		ImageWarning("Fail to create file [%s]\n", PPMFileName);
		return -1;
	}

	if (img->color == 1){
		fprintf(fp, "P5\n%d %d\n255\n", img->XSize[0], img->YSize[0]);
		if ((int)fwrite(img->Pixel[0], sizeof(unsigned char), 
         img->XSize[0]*img->YSize[0], fp) != img->XSize[0]*img->YSize[0]){

			ImageWarning("Fail to write image data\n");
			fclose(fp);
			return -1;
		}
	}
   else{
      for (c=1; c<3; c++){
         if (img->XSize[0] != img->XSize[c] || img->YSize[0] != img->YSize[c]){
            ImageWarning("PPM format must have same size components.\n");
            fclose(fp);
            return -1;
         }
      }   

		fprintf(fp,"P6\n%d %d\n255\n", img->XSize[0], img->YSize[0]);		
      for (i=0; i<img->XSize[0]*img->YSize[0]; i++){
         for (c=0; c<3; c++){
            if (fputc(img->Pixel[c][i], fp) == EOF){
               ImageWarning("Fail to write PPM image data\n");
               fclose(fp);
               return -1;
            }
         }
      }		
	}

   fclose(fp);
	return 0;
}

/*--------------------------------------------------------------------------------------------*/
int LCoImageWritePPM(LCoImage *img, char *PPMFileName)
{
   int i;
	CoImage *temp = CoImageAlloc(img->XSize, img->YSize, img->color);

   if (temp!=NULL){
      LCoImageToCoImage(temp, img);
      i = CoImageWritePPM(temp, PPMFileName);
      CoImageDealloc(temp);
      return i;
   }
   else{
      return 0;
   }
}

//--------------------------------------------------------------------------------------------
int FCoImageWritePPM(FCoImage *img, char *PPMFileName)
{
   int i;
	CoImage *temp = CoImageAlloc(img->XSize, img->YSize, img->color);

   if (temp!=NULL){
      FCoImageToCoImage(temp, img);
      i = CoImageWritePPM(temp, PPMFileName);
      CoImageDealloc(temp);
      return i;
   }
   else{
      return 0;
   }
}

//--------------------------------------------------------------------------------------------
// Copy image components
//--------------------------------------------------------------------------------------------
void CoImageCopy(CoImage *pDest, CoImage *pSrc)
{
   int c;
	assert(pDest);	assert(pSrc);
   if (pDest->color!=pSrc->color){
      ImageError("Different number of color components.\n");
   }
   for (c=0; c<pDest->color; c++){
      memcpy(pDest->Pixel[c], pSrc->Pixel[c], 
         sizeof(unsigned char)*pSrc->XSize[c]*pSrc->YSize[c]);
   }
}

//--------------------------------------------------------------------------------------------
void LCoImageCopy(LCoImage *pDest, LCoImage *pSrc)
{
   int c;
	assert(pDest);	assert(pSrc);
   if (pDest->color!=pSrc->color){
      ImageError("Different number of color components.\n");
   }
   for (c=0; c<pDest->color; c++){
      memcpy(pDest->Pixel[c], pSrc->Pixel[c], 
         sizeof(int)*pSrc->XSize[c]*pSrc->YSize[c]);
   }
}

//--------------------------------------------------------------------------------------------
void FCoImageCopy(FCoImage *pDest, FCoImage *pSrc)
{
   int c;
	assert(pDest);	assert(pSrc);
   if (pDest->color!=pSrc->color){
      ImageError("Different number of color components.\n");
   }
   for (c=0; c<pDest->color; c++){
      memcpy(pDest->Pixel[c], pSrc->Pixel[c], 
         sizeof(REAL)*pSrc->XSize[c]*pSrc->YSize[c]);
   }
}

//--------------------------------------------------------------------------------------------
// Convert from 8 bits components 
//--------------------------------------------------------------------------------------------
void CoImageToLCoImage(LCoImage *pDest, CoImage *pSrc)
{
	int i, c;
	assert(pDest);	assert(pSrc);
   for (c=0; c<pDest->color; c++){
      for (i=0; i<pDest->XSize[c]*pDest->YSize[c]; i++){
         pDest->Pixel[c][i] = (int) pSrc->Pixel[c][i];						  
      }
	}
}

//--------------------------------------------------------------------------------------------
void CoImageToFCoImage(FCoImage *pDest, CoImage *pSrc)
{
	int i, c;
	assert(pDest);	assert(pSrc);
   for (c=0; c<pDest->color; c++){
      for (i=0; i<pDest->XSize[c]*pDest->YSize[c]; i++){
         pDest->Pixel[c][i] = (REAL) pSrc->Pixel[c][i];						  
      }
	}
}

//--------------------------------------------------------------------------------------------
// Convert to 8 bits image components
//--------------------------------------------------------------------------------------------
void LCoImageToCoImage(CoImage *pDest, LCoImage *pSrc)
{
	int i, c;
   int temp;
	assert(pDest);	assert(pSrc);

   for (c=0; c<pDest->color; c++){
      for (i=0; i<pDest->XSize[c]*pDest->YSize[c]; i++){
         temp = (int) (pSrc->Pixel[c][i]);
         pDest->Pixel[c][i] 
            = (unsigned char) (temp<0 ? 0 : (temp>255 ? 255 : temp));
      }      
	}
}

//--------------------------------------------------------------------------------------------
void FCoImageToCoImage(CoImage *pDest, FCoImage *pSrc)
{
	int i, c;
   REAL temp;
	assert(pDest);	assert(pSrc);

   for (c=0; c<pDest->color; c++){
      for (i=0; i<pDest->XSize[c]*pDest->YSize[c]; i++){
         temp = (REAL) (pSrc->Pixel[c][i] + 0.499);
         pDest->Pixel[c][i] 
            = (unsigned char) (temp < 0.0 ? 0 : (temp > 255.0 ? 255 : temp));
      }      
	}
}


/*--------------------------------------------------------------------------------------------*/
/*--------------------------------------------------------------------------------------------*/
/* The following codes for computing the inter-component KLT matrix is obtained
 * from the source code of SPIHT, copyrighted by Amir Said and W.A. Pearlman.
 * Please refer to the source code for the exact codes.
 */
/*--------------------------------------------------------------------------------------------*/
/*--------------------------------------------------------------------------------------------*/

/*--------------------------------------------------------------------------------------------*/
#define ROTATION(M,i,j,k,l) g = M[i][j];  h = M[k][l];\
         M[i][j] = g-s*(h+g*tau);  M[k][l] = h+s*(g-h*tau)

/*--------------------------------------------------------------------------------------------*/
void Jacobi(double m[3], double A[3][3], double R[3][3])
{
   int i, j, k, p, q, itr = 0;
   
   double v, tau, s, t, h, g, c;
   
   double b[3], d[3], X[3][3];

   for (i = 0; i < 3; i++) {
      for (j = i; j < 3; j++) {
         X[i][j] = X[j][i] = 0;         
      }      
      m[i] = 0;
      b[i] = d[i] = A[i][i];  
      X[i][i] = 1; 
   }

   for (;;) {
      
      for (s = i = 0; i < 2; i++){
         for (j = i + 1; j < 3; j++){
            
            if (s < (t = fabs(A[i][j]))) {
               s = t;  
               p = i;  
               q = j; 
            }

         }
      }

      if (s == 0){
         break;
      }

      if (++itr > 50){
         Error("Too many iterations in function Jacobi.\n");
      }

      g = 100 * s;  
      h = d[q] - d[p];

      if (g + fabs(h) == fabs(h)){
         t = A[p][q] / h;
      }
      else {
         v = 0.5 * h / A[p][q];  
         t = 1 / (fabs(v) + sqrt(1 + v * v));
         if (v < 0){ 
            t = -t; 
         }
      }

      c = 1 / sqrt(1 + t * t);  
      s = t * c;
      tau = s / (1 + c); 
      h = t * A[p][q];
      A[p][q] = 0;
      m[p] -= h; 
      m[q] += h;
      d[p] -= h;
      d[q] += h;

      for (k = 0; k < p; k++)     { ROTATION(A,k,p,k,q); }
      for (k = p + 1; k < q; k++) { ROTATION(A,p,k,k,q); }
      for (k = q + 1; k < 3; k++) { ROTATION(A,p,k,q,k); }
      for (k = 0; k < 3; k++)     { ROTATION(X,k,p,k,q); }
      
      for (k = 0; k < 3; k++) {
         b[k] += m[k]; 
         d[k] = b[k]; 
         m[k] = 0; 
      }
   }
   
   for (i = 0; i < 3; i++) {
      for (s = k = -1, j = 0; j < 3; j++){
         if (d[j] > s) {
            s = d[k=j];
         }
      }
      if (k < 0){ 
         Error("Incorrect decomposition");
      }
      for (d[k] = -2, j = 0; j < 3;j++){ 
         //R[i][j] = 1e-4 * ROUND(1e4 * X[j][k]); 
         R[i][j] = X[j][k]; 
      }
   }
}

/*--------------------------------------------------------------------------------------------*/
void FCoImageKLTMatrix(FCoImage *img, int M[3], int T[3][3], int precision)
{
   double p[3], mean[3], R[3][3], S[3][3];
   int k, l, i, j;

   // form covariance matrix
   for (k = 0; k < 3; k++){      
      mean[k] = 0;
      for (l = 0; l < 3; l++){
         S[k][l] = 0;
      }
   }
   
   // find cross product terms
   for (j=0; j<img->YSize[0]; j++){
      for (i=0; i<img->XSize[0]; i++){
         
         for (k = 0; k < 3; k++) {
            p[k] = img->Pixel[k][j*img->XSize[0]+i];  
            mean[k] += p[k]; 
         }
         
         for (k = 0; k < 3; k++){
            for (l = k; l < 3; l++) {
               S[k][l] += p[k] * p[l];