fft_image.c

这是一个关于fft的算法

/* -------------------------------------- */
/*    FFT transform of the input image    */
/* -------------------------------------- */

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

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

#define DIM    256               /* the dimension of the input image */       
#define Pi     (acos(0.0)*2.0) 
#define level  255               /* the level of gray */

/* ---------------------- */
/*    Global variables    */
/* ---------------------- */

double    in_image[DIM][DIM];  /* the input image */
double    out_image[DIM][DIM]; /* the reconstructed image*/
double    out[DIM][DIM];       /* the power spectrum of the input image */
double    ZR[DIM][DIM];        /* the real part of the fft of the input image*/
double    ZI[DIM][DIM];

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

main(int argc, char *argv[])
{

	int    i, j;
	double *temp;

	if(argc!=4) {
		printf("\n Usage: --- fft2 in_image fft_image out_image\n");
		exit(-1);
	}

	temp=(double *)malloc(sizeof(double)*(DIM*DIM));
	if(temp==NULL) {
		printf("\nAllocation error in main() of temp!\n");
		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];

	fft2_image();

	for(i=0;i<DIM;i++)
		for(j=0;j<DIM;j++)
			temp[i*DIM+j]=out[i][j];

	write_bmp_image(temp, DIM, argv[1], argv[2]);

	ifft2_image();

	for(i=0;i<DIM;i++)
		for(j=0;j<DIM;j++)
			temp[i*DIM+j]=out_image[i][j];

	write_bmp_image(temp, DIM, argv[1], argv[3]);
}

/* -------------------- */
/*   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);

}

/* ------------------ */
/*    fft2_image()    */
/* ------------------ */

fft2_image()
{
	int    i, j;
	double *A, *B;
	double **temp;            /* used in computing power spectrum */
	double gmax, tmp;

	/* allocation memory for temp vector A and B */

	A=(double *)malloc(sizeof(double)*DIM);
	B=(double *)malloc(sizeof(double)*DIM);
	if((A==NULL)||(B==NULL)) {
		printf("\n Allocation error in fft2_image() of A or B \n");
		exit(-1);
	}

	/* allocation memory for temp matrix */

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

	/* initialization */

	for(i=0;i<DIM;i++)
		for(j=0;j<DIM;j++)
		{
			ZR[i][j]=in_image[i][j];
			ZI[i][j]=0.0;
		}

	/* column fft */

	for(i=0;i<DIM;i++)
	{
		for(j=0;j<DIM;j++)
		{
			A[j]=ZR[i][j];
			B[j]=ZI[i][j];
		}

		fft(A, B);          /* one dimension fft */

		for(j=0;j<DIM;j++)
		{
			ZR[i][j]=A[j];
			ZI[i][j]=B[j];
		}
	}

	/* row fft */

	for(j=0;j<DIM;j++)
	{
		for(i=0;i<DIM;i++)
		{
			A[i]=ZR[i][j];
			B[i]=ZI[i][j];
		}

		fft(A, B);          /* one dimension fft */

		for(i=0;i<DIM;i++)
		{
			ZR[i][j]=A[i];
			ZI[i][j]=B[i];
		}
	}

	/* compute the power spectrum and normalization */

	gmax=0.0;

	for(i=0;i<DIM;i++)
		for(j=0;j<DIM;j++)
		{
			tmp=temp[i][j]=sqrt(ZR[i][j]*ZR[i][j]+ZI[i][j]*ZI[i][j]);
			if(tmp > gmax)    gmax=tmp;
		}

	for(i=0;i<DIM;i++)
		for(j=0;j<DIM;j++)
			out[i][j]=level*(temp[i][j]/gmax);

	/* free the allocated memory */

	for(i=DIM-1;i>=0;i--)
		free((char *)temp[i]);
	free(temp);
	
	free(A); free(B);
	
}

/* ----------- */
/*    fft()    */
/* ----------- */

fft(double *A, double *B)                  /* one dimension fft */
{
	int     k, ii, jj, md2, m1, j1, j2, lmx, lo, lix, li, lm;

	double  darg=2*Pi/DIM;
	double  arg, temp1;
	double  c, s, t1, t2;
	int     l = 8;

	lmx=DIM;

	for(lo=0; lo < l; lo++)
	{
		lix=lmx;
		lmx=lmx/2.0;
		arg=0.0;

		for(lm=0; lm<lmx; lm++)
		{
			c=cos(arg);
			s=-sin(arg);
			arg+=darg;

			for(li=lix; li<DIM+1; li+=lix)
			{
				j1=li-lix+lm;
				j2=j1+lmx;
				t1=A[j1]-A[j2];
				t2=B[j1]-B[j2];
				A[j2]=c*t1+s*t2;
				B[j2]=c*t2-s*t1;
			}
		}
		darg=2.0*darg;
	}

	jj=0; md2=DIM/2; m1=DIM-1;

	for(ii=0; ii<m1; ii++)
	{
		if( (ii-jj) < 0 )
		{
			/* swap(A[ii], A[jj]) */
			temp1=A[ii];
			A[ii]=A[jj];
			A[jj]=temp1;

			/* swap(B[ii], B[jj]) */
			temp1=B[ii];
			B[ii]=B[jj];
			B[jj]=temp1;
		}

		k=md2;

		while( (k-jj) < 1 )
		{ 
			jj=jj-k; k=k/2;
		}

		jj=jj+k;
	}
	
	for(ii=0; ii<DIM; ii++)
	{
		A[ii]=A[ii]/DIM;
		B[ii]=B[ii]/DIM;
	}
}

/* ------------------- */
/*    ifft2_image()    */
/* ------------------- */

ifft2_image()
{
	int i, j, l, m;
	double *A, *B, **temp;
	double gmax=0.0, tmp;

	/* allocate the memory for invariables */
	
	A=(double *)malloc(sizeof(double)*DIM);
	B=(double *)malloc(sizeof(double)*DIM);

	if((A==NULL)||(B==NULL)) {
		printf("\n Allocation error in ifft2_image of A or B \n");
		exit(-1);
	}

	temp=(double **)malloc(sizeof(double *)*DIM);
	if(temp==NULL) {
		printf("\n Allocation error in ifft2_image of temp\n");
		exit(-1);
	}

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

	for(i=0;i<DIM;i++)	/* row inverse fft */
	{
		for(j=0;j<DIM;j++)
		{
			A[j]=ZR[i][j];
			B[j]=ZI[i][j];
		}

		ifft(A, B);       /* one_dimension inverse fft transform */

		for(j=0;j<DIM;j++)
		{
			ZR[i][j]=A[j];
			ZI[i][j]=B[j];
		}
	}

	for(j=0;j<DIM;j++)        /* column inverse fft */
	{
		for(i=0;i<DIM;i++)
		{
			A[i]=ZR[i][j];
			B[i]=ZI[i][j];
		}

		ifft(A, B);       /* one_dimension inverse fft transform */

		for(i=0;i<DIM;i++)
		{
			ZR[i][j]=A[i];
			ZI[i][j]=B[i];
		}
	}

	for(i=0;i<DIM;i++)
		for(j=0;j<DIM;j++)
		{
			tmp=temp[i][j]=sqrt(ZR[i][j]*ZR[i][j]+ZI[i][j]*ZI[i][j]);
			if(tmp > gmax) gmax=tmp;
		}

	/* reconstruct the image */

	for(i=0;i<DIM;i++)
		for(j=0;j<DIM;j++)
			out_image[i][j]=255*(temp[i][j]/gmax);

	/*for(i=0;i<DIM;i++)
		for(j=0;j<DIM/2;j++)
		{
			tmp=out_image[i][j];
			out_image[i][j]=out_image[i][j+DIM/2];
			out_image[i][j+DIM/2]=tmp;
		}*/	

	/* free the allocated memory */

	for(i=DIM-1;i>=0;i--)
		free((char *)temp[i]);
	free(temp);
	
	free(A); free(B);
}

/* ------------ */
/*    ifft()    */
/* ------------ */

ifft(double *A, double *B)                  /* one dimension fft */
{
	int     k, ii, jj, md2, m1, j1, j2, lmx, lo, lix, li, lm;

	double  darg=2*Pi/DIM;
	double  arg, temp1;
	double  c, s, t1, t2;
	int     l = 8;

	lmx=DIM;

	for(lo=0; lo < l; lo++)
	{
		lix=lmx;
		lmx=lmx/2.0;
		arg=0.0;

		for(lm=0; lm<lmx; lm++)
		{
			c=cos(arg);
			s=sin(arg);
			arg+=darg;

			for(li=lix; li<DIM+1; li+=lix)
			{
				j1=li-lix+lm;
				j2=j1+lmx;
				t1=A[j1]-A[j2];
				t2=B[j1]-B[j2];
				A[j2]=c*t1+s*t2;
				B[j2]=c*t2-s*t1;
			}
		}
		darg=2.0*darg;
	}

	jj=0; md2=DIM/2; m1=DIM-1;

	for(ii=0; ii<m1; ii++)
	{
		if( ii-jj < 0 )
		{
			/* swap(A[ii], A[jj]) */
			temp1=A[ii];
			A[ii]=A[jj];
			A[jj]=temp1;

			/* swap(B[ii], B[jj]) */
			temp1=B[ii];
			B[ii]=B[jj];
			B[jj]=temp1;
		}

		k=md2;

		while( (k-jj) < 1 )
		{ 
			jj=jj-k; k=k/2;
		}

		jj=jj+k;
	}
	
}

/* ----------------------- */
/*    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);

}