ezw.c

该程序把数字图像处理与小波变换结合起来

/*----------------------------------------------------------------------------*/	
/*----------------------------------------------------------------------------*/	
/* EZW2 - Embedded Zerotree Wavelet Coder
 * 
 * Author : Mow-Song, Ng
 * Data   : 20-07-2002
 *
 * Last update : 31-08-2002
 *
 * A large portion/idea of the codes come from:
 * C. Valens - Embedded Zerotree Wavelet Encoder Tutorial.
 * G. Davis - Baseline Wavelet Transform Coder Construction Kit
 *
 * The authors copyrights to the codes are acknowledged.
 * The EZW coder is patented by Jerome Shapiro. As far as I know, these are
 * the patent numbers 
 * 5315670 (Issued May 24, 1994), 5321776 (Issued June 14, 1994) and 
 * 5412741 (Issued May 2, 1995).
 *
 *
 * Well, I wrote the coder with some one else codes, add some lines, remove 
 * some lines, sometimes just verbatim. Any way, if you feel that I am worth
 * the acknowledgement, please do so, or perhaps drop me an e-mail.
 *
 * Please refer to my website for more info.
 *
 * My contact:
 * [e-mail]
 * http://www.pesona.mmu.edu.my/~msng
 *
 */
/*----------------------------------------------------------------------------*/	
/*----------------------------------------------------------------------------*/	

#include "ezw.h"

/*----------------------------------------------------------------------------*/	
/*----------------------------------------------------------------------------*/	
void EZW2Encode(WTRANSFORM *transform, ArithEncoder *encoder, 
					 int threshold, int HistoCount, int ByteBudget, 
					 int ScanOrder);

void EZW2BuildZerotreeMap(WTRANSFORM *transform, SubbandSymbol *ztmap);

void EncodeDominantPassMortonScan(DLNode *DomList, SLNode *SubList,
											 int *SubListSize, WTRANSFORM *transform, 
											 MAP *map, int threshold, 
											 ArithEncoder *encoder, 
											 BasicCoder **DominantPassCoder, 
											 SubbandSymbol *ztmap, int ByteBudget);

void EncodeNode(SLNode *SubList, int *SubListSize, WTRANSFORM *transform, 
					 DLNode *node, SubbandSymbol *ztmap, int threshold);

void EncodeSubordinatePass(SLNode *SubList, int SubListSize, int threshold, 
							ArithEncoder *encoder, BasicCoder *SubordinateListCoder, 
							int ByteBudget);

void EZW2Error(char *fmt, ...);
void EZW2Warning(char *fmt, ...);

/* termination flag */
Boolean EndEncoding;

/*----------------------------------------------------------------------------*/	
/*----------------------------------------------------------------------------*/	
void main(int argc, char **argv)
{
	FIMAGE *InputImage, *DecomposeImage;
   WAVELET *wavelet;
	WTRANSFORM *transform;
	char InputImageFileName[_MAX_PATH];
	char DecomposeImageFileName[_MAX_PATH];
	char EncodeImageFileName[_MAX_PATH];
	int NScale;
	int ByteBudget;
	Boolean KeepDecompose = FALSE;
	Boolean EncodeFileNameSpecify = FALSE;
	ArithEncoder *encoder;
	BIT_FILE *EncodeStream;
	float AbsMaxCoeff;
	float ImageMean;
	int Threshold;
	int i, j;
	int MaxHistogramCount;
	int *bin, nbin;
	int total;
	double TargetBitRate;
	char *ProgramName, *InputImageName;
	int WaveletIndex;
	FILTERSET *Filter[16];
	int ScanOrder;
	double t1;
	
	EndEncoding = FALSE;

	fprintf(stdout, "Embedded Zerotree Wavelet Coder - EZW2 (Encoder)\n");
	
	InitializeFilterSets();
	/* Initialize the filters used */
	Filter[0] = Antonini;
	Filter[1] = Adelson;
	Filter[2] = Odegard;
	Filter[3] = Villa1810;
	Filter[4] = Brislawn;
	Filter[5] = Brislawn2;
	Filter[6] = Haar;
	Filter[7] = Daub4;
	Filter[8] = Daub6;
	Filter[9] = Daub8;
	Filter[10] = Villa1;
	Filter[11] = Villa2;
	Filter[12] = Villa3;
	Filter[13] = Villa4;
	Filter[14] = Villa5;
	Filter[15] = Villa6;
	
	ProgramName = argv[0];
	/* default settings */
	NScale = 6;
	TargetBitRate = 0.25;
	MaxHistogramCount = 256;
	WaveletIndex = 0;
	ScanOrder = 0;

	if (argc<3){
      fprintf(stderr, "Usage: %s [image] [bitrate]\n", ProgramName);
		fprintf(stderr, "\t -l <number of levels>\n");
		fprintf(stderr, "\t -w <wavelet>\n");
		fprintf(stderr, "\t     0 =Antonini 1 =Adelson   2 =Odegard 3 =Villa1810\n");
		fprintf(stderr, "\t     4 =Brislawn 5 =Brislawn2 6 =Haar    7 =Daub4\n");
		fprintf(stderr, "\t     8 =Daub6    9 =Daub8     10=Villa1  11=Villa2\n");
		fprintf(stderr, "\t     12=Villa3   13=Villa4    14=Villa5  15=Villa6\n");
		fprintf(stderr, "\t -h <maximum histogram count>\n");
		fprintf(stderr, "\t -e <encode filename>\n");
		fprintf(stderr, "\t -k  keep decompose image\n");
		//fprintf(stderr, "\t -s <scan order>\n");
		//fprintf(stderr, "\t     0=Morton/1=Raster\n");
      return;
	}

	InputImageName = argv[1];
   TargetBitRate  = atof(argv[2]);
   argv+=3; argc-=3;

	while (argc > 0){
      if (!strcmp("-l", *argv)){
			argv++; argc--;
			NScale = atoi(*argv);
			argv++; argc--;
      } 
		else if (!strcmp("-h", *argv)){
			argv++; argc--;
			MaxHistogramCount = atoi(*argv);
			argv++; argc--;
		}
		else if (!strcmp("-w", *argv)){
			argv++; argc--;
			WaveletIndex = atoi(*argv);
			if (WaveletIndex<0 || WaveletIndex > 15){
				WaveletIndex = 0;
			}
			argv++; argc--;
		}
		else if (!strcmp("-e", *argv)){
			argv++; argc--;
			sprintf(EncodeImageFileName, "%s.ezw", *argv);
			EncodeFileNameSpecify = TRUE;
			argv++; argc--; 
		}
		else if (!strcmp("-k", *argv)){
			argv++; argc--;
			KeepDecompose = TRUE;
		}
		/*
		else if (!strcmp("-s", *argv)){
			argv++; argc--;
			ScanOrder = atoi(*argv);
			if (ScanOrder<0 || ScanOrder > 1){
				ScanOrder = 0;
			}
			argv++; argc--;
		}
		*/
		else{
			fprintf(stderr, "Invalid arguement %s\n", *argv);
			fprintf(stderr, "Usage: %s [image] [bitrate]\n", ProgramName);
			fprintf(stderr, "\t -l <number of levels>\n");
			fprintf(stderr, "\t -w <wavelet>\n");
			fprintf(stderr, "\t     0 =Antonini 1 =Adelson   2 =Odegard 3 =Villa1810\n");
			fprintf(stderr, "\t     4 =Brislawn 5 =Brislawn2 6 =Haar    7 =Daub4\n");
			fprintf(stderr, "\t     8 =Daub6    9 =Daub8     10=Villa1  11=Villa2\n");
			fprintf(stderr, "\t     12=Villa3   13=Villa4    14=Villa5  15=Villa6\n");
			fprintf(stderr, "\t -h <maximum histogram count>\n");
			fprintf(stderr, "\t -e <encoded filename>\n");
			fprintf(stderr, "\t -k  keep decompose image\n");
			//fprintf(stderr, "\t -s <scan order>\n");
			//fprintf(stderr, "\t     0=Morton/1=Raster\n");
			return;
		}
	}
      	
	/* build filenames */
	sprintf(InputImageFileName, "%s.pgm", InputImageName);
	sprintf(DecomposeImageFileName, "%s_%d_%s.decompose.pgm", 
		InputImageName, NScale, FilterName[WaveletIndex]);
	
	if(!EncodeFileNameSpecify){
		sprintf(EncodeImageFileName, "%s_%.6f_%d_%s_%d_%d.ezw", InputImageName, 
			TargetBitRate, NScale, FilterName[WaveletIndex], 
			MaxHistogramCount, ScanOrder);
	}

	/* read data */
	if ((InputImage = ReadPGMToFloat(InputImageFileName))==NULL){
		fprintf(stderr, "Fail to read input image: %s\n", InputImageFileName);
		exit(-1);
	}
	
	if ((InputImage->xsize>>NScale)<<NScale != InputImage->xsize || 
		(InputImage->ysize>>NScale)<<NScale != InputImage->ysize){
		EZW2Error("Image size is not consistent with number of scales.\n");
	}

	/* get image mean */
	ImageMean = 0.0;
	for (i=0; i<InputImage->xsize*InputImage->ysize; i++){
		ImageMean += (float)InputImage->pixelLinear[i];
	}
	ImageMean /=(InputImage->xsize*InputImage->ysize);

	/* subtract mean */
	for (i=0; i<InputImage->xsize*InputImage->ysize; i++){
		InputImage->pixelLinear[i] -=ImageMean;
	}

	/* initialize wavelet filter */
	wavelet = WaveletAlloc(Filter[WaveletIndex]);

	/* transform */
	transform = WaveletTransformAlloc(wavelet, InputImage, NScale, -1);
	
	/* decompose image */
	if (KeepDecompose){
		DecomposeImage = WaveletTransformCodeCoeff(transform, 1);
		WriteFloatToPGM(DecomposeImage, DecomposeImageFileName);
		FImageFree(DecomposeImage);
	}

	fprintf(stdout, "Original Image: %s\nEncoded Image : %s\n", 
		InputImageFileName, EncodeImageFileName);
	fprintf(stdout, "Width: %d   Height: %d   ", 
		InputImage->xsize, InputImage->ysize);
	fprintf(stdout, "Mean : %.3f   ", ImageMean);
	
	/* compute the byte budget */
	ByteBudget = (int)(InputImage->xsize*InputImage->ysize*TargetBitRate / 8.0);
	fprintf(stdout, "Bytes: %d (%.5fbpp)\n", ByteBudget, TargetBitRate);
	
	/* get largest magnitude coefficient */
	AbsMaxCoeff = -1.0;
	for (i=0; i<transform->vsize*transform->hsize; i++){
		if (fabs(transform->value[i]) > AbsMaxCoeff){
			AbsMaxCoeff = (float)fabs(transform->value[i]);
		}
	}
	fprintf(stdout, "Abs Max Coeff : %f\n", AbsMaxCoeff);
	
	/* get threshold */
	Threshold = 1 << (int)(floor(log10(AbsMaxCoeff)/log10(2)));
	fprintf(stdout, "Initial threshold: %d   Levels: %d   Wavelet: %s   Scan: %d\n\n", 
		Threshold, NScale, FilterName[WaveletIndex], ScanOrder);
	
	/* - debug
	printf("Analysing data.\n");

	nbin=log2(Threshold)+2;
	bin = (int*)malloc(sizeof(int)*nbin);

	for (i=0; i<nbin; i++){
		bin[i]=0;
	}

	for (i=0; i<transform->vsize*transform->hsize; i++){
		for (j=0; j<nbin; j++){
			if (fabs(transform->value[i]) < (1<<j)){
				bin[j]++;
				break;
			}
		}
	}

	printf("Thres   Number of Coeffs\n");
	total = 0;
	for (i=0; i<nbin; i++){
		printf("%5d %8d\n", 1<<i, bin[i]);
		total += bin[i];
	}
	printf("Total number of coefficients: %d\n", total);
	*/

	/* initialize encode stream */
	EncodeStream = OpenOutputBitFile(EncodeImageFileName);
	encoder = ArithEncoderAlloc(EncodeStream);
	ArithEncoderStart(encoder);
	
	/* write header */
	BasicCoderWriteNBits(encoder, InputImage->xsize, ImageXSizeBits);					
	BasicCoderWriteNBits(encoder, InputImage->ysize, ImageYSizeBits);			
	BasicCoderWriteNBits(encoder, ScanOrder, ScanOrderBits);
	BasicCoderWriteNBits(encoder, NScale, ScaleBits);
	BasicCoderWriteNBits(encoder, WaveletIndex, WaveletIndexBits);
	BasicCoderWriteNBits(encoder, *((unsigned int *)&ImageMean), ImageMeanBits);	
	BasicCoderWriteNBits(encoder, log2(Threshold), ThresholdBits);					
	BasicCoderWriteNBits(encoder, MaxHistogramCount, MaxHistoCountBits);									

	/* EZW coding */
	EZW2Encode(transform, encoder, Threshold, MaxHistogramCount, 
		ByteBudget, ScanOrder);

	/* Clean up */
	ArithEncoderDone(encoder);
	ArithEncoderDealloc(encoder);
	CloseOutputBitFile(EncodeStream);
	FImageFree(InputImage);
	WaveletTransformDealloc(transform);
	WaveletDealloc(wavelet);
	RemoveFilterSets();
	fprintf(stderr, "Done\n");

	PrintLeaks();

	return;
}

/*----------------------------------------------------------------------------*/	
/*----------------------------------------------------------------------------*/	
void EZW2Encode(WTRANSFORM *transform, ArithEncoder *encoder, 
					int threshold, int HistoCount, int ByteBudget, int ScanOrder)
{
	//定义4个主编码扫描和一个从属编码扫描
	BasicCoder *DominantPassCoder[4];
	BasicCoder *SubordinateListCoder;
	int i, nPass=0;
	int BitsCount, TempCount;
	MAP *map;
	SubbandSymbol *ztmap;
	DLNode *DomList, *DomListTemp1, *DomListTemp2;
	SLNode *SubList;
	int SubListSize;
	
	//建立图像小波系数的零树图
	if ((ztmap=SubbandSymbolAlloc(transform->hsize,
					transform->vsize, transform->nsteps)) == NULL)
	{
		EZW2Error("Fail to allocate subband symbol.\n");
	}	
	EZW2BuildZerotreeMap(transform, ztmap);
	
	//建立主编码扫描的小波系数表
	if ((DomList = (DLNode *)malloc(sizeof(DLNode)
						*transform->vsize*transform->hsize))==NULL)
	{
		EZW2Error("Fail to allocate dominant list.\n");
	}
		
	//建立从属编码扫描的小波系数表
	if ((SubList = (SLNode *)malloc(sizeof(SLNode)
						*transform->vsize*transform->hsize))==NULL)
	{
		EZW2Error("Fail to allocate subordinate list.\n");
	}

	//在主编码扫描未完成的时候，从属过程的系数表为空
	SubListSize = 0;
	
	//主编码扫描过程建立4张系数表
	for (i=0; i<4; i++)
	{
		if ((DominantPassCoder[i] = BasicCoderAlloc(4, HistoCount)) == NULL)
		{
			EZW2Error("Fail to allocate dominanat pass coder.\n");
		}
	}

	//建立一个从属过程的小波系数上下文表
	if ((SubordinateListCoder = BasicCoderAlloc(2, HistoCount)) == NULL){
		EZW2Error("Fail to allocate subordinate list coder.\n");
	}

	//记录编码为POS和NEG的小波系数信息
	if ((map=MapAlloc(transform->hsize, transform->vsize, 
			transform->nsteps))==NULL){
		EZW2Error("Fail to allocate map.\n");
	}

	//统计目前为止编码器的输出数据流量
	BitsCount = ArithEncoderNBitsOutput(encoder);

	while((threshold!=0) && (EndEncoding!=TRUE))
	{
		printf("Pass #%2d - ", ++nPass);
		
		if (ScanOrder==0){
			//主编码扫描过程
			EncodeDominantPassMortonScan(DomList, SubList, &SubListSize, 
				transform, map, threshold, encoder,
				DominantPassCoder, ztmap, ByteBudget);
		}
		else{
			//否则，空操作
		}
		//统计当前编码输出流量
		TempCount = ArithEncoderNBitsOutput(encoder);
		//从属编码过程
		EncodeSubordinatePass(SubList, SubListSize, threshold>>1, encoder, 
							SubordinateListCoder, ByteBudget);
		//统计输出流量
		BitsCount = ArithEncoderNBitsOutput(encoder);

		printf("Sign Coeff: %6d  (%5d)  Bytes: %-6d (SB:%d)\n",
			SubListSize, threshold, ROUND((BitsCount/8.0)), BitsCount-TempCount);
		//完成编码后，将编码器复位
		if (EndEncoding != TRUE){
			threshold >>= 1;
			for (i=0; i<4; i++){
				BasicCoderReset(DominantPassCoder[i]);
			}
			BasicCoderReset(SubordinateListCoder);
		}
	}