fft.cpp

compuete the fast fourier transformation of the image. it is very fast compare than fourier transfor

#include <math.h>
#include "stdafx.h"
#include "drishti.h"
#include "MemUtil.h"
#include "FFT.h"

/*
void FFT::ReadFile(char* fileName, Complex** imageArray, int N){
	CFile FileObjectIn;
	{
		if(!FileObjectIn.Open(fileName+CString(".raw"), CFile::modeRead)){
			char a[100];
			sprintf(a, "Unable to open file: %s", fileName);
			throw a;
		}
	}

	unsigned char *buff = new unsigned char[N];
	for(unsigned int i=0;i<N;i++){
		UINT readReturn = FileObjectIn.Read(buff, N);
		if( readReturn != N){
			throw "Error while reading file."; 
		}
		for(int j=0;j<N;j++){
			imageArray[i][j].real( (double) buff[j]);
		}
	}

	FileObjectIn.Close();
	delete [] buff;
}

void FFT::ReadFile(char* fileName, Complex ***imageArray, int N){
	CFile FileObjectIn;
	{
		if(!FileObjectIn.Open(fileName+CString(".raw"), CFile::modeRead)){
			char a[100];
			sprintf(a, "Unable to open file: %s.raw", fileName);
			throw a;
		}
	}

	unsigned char *buff = new unsigned char[3*N];
	for(unsigned int i=0;i<N;i++){
		UINT readReturn = FileObjectIn.Read(buff, 3*N);
		if( readReturn != 3*N){
			char a[100];
			sprintf(a, "Error while reading file at row %d.", i);
			throw a;
		}
		for(int j=0;j<N;j++){
			for(int k=0;k<3;k++){
				imageArray[k][i][j].real( (double) buff[j*3+k]);
			}
		}
	}

	FileObjectIn.Close();
	delete []buff;
}

void FFT::SaveFile(char *fileName, char *ext, Complex **imageArray, int N){
	CFile file;
	{
		if(!file.Open(fileName+CString(ext), CFile::modeWrite|CFile::modeCreate)){
			char a[100];
			sprintf(a, "Unable to save file: %s", fileName);
			throw a;
		}
	}

	unsigned char *buf = new unsigned char[N];
	for(int i=0;i<N;i++){
		for(int j=0;j<N;j++){
			buf[j] = (unsigned char)(round(imageArray[i][j].real()));
		}
		file.Write(buf, N);
	}
	
	file.Close();
	delete []buf;
}

void FFT::SaveFile(char *fileName, char *ext, Complex ***imageArray, int N){
	CFile file;
	{
		if(!file.Open(fileName+CString(ext), CFile::modeWrite|CFile::modeCreate)){
			char a[100];
			sprintf(a, "Unable to save file: %s", fileName+CString(ext));
			throw a;
		}
	}

	unsigned char *buf = new unsigned char[3*N];
	for(int i=0;i<N;i++){
		for(int j=0;j<N;j++){
			for(int k=0;k<3;k++){
				buf[j*3+k] = (unsigned char)(round(imageArray[k][i][j].real()));
			}
		}
		file.Write(buf, 3*N);
	}
	
	file.Close();
	delete []buf;
}

Complex** FFT::Magnitude(Complex** f, int N, int strech){
	Complex** imageArray;
	{
		imageArray = (Complex **)new Complex*[N];
		for(int k1=0;k1<N;k1++)
			imageArray[k1] = new Complex[N];
	}

	{//Initializing different array for image enhancement
		for(int i=0;i<N;i++)
			for(int j=0;j<N;j++)
				imageArray[i][j]=f[i][j];
	}

	for(int i=0;i<N;i++){
		for(int j=0;j<N;j++){
			imageArray[i][j].real(sqrt(pow(imageArray[i][j].real(), 2)+pow(imageArray[i][j].imag(), 2)));
			imageArray[i][j].imag(atan(imageArray[i][j].imag()/imageArray[i][j].real()));
		}
	}

	if(strech)
		Strech(imageArray, N, FFT_LOGARITHM);

	return imageArray;
}

void FFT::WalshTransform(Complex *F, int N, int dir) {
	Complex T;
	int LN=round(log((double)N)/log(2.0));
	int J, I, L, LE, LE1, IP;

	BitReorder(F, N);

	for(L=1;L<=LN;L++){
		LE=pow(2, L);
		LE1=LE/2;
		for(J=1;J<=LE1;J++){
			for(I=J;I<=N;I+=LE){
				IP=I+LE1;
				T=F[IP-1];
				F[IP-1]=F[I-1]-T;
				F[I-1]=F[I-1]+T;
			}
		}
	}

	if(dir==FORWARD){
		for(I=0;I<N;I++){
			F[I]/=(double)N;
		}
	}
}
*/
void FFT::BitReorder(Complex *F, int N) {
	int i, j, k;

	j=1;
	for(i=1;i<=N;i++)
	{
		if(i<j)	{
			Complex temp = F[j-1];
			F[j-1] = F[i-1];
			F[i-1] = temp;
		}
		k=N/2;
		while(j>k){
			j-=k;
			k=(k+1)/2;
		}
		j+=k;
	}
}

void FFT::Four1D(Complex *F, int N, int dir) {
	Complex T;
	int LN=round(log((double)N)/log(2.0));
	int J, I, L, LE, LE1, IP;

	BitReorder(F, N);

	for(L=1;L<=LN;L++){
		LE=pow(2, L);
		LE1=LE/2;
		Complex U(1.0, 0.0);
		Complex W(cos(pi/LE1), dir*sin(pi/LE1));
		for(J=1;J<=LE1;J++){
			for(I=J;I<=N;I+=LE){
				IP=I+LE1;
				T=F[IP-1]*U;
				F[IP-1]=F[I-1]-T;
				F[I-1]=F[I-1]+T;
			}
			U*=W;
		}
	}

	if(dir==FORWARD){
		for(I=0;I<N;I++){
			F[I]/=(double)N;
		}
	}
}

void FFT::Four2D(Complex **f, int N, int dir)
{
	int i,j;

	// Centeralizing
	if(dir==FORWARD){
		//printf("Centralizing...\n");
		for(i=0;i<N;i++){
			for(j=0;j<N;j++){
				float d = pow(-1.0, ((i+j)&1));
				f[i][j].real(f[i][j].real()*d);
				//f[i][j].imag(f[i][j].imag()*d);
			}
		}
	}
	
	// Transform the columns
	Complex *row = new Complex[N];
	for (j=0;j<N;j++) {
		for (i=0;i<N;i++) row[i] = f[i][j]; 
		Four1D(row, N, dir);
		for (i=0;i<N;i++) f[i][j] = row[i];
	}
	delete [] row;
	
	for (i=0;i<N;i++) {
		Four1D(f[i], N, dir);
	}

	// Centeralizing
	if(dir==REVERSE){
		for(int i=0;i<N;i++){
			for(int j=0;j<N;j++){
				f[i][j].real( f[i][j].real()*pow(-1, ((i+j)&1)) );
				//f[i][j].imag( pow(-1, i+j) * f[i][j].imag() );
				f[i][j].imag( 0.0 );
			}
		}
	}
}
/*
void FFT::Four3D(Complex ***c, int N, int dir, int bands){
		for(int i=0;i<bands;i++)
			Four2D(c[i], N, dir);
}

void FFT::Strech(Complex **f, int N, char type){
	Complex max, min;
	int i,j;
	max = max * 0.0;
	min = min / min;
	for(i=0;i<N;i++){
		for(j=0;j<N;j++){
				if(max.real()<f[i][j].real())
					max.real(f[i][j].real());
				if(max.imag()<f[i][j].imag())
					max.imag(f[i][j].imag());
				
				if(min.real()>f[i][j].real())
					min.real(f[i][j].real());
				if(min.imag()>f[i][j].imag())
					min.imag(f[i][j].imag());
		}
	}
	//Contrast Stretching (f(u)-min)*(255/(max-min))
	if(type==FFT_LINEAR){
			for(i=0;i<N;i++){
				for(j=0;j<N;j++){
					f[i][j].real( (f[i][j].real()-min.real()) * (255.0/(max.real()-min.real())) );
					f[i][j].imag( (f[i][j].imag()-min.imag()) * (255.0/(max.imag()-min.imag())) );
				}
			}
	} else if (type==FFT_LOGARITHM){
		//Logarithm Stretching c*log( 1+abs(f(u)) ) Where c = 255 / (log(1+abs(max)))
		double cReal, cImag;
		cReal = 255.0 / ( log( 1 + fabs(max.real()) ));
		cImag = 255.0 / ( log( 1 + fabs(max.imag()) ));
			for(i=0;i<N;i++){
				for(j=0;j<N;j++){
					double lg = log( 1 + fabs(f[i][j].real()) );
					f[i][j].real( cReal * lg);
					lg = log( 1 + fabs(f[i][j].imag()) );
					f[i][j].imag( cImag * lg);
				}
		}
	}
}

void FFT::Strech(Complex ***f, int N, int bands, char type){
	Complex max, min;
	int i,j,k;
	for(i=0;i<N;i++){
		for(j=0;j<N;j++){
			for(k=0;k<bands;k++){
				if(max.real()<f[k][i][j].real())
					max.real(f[k][i][j].real());
				if(max.imag()<f[k][i][j].imag())
					max.imag(f[k][i][j].imag());
				
				if(min.real()>f[k][i][j].real())
					min.real(f[k][i][j].real());
				if(min.imag()>f[k][i][j].imag())
					min.imag(f[k][i][j].imag());
			}
		}
	}
	//Contrast Stretching (f(u)-min)*(255/(max-min))
	if(type==FFT_LINEAR){
		for(k=0;k<bands;k++){
			for(i=0;i<N;i++){
				for(j=0;j<N;j++){
					f[k][i][j].real( (f[k][i][j].real()-min.real()) * (255/(max.real()-min.real())) );
					f[k][i][j].imag( (f[k][i][j].imag()-min.imag()) * (255/(max.imag()-min.imag())) );
				}
			}
		}
	} else if (type==FFT_LOGARITHM){
		//Logarithm Stretching c*log( 1+abs(f(u)) ) Where c = 255 / (log(1+abs(max)))
		double cReal, cImag;
		cReal = 255.0 / ( log( 1 + fabs(max.real()) ));
		cImag = 255.0 / ( log( 1 + fabs(max.imag()) ));
		for(k=0;k<bands;k++){
			for(i=0;i<N;i++){
				for(j=0;j<N;j++){
					double lg = log( 1 + fabs(f[k][i][j].real()) );
					f[k][i][j].real( cReal * lg);
					lg = log( 1 + fabs(f[k][i][j].imag()) );
					f[k][i][j].imag( cImag * lg);
				}
			}
		}
	}
}
*/
void FFT::BitReorder(Complex **F, int N) {
	int i, j, k;

	j=1;
	for(i=1;i<=N;i++)
	{
		if(i<j)	{
			swap(F[0][i-1], F[0][j-1]);
			swap(F[1][i-1], F[1][j-1]);
			swap(F[2][i-1], F[2][j-1]);
		}
		k=N/2;
		while(j>k){
			j-=k;
			k=(k+1)/2;
		}
		j+=k;
	}
}

void FFT::Four1D(Complex **F, int N, int dir) {
	Complex T;
	int LN=round(log((double)N)/log(2.0));
	int J, I, L, LE, LE1, IP;

	BitReorder(F, N);

	for(L=1;L<=LN;L++){
		LE=pow(2, L);
		LE1=LE/2;
		Complex U(1.0, 0.0);
		Complex W(cos(pi/LE1), dir*sin(pi/LE1));
		for(J=1;J<=LE1;J++){
			for(I=J;I<=N;I+=LE){
				IP=I+LE1;

				T=F[0][IP-1]*U;
				F[0][IP-1]=F[0][I-1]-T;
				F[0][I-1]=F[0][I-1]+T;

				T=F[1][IP-1]*U;
				F[1][IP-1]=F[1][I-1]-T;
				F[1][I-1]=F[1][I-1]+T;

				T=F[2][IP-1]*U;
				F[2][IP-1]=F[2][I-1]-T;
				F[2][I-1]=F[2][I-1]+T;
			}
			U*=W;
		}
	}

	if(dir==FORWARD){
		for(I=0;I<N;I++){
			F[0][I]/=(double)N;
			F[1][I]/=(double)N;
			F[2][I]/=(double)N;
		}
	}
}

void FFT::Four2D(Complex ***f, int N, int dir)
{
	int i,j;
	MemUtil u;

	// Centeralizing
	if(dir==FORWARD){
		//printf("Centralizing...\n");
		for(i=0;i<N;i++){
			for(j=0;j<N;j++){
				double d = pow(-1.0, ((i+j)&1));
				f[0][i][j].real(f[0][i][j].real()*d);
				f[1][i][j].real(f[1][i][j].real()*d);
				f[2][i][j].real(f[2][i][j].real()*d);
				//f[i][j].imag(f[i][j].imag()*d);
			}
		}
	}
		
	Complex **row = (Complex **) new Complex*[0];
	u.GetMemory(row, 4, N);
	// Transform the columns
	for (j=0;j<N;j++) {
		for (i=0;i<N;i++) { 
			row[0][i] = f[0][i][j]; 
			row[1][i] = f[1][i][j]; 
			row[2][i] = f[2][i][j]; 
		}
		Four1D(row, N, dir);
		for (i=0;i<N;i++) { 
			f[0][i][j] = row[0][i];
			f[1][i][j] = row[1][i];
			f[2][i][j] = row[2][i];
		}
	}
	u.FreeMemory(row, 4, N);
	
	// Transform the rows
	for (i=0;i<N;i++) {
		Four1D(f[0][i], N, dir);
		Four1D(f[1][i], N, dir);
		Four1D(f[2][i], N, dir);
	}
	
	// Centeralizing
	if(dir==REVERSE){
		for(int i=0;i<N;i++){
			for(int j=0;j<N;j++){
				float num = pow(-1, ((i+j)&1));
				f[0][i][j].real( f[0][i][j].real()*num );
				f[0][i][j].imag( 0.0 );
				f[1][i][j].real( f[1][i][j].real()*num );
				f[1][i][j].imag( 0.0 );
				f[2][i][j].real( f[2][i][j].real()*num );
				f[2][i][j].imag( 0.0 );
			}
		}
	}
}