//=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*
//   Program fraunhofer
//=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*

#include "cip.h"
#include "complex.h"

//---- definition of fundermental constants
//---- experimental settings
const int N = 256;                 // size of image

//---- prototypes of functions
void Load( char* fname, Complex h[N][N] );
void Save( char* fname, Complex H[N][N], double amp );
void FFT1( int sign, Complex src[N], Complex dst[N] );
void FFT2( int sign, Complex src[N][N], Complex dst[N][N] );

//---- tables
class ExpTable
{
  Complex expo[2*N];
public:
  ExpTable( void ){
    for( int p=-N+1 ; p<=N-1 ; p++ ){
      expo[p+N] = Complex( cos(2*M_PI/N*p), sin(2*M_PI/N*p) );
    }
  }
  inline Complex operator [] ( int p ){
    return expo[p+N];
  }
} Exp;

class RevTable
{
  int reve[N];
public:
  RevTable( void ){
    for( int i=0 ; i<N ; i++ ){
      int rev=0;
      for( u_int mask1=1, mask2=N/2 ; mask2 ; mask1<<=1, mask2>>=1 ){
	if( i&mask1 ) rev|=mask2;
      }
      reve[i] = rev;
    }
  }
  inline int operator [] ( int i ){
    return reve[i];
  }
} Rev;

//---- main function
int main( void )
{
  static Complex h[N][N], H[N][N];
  char buf[32], film[32], screen[32];

  printf("The input file name of the film image    = ");
  fgets( buf, 32, stdin );
  if( 1 != sscanf( buf, "%s", film ) ) return 1;

  printf("The output file name of the screen image = ");
  fgets( buf, 32, stdin );
  if( 1 != sscanf( buf, "%s", screen ) ) return 1;

  printf("Lighting up the film....\n");
  Load( film, H );
  printf("Propagating wave from the film to the screen....\n");
  FFT2( -1, H, h );
  printf("Taking this picture to the another film\n");
  Save( screen, h, 1.0 );

  return 0;
}

//---- load pgm file(raw) file to the given complex array  
void Load( char* fname, Complex h[N][N] )
{
  FILE* fptr = fopen( fname, "rb" );
  if( fptr == NULL ){
    fprintf( stderr, "File not found!\n" );
    exit(1);
  }
  char header[64];
  for( int significant=0 ; significant<3 ; ){
    fgets( header, 64, fptr );
    if( header[0] != '#' ) significant++;
  }

  u_char buf[N*N];
  fread( buf, sizeof(u_char), N*N, fptr );

  u_char* ptr = buf;
  for( int i=0 ; i<N ; i++ ){
    for( int j=0 ; j<N ; j++ ){
      h[i][j].real = double(*ptr++)/255.0;
      h[i][j].imag = 0.0;
    }
  }
  fclose( fptr );
}

//---- save pgm file(raw) file from the given complex array  
void Save( char* fname, Complex H[N][N], double amp )
{
  FILE* fptr = fopen( fname, "wb" );
  if( fptr == NULL ){
    fprintf( stderr, "File not found!\n" );
    exit(1);
  }

  char header[64];
  sprintf( header, "P5\n%d %d\n255\n", N, N );
  fputs( header, fptr );

  u_char buf[N*N];
  u_char* ptr = buf;

  double bright;
  for( int i=0 ; i<N ; i++ ){
    for( int j=0 ; j<N ; j++ ){
      bright = amp*norm( H[i][j] );
      if( bright > 1.0 ) bright = 1.0;
      if( bright < 0.0 ) bright =   0;
      *ptr++ = (u_char)(255.0*bright);
    }
  }

  fwrite( buf, sizeof(u_char), N*N, fptr );
  fclose( fptr );
}




//---- 1-Dimensional Fast Fourier Transform (not used in this program)
void FFT1( int sign, Complex src[N], Complex dst[N] )
{
  int k;
    // copy from source array at reversed index to destination array
  for( int j=0 ; j<N ; j++ ){
    dst[Rev[j]] = src[j];
  }
    // x=0, k=0 are treated as the center of the arrays.
  for( k=N/2 ; k<N ; k++ ) dst[k] *= -1;

    // Danielson-Lanczos's algorithm
  int k0, k1, k2, dk, p, dp;
  Complex ex, work;
  for( dk=0x01, dp=N/2 ; dp ; dk<<=1, dp>>=1 ){
    for( k0=0, p=0 ; k0<dk ; k0++, p+=sign*dp ){
      for( k1=k0, ex=Exp[p] ; k1<N ; k1=k2+dk ){
        k2=k1+dk;
        work     = dst[k2];
        work    *= ex;
        dst[k2]  = dst[k1];
        dst[k2] -= work;
        dst[k1] += work;
      }
    }
  }
    // x=0, k=0 are treated as the center of the arrays.
  for( k=1 ; k<N ; k+=2 ) dst[k] *= -1;
}

//---- 2-Dimensional Fast Fourier Transform
void FFT2( int sign, Complex src[N][N], Complex dst[N][N] )
{
  int kx, ky;
    // copy from source array at reversed index to destination array
  for( int jx=0 ; jx<N ; jx++ ){
    for( int jy=0 ; jy<N ; jy++ ){
      dst[Rev[jx]][Rev[jy]] = (1.0/N)*src[jx][jy];
    }
  }
    // x=0, k=0 are treated as the center of the arrays.
  for( kx=N/2 ; kx<N ; kx++ ){
    for( ky=0 ; ky<N/2 ; ky++ ){
      dst[kx][ky] *= -1.0;
      dst[ky][kx] *= -1.0;
    }
  }
    // Danielson-Lanczos's algorithm
  int k0, k1, k2, dk, p, dp;
  Complex ex, work;
  for( dk=0x01, dp=N/2 ; dp ; dk<<=1, dp>>=1 ){
    for( k0=0, p=0 ; k0<dk ; k0++, p+=sign*dp ){
      for( k1=k0, ex=Exp[p] ; k1<N ; k1=k2+dk ){
        for( ky=0, k2=k1+dk ; ky<N ; ky++ ){      // transform columns
          work         = dst[k2][ky];
          work        *= ex;
          dst[k2][ky]  = dst[k1][ky];
          dst[k1][ky] += work;
          dst[k2][ky] -= work;
        }
        for( kx=0, k2=k1+dk ; kx<N ; kx++ ){      // transform rows
          work         = dst[kx][k2];
          work        *= ex;
          dst[kx][k2]  = dst[kx][k1];
          dst[kx][k1] += work;
          dst[kx][k2] -= work;
        }
      }
    }
  }
    // x=0, k=0 are treated as the center of the arrays.
  for( kx=0 ; kx<N ; kx++ ){
    for( ky=kx+1 ; ky<N ; ky+=2 ){
      dst[kx][ky] *= -1.0;
      dst[ky][kx] *= -1.0;
    }
  }
}