wavelet.c

这是一个关于wavlet的算法

/* --------------------------------------------------------- */
/*    This program is to fulfil the wavelet decomposition    */
/* --------------------------------------------------------- */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

/* ----------- */
/*    Macro    */
/* ----------- */

#define    DIM       256   /* the size of the images */ 
#define    level     3     /* the decomposition level */

/* ---------------------- */
/*    Global Variables    */
/* ---------------------- */

int    L;                     /* the length of filters must be even */
double *h, *g;                /* the coefficient vectors of filter bank */

char   exp_style;             /* the expanding style of image signal */

double in_image[DIM][DIM];    /* the input image */
double r_image[DIM][DIM];     /* the reconstructed image */

/* ------------------- */
/*    filter_bank()    */
/* ------------------- */

void filter_bank()
{
	int i;

	printf("\n Please input the length of filters: ");
	scanf("%d",&L); //getchar();

	/* allocate the memory for h and g */
	h=(double *)malloc(sizeof(double)*L);  /* low frequency filter */
	g=(double *)malloc(sizeof(double)*L);  /* high frequency filter */

	if((h==NULL)||(g==NULL)) {
		printf("\n Allocation error for h or g in filter_bank() \n");
		exit(1);
	}

	switch(L) {

		case 4:
			h[0]= 0.482962913145;
			h[1]= 0.836516303738;
			h[2]= 0.224143868042;
			h[3]=-0.129409522551;
				g[3]= 0.482962913145;   	//g[1]
				g[2]=-0.836516303738;          // g[0]
				g[1]= 0.224143868042;          // g[-1]
				g[0]= 0.129409522551;          // g[-2]
		break;

		case 6:
			h[0]= 0.332670552950;
			h[1]= 0.806891509311;
			h[2]= 0.459877502118;
			h[3]=-0.135011020010;
			h[4]=-0.085441273882;
			h[5]= 0.035226291882;
				g[5]= 0.332670552950;		//g[1]
				g[4]=-0.806891509311;		//g[0]
				g[3]= 0.459877502118;		//g[-1]
				g[2]= 0.135011020010;   	//g[-2]
				g[1]=-0.085441273882;          // g[-3]
				g[0]=-0.035226291882;          // g[-4]
		break;

		case 8:
			h[0]= 0.230377813309;
			h[1]= 0.714846570553;
			h[2]= 0.630880767930;
			h[3]=-0.027983769417;
			h[4]=-0.187034811719;
			h[5]= 0.030841381836;
			h[6]= 0.032883011667;
			h[7]=-0.010597401785;
				g[7]= 0.230377813309;		//g[1]
				g[6]=-0.714846570553;		//g[0]
				g[5]= 0.630880767930;		//g[-1]
				g[4]= 0.027983769417;		//g[-2]
				g[3]=-0.187034811719;		//g[-3]
				g[2]=-0.030841381836;		//g[-4]
				g[1]= 0.032883011667;		//g[-5]
				g[0]= 0.010597401785;           // g[-6]
		break;

		case 10:
			h[0]= 0.160102397974;
			h[1]= 0.603829269797;
			h[2]= 0.724308528438;
			h[3]= 0.138428145901;
			h[4]=-0.242294887066;
			h[5]=-0.032244869585;
			h[6]= 0.077571493840;
			h[7]=-0.006241490213;
			h[8]=-0.012580751999;
			h[9]= 0.003335725285;
				g[9]= 0.160102397974;		//g[1]
				g[8]=-0.603829269797;		//g[0]
				g[7]= 0.724308528438;		//g[-1]
				g[6]=-0.138428145901;		//g[-2]
				g[5]=-0.242294887066;		//g[-3]
				g[4]= 0.032244869585;		//g[-4]
				g[3]= 0.077571493840;		//g[-5]
				g[2]= 0.006241490213;		//g[-6]
				g[1]=-0.012580751999;		//g[-7]
				g[0]=-0.003335725285;		//g[-8]
		break;

		case 12:
			h[0]= 0.111540743350;
			h[1]= 0.494623890398;
			h[2]= 0.751133908021;
			h[3]= 0.315250350709;
			h[4]=-0.226264693965;
			h[5]=-0.129766867567;
			h[6]= 0.097501605587;
			h[7]= 0.027522865530;
			h[8]=-0.031582039318;
			h[9]= 0.000553842201;
			h[10]= 0.004777257511;
			h[11]=-0.001077301085;
				g[11]= 0.111540743350;		//g[1]
				g[10]=-0.494623890398;		//g[0]
				g[9]= 0.751133908021;		//g[-1]
				g[8]=-0.315250350709;		//g[-2]
				g[7]=-0.226264693965;		//g[-3]
				g[6]= 0.129766867567;		//g[-4]
				g[5]= 0.097501605587;		//g[-5]
				g[4]=-0.027522865530;		//g[-6]
				g[3]=-0.031582039318;		//g[-7]
				g[2]=-0.000553842201;		//g[-8]
				g[1]= 0.004777257511;		//g[-9]
				g[0]= 0.001077301085;		//g[-10]
		break;

		case 14:
			h[0]= 0.077852054085;
			h[1]= 0.396539319482;
			h[2]= 0.729132090846;
			h[3]= 0.469782287405;
			h[4]=-0.143906003929;
			h[5]=-0.224036184994;
			h[6]= 0.071309219267;
			h[7]= 0.080612609151;
			h[8]=-0.038029936935;
			h[9]=-0.016574541631;
			h[10]= 0.012550998556;
			h[11]= 0.000429577973;
			h[12]=-0.001801640704;
			h[13]= 0.000353713800;
				 g[13]= 0.077852054085;		//g[1]
				 g[12]=-0.396539319482;		//g[0]
				 g[11]= 0.729132090846;		//g[-1]
				 g[10]=-0.469782287405;		//g[-2]
				 g[9]=-0.143906003929;		//g[-3]
				 g[8]= 0.224036184994;		//g[-4]
				 g[7]= 0.071309219267;		//g[-5]
				 g[6]=-0.080612609151;		//g[-6]
				 g[5]=-0.038029936935;		//g[-7]
				 g[4]= 0.016574541631;		//g[-8]
				 g[3]= 0.012550998556;		//g[-9]
				 g[2]=-0.000429577973;		//g[-10]
				 g[1]=-0.001801640704;		//g[-11]
				 g[0]=-0.000353713800;          //g[-12]
		break;

		case 16:
			h[0]= 0.054415842243;
			h[1]= 0.312871590914;
			h[2]= 0.675630736297;
			h[3]= 0.585354683654;
			h[4]=-0.015829105256;
			h[5]=-0.284015542962;
			h[6]= 0.000472484574;
			h[7]= 0.128747426620;
			h[8]=-0.017369301002;
			h[9]=-0.044088253931;
			h[10]= 0.013981027917;
			h[11]= 0.008746094047;
			h[12]=-0.004870352993;
			h[13]=-0.000391740373;
			h[14]= 0.000675449406;
			h[15]=-0.000117476784;
				 g[15]= 0.054415842243;         //g[1]
				 g[14]=-0.312871590914;		//g[0]
				 g[13]= 0.675630736297;		//g[-1]
				 g[12]=-0.585354683654;		//g[-2]
				 g[11]=-0.015829105256;		//g[-3]
				 g[10]= 0.284015542962;         //g[-4]
				 g[9]= 0.000472484574;		//g[-5]
				 g[8]=-0.128747426620;		//g[-6]
				 g[7]=-0.017369301002;		//g[-7]
				 g[6]= 0.044088253931;		//g[-8]
				 g[5]= 0.013981027917;		//g[-9]
				 g[4]=-0.008746094047;		//g[-10]
				 g[3]=-0.004870352993;		//g[-11]
				 g[2]= 0.000391740373;		//g[-12]
				 g[1]= 0.000675449406;		//g[-13]
				 g[0]= 0.000117476784;		//g[-14]
		break;

		default:printf("Not prepare !");
		exit(1);
	}

	/*for(i=0;i<L;i+2)
		g[i]=-h[L-i-1];
	for(i=1;i<L;i+2)
		g[i]=h[L-i-1];*/
}

/* --------------- */
/*    wavelet()    */
/* --------------- */

void wavelet(double *c, int n)
{
	int i,j;
	double *tempc;
	
	/* allocate the memory for tempc */
	tempc=(double *)malloc(sizeof(double)*(n+L+L-4));
	if(tempc==NULL) {
		printf("\n Allocation error for tempc in wavelet() \n");
		exit(1);
	}

	for(i=0;i<n;i++) tempc[i+L-2]=c[i];

	for(i=L-3;i>=0;i--)
	{
		switch(exp_style) {
			case 's' :  
			tempc[i]=c[L-3-i];
			tempc[n+L+L-i-5]=c[n-(L-2)+i];
			break;
			
			case 'p':
			tempc[i]=c[n-(L-2)+i];
			tempc[n+L+L-i-5]=c[L-3-i];
			break;

			case 'd':
			tempc[i]=tempc[i+1]*2-tempc[i+2];
			tempc[n+L+L-i-5]=tempc[n+L+L-i-6]*2-tempc[n+L+L-i-7];
			break;

			default:
			printf("\n You should choose s, p or d \n");
			exit(1);
		}
	}

	for(i=0;i<n/2;i++)
	{
		c[i]=0;
		for(j=0;j<L;j++)
			c[i]+=h[j]*tempc[L-2+i+i+j];
		c[i+n/2]=0;
		for(j=0;j<L;j++)
			c[i+n/2]+=g[j]*tempc[i+i+j];
	}
	
	free(tempc);
}

/* ------------ */
/*    wave()    */
/* ------------ */

void wave(double **buf, int n)
{
	int i,j;
	double *c;

	/* allocate the memory for c */
	c=(double *)malloc(sizeof(double)*n);
	if(c==NULL) {
		printf("\n allocation error for c in wave()\n");
		exit(1);
	}

	/* for row */
	for(i=0;i<n;i++)
	{
		for(j=0;j<n;j++)
			c[j]=buf[i][j];	
		wavelet(c,n);
		for(j=0;j<n;j++)
			buf[i][j]=c[j];
	}

	/* for column */
	for(j=0;j<n;j++)
	{
		for(i=0;i<n;i++)
			c[i]=buf[i][j];
		wavelet(c,n);
		for(i=0;i<n;i++)
			buf[i][j]=c[i];
	}

	free(c);
}

/* ---------------- */
/*    iwavelet()    */
/* ---------------- */

void iwavelet(double *c, int n)
{
	int i,j;
	double *tempc;
	
	/* allocate the memory for tempc */
	tempc=(double *)malloc(sizeof(double)*(n+L-2));
	if(tempc==NULL) {
		printf("\n Allocation error for tempc in wavelet() \n");
		exit(1);
	}

	for(i=0;i<n;i++) tempc[i+L/2-1]=c[i];

	for(i=L/2-2;i>=0;i--)
	{
		switch(exp_style) {
			case 's' :  
			tempc[i]=c[L/2-2-i];
			tempc[n+L-i-3]=c[n-(L/2-1)+i];
			break;
			
			case 'p':
			tempc[i]=c[n-(L/2-1)+i];
			tempc[n+L-i-3]=c[L/2-2-i];
			break;

			case 'd':
			tempc[i]=tempc[i+1]*2-tempc[i+2];
			tempc[n+L-i-3]=tempc[n+L-i-4]*2-tempc[n+L-i-5];
			break;

			default:
			printf("\n You should choose s, p or d \n");
			exit(1);
		}
	}

	for(i=0;i<n/2;i++)
	{
		c[i+i+1]=0;
		c[i+i]=0;
		for(j=0;j<L/2;j++)
		{
			c[i+i]+=h[j+j]*tempc[L/2-1+i-j]+g[j+j]*tempc[n/2+L-2+i-j];
			c[i+i+1]+=h[j+j+1]*tempc[L/2-1+i-j]+g[j+j+1]*tempc[n/2+L-2+i-j];
		}	
	}
	
	free(tempc);
}

/* ------------- */
/*    iwave()    */
/* ------------- */

void iwave(double **buf, int n)
{
	int i,j;
	double *c;

	/* allocate the memory for c */
	c=(double *)malloc(sizeof(double)*n);
	if(c==NULL) {
		printf("\n allocation error for c in wave()\n");
		exit(1);
	}

	/* for row */
	for(i=0;i<n;i++)
	{
		for(j=0;j<n;j++)
			c[j]=buf[i][j];	
		iwavelet(c,n);
		for(j=0;j<n;j++)
			buf[i][j]=c[j];
	}

	/* for column */
	for(j=0;j<n;j++)
	{
		for(i=0;i<n;i++)
			c[i]=buf[i][j];
		iwavelet(c,n);
		for(i=0;i<n;i++)
			buf[i][j]=c[i];
	}

	free(c);
}

/* -------------------- */
/*   read_bmp_image()   */
/* -------------------- */

read_bmp_image(double *A, int size, char fn[])
/* A -- image matrix */
/* size --image size */
/* fn -- image file */
{

	int i;
	unsigned char *temp,  /* tempory unsigned char version of A */
		      *ptemp; /* point of temp */	
	double *pA;           /* point of A*/
	FILE *fp;
	long length;          /* the length of head file in bmp image */

	temp=(unsigned char *)malloc(sizeof(unsigned char)*(size*size));
	if(temp==NULL) {
		printf("\n Memory allocation error in read_image_bmp! \n");
		exit(1);
	}

	if((fp=fopen(fn,"rb"))==NULL) {
		printf("\n Cannot open image file %s ! \n",fn);
		exit(1);
	}

	fseek(fp,10,SEEK_SET);
	fread(&length,4,1,fp);

	fseek(fp,length,SEEK_SET);
	if(fread(temp,sizeof(unsigned char),size*size,fp)!=size*size) {
		if(feof(fp)) printf("\n Premature end of file %s \n",fn);
		else         printf("\n File read error %s\n",fn);
		exit(1);
	}
	fclose(fp);

	pA=A; ptemp=temp;
	for(i=0;i<size*size;i++)
	{
		*pA = (double) *ptemp;
		pA++; ptemp++;
	}
	free(temp);

}

/* ----------------------- */
/*    Write_bmp_image()    */
/* ----------------------- */

write_bmp_image(double *A, int size, char fn1[], char fn2[])
{

	int i;
	unsigned char *temp,   /* tempory unsigned char version of A */
		      *ptemp;  /* point of temp */
	double *pA;     /* point of A */
	FILE *fpi, *fpo;
	long length;           /* the length of head file in bmp image*/
	unsigned char *buf;

	if((fpi=fopen(fn1,"rb"))==NULL) {
		printf("\n Cannot open the input file %s \n",fn1);
		exit(1);
	}

	if((fpo=fopen(fn2,"wb"))==NULL) {
		printf("\n Cannot open the output file  %s \n",fn2);
		exit(1);
	}

	fseek(fpi,10,SEEK_SET);
	fread(&length,4,1,fpi);

	buf=(unsigned char *)malloc(sizeof(unsigned char)*length);
	if(buf==NULL) {
		printf("\n Mallocation error in write_image_bmp! \n");
		exit(1);
	}

	fseek(fpi,0,SEEK_SET);
	fread(buf,length,1,fpi);
	fwrite(buf,1,length,fpo);

	temp=(unsigned char *)malloc(sizeof(unsigned char)*size*size);
	if(temp==NULL) {
		printf("\n Mallocation error in write_image_bmp! \n");
		exit(1);
	}

	pA=A;  ptemp=temp;
	for(i=0;i<size*size;i++)
	{
		if(*pA<0) *ptemp=0;
		else {
			if(*pA>255) *ptemp=255;
			else *ptemp = (unsigned char) *pA;
		}
		pA++;  ptemp++;
	}

	fseek(fpo,length,SEEK_SET);
	fwrite(temp, sizeof(unsigned char), size*size, fpo);

	fclose(fpi);  fclose(fpo);  free(temp); free(buf);
}

/* ------------ */
/*    Main()    */
/* ------------ */

void main(int argc, char *argv[])
{
	int i, j, iter, dim=DIM;
	double *temp;
	double **w_image;

	if(argc!=4) {
		printf("Usage: --- wavelet in_image w_image r_image");
		exit(1);
	}

	/* allocate the memory for temp */
	temp=(double *)malloc(sizeof(double)*DIM*DIM);
	if(temp==NULL) {
		printf("\n Allocation error for temp in main() \n");
		exit(1);
	}

	w_image=(double **)malloc(sizeof(double *)*DIM);
	if(w_image==NULL) {
		printf("\n Allocation error for w_image in main() \n");
		exit(1);
	}
	for(i=0;i<DIM;i++)
	{
		w_image[i]=(double *)malloc(sizeof(double)*DIM);
		if(w_image[i]==NULL) {
			printf("\n Allocation error for w_image[%d] in main() \n",i);
			exit(1);
		}
	}

	read_bmp_image(temp, DIM, argv[1]);

	for(i=0;i<DIM;i++)
		for(j=0;j<DIM;j++)
		{
			in_image[i][j]=temp[i*DIM+j];
			w_image[i][j]=temp[i*DIM+j];
		}

	filter_bank();

	printf("\n s --- symmetry period expanded \n");
	printf("\n p --- period expanded \n");
	printf("\n d --- double symmetry period expanded \n");

	printf("\n Please input the expanding style: ");
	scanf("%s",&exp_style);

	/* wavelet analysis */
	for(iter=0;iter<level;iter++)
	{
		wave(w_image, dim);
		dim=dim/2;
	}

	dim=dim*2;
	
	/* output the wavelet decomposition result */
	for(i=0;i<DIM;i++)
		for(j=0;j<DIM;j++)
			temp[i*DIM+j]=w_image[i][j];
	write_bmp_image(temp, DIM, argv[1], argv[2]);

	/* wavelet reconstruct */
	for(iter=0;iter<level;iter++)
	{
		iwave(w_image, dim);
		dim=dim*2;
	}

	for(i=0;i<DIM;i++)
		for(j=0;j<DIM;j++)
		{
			r_image[i][j]=w_image[i][j];
			temp[i*DIM+j]=w_image[i][j];
		}
	write_bmp_image(temp, DIM, argv[1], argv[3]);
	
	for(i=DIM;i>=0;i--)
		free((char *)w_image[i]);
	free(w_image);
	free(temp);
}