main.c

用C++语言实现的基于小波分析的源代码,实现了小波分析的诸多算法

/******
todos :

	. lots of <> marks scattered around in the many modules

	. at "stop" functionality to ZT coders by including a "stop_subband"
		and setting the ZT "bottom" flag in that subband.  
		This would be something like the !propZTs of golfz

	. lifting transform

	. quantizers.c is poorly estimating distortion & entropy
		chooseQuantizers needs lots of work
			(stopAtQuantizer makes this a low-priority problem)
	. DCT works great with uniform quantizers, but wavelets (especially the S+P)
		are killed by large quantization of the top LL band
	. my old adaptive quantizers idea

	. the ultimate goal is for the encoders to have tune-able parameters
		that are optimized *per band* (such as context sizes).  The
		sever form of this is to actually choose different coders for
		different bands, but that is more tricky, because it is critical
		that coders are given passes to adapt; (this could be done by
		manually conditioning the new coder each time there's a change,
		but that would be slow for the decoder)

*********/

//#define LL_WAVE_CNTX

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <crblib/inc.h>
#include <crblib/timer.h>
#include <crblib/cindcatr.h>
#include <crblib/fileutil.h>
#include <wtlib/image.h>
#include <wtlib/coder.h>
#include <wtlib/codeimage.h>
#include <wtlib/quantim.h>
#include <crblib/tsc.h>
#include <wtlib/subbands.h>
#include <wtlib/wpacket.h>

#define  readImageFile	readImageFileBytes
#define writeImageFile writeImageFileBytes
#define 	dataBytes	1

/** some globals set by the mainline, for reference by the coders ;
/*** the various coder templates available : ***/

extern int tune_param;

const static int default_lossless_coder = 11; // an order-1

int main(int argc,char *argv[])
{
image *im=NULL,*save=NULL,*im_work=NULL;	// allocs
wavelet *wavelet = NULL;
char *sourceName;
int i,qtype,qflags;
int complen,stop_comp_len;
int width,height,planes,ll_width,ll_height;
int levels,coderN,losslessN,transformN;
double quantizer,stop_comp_rate;
bool writeWave,writeDec,doLossless,writeDeltas,doQuant,doYUV,interleaved,do_pause,oldSchool,wavePacket;
int autotune_min,autotune_max,autotune_step;
const coder * coder_template;
double forw_time,enc_time,dec_time,inv_time;

	errputs("WaveCoder v1.7, (c) 1998 by cBloom");

	if ( argc < 4 ) {
		errputs("wavecoder : tree-structured transforms and coders");
		errputs("usage: wt <source> <width> <height> -[flag] ...");
		errputs("flags:");
		errputs("	o	: old Skool , yo, no subband waveletpackets");
		errputs("	c#	: choose coder #");
		errputs("	l#	: choose # of levels (default 4)");
		errputs("	p#	: set # of planes > 1");
		errputs("	x#	: set tune param");
		errputs("	a#,#,# : auto-tune param <min,max,step>");
		errputs("	s#	: stop at bitrate #");
		errputs("	y 	: do YUV (also sets planes to 3)");
		errputs("	i 	: interleaved multi-planes");
		errputs("	ll	: code deltas to make lossless");
		errputs("	ww	: write wavelet coefficients to wave.raw");
		errputs("	wr	: write decompressed image to dec.raw");
		errputs("	wd	: write deltas to deltas.raw");
		errputs("	td	: use tree-structured dct");
		errputs("	ts	: use S+P transform (lifting)");
		errputs("	th	: use Haar transform (lifting)");
		errputs("	tw	: Daubechies 9/7 wavelet transform [default]");
		errputs("	tc2	: CDF (2,2) wavelet transform (lifting)");
		errputs("	tc4	: CDF (2,4) wavelet transform (lifting)");
		errputs("	tb	: BCW3 wavelet transform (lifting)");
		errputs("	t4	: D4 wavelet transform (lifting)");
		errputs("	t9f	: Floating (9,7) wavelet transform");
		errputs("	t9b	: Binary (9,7) wavelet transform");
		errputs("	t9l	: Lifting (9,7) wavelet transform (lifting)");
		errputs("	qu#	: quantize uniformly by (floating) factor");
		errputs("	qs# : quantize to target rate # & truncate to exact");
		errputs("	qf  : fiddle with quantizers");
		errputs("	qt# : set quantizer type");
		errputs("	g	: getch() at end (pause)");
		errputs("----------");
		fprintf(stderr,"(press a key)\r"); fflush(stderr); getchar();
		errputs("coders:                    ");
		for(i=0;i<num_coders;i++) fprintf(stderr,"%d : %s\n",i,coder_list[i]->name);
		errputs("");
		exit(0);
	}

	sourceName = argv[1];
	width = atol(argv[2]);
	height = atol(argv[3]);

	planes = 1;	levels = 4;	
	transformN = TRANS_DWT; coderN = -1; losslessN = default_lossless_coder;
	interleaved = doYUV = do_pause = doQuant = doLossless = writeWave = writeDec = writeDeltas = oldSchool = false;
	wavePacket = true;
	qtype = qflags = 0; quantizer = 1.0;
	stop_comp_len = LONG_MAX; stop_comp_rate = 8.0;
	autotune_step = 0;

	for(i=4;i<argc;i++) {
		char *s = argv[i];
		if ( *s == '-' || *s == '/' ) s++;
		s++;
		switch(toupper(s[-1])) {
			case 'O':
				oldSchool = true;
				errputs("got option: Old School : no wavelet packets or subband structs");
				break;
			case 'G':
				do_pause = true;
				errputs("got option: pause");
				break;
			case 'P':
				planes = atol(s); if ( planes < 1 ) planes = 1;
				fprintf(stderr,"got option: planes = %d\n",planes);
				break;
			case 'Y':
				doYUV = true; planes = 3;
				fprintf(stderr,"got option: do 3-planes with YUV\n");
				break;
			case 'I':
				interleaved = true;
				fprintf(stderr,"got option: interleaved planes\n");
				break;
			case 'X':
				tune_param = atol(s);
				fprintf(stderr,"got option: tune_param = %d\n",tune_param);
				break;
			case 'A':
				autotune_min = strtol(s,&s,10); s++;
				autotune_max = strtol(s,&s,10); s++;
				autotune_step = strtol(s,&s,10);
				fprintf(stderr,"got option: autotune : min = %d, max = %d, step = %d\n",autotune_min,autotune_max,autotune_step);
				break;
			case 'S':
				stop_comp_rate = atof(s); stop_comp_len = 0;
				fprintf(stderr,"got option: stop rate = %1.3f\n",stop_comp_rate);
				break;
			case 'L':
				if ( *s >= '0' && *s <= '9' ) {
					levels = atol(s); if ( levels < 0 ) levels = 0;
					fprintf(stderr,"got option: levels = %d\n",levels);
				} else {
					s++;
					switch(toupper(s[-1])) { 
						case 'L':
							doLossless = true;
							losslessN = atol(s);
							fprintf(stderr,"got option: do lossless, %d = %s\n",losslessN,coder_list[losslessN]->name);
							break;
						default:
							fprintf(stderr,"unknown option: %s\n",s-1);
							break;
					}
				}
				break;
			case 'C':
				coderN = atol(s);
				if ( coderN < num_coders && coderN >= 0 ) {
					fprintf(stderr,"got option: coderN = %d = %s\n",coderN,coder_list[coderN]->name);
				} else errputs("invalid coder number on command line\n");
				break;
			case 'T': 
				s++;
				switch(toupper(s[-1])) {
					case 'D':
						transformN = TRANS_DCT;
						fprintf(stderr,"got option: use DCT\n");
						break;
					case 'S':
						transformN = TRANS_SPT;
						fprintf(stderr,"got option: use S+P transform\n");
						break;
					case 'C':
						if ( atol(s) == 4 ) {
							transformN = TRANS_CDF24;
							fprintf(stderr,"got option: use CDF (2,4) transform\n");
						} else {
							transformN = TRANS_CDF22;
							fprintf(stderr,"got option: use CDF (2,2) transform\n");
						}
						break;
					case '9':
						switch(toupper(s[0])) {
							case 'F':
								transformN = TRANS_F97;
								fprintf(stderr,"got option: use F (9,7) transform\n");
								break;
							case 'B':
								transformN = TRANS_B97;
								fprintf(stderr,"got option: use B (9,7) transform\n");
								break;
							case 'L':
							default:
								transformN = TRANS_L97;
								fprintf(stderr,"got option: use L (9,7) transform\n");
								break;
						}
						break;
					case 'B':
						transformN = TRANS_BCW3;
						fprintf(stderr,"got option: use BCW3 transform\n");
						break;
					case '4':
						transformN = TRANS_D4;
						fprintf(stderr,"got option: use D4 transform\n");
						break;
					case 'H':
						transformN = TRANS_HAAR;
						fprintf(stderr,"got option: use Haar (S) transform\n");
						break;
					case 'W':
						transformN = TRANS_DWT;
						fprintf(stderr,"got option: use Daub Wavelet transform\n");
						break;
					default:
						fprintf(stderr,"unknown option: %s\n",s-1);
						break;
				}
				break;
			case 'W': s++;
				switch(toupper(s[-1])) {
					case 'W':
						writeWave = true;
						fprintf(stderr,"got option: write Wave to wave.raw\n");
						break;
					case 'R':
						writeDec = true;
						fprintf(stderr,"got option: write Raw Output to dec.raw\n");
						break;
					case 'D':
						writeDeltas = true;
						if ( ! doLossless ) { doLossless = true; losslessN = default_lossless_coder; }
						fprintf(stderr,"got option: write deltas to deltas.raw\n");
						break;
					default:
						fprintf(stderr,"unknown option: %s\n",s-1);
						break;
				}
				break;
			case 'Q':
				s++;
				switch(toupper(s[-1])) {
					case 'U':
						quantizer = atof(s); doQuant = true;
						fprintf(stderr,"got option: quantizer = %f\n",quantizer);
						break;
					case 'S': case 'R':
						stop_comp_rate = atof(s); stop_comp_len = 0; doQuant = true;
						qflags |= QFLAG_SEEKSTOP;
						fprintf(stderr,"got option: rate target stopping = %f\n",stop_comp_rate);
						break;
					case 'F':
						doQuant = true;
						qflags |= QFLAG_FIDDLE;
						errputs("got option: fiddle quantizers");
						break;
					case 'T':
						qtype = atol(s); doQuant = true;
						if ( qtype > QTYPE_MAX ) qtype = QTYPE_MAX;
						fprintf(stderr,"got option: quantizer type = %d = %s\n",qtype,qtype_names[qtype]);
						break;
					default:
						fprintf(stderr,"unknown option: %s\n",s-1);
						break;
				}
				break;
			default:
				fprintf(stderr,"unknown option: %s\n",s-1);
				break;
		}
	}

	/** fix up options received **/

	if (! (transformN & TRANS_PACKETS) ) oldSchool = true;
	if ( oldSchool ) wavePacket = false;

	if ( transformN == TRANS_DCT ) levels = 3;

	if ( coderN >= num_coders || coderN < 0 ) {
		fprintf(stderr,"using default coder : %s\n",coder_list[0]->name);
		coderN = 0;
	}
	coder_template = coder_list[coderN];

	if ( coder_template->encodeBandZT ) { wavePacket = false; }
	if ( coder_template->encodeSubbandBP ) { oldSchool = false; }

	if ( wavePacket ) errputs("using Wavelet-Packets");
	else errputs("forced to use wavelets without packets");

	ll_width = width>>levels;
	ll_height = height>>levels;

	if ( ((ll_width>>1)<<(levels+1)) != width )	
		errexit("width must be a factor of levels+1");
	if ( ((ll_height>>1)<<(levels+1)) != height )	
		errexit("width must be a factor of levels+1");
	/** <> don't abort, just make the extension **/

	/** start allocs and go to it **/

	if ( (im = newImage(width,height,planes)) == NULL )
		errexit("newImage failed");
	if ( (im_work = newImageFrom(im)) == NULL )
		errexit("newImage failed");

	if ( ! readImageFile(sourceName,im,interleaved) )
		errexit("readImage failed");

	if ( (save = copyImage(im)) == NULL )
		errexit("copyimage failed");

	stop_comp_len = (int)((dataBytes * stop_comp_rate * im->tot_size + 4) / 8);

	errputs("forward transform...");

	pushTSC();

	if ( doYUV ) RGBtoYUV(im);

	if ( oldSchool ) {
		wavelet = makeTransQuantWavelet(im,levels,transformN,
			qtype,qflags,quantizer,
			stop_comp_len,coder_template);
	} else {
		wavelet = makeSubband(im,levels,quantizer,
			stop_comp_len,coder_template,
			wavePacket,transformN);
	}

	forw_time = popTSC();

	/** im now contains the quantized transformed data **/

	if ( writeWave ) {
		patchImage(im,im_work,0,0,width,height,0,0);

		if ( ! writeImageFile("wave.raw",im_work,interleaved) )
			errexit("writeImage failed");
	}

	if ( autotune_step != 0 ) {
		int best_size,best_param,size;

		errputs("doing the autotune:");

		best_size = LONG_MAX; best_param = autotune_min;
		for(tune_param=autotune_min;tune_param<=autotune_max;tune_param += autotune_step) {
			size = packedImageSize(im,coder_template,levels);
			fprintf(stderr,"%d : %d\n",tune_param,size);
			if ( size < best_size ) {
				best_size = size;
				best_param = tune_param;
			}
		}
		tune_param = best_param;
		printf("autotuned to : param = %d at size = %d\n",tune_param,best_size);
	}

	errputs("encoding...");	/********************************/

	pushTSC();

	if ( oldSchool ) {
		encodeWaveletAndLL(im,wavelet,coder_template,stop_comp_len);
	} else {
		encodeWaveletSubbands(wavelet);
	}

	enc_time = popTSC();

	complen = wavelet->complen;

	errputs("decoding...");	/********************************/

		/** secret side information is now kept in the wavelet structure
		**/

	zeroImage(im);	
		/** must do this for decoders that terminate early ; 
		***		also makes sure I don't cheat 
		**/			

	pushTSC();

	if ( oldSchool ) {
		decodeWaveletAndLL(im,wavelet);
	} else {
		decodeWaveletSubbands(wavelet);
	}

	dec_time = popTSC();

	errputs("reverse transform..."); 	/********************************/

	if ( doLossless ) {
		/** save the quantized wavelet in im_work **/
		patchImage(im,im_work,0,0,width,height,0,0);
	}

	pushTSC();

	if ( oldSchool ) {
		deTransQuantWavelet(wavelet);
	} else {
		unSubband(wavelet);
	}

	if ( doYUV ) YUVtoRGB(im);

	inv_time = popTSC();

	/********************************/

	printf("%-18s:",FilePart(sourceName)); 
	TheCompressionIndicator(dataBytes*im->tot_size,complen,stdout); puts("");

	imageCompare(save,im,stdout);

	showTSC("forward",stdout,forw_time);
	showTSC("encode",stdout,enc_time);
	showTSC("decode",stdout,dec_time);
	showTSC("inverse",stdout,inv_time);
	puts("");

	if ( writeDec )
		if ( ! writeImageFile("dec.raw",im,interleaved) ) errexit("writeImage failed");

	if ( doLossless ) {
		int complen_deltas;
		coder *encoder;

		/** im contains the reconstruction, save the original, and im_work the wavelet ***/

		if ( (encoder = coder_create_write(coder_list[losslessN],wavelet,LONG_MAX)) == NULL )
			errexit("coder_create deltas failed!");

		encoder->init(encoder);

		subtractImage(im,save);

		encodeDeltaIm(encoder,im,im_work);

		coder_flush_write(encoder);	
		encoder->free(encoder);
		coder_destroy(encoder); 

		if ( writeDeltas ) {
			if ( ! writeImageFile("deltas.raw",im,interleaved) )
				errexit("writeImage failed");
			if ( ! writeImageFile("del_cntx.raw",im_work,interleaved) )
				errexit("writeImage failed");
		}

		complen_deltas = wavelet->complen;

		printf("%-18s:","deltas"); 
		TheCompressionIndicator(dataBytes*im->tot_size,complen_deltas,stdout); puts("");
		printf("%-18s:","lossless"); 
		TheCompressionIndicator(dataBytes*im->tot_size,complen_deltas+complen,stdout); puts("");
	}

	freeWavelet(wavelet);
	freeImage(im); im = NULL;
	freeImage(im_work); im_work = NULL;
	freeImage(save); save = NULL;

	if ( do_pause ) { errputs("press a key to close WT"); getchar(); }

return 0;
}