www.pudn.com > wavecode.rar > dwt.c

 
 
#include  
#include  
#include  
#include  
#include  
#include  
 
 
#define WT_ROUND_IN	0 
#define WT_ROUND_OUT	0.5 
 
#define F_EXTPAD 4 
#define D_EXTPAD 2 
 
static double *x_alloc = NULL;	//** temp work 
 
void forwardDWT(double *x_in, int N) 
{ 
int i, n, half; 
double *x,*r,*d; 
 
	x = x_alloc + F_EXTPAD; 
	memcpy(x,x_in,sizeof(double)*N); 
	for(i=1;i<=F_EXTPAD;i++)  
	  { 
		x[-i] = x[i]; 
		x[(N-1) + i] = x[(N-1) - i]; 
	  } 
 
	half = N>>1;   
 
	r = x_in; d = x_in + half; 
	for (n=half;n--;)  
	  { 
 
		*r++ = 0.0378285 * (x[ 4] + x[-4]) - 0.0238495 * (x[ 3] + x[-3])  
		      - 0.110624 * (x[ 2] + x[-2]) +0.377403 * (x[ 1] + x[-1])  
		      + 0.852699 * x[0]; 
 
		x++; 
 
		*d++ =  - 0.064539 * (x[ 3] + x[-3]) + 0.040689 * (x[ 2] + x[-2]) +  
				0.418092 * (x[ 1] + x[-1]) - 0.788486 * x[ 0] ; 
 
		x++; 
    } 
} 
 
void inverseDWT(double *x, int N) 
{ 
int i, n, half; 
double *r, *d; 
 
	half = N/2; 
 
	r = x_alloc + D_EXTPAD;   
	d = x_alloc + D_EXTPAD + half+D_EXTPAD+D_EXTPAD;   
	memcpy(r,x,half*sizeof(double)); 
	memcpy(d,x+half,half*sizeof(double)); 
 
	for(i=1;i<=D_EXTPAD;i++)  
	 { 
		r[-i] = r[i]; 
		r[(half-1)+i] = r[half - i]; 
		d[-i] = d[i-1]; 
		d[(half-1)+i] = d[(half-1) - i]; 
	 } 
  
	for (n = half; n--;) 
	  { 
 
		*x++ = 0.788486 * r[0] - 0.0406894 * ( r[1] + r[-1] ) 
		     - 0.023849 * (d[1] + d[-2]) + 0.377403 * (d[0] + d[-1]); 
 
		*x++ = 0.418092 * (r[1] + r[0]) - 0.0645389 * (r[2] + r[-1]) 
		     - 0.037829 * (d[2] + d[-2]) + 0.110624 * (d[1] + d[-1]) - 0.852699 * d[0]; 
 
		d++; r++; 
	  } 
} 
 
void waveletTransform2D(int **raw, 
                        double **wave, 
                        int width, 
                        int height, 
                        int levels, 
                        bool inverse) 
{ 
int x, y, w, yw, h, l; 
int *raw_ptr; 
double *buffer,*wave_ptr; 
 
    /* Check the dimensions for compatability. */ 
 
    if (width%(1 << levels) || height%(1 << levels)) 
     { 
		errputs("width and height must be divisible by 2^levels"); 
		exit(10); 
	  } 
 
	if ( (x_alloc = malloc(sizeof(double)*(width+height+F_EXTPAD+F_EXTPAD))) == NULL )  
	 { 
		errputs("malloc failed"); exit(10); 
	 } 
   
    /* Allocate a work array (for transposing columns) */ 
     
   	if ( (buffer = newarray(double,height)) == NULL )  
   	{ 
		 errputs("malloc failed"); exit(10); 
	   } 
 
    /* Compute the wavelet transform. */ 
   
	if ( !inverse )  
	 { /* forward transform. */ 
 
		/** raw -> wave **/ 
 
		wave_ptr = wave[0]; raw_ptr = raw[0]; 
		for(x=width*height;x--;) *wave_ptr++ = (double)(*raw_ptr++) + WT_ROUND_IN; 
		 
		for (l = 0; l < levels; l++)  
		  { 
			w = width >> l; 
			h = height >> l; 
       
			/* Rows. */ 
	 
			for (y = 0; y < h; y++) 
				forwardDWT(wave[y], w); 
     
			/* Columns. */ 
	 
			wave_ptr = wave[0]; 
			for (x = 0; x < w; x++)  
			{ 
				for (y = 0,yw=0; y < h; y++,yw+=width) buffer[y] = wave_ptr[yw]; 
				forwardDWT(buffer, h); 
				for (y = 0,yw=0; y < h; y++,yw+=width) wave_ptr[yw] = buffer[y]; 
				wave_ptr ++; 
			} 
		} 
 
   }  
   else { 
 
		for (l = levels-1; l >= 0; l--) 
		 { 
			w = width >> l; 
			h = height >> l; 
 
			/* Columns. */ 
	 
			wave_ptr = wave[0]; 
			for (x = 0; x < w; x++)  
			{ 
				for (y = 0,yw=0; y < h; y++,yw+=width) buffer[y] = wave_ptr[yw]; 
				inverseDWT(buffer, h); 
				for (y = 0,yw=0; y < h; y++,yw+=width) wave_ptr[yw] = buffer[y]; 
				wave_ptr ++; 
			} 
 
			/* Rows. */ 
	 
			for (y = 0; y < h; y++) 
				inverseDWT(wave[y], w); 
    	} 
 
		wave_ptr = wave[0]; 
		raw_ptr = raw[0]; 
		for(x=width*height;x--;) 
			*raw_ptr++ = (int) ( (*wave_ptr++)  + WT_ROUND_OUT ); 
	} 
 
	free(x_alloc); x_alloc = NULL; 
	free(buffer); 
}