// 2d_fourier.c
//
// Karl's Nifty 2D FFT thingy
// Nov 4, 2000
#include <stdio.h> 
#include <stdlib.h>
#include <string>
#include <fstream.h>
#include <math.h>
#include <iostream.h>

#define MAXSIGNALSIZE		520
#define INSIGNALSOURCE		argv[1]	//"2D_White_Box.pgm"
#define OUTSIGNALDEST		argv[2] //"2D_Square_Square.pgm"

struct complex {
	double R,I;
};

int MAX(int x, int y);
int MIN(int x, int y);
int getImageFromFile(char *filename, int signalarray[MAXSIGNALSIZE][MAXSIGNALSIZE],
				 int *sizex, int *sizey, int *sdepth);
int writeImageToFile(char *filename, int signalarray[MAXSIGNALSIZE][MAXSIGNALSIZE], 
			int sizex, int sizey, int sdepth);
void TwoDUFFT(complex f[MAXSIGNALSIZE][MAXSIGNALSIZE], int sizex, int sizey, int depth);
void TwoDFFT(complex f[MAXSIGNALSIZE][MAXSIGNALSIZE], int sizex, int sizey, int depth);
void FFT(complex signalarray[], int size);
void uFFT(complex signalarray[], int size);
void FFTs(complex signalarray[], int size);
void uFFTs(complex signalarray[], int size);
int maxpix;


//Main Program **************************************************
int main(int argc, char *argv[]) {
	
	//cout << "Hello\n";
	maxpix=0;
	if (argc<3) {
		cout << "\nUsage: " << argv[0] << " signal.dat out.dat\n\n";
		exit(0);
	}
	
	int insignal[MAXSIGNALSIZE][MAXSIGNALSIZE];
	complex sig[MAXSIGNALSIZE][MAXSIGNALSIZE];
	//cout << "Allocated\n";
	int insizex, insizey, depth;
	//char *out[5];
	
	if (!getImageFromFile(INSIGNALSOURCE, insignal, &insizex, &insizey, &depth)) {
		cout << "\nError, could not open file " << INSIGNALSOURCE << "\n\n";
		exit(0);
	}

	//cout << "Read image " << insizex << " " << insizey << "\n";
	
	for (int x=0; x<insizex; x++) {
		for (int y=0; y<insizey; y++) {
			sig[x][y].R = insignal[x][y];
			sig[x][y].I = 0;
		}
	}

	TwoDFFT(sig,insizex,insizey,depth);

	int hix = 0, hiy = 0, ix = 0, iy = 0;
	float hival = 0;
	for (int x=0; x<insizex; x++) {
		for (int y=0; y<insizey/2; y++) {
			double sum=0;
			for (int x1=x-1; x1<x+2; x1++) {
				for (int y1=y-1; y1<y+2; y1++) {
					int x2=x1,y2=y1;
					if (x1<0) x2=insizex-1;
					if (y1<0) y2=insizey-1;
					if (y1>=insizey) y2=0;
					if (x1>=insizex) x2=0;
					sum+=(sig[x2][y2].R*sig[x2][y2].R)+
						 (sig[x2][y2].I*sig[x2][y2].I);
				}
			}
			float tval = (((sig[x][y].R*sig[x][y].R)+
				(sig[x][y].I*sig[x][y].I)+sum)/sum);
			if (tval > hival) {
				hival=tval;
				hix=x;
				hiy=y;
			}
		}
	}
	
	ix=512-hix;
	iy=512-hiy;

	//cout << "Removed Interference at: " << hix << "," << hiy <<
	//	"	& " << ix << "," << iy << "\n";

	cout << "(" << sig[hix][hiy].R << "," << sig[hix][hiy].I << ")	";
	cout << "(" << sig[ix][iy].R << "," << sig[ix][iy].I << ")\n";

	sig[hix][hiy].R=0;
//(sig[1][1].R+sig[0][1].R+sig[0][511].R+sig[511][0].R)/4;
	sig[hix][hiy].I=0;
//(sig[1][1].I+sig[0][1].I+sig[0][511].I+sig[511][0].I)/4;
	sig[ix][iy].R=0;
	sig[ix][iy].I=0;



	//sig[0][0].R = 0;
	//sig[0][0].I = 0;
	//sig[0][511].R = 0;
	//sig[0][511].I = 0;

	TwoDUFFT(sig,insizex,insizey,depth);

/*	out[0] = (char *)malloc(sizeof(char)*20);
	strcpy(out[0], OUTSIGNALDEST);
	strcat(out[0],"R.pgm");
	out[1] = (char *)malloc(sizeof(char)*20);
	strcpy(out[1], OUTSIGNALDEST);
	strcat(out[1],"I.pgm");
	out[2] = (char *)malloc(sizeof(char)*20);
	strcpy(out[2], OUTSIGNALDEST);
	strcat(out[2],"P.pgm");
	out[3] = (char *)malloc(sizeof(char)*20);
	strcpy(out[3], OUTSIGNALDEST);
	strcat(out[3],"M.pgm");
	out[4] = (char *)malloc(sizeof(char)*20);
	strcpy(out[4], OUTSIGNALDEST);
	strcat(out[4],"PO.pgm");
*/
	
	for (int y=0; y<insizey; y++) {
		for (int x=0; x<insizex; x++) {
			insignal[x][y] = (int)sig[x][y].R;
			//insignal[x][y] = (int)sig[x][y].I;
			//insignal[x][y] = (int)atan(sig[x][y].I/sig[x][y].R);
			//insignal[x][y] = (int)sqrt(sig[x][y].R*sig[x][y].R + sig[x][y].I*sig[x][y].I);
			//insignal[x][y] = (int)(insignal[x][y] * insignal[x][y]);
		}
	}
	
	//cout << "maxpix " << maxpix << "\n";
	
//	for (int loop=0; loop<5; loop++) {
			if (!writeImageToFile(OUTSIGNALDEST, insignal, insizex, insizey, depth)) {
			cout << "\nError, could not write to file " << OUTSIGNALDEST << "\n";
		}
//	}
	
/*
	cout << "Real written:	" << out[0] << "\n";
	cout << "Imaginary written:	" << out[1] << "\n";
	cout << "Phase written:	" << out[2] << "\n";
	cout << "Phase written:	" << out[3] << "\n";
	cout << "Magnitude written:	" << out[4] << "\n\n";	
*/
}

//Dump signal file to an array ***********************************
int getImageFromFile(char *filename, int signalarray[MAXSIGNALSIZE][MAXSIGNALSIZE],
			 int *sizex, int *sizey, int *sdepth) {
	ifstream Infile;
	Infile.open(filename, ios::in);
	if (!Infile) return 0;
	else {
		string image_type;
		int height, width, depth, xin, yin;

		Infile >> image_type >> width >> height >> depth;

		//cout << "Reading ... size = " << height << ", " << width << "\n";

		*sizex = height;
		*sizey = width;
		*sdepth = depth;
		
		if (image_type!="P5") {
			cout << "\nError: " << image_type << " is not a valid image type.\n\n";
			exit(1);
		}

		xin = 0; yin = 0;
		unsigned char greyscale;
		while (!Infile.eof()) {
			Infile.read(&greyscale,1);
			if (greyscale>depth) {
				cout << "Greyscale Out of Bounds! = " << greyscale << "\n";
			}
			else {
				signalarray[xin][yin]=(int)greyscale;
			}
			if (xin<(width-1)) {
				xin++;
			}
			else {
				xin=0;
				yin++;
			}
		}
		Infile.close();
		return 1;
	}
}

//Dump array to signal file ***************************************
int writeImageToFile(char *filename, int signalarray[MAXSIGNALSIZE][MAXSIGNALSIZE],
				int sizex, int sizey, int sdepth) {
	maxpix=0;
	for (int x=0; x<sizex; x++) {
		for (int y=0; y<sizey; y++) {
			if (signalarray[x][y]>maxpix) { maxpix=signalarray[x][y]; }
		}
	}
	//cout << "MAxpix image " << maxpix << "\n";
	ofstream outfile;
	outfile.open(filename, ios::out);
	outfile << "P5\n" << sizex << " " << sizey << "\n"
		<< sdepth << "\n";
	if (!outfile) {
		return 0;
	}
	else {
		outfile.setf(ios::fixed);
		for (int y=0; y<sizey; y++) {
			for (int x=0; x<sizex; x++) {
				//cout << signalarray[x][y] << "\t";
				if (signalarray[x][y]<0) signalarray[x][y]=0;
				outfile << (unsigned char)
					(((double)signalarray[x][y]/(double)maxpix)*(sdepth));
			}
		}
		outfile.close();
		return 1;
	}
}

void TwoDFFT(complex f[MAXSIGNALSIZE][MAXSIGNALSIZE], int sizex, int sizey, int depth) {
	for (int y=0; y<sizey; y++) {
		complex t1[MAXSIGNALSIZE];
		for (int x=0; x<sizex; x++) {
			t1[x].R = f[x][y].R;
			t1[x].I = f[x][y].I;
		}
		FFTs(t1,sizex);
		for (int x=0; x<sizex; x++) {
			f[x][y].R = t1[x].R;
			f[x][y].I = t1[x].I;
		}
	}
	for (int x=0; x<sizex; x++) {
		complex t1[MAXSIGNALSIZE];
		for (int y=0; y<sizey; y++) {
			t1[y].R = f[x][y].R;
			t1[y].I = f[x][y].I;
		}
		FFTs(t1,sizey);
		for (int y=0; y<sizey; y++) {
			f[x][y].R = t1[y].R;
			f[x][y].I = t1[y].I;
		}
	}
}

//FFT starting function ********************************************
void FFTs(complex signalarray[], int size) {
	FFT(signalarray,size);
	//size=1;
	for (int loop = 0; loop<size; loop++) {
		signalarray[loop].R = signalarray[loop].R/sqrt(size);
		signalarray[loop].I = signalarray[loop].I/sqrt(size);
	}
}

void FFT(complex signalarray[], int size) {
	if (size==2) {
		complex temp = signalarray[0];
		signalarray[0].R = signalarray[0].R + signalarray[1].R;
		signalarray[0].I = signalarray[0].I + signalarray[1].I;
		signalarray[1].R = temp.R - signalarray[1].R;
		signalarray[1].I = temp.I - signalarray[1].I;
		//cout << "SigArray (Base) -> (" << signalarray[0].R << "," << signalarray[0].I 
		//	<< ") (" << signalarray[1].R << "," << signalarray[1].I << ")\n";
	}
	else {
		complex g[size/2];
		complex h[size/2];
		for (int loop=0; loop<(size/2); loop++) {
			g[loop]=signalarray[2*loop];
			h[loop]=signalarray[2*loop + 1];
		}
		FFT(g, size/2);
		FFT(h, size/2);
		complex gsum,hsum;
		//cout << "Sig array tally -> ";
		for (int loop=0; loop<size; loop++) {
			signalarray[loop].R = g[loop % (size/2)].R 
									+ cos(2*3.1415*loop/size)*(h[loop % (size/2)].R)
									+ sin(2*3.1415*loop/size)*h[loop % (size/2)].I;
			signalarray[loop].I = g[loop % (size/2)].I 
									- sin(2*3.1415*loop/size)*(h[loop % (size/2)].R)
									+ cos(2*3.1415*loop/size)*h[loop % (size/2)].I;
									
			//cout <<"(" << signalarray[loop].R << "," << signalarray[loop].I << ") ";
		}
		//cout << "\n";
	}
}


void uFFTs(complex signalarray[], int size) {
	uFFT(signalarray,size);
	for (int loop = 0; loop<size; loop++) {
		signalarray[loop].R = signalarray[loop].R/sqrt(size);
		signalarray[loop].I = signalarray[loop].I/sqrt(size);
		//cout <<"(" << signalarray[loop].R << "," << signalarray[loop].I << ") ";
	}
	//cout << "\n";
}

void uFFT(complex signalarray[], int size) {
	if (size==2) {
		complex temp = signalarray[0];
		signalarray[0].R = signalarray[0].R + signalarray[1].R;
		signalarray[0].I = signalarray[0].I + signalarray[1].I;
		signalarray[1].R = temp.R - signalarray[1].R;
		signalarray[1].I = temp.I - signalarray[1].I;
		//cout << "SigArray (Base) -> (" << signalarray[0].R << "," << signalarray[0].I 
		//	<< ") (" << signalarray[1].R << "," << signalarray[1].I << ")\n";
	}
	else {
		complex g[size/2];
		complex h[size/2];
		for (int loop=0; loop<(size/2); loop++) {
			g[loop]=signalarray[2*loop];
			h[loop]=signalarray[2*loop + 1];
		}
		uFFT(g, size/2);
		uFFT(h, size/2);
		complex gsum,hsum;
		//cout << "Sig array tally -> ";
		for (int loop=0; loop<size; loop++) {
			signalarray[loop].R = g[loop % (size/2)].R 
									+ cos(2*3.1415*loop/size)*(h[loop % (size/2)].R)
									- sin(2*3.1415*loop/size)*h[loop % (size/2)].I;
			signalarray[loop].I = g[loop % (size/2)].I 
									+ sin(2*3.1415*loop/size)*(h[loop % (size/2)].R)
									+ cos(2*3.1415*loop/size)*h[loop % (size/2)].I;
									
			//cout <<"(" << signalarray[loop].R << "," << signalarray[loop].I << ") ";
		}
		//cout << "\n";
	}
}


void TwoDUFFT(complex f[MAXSIGNALSIZE][MAXSIGNALSIZE], int sizex, int sizey, int depth) {
	for (int y=0; y<sizey; y++) {
		complex t1[MAXSIGNALSIZE];
		for (int x=0; x<sizex; x++) {
			t1[x].R = f[x][y].R;
			t1[x].I = f[x][y].I;
		}
		uFFTs(t1,sizex);
		for (int x=0; x<sizex; x++) {
			f[x][y].R = t1[x].R;
			f[x][y].I = t1[x].I;
		}
	}
	for (int x=0; x<sizex; x++) {
		complex t1[MAXSIGNALSIZE];
		for (int y=0; y<sizey; y++) {
			t1[y].R = f[x][y].R;
			t1[y].I = f[x][y].I;
		}
		uFFTs(t1,sizey);
		for (int y=0; y<sizey; y++) {
			f[x][y].R = t1[y].R;
			f[x][y].I = t1[y].I;
		}
	}
}

int MAX(int x,int y) {
	if (x>y) return x;
	else return y;
}

int MIN(int x,int y) {
	if (x>y) return y;
	else return x;
}