subbands.c

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

/**********

the quantim.c version of quantizers per band is archaic; we no longer
have an gauranteed band structure.  Any quantizer who even uses "levels"
is archaic.

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

target structure :

	image -> [quant & trans] -> (subbands + wavelet struct + quantinfo)
	subbands -> [sender]
		sends activities & bitmap of the subband tree structure
		sorts by activity
		sets optimal cntx pointers
	[sender] -> subband_leafs to coders
		coders which care about contexts must look at cntx widths now

----------

choosing cntx1 & cntx2 is not trivial ; for one thing, the end-coder has
some preferences.  The great "bpb" and "bpb2" coders want sister contexts
(same size as current); in this case we should use the two similar-activity
contexts which are closest *and of the same size*

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

#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <crblib/inc.h>
#include <crblib/intmath.h>

#include "image.h"
#include "subbands.h"
#include "coder.h"
#include "codezt.h"
#include <crblib/arithc.h>
#include <crblib/codeutil.h>

void freeSubbands(subband *root);

void arithRLEinit(void);
void arithBitRLE(arithInfo * ari,bool bit);
bool arithGetBitRLE(arithInfo * ari);

subband * imageToSubbands(image *im,int levels);
subband * makeSubbandQuad(int *band,int w,int h,int fullw,int levels,subband_node *higher);

subband_node *new_subband_node();
subband_node *new_subband_node_foliated();
subband_leaf *new_subband_leaf();

void encodeSubbands(coder * encoder,subband * root);
void decodeSubbands(coder * encoder,subband * root);

void encodeWaveletSubbands(wavelet *wavelet);
void decodeWaveletSubbands(wavelet *wavelet);

void encodeLowestSubband(coder *encoder,subband_leaf *sbl);
void decodeLowestSubband(coder *decoder,subband_leaf *sbl);

void activitySubbands(subband * sb);
void getActivities(coder * decoder,subband * sb);
void sendActivities(coder * encoder,subband * sb);
void initActivityRLE(void);

subband_leaf * LLleaf(subband *sb);
subband_leaf * seekLeaf(subband_leaf **leaves,int num,int w,int type);
subband_leaf * seekLeafByType(subband_leaf **leaves,int num,int w,int type);
subband * seekNode(subband_leaf **leaves,int num,int w,int type);
subband * seekNodeByType(subband_leaf **leaves,int num,int w,int type);

#define TOP_BITMASK (CODE_MAX_VAL>>1)

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

/** ! we send the LL band here ! **/

void encodeWaveletSubbands(wavelet *wavelet)
{
coder *encoder;
 
	if ( (encoder = coder_create_write(wavelet->coder_template,wavelet,wavelet->stoplen)) == NULL )
		errexit("coder_create_write failed!");

	encoder->init(encoder);

	if ( encoder->encodeBandZT ) {
		int w,h,p;
		w= (wavelet->im->width)>>(wavelet->levels);
		h= (wavelet->im->height)>>(wavelet->levels);
		for(p=0;p<(wavelet->im->planes);p++)
			coder_encodeDPCM(encoder,wavelet->im->data[p][0],w,h,(wavelet->im->width)-w);

		encodeImageZT(encoder,wavelet->im,wavelet->levels);
	} else {
//		sendBandMap(encoder,wavelet->subband_root); , easy
		encodeSubbands(encoder,wavelet->subband_root);
	}

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

void decodeWaveletSubbands(wavelet *wavelet)
{
coder *decoder;

	if ( (decoder = coder_create_read(wavelet)) == NULL )
		errexit("decoder_create main failed!");

	decoder->init(decoder);

	if ( decoder->decodeBandZT ) {
		int w,h,p;
		w= (wavelet->im->width)>>(wavelet->levels);
		h= (wavelet->im->height)>>(wavelet->levels);
		for(p=0;p<(wavelet->im->planes);p++)
			coder_decodeDPCM(decoder,wavelet->im->data[p][0],w,h,(wavelet->im->width)-w);

		decodeImageZT(decoder,wavelet->im,wavelet->levels);
	} else {
//		if ( wavelet->subband_root ) freeSubbands(wavelet->subband_root);
//		wavelet->subband_root = getBandMapImage(decoder,im);
		decodeSubbands(decoder,wavelet->subband_root);
	}

	decoder->free(decoder);
	coder_flush_read(decoder);
	coder_destroy(decoder);
}

void encodeSubbandLeaf(coder * encoder,subband_leaf * sbl)
{
	if ( !sbl ) return;
	if ( coder_timetostop(encoder) ) return;

	if ( sbl->parent == NULL ) {
		encodeLowestSubband(encoder,sbl);
	} else {
		encoder->encodeBand(encoder,sbl->band,sbl->width,sbl->height,
			sbl->width + sbl->rowpad,sbl->parent->band);
	}
}

void decodeSubbandLeaf(coder * decoder,subband_leaf * sbl)
{
	if ( !sbl ) return;
	if ( coder_timetostopd(decoder,0) ) return;

	if ( sbl->parent == NULL ) {
		decodeLowestSubband(decoder,sbl);
	} else {
		decoder->decodeBand(decoder,sbl->band,sbl->width,sbl->height,
			sbl->width + sbl->rowpad,sbl->parent->band);
	}
}

void encodeSubbandLeafBP(coder * encoder,subband_leaf * sbl,int bitmask)
{
	if ( !sbl ) return;
	if ( coder_timetostop(encoder) ) return;

	if ( sbl->parent == NULL ) {
		if ( bitmask == TOP_BITMASK )
			encodeLowestSubband(encoder,sbl);
		return;
	}

	if ( bitmask > sbl->max_value ) {
		arithBitRLE(encoder->arith,0);
		return;
	} else {
		arithBitRLE(encoder->arith,1);
	}

	if ( encoder->encodeBandBP ) {
		encoder->encodeBandBP(encoder,sbl->band,sbl->width,sbl->height,
			sbl->width + sbl->rowpad,sbl->parent->band,bitmask);
	} else {
		encoder->encodeSubbandBP(encoder,sbl,bitmask);
	}
}

void decodeSubbandLeafBP(coder * decoder,subband_leaf * sbl,int bitmask)
{
	if ( !sbl ) return;
	if ( coder_timetostopd(decoder,0) ) return;

	if ( sbl->parent == NULL ) {
		if ( bitmask == TOP_BITMASK )
			decodeLowestSubband(decoder,sbl);
		return;
	}

	if ( ! arithGetBitRLE(decoder->arith) )
		return;

	if ( decoder->decodeBandBP ) {
		decoder->decodeBandBP(decoder,sbl->band,sbl->width,sbl->height,
			sbl->width + sbl->rowpad,sbl->parent->band,bitmask);
	} else {
		decoder->decodeSubbandBP(decoder,sbl,bitmask);
	}
}

void encodeLowestSubband(coder *encoder,subband_leaf *sbl)
{
	coder_encodeDPCM(encoder,sbl->band,sbl->width,sbl->height,sbl->rowpad);
}

void decodeLowestSubband(coder *decoder,subband_leaf *sbl)
{
	coder_decodeDPCM(decoder,sbl->band,sbl->width,sbl->height,sbl->rowpad);
}

int countLeaves(subband *sb)
{
	if ( ! sb ) return 0;
	if ( sb->leaf ) return 1;
	else {
		subband_node * sbn;
		sbn = sb;
		return countLeaves(sbn->LL) +countLeaves(sbn->LH) +countLeaves(sbn->HL) +countLeaves(sbn->HH);
	}
}

void readInLeaves(subband *sb,subband_leaf ** leaves,int *index)
{
	if ( !sb ) return;
	if ( sb->leaf ) {
		leaves[ *index ] = (subband_leaf *) sb;
		(*index) ++;
		return;
	} else {
		subband_node * sbn = sb;
		readInLeaves(sbn->LL,leaves,index);
		readInLeaves(sbn->LH,leaves,index);
		readInLeaves(sbn->HL,leaves,index);
		readInLeaves(sbn->HH,leaves,index);
	}
}

void setLeaveContexts(subband_leaf ** leaves,int nLeaves)
{
int i,t,j;
subband_leaf *cur,*c1,*c2;

	for(i=3;i<nLeaves;i++) {
		cur = leaves[i]; 
		for(j = i-1; j>= 0;j--) {
			c1  = leaves[j];
			if ( c1->width == cur->width )
				break;
		}
		if ( j < 0 ) c1 = leaves[i-1];
		for(j = j-1; j>= 0;j--) {
			c2  = leaves[j];
			if ( c2->width == cur->width )
				break;
		}
		if ( j < 0 ) { c2 = leaves[i-2];
			if ( c2 == c1 ) c2 = leaves[i-1];
		}

		cur->cntx1 = c1;	cur->cntx2 = c2;
		cur->cntx1shift = cur->cntx2shift = 0;

		for(t=cur->width; t > c1->width ; t>>= 1) cur->cntx1shift++;
		for(t=cur->width; t < c1->width ; t<<= 1) cur->cntx1shift--;
		for(t=cur->width; t > c2->width ; t>>= 1) cur->cntx2shift++;
		for(t=cur->width; t < c2->width ; t<<= 1) cur->cntx2shift--;
	}	
}

/*}{********* sort & find parents *******/
#if 1

void sortLeaves(subband_leaf ** leaves,int nLeaves)
{
subband_leaf *cur;
int i,p,j;
bool didStuff;

	// the idea here is keep moving a leaf up the list while its
	//	activity is higher than its predecessor, and its predecessor isnt its parent
	//
	// btw we should really use D/R for the activity

	do {
		didStuff = false;
		for(i=0;i<nLeaves;i++) {
			cur = leaves[i];
			if ( cur->parent == NULL ) p = 0;
			else {
				p = 0;
				while ( leaves[p] != cur->parent ) p++;
				p++;
				if ( p >= i ) continue;
			}
			// try to swap cur up to p
			while(p < i ) {
				if ( cur->activity > leaves[p]->activity ) {
					// do it
					for (j=i-1;j>=p;j--) {
						leaves[j+1] = leaves[j];
					}
					leaves[p] = cur;
					didStuff = true;
					break;
				}
				p++;
			}
		}
	} while( didStuff);

	setLeaveContexts(leaves,nLeaves);
}

#else
/*}{********* sort & find parents : sort-favoring way *******/

int sortLeavesCompare(const void *p1,const void *p2)
{
subband_leaf *l1,*l2;
	l1 = *((subband_leaf **)p1);
	l2 = *((subband_leaf **)p2);
	if ( l2->width != l1->width ) 
		return (l1->width) - (l2->width);
return (l2->activity) - (l1->activity);
}

void sortLeaves(subband_leaf ** leaves,int nLeaves)
{
subband_leaf *cur;
int i,w,type;

	// sort by activity

	qsort(leaves,nLeaves,sizeofpointer,sortLeavesCompare);

	/* now pick up parents	(unlike contexts, we will move
	*	leaves around in order to get good parents
	*
	*	try to find a parent of half the size, 
	*		and the same type
	*	we can use a node, instead of a leaf, as the parent
	*		in this case, all children of the node must be
	*		sent before us.
	*
	*/

	for(i=0;i<nLeaves;i++) {
		cur = leaves[i];	cur->parent = NULL;
		type = cur->type; w = (cur->width)>>1;

		if ( cur->parent = seekLeafByType(leaves,i,(w<<0),type) ) goto found;
		if ( cur->parent = seekLeafByType(leaves,i,(w<<1),type) ) goto found;
		if ( cur->parent = seekLeafByType(leaves,i,(w<<2),type) ) goto found;
		if ( cur->parent = seekLeafByType(leaves,i,(w<<3),type) ) goto found;

		if ( cur->parent = seekNodeByType(leaves,i,(w<<0),type) ) goto found;
		if ( cur->parent = seekNodeByType(leaves,i,(w<<1),type) ) goto found;
		if ( cur->parent = seekNodeByType(leaves,i,(w<<2),type) ) goto found;
		if ( cur->parent = seekNodeByType(leaves,i,(w<<3),type) ) goto found;

		// not found any!
		continue;

		found:
		cur->parent = LLleaf(cur->parent);
		// <> this is kind of sloppy, and means the parent won't be
		//	quite right
	}

	setLeaveContexts(leaves,nLeaves);
}

#endif
/*}{******* sort & find parents *************/

void sendBandMap(coder *encoder,subband *sb)
{
	if ( sb->leaf ) {
		arithBit(encoder->arith,0);
	} else {
		arithBit(encoder->arith,1);
		sendBandMap(encoder,sb->LL);
		sendBandMap(encoder,sb->LH);
		sendBandMap(encoder,sb->HL);
		sendBandMap(encoder,sb->HH);
	}
}

void encodeSubbands(coder * encoder,subband * root)
{
subband_leaf ** leaves;
int nLeaves,i;

	initActivityRLE();
	sendActivities(encoder,root);

	encoder->w->stopline = - 1;

	nLeaves = countLeaves(root);

	if ( (leaves = malloc(sizeofpointer*nLeaves)) == NULL )
		errexit("can't alloc leaves");

	i = 0;
	readInLeaves(root,leaves,&i);

	/**
	*	sorting seems to be a nice little win
	**/

	sortLeaves(leaves,nLeaves);

// printSubbandInfo(root,stdout,0);

	arithRLEinit();

	if ( encoder->encodeBand ) {
		for(i=0;i<nLeaves;i++) {
			encodeSubbandLeaf(encoder,leaves[i]);
		}
	} else if ( encoder->encodeBandBP || encoder->encodeSubbandBP ) {
		int bitmask;

		for(bitmask = TOP_BITMASK;bitmask>=1;bitmask>>=1) {
			for(i=0;i<nLeaves;i++) {
				encodeSubbandLeafBP(encoder,leaves[i],bitmask);
			}
		}

	} else {
		errexit("only encodeBand() & BP() can do subbands right now");
	}

	if ( encoder->w->stopline < 0 ) encoder->w->stopline = INT_MAX;

	free(leaves);
}

void decodeSubbands(coder * decoder,subband * root)
{
subband_leaf ** leaves;
int nLeaves,i;

	initActivityRLE();
	getActivities(decoder,root);

	nLeaves = countLeaves(root);

	if ( (leaves = malloc(sizeofpointer*nLeaves)) == NULL )
		errexit("can't alloc leaves");

	i = 0;
	readInLeaves(root,leaves,&i);

	sortLeaves(leaves,nLeaves);

	arithRLEinit();

	if ( decoder->decodeBand ) {
		for(i=0;i<nLeaves;i++) {
			decodeSubbandLeaf(decoder,leaves[i]);
		}
	} else if ( decoder->decodeBandBP || decoder->decodeSubbandBP ) {
		int bitmask;

		for(bitmask = TOP_BITMASK;bitmask>=1;bitmask>>=1) {
			for(i=0;i<nLeaves;i++) {
				decodeSubbandLeafBP(decoder,leaves[i],bitmask);
			}
		}

	} else {
		errexit("only decodeBand() & BP() can do subbands right now");
	}

	free(leaves);