wav_trf.c

matlab 例程

#include "wav_basic.h"
#include "wav_filters_extern.h"
#include "alloc.h"
#include "wav_gen.h"


#define set_filter_param(fname) \
{ \
	MFLP= fname ## _FLP; \
	MFHP= fname ## _FHP; \
	MILP= fname ## _ILP; \
	MIHP= fname ## _IHP; \
	\
	Nflp= fname ## _Nflp; \
	Nfhp= fname ## _Nfhp; \
	Nilp= fname ## _Nilp; \
	Nihp= fname ## _Nihp; \
	\
	PS= fname ## _PS; \
	\
	begflp=mfl= fname ## _begflp; \
	begfhp=mfh= fname ## _begfhp; \
	begilp=mil= fname ## _begilp; \
	begihp=mih= fname ## _begihp; \
}


void choose_filter(char name,int tap)

{
	switch (name) {
	   	case 'A':
			switch (tap) {
				default:
					set_filter_param(Antonini);
			}
			break;

		case 'H':
			switch (tap) {
				default:
					set_filter_param(Haar);
			}
			break;

	
        
		case 'B':
			switch (tap) {
				default:
					set_filter_param(DB2);
			}
			break;
	case 'D':
			switch (tap) {
				default:
					set_filter_param(DB3);
			}
			break;
		
		default:
			switch (tap) {
				case 9:
					set_filter_param(D79);
					break;

				default:
					set_filter_param(D8);
			}
	}
}


//
void wav2d_inpl(float **buffer_t,int Ni,int Nj,int levs,float *lp,int Nl,
						float *hp,int Nh,char forw,int *shift_arr_row,int *shift_arr_col,int k)

{
	int i,dimi,dimj;

	dimi=(forw)?Ni:(Ni>>(levs-1)); //forw为1则为正小波变换,为0为反小波变换
	dimj=(forw)?Nj:(Nj>>(levs-1));
if(k==1)
{
	for(i=0;i<levs;i++) {

		if(forw) {


			begflp=mfl+shift_arr_row[i];   
			begfhp=mfh-shift_arr_row[i];

			filt_n_dec_all_rows(buffer_t,dimi,dimj,lp,Nl,hp,Nh,k);


			begflp=mfl+shift_arr_col[i];   
			begfhp=mfh-shift_arr_col[i];

			filt_n_dec_all_cols(buffer_t,dimi,dimj,lp,Nl,hp,Nh,k);
		}
		else {

			begilp=mil+shift_arr_row[levs-1-i];
			begihp=mih-shift_arr_row[levs-1-i];

			ups_n_filt_all_rows(buffer_t,dimi,dimj,lp,Nl,hp,Nh,levs-1-i,shift_arr_row);

			begilp=mil+shift_arr_col[levs-1-i];
			begihp=mih-shift_arr_col[levs-1-i];

			ups_n_filt_all_cols(buffer_t,dimi,dimj,lp,Nl,hp,Nh,levs-1-i,shift_arr_col);
		}

		dimi=(forw)?(dimi>>1):(dimi<<1);
		dimj=(forw)?(dimj>>1):(dimj<<1);
	}
}
if(k==2)
{
	for(i=0;i<levs;i++) {

		if(forw) {


			begflp=mfl+shift_arr_row[i];   
			begfhp=mfh-shift_arr_row[i];

			filt_n_dec_all_rows(buffer_t,dimi,dimj,lp,Nl,hp,Nh,k+1);


			begflp=mfl+shift_arr_col[i];   
			begfhp=mfh-shift_arr_col[i];

			filt_n_dec_all_cols(buffer_t,dimi,dimj,lp,Nl,hp,Nh,k+1);
            

		}
	

		dimi=(forw)?(dimi>>1):(dimi<<1);
		dimj=(forw)?(dimj>>1):(dimj<<1);
	}
}
if(k==3)
{
	for(i=0;i<levs;i++) {

		if(forw) {


			begflp=mfl+shift_arr_row[i];   
			begfhp=mfh-shift_arr_row[i];

			filt_n_dec_all_rows(buffer_t,dimi,dimj,lp,Nl,hp,Nh,k-1);


			begflp=mfl+shift_arr_col[i];   
			begfhp=mfh-shift_arr_col[i];

			filt_n_dec_all_cols(buffer_t,dimi,dimj,lp,Nl,hp,Nh,k-1);
            begflp=mfl+shift_arr_row[i];   
			begfhp=mfh-shift_arr_row[i];

			filt_n_dec_all_rows(buffer_t,dimi,dimj,lp,Nl,hp,Nh,k);
	        begflp=mfl+shift_arr_col[i];   
			begfhp=mfh-shift_arr_col[i];

			filt_n_dec_all_cols(buffer_t,dimi,dimj,lp,Nl,hp,Nh,k);
            

		}
	

		dimi=(forw)?(dimi>>1):(dimi<<1);
		dimj=(forw)?(dimj>>1):(dimj<<1);
	}
}
if(k==4)
{
	for(i=0;i<levs;i++) {

		if(forw) {


			begflp=mfl+shift_arr_row[i];   
			begfhp=mfh-shift_arr_row[i];

			filt_n_dec_all_rows(buffer_t,dimi,dimj,lp,Nl,hp,Nh,k-2);


			begflp=mfl+shift_arr_col[i];   
			begfhp=mfh-shift_arr_col[i];

			filt_n_dec_all_cols(buffer_t,dimi,dimj,lp,Nl,hp,Nh,k-2);
            begflp=mfl+shift_arr_row[i];   
			begfhp=mfh-shift_arr_row[i];

			filt_n_dec_all_rows(buffer_t,dimi,dimj,lp,Nl,hp,Nh,k-2);
	        begflp=mfl+shift_arr_col[i];   
			begfhp=mfh-shift_arr_col[i];

			filt_n_dec_all_cols(buffer_t,dimi,dimj,lp,Nl,hp,Nh,k-2);
            begflp=mfl+shift_arr_row[i];   
			begfhp=mfh-shift_arr_row[i];

			filt_n_dec_all_rows(buffer_t,dimi,dimj,lp,Nl,hp,Nh,k-1);
            begflp=mfl+shift_arr_col[i];   
			begfhp=mfh-shift_arr_col[i];

			filt_n_dec_all_cols(buffer_t,dimi,dimj,lp,Nl,hp,Nh,k-1);

		}
	

		dimi=(forw)?(dimi>>1):(dimi<<1);
		dimj=(forw)?(dimj>>1):(dimj<<1);
	}
}
}

//
void wav2d_inpl_edge(float **buffer_t,int Ni,int Nj,int levs,float *lp,int Nl,
						float *hp,int Nh,char forw,int *shift_arr_row,int *shift_arr_col)

{
	int i,dimi,dimj;

	dimi=(forw)?Ni:(Ni>>(levs-1)); //forw为1则为正小波变换,为0为反小波变换
	dimj=(forw)?Nj:(Nj>>(levs-1));

	for(i=0;i<levs;i++) {

		if(forw) {


			begflp=mfl+shift_arr_row[i];   
			begfhp=mfh-shift_arr_row[i];

			filt_n_dec_all_rows_edge(buffer_t,dimi,dimj,lp,Nl,hp,Nh);


			begflp=mfl+shift_arr_col[i];   
			begfhp=mfh-shift_arr_col[i];

			filt_n_dec_all_cols_edge(buffer_t,dimi,dimj,hp,Nh,lp,Nl);
		}
		else {

			begilp=mil+shift_arr_row[levs-1-i];
			begihp=mih-shift_arr_row[levs-1-i];

			ups_n_filt_all_rows(buffer_t,dimi,dimj,lp,Nl,hp,Nh,levs-1-i,shift_arr_row);

			begilp=mil+shift_arr_col[levs-1-i];
			begihp=mih-shift_arr_col[levs-1-i];

			ups_n_filt_all_cols(buffer_t,dimi,dimj,lp,Nl,hp,Nh,levs-1-i,shift_arr_col);
		}

		dimi=(forw)?(dimi>>1):(dimi<<1);
		dimj=(forw)?(dimj>>1):(dimj<<1);
	}
}




//
void wavpack2d_inpl(float **image,int Ni,int Nj,int levs,float *lp,int Nl,
						float *hp,int Nh,char forw,int *shift_arr_row,int *shift_arr_col)

{
	int i,k,l,dimi,dimj;
	int p,q,rcnt,ccnt;
	float **buffer;

	buffer=allocate_2d_float(Ni,Nj,0);

	dimi=(forw)?Ni:(Ni>>(levs-1));
	dimj=(forw)?Nj:(Nj>>(levs-1));

	for(i=0;i<levs;i++) {

		rcnt=Ni/dimi;

		ccnt=Nj/dimj;

		if(forw&&(i>=LL_ONLY_AFTER_LEV)) {
			rcnt=ccnt=1;
		}
		else if((!forw)&&(levs-1-i>=LL_ONLY_AFTER_LEV)) {
			rcnt=ccnt=1;
		}

		for(p=0;p<rcnt;p++)
			for(q=0;q<ccnt;q++) {

				for(k=0;k<dimi;k++)
					for(l=0;l<dimj;l++)
						buffer[k][l]=image[p*dimi+k][q*dimj+l];

				if(forw) {


					begflp=mfl+shift_arr_row[i];   
					begfhp=mfh-shift_arr_row[i];

					filt_n_dec_all_rows(buffer,dimi,dimj,lp,Nl,hp,Nh,1);

					begflp=mfl+shift_arr_col[i];   
					begfhp=mfh-shift_arr_col[i];

					filt_n_dec_all_cols(buffer,dimi,dimj,lp,Nl,hp,Nh,1);
				}
				else {

					begilp=mil+shift_arr_row[levs-1-i];
					begihp=mih-shift_arr_row[levs-1-i];

					ups_n_filt_all_rows(buffer,dimi,dimj,lp,Nl,hp,Nh,levs-1-i,shift_arr_row);

					begilp=mil+shift_arr_col[levs-1-i];
					begihp=mih-shift_arr_col[levs-1-i];

					ups_n_filt_all_cols(buffer,dimi,dimj,lp,Nl,hp,Nh,levs-1-i,shift_arr_col);
				}

				for(k=0;k<dimi;k++)
					for(l=0;l<dimj;l++)
						image[p*dimi+k][q*dimj+l]=buffer[k][l];
			}

		dimi=(forw)?(dimi>>1):(dimi<<1);
		dimj=(forw)?(dimj>>1):(dimj<<1);
	}

	free_2d_float(buffer,Ni);
}



void wavelet_transform_test(float **original,float **original1,int Ni,int Nj,int k)

{  	int i,j;
	float **buffer_t;
	int shift_arr_r[MAX_ARR_SIZE],shift_arr_c[MAX_ARR_SIZE];
    int PACKET=0;
	
	const int levs=1;
	int Nl,Nh;
	float *lp,*hp;


	
	choose_filter('D',1);

	
	buffer_t=allocate_2d_float(Ni,Nj,0);

	for(i=levs-1;i>=0;i--) {

		shift_arr_r[i]=shift_arr_c[i]=1;
	}

			

	for(i=0;i<Ni;i++) {
		
		for(j=0;j<Nj;j++) {
			
			buffer_t[i][j]=original[i][j];
		}
	}
				
	lp=MFLP;Nl=Nflp;  
	hp=MFHP;Nh=Nfhp;
				
	if(PACKET) {
		wavpack2d_inpl(buffer_t,Ni,Nj,levs,lp,Nl,hp,Nh,1,shift_arr_r,shift_arr_c);
	}
	else {
		
		wav2d_inpl(buffer_t,Ni,Nj,levs,lp,Nl,hp,Nh,1,shift_arr_r,shift_arr_c,k);
	}
	

					
				
	
			
for(i=0;i<Ni;i++) {
		
		for(j=0;j<Nj;j++) {
			
		original1[i][j]=buffer_t[i][j];
		
		}
	}
for(i=0;i<Ni;i++) {
		
		for(j=0;j<Nj;j++) {
			if(original1[i][j]>255)
				original1[i][j]=255;
		    if(original1[i][j]<0)
				original1[i][j]=0;
		original1[i][j]=original1[i][j];	
		}
}

	free_2d_float(buffer_t,Ni);
}

void wavelet_transform_edge(float **original,float **original1,int Ni,int Nj,int n)

{  	int i,j,k;
	float **buffer_t;
	int shift_arr_r[MAX_ARR_SIZE],shift_arr_c[MAX_ARR_SIZE];
    int PACKET=0;
	
	const int levs=2;
	int Nl,Nh;
	float *lp,*hp;


	
	choose_filter('D',1);

	
	buffer_t=allocate_2d_float(Ni,Nj,0);

	for(i=levs-1;i>=0;i--) {

		shift_arr_r[i]=shift_arr_c[i]=1;
	}

			

	for(i=0;i<Ni;i++) {
		
		for(j=0;j<Nj;j++) {
			
			buffer_t[i][j]=original[i][j];
		}
	}
	for(k=levs-1;k>0;k--)				
{	lp=MFLP;Nl=Nflp;  
	hp=MFHP;Nh=Nfhp;
				
	if(PACKET) {
		wavpack2d_inpl(buffer_t,Ni,Nj,levs,lp,Nl,hp,Nh,1,shift_arr_r,shift_arr_c);
	}
	else {
	    wav2d_inpl_edge(buffer_t,Ni,Nj,levs,lp,Nl,lp,Nl,1,shift_arr_r,shift_arr_c);
	}	    
}			
	lp=MFLP;Nl=Nflp;  
	hp=MFHP;Nh=Nfhp;
				
	if(PACKET) {
		wavpack2d_inpl(buffer_t,Ni,Nj,levs,lp,Nl,hp,Nh,1,shift_arr_r,shift_arr_c);
	}
	else {
		if(n==0) wav2d_inpl_edge(buffer_t,Ni,Nj,levs,lp,Nl,hp,Nh,1,shift_arr_r,shift_arr_c);
	     else    wav2d_inpl_edge(buffer_t,Ni,Nj,levs,hp,Nh,lp,Nl,1,shift_arr_r,shift_arr_c);
	}
	

for(i=0;i<Ni;i++) {
		
		for(j=0;j<Nj;j++) {
			
		original1[i][j]=buffer_t[i][j];
		
		}
	}
for(i=0;i<Ni;i++) {
		
		for(j=0;j<Nj;j++) {
			if(original1[i][j]>255)
				original1[i][j]=255;
		    if(original1[i][j]<0)
				original1[i][j]=0;
		original1[i][j]=original1[i][j];	
		}
}

	free_2d_float(buffer_t,Ni);
}

void wavelet_inverse_test(float **original,float **original1,int Ni,int Nj)

{  	int i,j;
	float **buffer_t;
	int shift_arr_r[MAX_ARR_SIZE],shift_arr_c[MAX_ARR_SIZE];
    int PACKET=0;
	const int levs=2;
	int Nl,Nh;
	float *lp,*hp;
	float mse=0;

	choose_filter('B',1);

	buffer_t=allocate_2d_float(Ni,Nj,0);

	for(i=levs-1;i>=0;i--) {

		shift_arr_r[i]=shift_arr_c[i]=0;
	}

			

	for(i=0;i<Ni;i++) {
		
		for(j=0;j<Nj;j++) {
			
			buffer_t[i][j]=original[i][j];
		}
	}

	lp=MFLP;Nl=Nflp;  
	hp=MFHP;Nh=Nfhp;
				
	if(PACKET) {
		
		// Packet transform.
		wavpack2d_inpl(buffer_t,Ni,Nj,levs,lp,Nl,hp,Nh,1,shift_arr_r,shift_arr_c);
	}
	else {
		
		// Regular wavelets.
		wav2d_inpl(buffer_t,Ni,Nj,levs,lp,Nl,hp,Nh,1,shift_arr_r,shift_arr_c,1);
	}
		
	lp=MILP;Nl=Nilp;  
	hp=MIHP;Nh=Nihp;
				
	if(PACKET)
		wavpack2d_inpl(buffer_t,Ni,Nj,levs,lp,Nl,hp,Nh,0,shift_arr_r,shift_arr_c);
	else
		wav2d_inpl(buffer_t,Ni,Nj,levs,lp,Nl,hp,Nh,0,shift_arr_r,shift_arr_c,1);
				
					
				
	
			
for(i=0;i<Ni;i++) {
		
		for(j=0;j<Nj;j++) {
			
		original1[i][j]=buffer_t[i][j];
		
		}
	}
for(i=0;i<Ni;i++) {
		
		for(j=0;j<Nj;j++) {
			if(original1[i][j]>256)
				original1[i][j]=255;
		    if(original1[i][j]<0)
				original1[i][j]=0;
		original1[i][j]=(float)((int)(original1[i][j]+.5));	
		}
}
	free_2d_float(buffer_t,Ni);
}