ezw3.c

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

/*--------------------------------------------------------------------------------------------*/
#include "ezw3.h"


/*--------------------------------------------------------------------------------------------*/
#define _CODE_MEAN_                 1
#define _REFINEMENT_SAME_PASS_      0
#define _MERGE_NODE_CONTEXT_        1
#define _MERGE_REFI_CONTEXT_        1
#define _RESET_MODELS_              1

/*--------------------------------------------------------------------------------------------*/
#define DEFAULT_MEAN 128

/*--------------------------------------------------------------------------------------------*/
#define SSIG_MAX	256
#define REFI_MAX	256

#define nSSigModel   4
#define nRefiModel   1

/*--------------------------------------------------------------------------------------------*/
Timer CodeTime, TotalTime;

/*--------------------------------------------------------------------------------------------*/
/* Code alphabet */
#define _PS_ 0x0
#define _NG_ 0x1
#define _ZS_ 0x2
#define _IZ_ 0x3

/*--------------------------------------------------------------------------------------------*/
/* Hilbert scan paths */
#define RS  0
#define LS  1
#define US  2
#define DS  3

#define MoveUp(px, py)     ((*py)--)
#define MoveDown(px, py)   ((*py)++)
#define MoveLeft(px, py)   ((*px)--)
#define MoveRight(px, py)  ((*px)++)

/*--------------------------------------------------------------------------------------------*/
#define CHKBIT(x, y) (((x) & (y))? 1 : 0)
#define SETBIT(x, y) ((x) | (y))
#define CLRBIT(x, y) ((x) & (~(y)))

/*--------------------------------------------------------------------------------------------*/
/*--------------------------------------------------------------------------------------------*/
void ChildrenPareCtxBitUpdate(Coder *coder, int n, int x, int y);
void SSigContextUpdate(Coder *coder, int i, int x, int y, int n);
void CSigContextUpdate(Coder *coder, int i, int x, int y, int n);
int  CSigContext(Coder *coder, int i, int n, int LastScanSig);
/*--------------------------------------------------------------------------------------------*/

/*--------------------------------------------------------------------------------------------*/
int Encode(Coder *coder)
{
	int i;
   REAL mean;

	coder->Mode = ENCODE;
	coder->CodedSymbols = 0;

	sprintf(coder->InputImageFileName, "%s.pgm", coder->InputImageName);
	if(!coder->EncodeFileNameSpecify){
		sprintf(coder->EncodeImageFileName, "%s_%.6f_%d_%s.%s", coder->InputImageName, 
			coder->TargetBitRate, coder->NScale, SubbandFilterName[coder->WaveletIndex],
			ENC_EXT);
	}
	
	if ((coder->InputImage = FImageReadPGM (coder->InputImageFileName))==NULL){
		Error("Encode: Fail to read input image: %s\n", coder->InputImageFileName);
	}

	coder->ImageXSize = coder->InputImage->XSize;
	coder->ImageYSize = coder->InputImage->YSize;

	/* check size constraints */
	if ((coder->ImageXSize>>(coder->NScale+1))<<(coder->NScale+1) != coder->ImageXSize || 
		 (coder->ImageYSize>>(coder->NScale+1))<<(coder->NScale+1) != coder->ImageYSize){
		Error("Encode: Image size is not consistent with number of scales.\n");
	}

#if _CODE_MEAN_

   mean = 0.0;
   for (i=0; i<coder->ImageXSize*coder->ImageYSize; i++){
		mean +=coder->InputImage->Pixel[i];
	}
   mean /= (REAL)(coder->ImageXSize*coder->ImageYSize);

   coder->mean = (int)(0.5 + mean);
   mean = (REAL)coder->mean;

#else
   
   mean = DEFAULT_MEAN;

#endif

	for (i=0; i<coder->ImageXSize*coder->ImageYSize; i++){
		coder->InputImage->Pixel[i] -= mean;
	}

	/* compute target bits, assume data is 8 bits per sample */
	coder->TargetBits = (int)(coder->TargetBitRate*coder->ImageXSize*coder->ImageYSize);
	coder->TargetBits = ((coder->TargetBits+7)>>3)<<3;
	
	fprintf(stdout, "Encoding (%s): %s (%d x %d) to %s\n", BUILD_STAMP, 
		coder->InputImageFileName,	coder->ImageXSize, coder->ImageYSize, 
		coder->EncodeImageFileName);

	fprintf(stdout, "Wavelet : %s, Levels : %d, Bitrate : %.4lf (%d bytes)\n",
		SubbandFilterName[coder->WaveletIndex], coder->NScale,	coder->TargetBitRate, 
		(int)((coder->TargetBits+7)>>3 ));
	
	/* transform */
	InitializeSubbandBuffers(coder->ImageYSize, coder->ImageXSize);

   Subband2DAnal(coder->InputImage->Pixel, coder->ImageYSize, 
		coder->ImageXSize, coder->NScale, coder->WaveletIndex);
	
	coder->dwt = DWTransformAlloc(coder->ImageXSize, coder->ImageYSize, coder->NScale);
	DWTransformInputFromMallat(coder->dwt, coder->InputImage->Pixel);
	FImageDealloc(coder->InputImage);

	/* coding structures */
	BuildCodingArrays(coder);
	ComputeMaxMagnitudeTree(coder);

	DWTransformDealloc(coder->dwt);
	FreeSubbandBuffers();

	/* setup encode stream */
	coder->BitsCount = 0;
	coder->Encoder = AREncoderAlloc();
	AREncoderOpenFile(coder->Encoder, coder->EncodeImageFileName);

   coder->TargetBits -= 8;

   AllocModels(coder);
	
	WriteHeader(coder);
	
	MainCodingPass(coder);
	
   DeallocModels(coder);	

	AREncoderCloseFile(coder->Encoder);
	AREncoderDealloc(coder->Encoder);

	DestroyCodingArrays(coder);

   return 0;
}

/*--------------------------------------------------------------------------------------------*/
int Decode(Coder *coder)
{
	int i, n, value;
	struct stat statistics; 
	int EncodedImageFileSize;
   REAL mean;

	coder->Mode = DECODE;
	coder->CodedSymbols = 0;

	/* get the encoded stream size - for maximum decode bits count */
 	if(stat(coder->EncodeImageFileName, &statistics) == -1){
		EncodedImageFileSize = 0;
	}
	else{
		EncodedImageFileSize = (int)statistics.st_size;
	}
	EncodedImageFileSize <<= 3;
	/* temporary for reading header */
	coder->TargetBits = EncodedImageFileSize;
	
	/* open encoded stream */
	coder->BitsCount = 0;
	coder->Decoder = ARDecoderAlloc();
	ARDecoderOpenFile(coder->Decoder, coder->EncodeImageFileName);
	
	ReadHeader(coder);

	coder->TargetBits = (int)(coder->TargetBitRate*coder->ImageXSize*coder->ImageYSize);
	coder->TargetBits = ((coder->TargetBits+7)>>3)<<3;
	
	if (coder->TargetBits == 0 || coder->TargetBits > EncodedImageFileSize){
		coder->TargetBits = EncodedImageFileSize;
		coder->TargetBitRate = (double)coder->TargetBits/(coder->ImageXSize*coder->ImageYSize);
	}
		
	/* build the decode filename */
	if (!coder->DecodeFileNameSpecify){
		if (coder->TargetBits==EncodedImageFileSize){
			sprintf(coder->DecodeImageFileName, "%s.%s.pgm", coder->EncodeImageName, ENC_EXT);	
		}
		else{
			sprintf(coder->DecodeImageFileName, "%s_%d.%s.pgm", coder->EncodeImageName,  
				coder->TargetBits>>3, ENC_EXT);
		}
	}
	
	fprintf(stdout, "Decoding (%s): %s to %s (%d x %d)\n", BUILD_STAMP, 
		coder->EncodeImageFileName, coder->DecodeImageFileName, 
		coder->ImageXSize, coder->ImageYSize);

	fprintf(stdout, "Wavelet: %s, Levels: %d, Bitrate: %.4lf (%d bytes)\n", 
		SubbandFilterName[coder->WaveletIndex], coder->NScale, coder->TargetBitRate, 
		(coder->TargetBits+7)>>3);

   coder->TargetBits += 16;

	coder->dwt = DWTransformAlloc(coder->ImageXSize, coder->ImageYSize, coder->NScale);
	
	/* Coding structures */
	BuildCodingArrays(coder);

   AllocModels(coder);

	MainCodingPass(coder);

   DeallocModels(coder);	

	ARDecoderCloseFile(coder->Decoder);
	ARDecoderDealloc(coder->Decoder);

	/* build the decoded image */
	InitializeSubbandBuffers(coder->ImageYSize, coder->ImageXSize);
	
	coder->DecodeImage = FImageAlloc(coder->ImageXSize, coder->ImageYSize);
	
	DWTransformOutputToMallat(coder->dwt, coder->DecodeImage->Pixel);

	Subband2DSynt(coder->DecodeImage->Pixel, coder->ImageYSize, 
		coder->ImageXSize, coder->NScale, coder->WaveletIndex);

	DWTransformDealloc(coder->dwt);
	FreeSubbandBuffers();

#if _CODE_MEAN_
   mean = (REAL)(coder->mean);
#else
   mean = DEFAULT_MEAN;
#endif

	for (i=0; i<coder->ImageXSize*coder->ImageYSize; i++){
		coder->DecodeImage->Pixel[i] += mean;
	}
	
	FImageWritePGM(coder->DecodeImage, coder->DecodeImageFileName);
	
	FImageDealloc(coder->DecodeImage);
	DestroyCodingArrays(coder);
	
	return 0;
}

/*--------------------------------------------------------------------------------------------*/
void MainCodingPass(Coder *coder)
{	
   int i, n;

   coder->nPS = 0;
   coder->nNG = 0;
   coder->nIZ = 0;
   coder->nZS = 0;
   coder->nRF = 0;

   DListInit(&coder->LSP, DeleteData, NULL, NULL);
   coder->RefinementMark[0] = NULL;

   fprintf(stdout, "BitP: %-2d LSP: %-8d\n", coder->CurrentBitPlane+1, DListSize(&coder->LSP));
   
   while (coder->CurrentThreshold > 0) {

      if (SortingPass(coder)){
         fprintf(stdout, "End at sorting pass.\n");
         break;
      }
      
      fprintf(stdout, "BitP: %-2d LSP: %-8d\n", coder->CurrentBitPlane, DListSize(&coder->LSP));

      if (RefinementPass(coder)){
         fprintf(stdout, "End at refinement pass.\n");
			break;
      }
     
      for (n=0; n<coder->nSubbands; n++){
         for (i=0; i<coder->SubbandSize[n]; i++){
            coder->State[n][i] = CLRBIT(coder->State[n][i], SCAN_CTX_BIT);
         }
      }

      if (coder->TestFullBitPlane && coder->CurrentBitPlane==coder->LastFullBitPlane){
			break;
		}
      
      coder->CurrentThreshold>>=1;
		coder->CurrentBitPlane--;

#if _RESET_MODELS_
      ResetModels(coder);
#endif

   }

   fprintf(stdout, "BitP: %-2d LSP: %-8d\n", coder->CurrentBitPlane, DListSize(&coder->LSP));

   DListDestroy(&coder->LSP);

	fprintf(stdout, "Coded Symbols = %2d\n", coder->CodedSymbols);
   fprintf(stdout, "PS=%d, NG=%d, IZ=%d, ZS=%d, RF=%d\n", 
      coder->nPS, coder->nNG, coder->nIZ, coder->nZS, coder->nRF);

}


/*--------------------------------------------------------------------------------------------*/
void SetResidualScanState(Coder *coder, int x, int y, int n)
{
   int bx, by, bsize, i;

   x <<=1;
   y <<=1;
   n+=3;
   bsize = 2;
   
   while (n<coder->nSubbands){
      
      for (by=y; by<y+bsize; by++){
         for (bx=x; bx<x+bsize; bx++){
            i = by*coder->SubbandXSize[n] + bx;
            coder->State[n][i] = SETBIT(coder->State[n][i], SCAN_CTX_BIT);
         }
      }
      
      bsize<<=1; 
      x<<=1;
      y<<=1;
      n+=3; // next child subband
      
   }
}

/*--------------------------------------------------------------------------------------------*/
int ProcessNode(Coder *coder, int n, int x, int y, int *LastScanSig)
{
   int i, s, sign, ctx;
   ListData *Data;
   int scan, sig;
   
   // check if out of bound 
   if (x>=coder->SubbandXSize[n] || y>=coder->SubbandYSize[n]){
      // no processing
      return 0;
   }

   i = y*coder->SubbandXSize[n]+x;

   scan = CHKBIT(coder->State[n][i], SCAN_CTX_BIT);
   sig  = CHKBIT(coder->State[n][i], CSIG_CTX_BIT);

   if (scan==0 && sig==0){
         
      // marked as scanned
      coder->State[n][i] = SETBIT(coder->State[n][i], SCAN_CTX_BIT);

      ctx = CSigContext(coder, i, n, *LastScanSig);

      if (coder->Mode==ENCODE){

         if (coder->MaxMagnitude[n][i] & coder->CurrentThreshold){
            // significant set
            
            if (coder->Magnitude[n][i] >= coder->CurrentThreshold){
               
               coder->State[n][i] = SETBIT(coder->State[n][i], CSIG_CTX_BIT);
               
               // significant root                         
               if (CHKBIT(coder->State[n][i], SIGN_CTX_BIT)){
                  // positive
                  s = _PS_;
               }
               else{
                  // negative
                  s = _NG_;	
               }
               
               CSigContextUpdate(coder, i, x, y, n);
               
            }
            else {
               // significant residual 
               s = _IZ_;
            }

            *LastScanSig = 1;

         }
         else{
            // zero set
            s = _ZS_;
            *LastScanSig = 0;
            
            // marked all residual coefficients as scanned
            if (n==0){
               coder->State[1][i] = SETBIT(coder->State[1][i], SCAN_CTX_BIT);
               coder->State[2][i] = SETBIT(coder->State[2][i], SCAN_CTX_BIT);
               coder->State[3][i] = SETBIT(coder->State[3][i], SCAN_CTX_BIT);
               SetResidualScanState(coder, x, y, 1);
               SetResidualScanState(coder, x, y, 2);
               SetResidualScanState(coder, x, y, 3);
            }
            else{
               SetResidualScanState(coder, x, y, n);
            }
         }
#if _MERGE_NODE_CONTEXT_
         if (WriteSymbol(coder, coder->SSigModel[0][ctx], s)){
            return 1;
         }
#else
         if (WriteSymbol(coder, coder->SSigModel[n][ctx], s)){
            return 1;
         }
#endif
         
      }
      else{
         
#if _MERGE_NODE_CONTEXT_
         if (ReadSymbol(coder, coder->SSigModel[0][ctx], &s)){
            return 1;
         }
#else
         if (ReadSymbol(coder, coder->SSigModel[n][ctx], &s)){
            return 1;
         }
#endif
         
         if (s==_ZS_){
            *LastScanSig = 0;
            // marked all residual coefficients as scanned
            if (n==0){
               coder->State[1][i] = SETBIT(coder->State[1][i], SCAN_CTX_BIT);
               coder->State[2][i] = SETBIT(coder->State[2][i], SCAN_CTX_BIT);
               coder->State[3][i] = SETBIT(coder->State[3][i], SCAN_CTX_BIT);
               SetResidualScanState(coder, x, y, 1);
               SetResidualScanState(coder, x, y, 2);
               SetResidualScanState(coder, x, y, 3);
            }
            else{
               SetResidualScanState(coder, x, y, n);
            }
         }
         else{
            if (s==_IZ_){               
            }
            else {
               // significant root
               coder->State[n][i] = SETBIT(coder->State[n][i], CSIG_CTX_BIT);               
               coder->Magnitude[n][i] = 1;
               
               coder->dwt->SubbandPtr[n][i]
                  = (Real)Dequantize(1, coder->CurrentThreshold, MSB_GAMMA);
               
               if (s==_PS_){
                  coder->State[n][i] = SETBIT(coder->State[n][i], SIGN_CTX_BIT);                          
               }
               else{
                  // inverse the sign
                  coder->dwt->SubbandPtr[n][i] *= -1;
               }