wpacket.c

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

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

	the activity measure for deciding whether to transform a level
	is critical.  The ideal measure is of course actually packing	
	it (that's a hassle cuz you've got to correctly assign parents
	and whatnot).

	In particular, we need to take the 'quantizer' into account.

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

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

#include <wtlib/subbands.h>
#include <wtlib/image.h>
#include <wtlib/coder.h>
#include <wtlib/quantim.h>

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

void transLeafToQuadMaybe(subband ** sb_ptr,int levels,int trans)
{
subband_leaf *sbl;
subband *sb;
int oldact;

	sbl = (subband_leaf *) *sb_ptr;
	if ( sbl->width < 4 ) return;

	activitySubbands((subband *)sbl);
	oldact = sbl->activity;

	transQuad(sbl->band,sbl->width,sbl->height,sbl->width+sbl->rowpad,0,trans);

	activitySubbands((subband *)sbl);
	/**
	*	the details of the "activity" measurement are critical for this decision
	*	mse & rmse seem quite poor
	*	the percentage of nonzeros seems to be best
	***/
#if 0
	if ( sbl->activity <= oldact ) {
		// on images where WP helps, this helps more, but on images
		//	where it hurts, this hurts even more
#else
	if ( sbl->activity <  oldact ) {
#endif
		// buono

		sb = makeSubbandQuad(sbl->band,sbl->width,sbl->height,sbl->width+sbl->rowpad,0,sbl->prev);
	
		free(sbl);

		*sb_ptr = sb;

		if ( levels ) {
			transLeafToQuadMaybe(&(sb->LL),levels-1,trans);
			transLeafToQuadMaybe(&(sb->LH),levels-1,trans);
			transLeafToQuadMaybe(&(sb->HL),levels-1,trans);
			transLeafToQuadMaybe(&(sb->HH),levels-1,trans);
		}
	} else {
		// undo it
		transQuad(sbl->band,sbl->width,sbl->height,sbl->width+sbl->rowpad,1,trans);
	}
}

subband * transLeafToQuad(subband_leaf *sbl,int levels,bool doWP,int trans)
{
subband *sb;

	transQuad(sbl->band,sbl->width,sbl->height,sbl->width+sbl->rowpad,0,trans);

	sb = makeSubbandQuad(sbl->band,sbl->width,sbl->height,sbl->width+sbl->rowpad,0,sbl->prev);
	
	free(sbl);

	if ( levels ) {
		sb->LL = transLeafToQuad((subband_leaf *)sb->LL,levels-1,doWP,trans);
		if ( doWP && levels > 1 ) {
			/** makes sure that wavelet packets have proper parents **/
			transLeafToQuadMaybe(&(sb->LH),levels-2,trans);
			transLeafToQuadMaybe(&(sb->HL),levels-2,trans);
			transLeafToQuadMaybe(&(sb->HH),levels-2,trans);
		}
	}

return sb;
}

void detransSubbands(subband *sb,int trans)
{
	if ( !sb || sb->leaf ) return;
	detransSubbands(sb->LL,trans);
	detransSubbands(sb->LH,trans);
	detransSubbands(sb->HL,trans);
	detransSubbands(sb->HH,trans);
	if ( sb->width ) {
		int w,h,fw;
		w = sb->width; fw = sb->width + sb->rowpad;
		h = sb->height;
		while( ! sb->leaf ) sb = sb->LL;
		
		transQuad(((subband_leaf *)sb)->band,w,h,fw,1,trans);
	}
}

subband * transPlaneSubbands(image *im,int plane,int levels,bool doWP,int trans)
{
subband_leaf * sbl;
	sbl = new_subband_leaf();
	sbl->band = im->data[plane][0];
	sbl->width = im->width;
	sbl->height= im->height;
	sbl->rowpad = 0;
	if ( levels )
		return transLeafToQuad(sbl,levels-1,doWP,trans);
	else
		return (subband *)sbl;
}


subband * transImageSubbands(image *im,int levels,bool doWP,int trans)
{
	if ( im->planes <= 1 ) return transPlaneSubbands(im,0,levels,doWP,trans);
	else {
		subband_node * sb;
		sb = new_subband_node();
		if ( im->planes >= 1 ) sb->LL = transPlaneSubbands(im,0,levels,doWP,trans);
		if ( im->planes >= 2 ) sb->LH = transPlaneSubbands(im,1,levels,doWP,trans);
		if ( im->planes >= 3 ) sb->HL = transPlaneSubbands(im,2,levels,doWP,trans);
		if ( im->planes >= 4 ) sb->HH = transPlaneSubbands(im,3,levels,doWP,trans);
		if ( im->planes >= 5 ) errexit("too many planes for my cheezy ass");
		return sb;
	}
}

wavelet * makeSubband(image *im,int levels,double q,
	int stoplen,const coder *coder_template,
	bool doWP,int transformN)
{
wavelet * w;

	if ( (w = newWavelet(im,levels)) == NULL )
		return NULL;

	w->stoplen = stoplen;
	w->coder_template = coder_template;
	w->transform = transformN;

	w->qi = findQuantizers(0,0,q, NULL,levels,0,NULL);

	w->subband_root = transImageSubbands(im,levels,doWP,transformN);

	transposeSubbandLHs(w->subband_root);

	if ( q != 1.0 ) quantizeImage(im,levels,w->qi);

	activitySubbands(w->subband_root);

return w;
}


void unSubband(wavelet * w)
{

	if ( w->qi && w->qi->nquants > 0 ) dequantizeImage(w->im,w->levels,w->qi);

	transposeSubbandLHs(w->subband_root);

	detransSubbands(w->subband_root,w->transform);
}