#include <Magick++.h>
#include <math.h>
#include <stdlib.h>
#include <fstream.h>
#include <complex.h>
#include <sys/time.h>
#include <unistd.h>


using namespace std;
using namespace Magick;

double pi=3.1415;

timeval t1,t2;

void startwatch()
{

  cout << "\n Stopwatch started..." << flush;
      gettimeofday(&t1,NULL);
      

}

void stopwatch()
{

  
  
      gettimeofday(&t2,NULL);
      
      long timediff=(t2.tv_sec-t1.tv_sec)*1000000 + (t2.tv_usec-t1.tv_usec);
      
      cout  << "\nTime Taken : " << timediff << flush;


      
}



class ImageFFT
{

private:
  
  // Bit reverse an integer
  int rev(int k,int b)
  {
    char a[b],i;
    for(i=0;i<b;i++)
      {
	a[i]=k%2;
	k/=2;
      }
    k=0;
    int s=1;
    for(i=b-1;i>=0;i--)
      {
	k+=a[i]*s;
	s*=2;
      
      }
    return k;
  }

  //____________Bit reverse an array - n is the size of array___
  void bitrev(complex <double> a[],complex <double> A[],int n)
  {
    int ln=(int)ceil(log(n)/log(2));
    for(int i=0;i<n;i++)
      A[rev(i,ln)]=a[i];
  }


  //_______Fast Fourier Transform ( also includes inverse )
  // mode = -1  ==> Inverse FFT
  void fft(complex <double> a[],complex <double> A[],int n,int mode=1)
  {
    bitrev(a,A,n);
    int ln=(int)ceil(log(n)/log(2));
    complex<double> t,u;
    int s,j,k;

    for(s=1;s<=ln;s++)
      {
	int m=(int)pow(2,s);
	complex<double> wm(cos(2*pi/m),sin(mode*2*pi/m));
	complex<double> w(1,0);
	for(j=0;j<=m/2-1;j++)
	  {

	    for(k=j;k<=n-1;k+=m)
	      {

		t=w*A[k+m/2];
		u=A[k];
		A[k]=u+t;
		A[k+m/2]=u-t;
	      }

	    w=w*wm;
	  }
      }
    if(mode==-1)
      {
	for(int i=0;i<n;i++)
	  A[i]/=n;
      }  
      
  }
  
  
public:
  
  
  // mode=-1 for inverse , =1 for forward FFT
  void DoImageFFT(complex<double>** buffer,int x,int y,int mode)
  {
    
    
    complex<double>* temp;
    complex<double>* temp1;
    
    
    int n;
    if (y>x) n=y; else n=x;
    
    temp = new complex<double>[n];
    temp1 = new complex<double>[n];
    
    //__Compute FFT for each column first___
    for(int i=0;i<x;i++)
      {
	for(int j=0;j<y;j++)
	  temp[j]=buffer[j][i];
	
	fft(temp,temp1,y,mode);
	
	for(int j=0;j<y;j++)
	  buffer[j][i]=temp1[j];
      }
   

    for(int i=0;i<y;i++)
      {
	for(int j=0;j<x;j++)
	  temp[j]=buffer[i][j];
	  
	fft(temp,temp1,x,mode);
	  
	for(int j=0;j<x;j++)
	  buffer[i][j]=temp1[j];
      }
   

  
  }  
  

};



void writeimage(char* filename,double** pixels,int r,int c,int ask=0)
{

  char buffer[80];
  if(ask)
    {

      char choice;
      cout << "\nWant to save the image(y/n)  : ";
      cin >> choice;
      if(choice!='y' && choice!='Y')
	return;

      cout << "Enter filename : " ;
      cin >> buffer;
    }

  
      
  Image fspectrum(Geometry(c,r),"white");
  fspectrum.classType(DirectClass);
  fspectrum.quantizeColorSpace(GRAYColorspace);
  Pixels fview(fspectrum);
  PixelPacket *fpixels = fview.get(0,0,c,r);

  for(int x=0;x<r;x++)
    for(int y=0;y<c;y++)


      {
	(*fpixels) = ColorGray(pixels[x][y]);
	fpixels++;
      }
    
  fview.sync();
  fspectrum.write((ask==1) ? buffer : filename);

}


void displayimage(double** pixels,int r,int c)
{
  Image fspectrum(Geometry(c,r),"white");
  fspectrum.classType(DirectClass);
  fspectrum.quantizeColorSpace(GRAYColorspace);
  Pixels fview(fspectrum);
  PixelPacket *fpixels = fview.get(0,0,c,r);

  for(int x=0;x<r;x++)
    for(int y=0;y<c;y++)
      {
	(*fpixels) = ColorGray(pixels[x][y]);
	fpixels++;
      }
    
  fview.sync();
  fspectrum.display();


}


PixelPacket* pixels;
int c,r;
complex<double>** inpixels;
int squaresize;
bool imageloaded=false;
double** fftout;

Image* image;
Pixels* view;


void loadimage(char* filename)
{

  image=new Image(filename);
  image->classType(DirectClass);
  view =new Pixels(*image);
  c=image->columns();
  r=image->rows();
  int i,j;


  // convert image to grayscale
  image->classType(DirectClass);
  image->quantizeColorSpace(GRAYColorspace);

  
  pixels=view->get(0,0,c,r);

  imageloaded=true;

}

void pixelsToComplex()
{

 
  int i,j; 

  for(i=0;i<squaresize;i++)
    {
      for(j=0;j<squaresize;j++)
	{

	  if(i<r && j<c)
	    inpixels[i][j]=complex<double>(pow(-1,i+j)* ( (ColorGray)*(pixels+i*c+j) ).shade(),0);
	  else
	    inpixels[i][j]=complex<double>(0,0);
	}
    }

}



void displayFourier(int squaresize)
{

  int i,j;
  
  double max=0;    
  
  for(int i=0; i<squaresize; i++)
    for(int j=0; j<squaresize; j++)
      {
	fftout[i][j] = abs(inpixels[i][j]);
	if (fftout[i][j] > max)
	  max = fftout[i][j];
      }
  //normalize from 0 to 1    

  
  for(int i=0; i<squaresize; i++)
    for(int j=0; j<squaresize; j++)
      {
	if (max!=0)
	  fftout[i][j] = fftout[i][j]/max;
      }

  // range compression    
  int c1;
  c1 = 1000000;
  max=0;
  
  for(int i=0; i<squaresize; i++)
    for(int j=0; j<squaresize; j++)
      {
	fftout[i][j] = log(1 + c1*fftout[i][j]);
	if (fftout[i][j] > max)
	  max = fftout[i][j];
      
      }

  //normalize
  for(int i=0; i<squaresize; i++)
    for(int j=0; j<squaresize; j++)
      {
	if (max!=0)
	  fftout[i][j] = fftout[i][j]/max;
      }
  
  
  displayimage(fftout,squaresize,squaresize);

  writeimage("a.jpg",fftout,squaresize,squaresize,1);
  

}

double distance(int x1,int y1,int x2,int y2)
{

  return sqrt(pow(x1-x2,2)+pow(y1-y2,2));
}



void measure_performance()
{

  int i=0;


  cout << "\nMeasuring performance of 2D FFT\n\n";

  int size;
  cout << "\nUpto what size ? ";
  cin >> size;

  
  int ln=(int)(log(size+1)/log(2));
	       

  size=pow(2,ln);
  
  fstream perf;
  
  
  perf.open("perf.txt",ios::out);

  complex<double>** matrix;
  matrix=new complex<double>*[size];
  for(i=0;i<size;i++)
    {
      matrix[i]=new complex<double>[size];
    }


  ImageFFT fft;
  for(i=0;i<=ln;i++)
    {
      int j=pow(2,i);
      timeval t1,t2;
      gettimeofday(&t1,NULL);
      
      
     
      fft.DoImageFFT(matrix,j,j,1);
      gettimeofday(&t2,NULL);
      
      long timediff=(t2.tv_sec-t1.tv_sec)*1000000 + (t2.tv_usec-t1.tv_usec);
      
      cout << "\n" <<  j << "\t" << timediff;

      perf << "\n" <<  j << "\t" << timediff;
    }

  perf.close();

      for(i=0;i<512;i++)
	delete[] matrix[i];
      delete[] matrix;
      
}




void freemem()
{

    // init memory for image
  int i;
  

  for(i=0;i<squaresize;i++)
    {
      delete[] inpixels[i];
      delete[] fftout[i];
    }

  delete[] inpixels;
  delete[] fftout;

  
}

void allocmem()
{

  int i;
  
  if(imageloaded)
    {

      // assume square image
      squaresize=r >c ? r : c;
      squaresize=(int)pow(2,(int)(log(squaresize)/log(2))+1);


      // init memory for image
      inpixels=new complex<double>*[squaresize];
      for(i=0;i<squaresize;i++)
	inpixels[i]=new complex<double>[squaresize];


    //alloc memory
      fftout=new double*[squaresize];
      for(i=0;i<squaresize;i++)

	{
	  //	  cout << i << flush;
	  fftout[i]=new double[squaresize];

	}

    }
  
}


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

  int choice;

  do

  {

    cout << "\n\n";
    cout << "             MAIN MENU             \n";
    cout << "------------------------------------";

    cout << "\n\n\n";
    cout << "1. Load image and compute Fourier Transform\n";

    cout << "2. Display image\n";
    cout << "3. Dynamic Range Compression for fourier transform\n";
    cout << "4. Ideal Low Pass\n";
    cout << "5. Ideal High Pass\n";
    cout << "6. Butterworth Low Pass\n";
    cout << "7. Butterworth High Pass\n";
    cout << "8. Do Translation for Fourier Transform\n";
    cout << "9. Do Rotation for Fourier Transform\n";
    cout << "10. Write Performance report for 2D FFT\n";
    cout << "15. Exit\n";

    cout << "\n\nEnter choice : ";

    cin >> choice;

    char filename[80];
    ImageFFT fft;
    double h;
    int mx,my;
    double f;
    int n;
    

    switch(choice)
      {

      case 1:


	cout << "\nEnter filename : ";
	cin >> filename;
	loadimage(filename);

	if(imageloaded)
	  freemem();
	else
	  imageloaded=true;


	allocmem();

	pixelsToComplex();

	//___Init data

	break;

      case 2:
	if(imageloaded)
	  image->display();
	else
	  cout << "\nNo image loaded...";
    

	break;
	
		

      case 3:
	
	cout << "Doing Fourier Transform..." << flush;

	startwatch();
	pixelsToComplex();
	fft.DoImageFFT(inpixels,squaresize,squaresize,1);

	stopwatch();
	
	//cout << "\nEnter filename : ";
	//cin >> filename;
	displayFourier(squaresize);
	break;
	

      case 4:

	
	cout << "Doing Fourier Transform..." << flush;

	

	
	cout << "\nLow Pass Filter\n";
	cout << "\nEnter cutoff frequency : ";
	cin >> f;
	

	startwatch();
	
	pixelsToComplex();
	
	fft.DoImageFFT(inpixels,squaresize,squaresize,1);
	
	mx=squaresize/2,my=squaresize/2;
  
  
	for(int i=0;i<squaresize;i++)
	  for(int j=0;j<squaresize;j++)
	    if (distance(j,i,mx,my) > f)
	      inpixels[i][j]= complex<double>(0,0);

	stopwatch();

	cout << "\nDisplaying Fourier..." << flush;
	displayFourier(squaresize);

	
	

	//__Inverse FFT
  
	fft.DoImageFFT(inpixels,squaresize,squaresize,-1);
	
	for(int i=0;i<squaresize;i++)
	  for(int j=0;j<squaresize;j++)
	    fftout[i][j] =  pow(-1,i+j) * inpixels[i][j].real();

	//fftout[i][j] = inpixels[i][j].real();

	    
	cout << "\nDisplaying image..." << flush;
	displayimage(fftout,r,c);
	  writeimage("",fftout,r,c,1);
	
	
	break;
     case 5:

	
	cout << "Doing Fourier Transform..." << flush;


	
	cout << "\nIdeal High Pass Filter\n";
	cout << "\nEnter cutoff frequency : ";
	cin >> f;
	

	startwatch();
	
	pixelsToComplex();
	
	fft.DoImageFFT(inpixels,squaresize,squaresize,1);

	mx=squaresize/2,my=squaresize/2;
  
  
	for(int i=0;i<squaresize;i++)
	  for(int j=0;j<squaresize;j++)
	    if (distance(j,i,mx,my) < f)
	      inpixels[i][j]= complex<double>(0,0);


	stopwatch();
	
	cout << "\nDisplaying Fourier..." << flush;
	displayFourier(squaresize);

	//__Inverse FFT
  
	fft.DoImageFFT(inpixels,squaresize,squaresize,-1);
	
	for(int i=0;i<squaresize;i++)
	  for(int j=0;j<squaresize;j++)
	    fftout[i][j] =  pow(-1,i+j) * inpixels[i][j].real();

	//fftout[i][j] = inpixels[i][j].real();

	    
	cout << "\nDisplaying image..." << flush;
	displayimage(fftout,r,c);
	  writeimage("",fftout,r,c,1);
	
	
	break;
     case 6:

	
	cout << "Doing Fourier Transform..." << flush;
	
	cout << "\nButterworth Low Pass Filter\n";
	cout << "\nEnter cutoff frequency : ";
	cin >> f;
	cout << "\nEnter order : ";
	cin >> n;

	startwatch();
	
	pixelsToComplex();
	
	fft.DoImageFFT(inpixels,squaresize,squaresize,1);


	
	
	
	mx=squaresize/2,my=squaresize/2;

	
  
	for(int i=0;i<squaresize;i++)
	  for(int j=0;j<squaresize;j++)

	    {
	      h = 1/(1 + 0.414*exp(2*n*log(distance(j,i,mx,my)/f)));
		inpixels[i][j]*= h;

	    }

	stopwatch();


	
	cout << "\nDisplaying Fourier..." << flush;
	displayFourier(squaresize);

	//__Inverse FFT
  
	fft.DoImageFFT(inpixels,squaresize,squaresize,-1);
	
	for(int i=0;i<squaresize;i++)
	  for(int j=0;j<squaresize;j++)
	    fftout[i][j] =  pow(-1,i+j) * inpixels[i][j].real();

	//fftout[i][j] = inpixels[i][j].real();

	    
	cout << "\nDisplaying image..." << flush;
	displayimage(fftout,r,c);
	  writeimage("",fftout,r,c,1);
	
	
	break;
     case 7:

	
	cout << "Doing Fourier Transform..." << flush;

	
	cout << "\nButterworth High Pass Filter\n";
	cout << "\nEnter cutoff frequency : ";
	cin >> f;
	cout << "\nEnter order : ";
	cin >> n;

	startwatch();
	pixelsToComplex();
	
	fft.DoImageFFT(inpixels,squaresize,squaresize,1);

	

	mx=squaresize/2,my=squaresize/2;
  
  
	for(int i=0;i<squaresize;i++)
	  for(int j=0;j<squaresize;j++)
	    {
	      h = 1/(1 + 0.414*exp(2*n*log(f/distance(j,i,mx,my))));
	      inpixels[i][j]*=h;
	    }

	stopwatch();
	
	cout << "\nDisplaying Fourier..." << flush;
	displayFourier(squaresize);

	//__Inverse FFT
  
	fft.DoImageFFT(inpixels,squaresize,squaresize,-1);
	
	for(int i=0;i<squaresize;i++)
	  for(int j=0;j<squaresize;j++)
	    fftout[i][j] =  pow(-1,i+j) * inpixels[i][j].real();

	//fftout[i][j] = inpixels[i][j].real();

	    
	cout << "\nDisplaying image..." << flush;
	displayimage(fftout,r,c);
	  writeimage("",fftout,r,c,1);
	
	  break;
      case 8:


	{

	  cout << "Doing Fourier Transform..." << flush;


	  int x,y;
	  complex<double> h;

	  cout << "\nTranslate using fourier transform";
	  cout << "\nEnter translation coord : ";
	  cin >> x >> y;

	  
	  pixelsToComplex();
	
	  fft.DoImageFFT(inpixels,squaresize,squaresize,1);

	
	

	
	  mx=squaresize/2,my=squaresize/2;
  
  
	  for(int i=0;i<squaresize;i++)
	    for(int j=0;j<squaresize;j++)
	      {
		h =complex<double>(cos(-2*pi*(j*x+i*y)/squaresize),sin(-2*pi*(j*x+i*y)/squaresize));
		inpixels[i][j]*=h;
	      }



	  stopwatch();
	  
	  cout << "\nDisplaying Fourier..." << flush;
	  displayFourier(squaresize);

	  //__Inverse FFT
  
	  fft.DoImageFFT(inpixels,squaresize,squaresize,-1);
	
	  for(int i=0;i<squaresize;i++)
	    for(int j=0;j<squaresize;j++)
	      fftout[i][j] =  pow(-1,i+j) * inpixels[i][j].real();

	  //fftout[i][j] = inpixels[i][j].real();

	    
	  cout << "\nDisplaying image..." << flush;
	  displayimage(fftout,r,c);
	  writeimage("",fftout,r,c,1);
	}
	break;
      case 9:

	{
	
	  double angle;
	  cout << "\nRotation of image and Its Fourier Transform\n";
	  cout << "\nEnter rotation angle for image : ";
	  cin >> angle;


	  cout << "\nInitial fourier transform..." << flush;

	  pixelsToComplex();
	
	  fft.DoImageFFT(inpixels,squaresize,squaresize,1);
	  displayFourier(squaresize);
	  

	  cout << "\nAfter rotation..." << flush;

	  cout << "\nImage..." << flush;

	  image->matteColor(ColorGray(0));
	  
	  image->matte(false);
	  
	  image->rotate(angle);
	  int tr=image->rows();
	  int tc=image->columns();
	  //cout << tr << "  " << tc << "  " << flush;
	  image->matteColor(ColorGray(0));

	  image->matte(false);

	    image->matteColor(ColorGray(0));
	  // crop it
	  char geom[80];
	  sprintf(geom,"%dx%d+%d+%d",c,r,(tc-c)/2,(tr-r)/2);
	  //cout << geom << flush;
	  image->crop(geom);
	 
	  image->display();

	  Pixels view(*image);
	  pixels=view.get(0,0,c,r);

	  cout << "\nFFT..." << flush;
	  
	  pixelsToComplex();
	
	  fft.DoImageFFT(inpixels,squaresize,squaresize,1);

	  displayFourier(squaresize);
	  
	 
	  break;
	}
	

      case 10:
	measure_performance();
	break;
	
      case 15:
	exit(0);
      default:
	break;
      }
    
  }
  while(1);

  
    
     
}