#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;

PixelPacket* pixels;
int c,r;
complex<double> inpixels[512][512];
complex<double> ifft[512][512];

int squaresize;
bool imageloaded=false;
double fftout[512][512];

Image* image;
Pixels* view;

complex<double>* t,t1;
double* s,s1;



timeval ts,te;


void startwatch()
{

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

}

void stopwatch()
{

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


      
}



class ImageFFT
{

  
private:

  complex<double> a[512];
  complex<double> A[512];
  
  // 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(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(int n,int mode=1)
  {
    bitrev(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=2 << (s-1);
	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(int x,int y,int mode)
  {
    
    //__Compute FFT for each column first___
    //256 times
    for(int i=0;i<x;i++)
      {

	for(int j=0;j<y;j++)
	  a[j]=ifft[j][i];
	
	fft(y,mode);
	
	for(int j=0;j<y;j++)
	  ifft[j][i]=A[j];
      }
   

    // 256 times
    for(int i=0;i<y;i++)
      {
	for(int j=0;j<x;j++)
	 a[j]=ifft[i][j];
	  
	fft(x,mode);
	  
	for(int j=0;j<x;j++)
	  ifft[i][j]=A[j];
      }

  }

  
  

  
 // mode=-1 for inverse , =1 for forward FFT
  void DoIImageFFT(int x,int y,int mode)
  {

   
    //__Compute FFT for each column first___
    //256 times
    for(int i=0;i<x;i++)
      {

	for(int j=0;j<y;j++)
	  a[j]=inpixels[j][i];
	
	fft(y,mode);
	
	for(int j=0;j<y;j++)
	  inpixels[j][i]=A[j];
      }
   

    // 256 times
    for(int i=0;i<y;i++)
      {
	for(int j=0;j<x;j++)
	 a[j]=inpixels[i][j];
	  
	fft(x,mode);
	  
	for(int j=0;j<x;j++)
	  inpixels[i][j]=A[j];
      }

  }


  
  
    void copyFFT(int x,int y)
      {
	int i,j;
	for(i=0;i<y;i++)
	  for(j=0;j<x;j++)
	    inpixels[i][j]=ifft[i][j];


	
      }
    
   

  
  
  

};



void writeimage(char* filename,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(fftout[x][y]);
	fpixels++;
      }
    
  fview.sync();
  fspectrum.write((ask==1) ? buffer : filename);

}


void displayimage(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(fftout[x][y]);
	fpixels++;
      }
    
  fview.sync();
  fspectrum.display();


}





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

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

      
  imageloaded=true;

}

void pixelsToComplex()
{

 
  int i,j; 

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

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

}



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(squaresize,squaresize);

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

}

inline 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);


  
  



  ImageFFT fft;
  for(i=0;i<=ln;i++)
    {
      int j=pow(2,i);
      timeval t1,t2;
      gettimeofday(&t1,NULL);
      
      
     
      fft.DoImageFFT(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();

    
      
}






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

	
	imageloaded=true;


	
	//pre compute fft
	
	pixelsToComplex();

	fft.DoImageFFT(squaresize,squaresize,1);
	

	

	//pixelsToComplex();

	//___Init data

	
	break;

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

	break;
	
		

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

	startwatch();
	

	

	fft.copyFFT(squaresize,squaresize);
	

	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();
	

	fft.copyFFT(squaresize,squaresize);
	

	
	
	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

	startwatch();
	
	fft.DoIImageFFT(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();


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

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

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


	startwatch();


	fft.copyFFT(squaresize,squaresize);

	

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

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


	startwatch();
	
	//__Inverse FFT

	
	fft.DoIImageFFT(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();

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

	    
	cout << "\nDisplaying image..." << flush;
	displayimage(r,c);
	  writeimage("",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();
	


	fft.copyFFT(squaresize,squaresize);
	
	
	
	
	
	
	
	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

	startwatch();
	
	fft.DoIImageFFT(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();

	stopwatch();
	
	    
	cout << "\nDisplaying image..." << flush;
	displayimage(r,c);
	  writeimage("",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();

	fft.copyFFT(squaresize,squaresize);

	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

	startwatch();
	
	fft.DoIImageFFT(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();


	stopwatch();
	
	cout << "\nDisplaying image..." << flush;
	displayimage(r,c);
	  writeimage("",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;


	  startwatch();
	  

	  fft.copyFFT(squaresize,squaresize);
	  
	  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
  
	  startwatch();
	  
	  fft.DoIImageFFT(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();

	  stopwatch();

	  
	    
	  cout << "\nDisplaying image..." << flush;
	  displayimage(r,c);
	  writeimage("",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;

	  fft.copyFFT(squaresize,squaresize);
	  

	  



	  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(squaresize,squaresize,1);

	  fft.copyFFT(squaresize,squaresize);
	  
	  displayFourier(squaresize);
	  
	 
	  break;
	}
	

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

  
    
     
}