subbands.c

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

}

int subbandO0Entropy(subband_leaf *sbl)
{
int *bp,x,y,z,ret;
int cnts[256];

	memclear(cnts,sizeof(int)*256);

	bp = sbl->band;
	for(y=0; y<(sbl->height);y++) {
		for(x=0; x<(sbl->width);x++) {
			z = *bp++; z=abs(z); z = min(z,0xFF);
			cnts[z]++;
		}
		bp += sbl->rowpad;
	}
	
	ret = 0; y = (sbl->width) * (sbl->height);
	for(x=0;x<256;x++) {
		z = cnts[x];
		if ( z && x ) {
			ret += z * intlog2x16(y/z);
			ret += intlog2x16(z);	//bits to send that prob.
		}
	}
	ret /= y;

return ret;
}

int subbandDeltaO0Entropy(subband_leaf *sbl)
{
int *bp,x,y,z,fullw,ret,pred;
int cnts[256];

	memclear(cnts,sizeof(int)*256);

	bp = sbl->band; fullw = sbl->width + sbl->rowpad;
	for(y=0; y<(sbl->height);y++) {
		for(x=0; x<(sbl->width);x++) {
			z = *bp;
			if ( x && y ) pred = (bp[-1] + bp[-fullw])>>1;
			else if ( x ) pred = bp[-1];
			else if ( y ) pred = bp[-fullw];
			else pred = 0;
			bp++;
			z = abs(z - pred); z = min(z,0xFF);
			cnts[z]++;
		}
		bp += sbl->rowpad;
	}
	
	ret = 0; y = (sbl->width) * (sbl->height);
	for(x=0;x<256;x++) {
		z = cnts[x];
		if ( z && x ) {	// doesn't count cost of sending zeros
			ret += z * intlog2x16(y/z);
			ret += intlog2x16(z);	//bits to send that prob.
		}
	}
	ret /= y;

return ret;
}

/** <> the transform seems *awfully* sensitive to these
*		parameters; on "tile1" we go from 3.432 to 3.8 to 4.5 bpp !
***/

#define MAX_QUANT	2
#define MAX_SHIFT 	6
#define MAX 		(1<<MAX_SHIFT)
#define MAX_SQR		(1<<(MAX_SHIFT+MAX_SHIFT))
#define MAX_MASK	(MAX-1)

int subbandO1Entropy(subband_leaf *sbl)
{
int *bp,i,x,y,z,fullw,ret,pred;
int * cnts;
int tots[MAX];
int escs[MAX];

	if ( ! (cnts = newarray(int,MAX_SQR)) ) return 0;
	for(i=0;i<MAX;i++) escs[i] = tots[i] = 1;

	ret = 0;
	bp = sbl->band; fullw = sbl->width + sbl->rowpad;
	for(y=0; y<(sbl->height);y++) {
		for(x=0; x<(sbl->width);x++) {
			z = abs(*bp);
			if ( x && y ) pred = (abs(bp[-1]) + abs(bp[-fullw]))>>1;
			else if ( x ) pred = abs(bp[-1]);
			else if ( y ) pred = abs(bp[-fullw]);
			else pred = 0;
			bp++;
			z >>= MAX_QUANT; pred >>= MAX_QUANT;
			z = min(z,MAX_MASK);	pred = min(pred,MAX_MASK);

			i = z + (pred<<MAX_SHIFT);

			ret += intlog2x16(tots[pred]);
			if ( cnts[i] ) {
				ret -= intlog2x16(cnts[i]);
				if ( cnts[i] == 1 && escs[pred] > 1 ) {
					escs[pred]--; tots[pred]--;
				}
			} else {
				ret -= intlog2x16(escs[pred]);
				ret += intlog2x16(z+1) + 1;
				escs[pred]++; tots[pred]++;
			}
			
			cnts[i]++;
			tots[pred]++;
		}
		bp += sbl->rowpad;
	}

	ret /= (sbl->height)*(sbl->width);

return ret;
}

void activitySubbands(subband * sb)
{
	if ( !sb ) return;
	if ( sb->leaf ) {
		subband_leaf *sbl;
		int *bp,x,y,z;
		sbl = (subband_leaf *) sb;

		bp = sbl->band;
		sbl->max_value = sbl->activity = 0;
		for(y=0; y<(sbl->height);y++) {
			for(x=0; x<(sbl->width);x++) {
				z = *bp++; if ( z < 0 ) z = -z;
				if ( z > sbl->max_value ) sbl->max_value = z;
				if ( z ) sbl->activity ++;
			}
			bp += sbl->rowpad;
		}

#if 0
		sbl->activity = (sbl->activity << 6) / ((sbl->height)*(sbl->width));
#else
		sbl->activity = subbandO1Entropy(sbl);
#endif

		return;
	} {
		subband_node *sbn;
		int num;
		sbn = (subband_node *) sb;
		sbn->activity = sbn->max_value = 0; num = 0;
		if ( sbn->LL ) { activitySubbands(sbn->LL); sbn->activity += sbn->LL->activity; sbn->max_value = max(sbn->max_value,sbn->LL->max_value); num++; } 
		if ( sbn->LH ) { activitySubbands(sbn->LH); sbn->activity += sbn->LH->activity; sbn->max_value = max(sbn->max_value,sbn->LH->max_value); num++; } 
		if ( sbn->HL ) { activitySubbands(sbn->HL); sbn->activity += sbn->HL->activity; sbn->max_value = max(sbn->max_value,sbn->HL->max_value); num++; } 
		if ( sbn->HH ) { activitySubbands(sbn->HH); sbn->activity += sbn->HH->activity; sbn->max_value = max(sbn->max_value,sbn->HH->max_value); num++; }
		if ( num ) sbn->activity /= num;
	}
}

static int last_activity;

void initActivityRLE(void) {
	last_activity = 0;
}

void sendActivities(coder * encoder,subband * sb)
{
	if ( !sb ) return;
	if ( sb->leaf ) {
		subband_leaf *sbl;
		sbl = (subband_leaf *) sb;

		cu_putExpandingSigned_ari(sbl->activity - last_activity,encoder->arith,30,30);
		last_activity = sbl->activity;

		return;
	} {
		subband_node *sbn;
		sbn = (subband_node *) sb;
		if ( sbn->LL ) sendActivities(encoder,sbn->LL);
		if ( sbn->LH ) sendActivities(encoder,sbn->LH);
		if ( sbn->HL ) sendActivities(encoder,sbn->HL);
		if ( sbn->HH ) sendActivities(encoder,sbn->HH);
	}
}

void getActivities(coder * decoder,subband * sb)
{
	if ( !sb ) return;
	if ( sb->leaf ) {
		subband_leaf *sbl;
		sbl = (subband_leaf *) sb;

		sbl->activity = cu_getExpandingSigned_ari(decoder->arith,30,30) + last_activity;
		last_activity = sbl->activity;

		sbl->max_value = 0;

		return;
	} {
		subband_node *sbn;
		sbn = (subband_node *) sb;
		sbn->activity = 0;
		if ( sbn->LL ) { getActivities(decoder,sbn->LL); sbn->activity += sbn->LL->activity; }
		if ( sbn->LH ) { getActivities(decoder,sbn->LH); sbn->activity += sbn->LH->activity; }
		if ( sbn->HL ) { getActivities(decoder,sbn->HL); sbn->activity += sbn->HL->activity; }
		if ( sbn->HH ) { getActivities(decoder,sbn->HH); sbn->activity += sbn->HH->activity; }
		sbn->max_value = 0;
	}
}

void transposeSubbandLHs(subband *sb)
{
	if ( !sb || sb->leaf ) return;
	transposeSubbandLHs(sb->LL);
	transposeSubbandLHs(sb->LH);
	transposeSubbandLHs(sb->HL);
	transposeSubbandLHs(sb->HH);
	if ( sb->LH && sb->LH->leaf ) {
		subband_leaf * sbl = (subband_leaf *) sb->LH;
		int *pa,*pb,x,y,z,fullw;
		fullw = sbl->width + sbl->rowpad;
		for(y=0;y< (sbl->height);y++) {
			pa = pb = (sbl->band) + y*(fullw + 1);
			for(x=y+1;x< (sbl->width);x++) {
				pa ++;	pb += fullw;
				z	= *pa;
				*pa = *pb;
				*pb = z;
			}
		}
		sbl->transposed = true;
	}
}

subband * imageToSubbands(image *im,int levels)
{
int p;
	assert( im->planes <= 4 );
	p = im->planes;
	if ( p == 1 ) 
		return makeSubbandQuad(im->data[0][0],im->width,im->height,im->width,levels-1,NULL);
	else {
		subband_node *sbn;
		sbn = new_subband_node();
		if ( p >= 1 ) sbn->LL = makeSubbandQuad(im->data[0][0],im->width,im->height,im->width,levels-1,NULL);
		if ( p >= 2 ) sbn->LH = makeSubbandQuad(im->data[1][0],im->width,im->height,im->width,levels-1,NULL);
		if ( p >= 3 ) sbn->HL = makeSubbandQuad(im->data[2][0],im->width,im->height,im->width,levels-1,NULL);
		if ( p >= 4 ) sbn->HH = makeSubbandQuad(im->data[3][0],im->width,im->height,im->width,levels-1,NULL);
		return (subband *)sbn;
	}
}

subband * makeSubbandQuad(int *band,int w,int h,int fullw,int levels,subband_node *higher)
{
subband_leaf *LH,*HL,*HH;
subband_node *sbn;

	sbn = new_subband_node_foliated();	
	sbn->prev = higher;
	sbn->width = w;
	sbn->height = h;
    sbn->rowpad = fullw - w;

	w >>= 1; h >>= 1;

	HL = (subband_leaf *)sbn->HL;
	LH = (subband_leaf *)sbn->LH;
	HH = (subband_leaf *)sbn->HH;

    HL->width = w;
    HL->height =h;
    HL->rowpad = fullw - w;
    HL->band = band + w;

    LH->width = w;
    LH->height =h;
    LH->rowpad = fullw - w;
    LH->band = band + h*fullw;

	LH->transposed = true;	// we called transposeLHs earlier
			/* <> should be set when the subbands are built,
			*		when we do the actual transpose
			*	there's a trouble here for the decoder, because he must know
			*		"transposed?" before actually doing it; this bit must be sent,
			*		unless we hard-code it to the LH (which isn't terrible)
			*/

    HH->width = w;
    HH->height =h;
    HH->rowpad = fullw - w;
    HH->band = LH->band + w;

	if ( levels ) {
		sbn->LL = makeSubbandQuad(band,w,h,fullw,levels-1,sbn);
	} else {
		subband_leaf *sbl;
		// this is the LL band :

		sbl = new_subband_leaf();
		sbn->LL = (subband *) sbl;
		sbl->prev = sbn;
		sbl->width = w;
		sbl->height =h;
		sbl->rowpad = fullw - w;
		sbl->band = band;
	}

	sbn->LL->type = TYPE_LL;
	sbn->LH->type = TYPE_LH;
	sbn->HL->type = TYPE_HL;
	sbn->HH->type = TYPE_HH;

return sbn;
}

subband_node *new_subband_node()
{
subband_node *ret;
	if ( (ret = new(subband_node)) == NULL ) errexit("malloc subband failed");
	ret->leaf = false;
return ret;
}

subband_node *new_subband_node_foliated()
{
subband_node *ret;
	if ( (ret = new(subband_node)) == NULL ) errexit("malloc subband failed");
	ret->leaf = false;
	ret->LH = (subband *) new_subband_leaf();
	ret->LH->prev = ret;
	ret->HL = (subband *) new_subband_leaf();
	ret->HL->prev = ret;
	ret->HH = (subband *) new_subband_leaf();
	ret->HH->prev = ret;
return ret;
}
subband_leaf *new_subband_leaf()
{
subband_leaf *ret;
	if ( (ret = new(subband_leaf)) == NULL ) errexit("malloc subband failed");
	ret->leaf = true;
return ret;
}

void freeSubbands(subband *root)
{
	if ( ! root ) return;
	if ( ! root->leaf ) {
		subband_node *node = root;
		freeSubbands(node->LL);
		freeSubbands(node->LH);
		freeSubbands(node->HL);
		freeSubbands(node->HH);
	}
	free(root);
}

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

static bool rle_last_bit;

#define RLE_BIT_TOT	8192
#define RLE_BIT_MID	358

void arithRLEinit(void)
{
rle_last_bit = 0;
}

void arithBitRLE(arithInfo * ari,bool bit)
{
	bit = bit?1:0;
	arithEncBit(ari,RLE_BIT_MID,RLE_BIT_TOT,( bit == rle_last_bit )?1:0);
	rle_last_bit = bit;
}

bool arithGetBitRLE(arithInfo * ari)
{
bool bit;
	bit = arithDecBit(ari,RLE_BIT_MID,RLE_BIT_TOT);
	if ( ! bit ) rle_last_bit ^= 1;

return rle_last_bit;
}

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

#define tabdepth(fp,depth) do { int i; for(i=0;i<depth;i++) fprintf(fp,"   "); } while(0)
#define depthputs(fp,str,depth) do { tabdepth(fp,depth); fputs(str,fp); } while(0)

void printSubbandInfo(subband *sb,FILE *fp,int depth)
{
	if ( !sb ) return;
	if ( sb->leaf ) {
		subband_leaf *sbl = (subband_leaf *) sb;
		//tabdepth(fp,depth); 
		fprintf(fp,"%08X ; prev = %08X , parent = %08X\n",sbl,sbl->prev,sbl->parent);
		tabdepth(fp,depth); fprintf(fp,"%dx%d",sbl->width,sbl->height);
		if ( sbl->prev ) fprintf(fp," prev = %dx%d ",sbl->prev->width,sbl->prev->height);
		if ( sbl->parent ) fprintf(fp," parent = %dx%d ",sbl->parent->width,sbl->parent->height);
		fprintf(fp,"\n");
		tabdepth(fp,depth); fprintf(fp,"act = %d, max = %d, type = %d",sbl->activity,sbl->max_value,sbl->type);
		if ( sbl->transposed ) fprintf(fp," transposed");	
		fprintf(fp,"\n");
		
	} else {
		fprintf(fp,"%08X ; prev = %08X\n",sb,sb->prev);
		if ( sb->LL ) { depthputs(fp,"LL:",depth); printSubbandInfo(sb->LL,fp,depth+1); }
		if ( sb->LH ) { depthputs(fp,"LH:",depth); printSubbandInfo(sb->LH,fp,depth+1); }
		if ( sb->HL ) { depthputs(fp,"HL:",depth); printSubbandInfo(sb->HL,fp,depth+1); }
		if ( sb->HH ) { depthputs(fp,"HH:",depth); printSubbandInfo(sb->HH,fp,depth+1); }
	}
}

subband_leaf * LLleaf(subband *sb)
{
	while ( ! sb->leaf ) sb = sb->LL;
	return (subband_leaf *) sb;
}


subband_leaf * seekLeaf(subband_leaf **leaves,int num,int w,int type)
{
subband_leaf * vs;
int j;
	for(j=num-1;j>=0;j--) {
		vs = leaves[j];
		if ( vs->width == w && vs->type == type )
			return vs;
	}
return NULL;
}
subband * seekNode(subband_leaf **leaves,int num,int w,int type)
{
subband * vs;
int j;
	for(j=num-1;j>=0;j--) {
		vs = (subband *) leaves[j];
		while(vs = vs->prev) {
			if ( vs->width == w && vs->type == type )
				return vs;
		}
	}
return NULL;
}

subband_leaf * seekLeafByType(subband_leaf **leaves,int num,int w,int type)
{
int j;
subband_leaf * ret;
	j = type;
	if ( j == TYPE_LL ) j = TYPE_HH;
	
	if ( ret = seekLeaf(leaves,num,w,j) ) return ret;
	if ( (++j) == 4 ) j = 0;
	if ( ret = seekLeaf(leaves,num,w,j) ) return ret;
	if ( (++j) == 4 ) j = 0;
	if ( ret = seekLeaf(leaves,num,w,j) ) return ret;
	if ( (++j) == 4 ) j = 0;
	if ( ret = seekLeaf(leaves,num,w,j) ) return ret;
return NULL;
}

subband * seekNodeByType(subband_leaf **leaves,int num,int w,int type)
{
int j;
subband * ret;
	j = type;
	if ( j == TYPE_LL ) j = TYPE_HH;
	
	if ( ret = seekNode(leaves,num,w,j) ) return ret;
	if ( (++j) == 4 ) j = 0;
	if ( ret = seekNode(leaves,num,w,j) ) return ret;
	if ( (++j) == 4 ) j = 0;
	if ( ret = seekNode(leaves,num,w,j) ) return ret;
	if ( (++j) == 4 ) j = 0;
	if ( ret = seekNode(leaves,num,w,j) ) return ret;
return NULL;
}