#include <fftw3.h>
#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include "project.h"

int main(int argc, char **argv) {
	int		n,m;
	fftw_plan	transform;
	fftw_complex 	*image, *out;
	int		smelly;
	char		c;
	int		i,j;
	double		max = 0;
	double		*value;

	//parse the input
	if (!(getchar() == 'P' && getchar() == '5')) {		//pgm files start with "P5"
		fprintf(stderr,"Need pgm data from stdin\n");
		return(EXIT_FAILURE);
	}
	NEXTLINE();
	while ((c = getchar()) == '#') NEXTLINE();		//ignore comments
  	ungetc(c,stdin);
	scanf("%i%i%i",&n,&m,&smelly);				//parse the important numbers
	NEXTLINE();
	fprintf(stderr,"%ix%i image read. %i - should be 255 :-)\n",n,m,smelly);

	//prepare the Row-major Format 2-D arrays
	image = fftw_malloc(sizeof(fftw_complex) * n*m);
	out = fftw_malloc(sizeof(fftw_complex) * n*m);
	value = (double *)malloc(sizeof(double)*n*m);

	//fill up the input array using the pgm picture
	for (j=0; j<m; j++) {
		for(i=0; i<n; i++) 
			image[j*n+i][0] = (double)getchar()*pow(-1,i+j);	//checkerboard to get fourier transform centered properly
	}

  
  	//do fourier transform (this is the easy bit ;-) )
	transform = fftw_plan_dft_2d(n,m,image,out, FFTW_FORWARD,FFTW_ESTIMATE);
	fftw_execute(transform);

	//take the size of the output array, and found out a suitible factor to scale by
	for (j=0; j<m; j++) {
		for (i=0; i<n; i++) {
			//note, FFTW scaled it's output by (n*m)
			value[j*n+i] = out[j*n+i][0]*out[j*n+i][0] + out[j*n+i][1]*out[j*n+i][1]/(n*m); 
			value[j*n+i] = log(value[j*n+i] +1.0);	//more detail if we do this
			max = value[j*n+i] > max ? value[j*n+i] : max;
		}
	}


	//output a pgm of the fourier transform
	PGM_HEADER(n,m)
	for (j=0; j<m; j++) {
		for (i=0; i<n; i++) {
			putchar((char)(value[j*n+i]/max*255));
		}
	}

	//fprintf(stderr,"Maximum value was %lf\n",max);
 

	//free memory
	fftw_destroy_plan(transform);
	fftw_free(image);
	fftw_free(out);
	free(value);
	
	return(EXIT_SUCCESS);
}