/*
 *  realce_espectro.c
 *
 *  Este programa ilustra as operacoes de processamento
 *  espectral para realcar uma imagem com dimensoes em potencia de 2.
 *  Utiliza a transformada de Fourier rapida para
 *  converter os dados entre dominio espacial e
 *  dominio da frequencia
 *
 *   'c':Alterna compressao da escala dinamica
 *   'a':Filtro Passa-Altas de Butterworth
 *   'b':Filtro Passa-Baixas de Butterworth
 *   'r':Reseta
 *   '+':Incrementa raio do filtro
 *   '-':Decrementa raio do filtro
 *   'i':Incrementa a ordem do filtro
 *   'd':Decrementa a ordem do filtro
 *
 *            Ting (01/11/08)
 */

#define _USE_MATH_DEFINES

#include <GL/glut.h>    // Header File For The GLUT Library
#include <stdio.h>      // Header file for standard file i/o.
#include <stdlib.h>     // Header file for malloc/free.
#include <math.h>       // Header file for math functions.

/* Include here the header of the image to be processed
   width and height are defined */
#include "child_profile_256.h"

/* ascii code for the escape key */
#define ESCAPE 27
#define TRUE 1
#define FALSE 0

typedef struct COMPLEX_struct {   /* pontos */
  double real, imag;
} COMPLEX;

int win;  // id da janela principal
int swin1, swin2, swin3, swin4; // id das subjanelas do glut

unsigned int c = 1; //compressao da escala do espectro
int radius = 200;
int n = 1;
GLubyte pixel[3];
char *data;
float *fourier;
GLubyte *imagem;
GLubyte **aux_imagem;
GLubyte *fft_imagem, *f_fft_imagem, *ifft_imagem;
COMPLEX **fft_aux, **ifft_fourier;

char opt='b';
/***********************************************************/
/* Processamento Espectral
   Perform a 2D FFT inplace given a complex 2D array
   The direction dir, 1 for forward, -1 for reverse
   The size of the array (nx,ny)
   Return false if there are memory problems or
      the dimensions are not powers of 2

Source: http://local.wasp.uwa.edu.au/~pbourke/other/dft/

*/

unsigned int Powerof2(int nx, int *m, int *twopm)
{
  *m=0;
  *twopm=1;

  while (*twopm < nx) {
    (*m)++;
    *twopm = (*twopm)*2;
  }

  if (*twopm == nx) return(TRUE);
  return(FALSE);
}

/*******************************************************
   This computes an in-place complex-to-complex FFT
   x and y are the real and imaginary arrays of 2^m points.
   dir =  1 gives forward transform
   dir = -1 gives reverse transform

     Formula: forward
                  N-1
                  ---
              1   \          - j k 2 pi n / N
      X(n) = ---   >   x(k) e                    = forward transform
              N   /                                n=0..N-1
                  ---
                  k=0

      Formula: reverse
                  N-1
                  ---
                  \          j k 2 pi n / N
      X(n) =       >   x(k) e                    = forward transform
                  /                                n=0..N-1
                  ---
                  k=0
*/
int FFT(int dir,int m,double *x,double *y)
{
   long nn,i,i1,j,k,i2,l,l1,l2;
   double c1,c2,tx,ty,t1,t2,u1,u2,z;

   /* Calculate the number of points */
   nn = 1;
   for (i=0;i<m;i++)
      nn *= 2;

   /* Do the bit reversal */
   i2 = nn >> 1;
   j = 0;
   for (i=0;i<nn-1;i++) {
      if (i < j) {
         tx = x[i];
         ty = y[i];
         x[i] = x[j];
         y[i] = y[j];
         x[j] = tx;
         y[j] = ty;
      }
      k = i2;
      while (k <= j) {
         j -= k;
         k >>= 1;
      }
      j += k;
   }

   /* Compute the FFT */
   c1 = -1.0;
   c2 = 0.0;
   l2 = 1;
   for (l=0;l<m;l++) {
      l1 = l2;
      l2 <<= 1;
      u1 = 1.0;
      u2 = 0.0;
      for (j=0;j<l1;j++) {
         for (i=j;i<nn;i+=l2) {
            i1 = i + l1;
            t1 = u1 * x[i1] - u2 * y[i1];
            t2 = u1 * y[i1] + u2 * x[i1];
            x[i1] = x[i] - t1;
            y[i1] = y[i] - t2;
            x[i] += t1;
            y[i] += t2;
         }
         z =  u1 * c1 - u2 * c2;
         u2 = u1 * c2 + u2 * c1;
         u1 = z;
      }
      c2 = sqrt((1.0 - c1) / 2.0);
      if (dir == 1)
         c2 = -c2;
      c1 = sqrt((1.0 + c1) / 2.0);
   }

   /* Scaling for forward transform */
   if (dir == 1) {
      for (i=0;i<nn;i++) {
         x[i] /= (double)nn;
         y[i] /= (double)nn;
      }
   }
   return(TRUE);
}

int FFT2D(COMPLEX **c,int nx,int ny,int dir)
{
   int i,j;
   int m,twopm;
   double *real,*imag;

   /* Transform the rows */
   real = (double *)malloc(nx * sizeof(double));
   imag = (double *)malloc(nx * sizeof(double));
   if (real == NULL || imag == NULL)
      return(FALSE);
   if (!Powerof2(nx,&m,&twopm) || twopm != nx) {
       printf("Imagem precisa ter dimensoes em potencia de 2!!\n");
       return(FALSE);
   }

   for (j=0;j<ny;j++) {
      for (i=0;i<nx;i++) {
         real[i] = c[i][j].real;
         imag[i] = c[i][j].imag;
      }
      FFT(dir,m,real,imag);
      for (i=0;i<nx;i++) {
         c[i][j].real = real[i];
         c[i][j].imag = imag[i];
      }
   }
   free(real);
   free(imag);

   /* Transform the columns */
   real = (double *)malloc(ny * sizeof(double));
   imag = (double *)malloc(ny * sizeof(double));
   if (real == NULL || imag == NULL)
      return(FALSE);
   if (!Powerof2(ny,&m,&twopm) || twopm != ny)
      return(FALSE);
   for (i=0;i<nx;i++) {
      for (j=0;j<ny;j++) {
         real[j] = c[i][j].real;
         imag[j] = c[i][j].imag;
      }
      FFT(dir,m,real,imag);
      for (j=0;j<ny;j++) {
         c[i][j].real = real[j];
         c[i][j].imag = imag[j];
      }
   }
   free(real);
   free(imag);

   return(TRUE);
}
/****************************************/
// Transformada de Fourier Rapida
void load_complex_image(void)
{
  int i, j;

   for (i = 0; i < height; i++) {
      for (j = 0; j < width; j++) {
	fft_aux[i][j].real = aux_imagem[i][j];
	fft_aux[i][j].imag = 0;
      }
   }
}

void fft()
{
  int u,v;

  load_complex_image();

  FFT2D(fft_aux, width, height, 1); // FFT direta

  float max=-1.E12, min=1E12;

  for (u=0; u<height; u++) {
    for (v=0; v<width; v++) {
       ifft_fourier[u][v].real = fft_aux[u][v].real;
       ifft_fourier[u][v].imag = fft_aux[u][v].imag;
       fourier[u*width+v] = sqrt(fft_aux[u][v].real*fft_aux[u][v].real + fft_aux[u][v].imag*fft_aux[u][v].imag);
      if (fourier[u*width+v] > max) max = fourier[u*width+v];
      if (fourier[u*width+v] < min) min = fourier[u*width+v];
    }
  }

  // Normalize and shift the values to display them symmetrically

  if (c == 0) {
    for (u=0; u<height; u++) {
      for (v=0; v<width; v++) {
	fft_imagem[((u+height/2)%height)*width+(v+width/2)%width] = (GLubyte)(255.*fourier[u*width+v]/max);
      }
    }
  } else {
    max = log(1+max);
    for (u=0; u<height; u++) {
      for (v=0; v<width; v++) {
	fft_imagem[((u+height/2)%height)*width+(v+width/2)%width] = (GLubyte)(255.*log(1+fourier[u*width+v])/max);
      }
    }
  }
}
/*********************************************************/
//Filtro Passa-Baixas de Butterworth (FPBB)
void FPBB (int raio)
{
  int u,v,R;
  double h;

  for (u = 0; u < height; u++) {
    for (v = 0; v < width; v++) {
	fft_aux[u][v].real = ifft_fourier[u][v].real;
	fft_aux[u][v].imag = ifft_fourier[u][v].imag;
      }
   }

  // Filtragem
  R = 2*n;
  float max=-1.E12, min=1E12;
  for (u = 0; u < height; u++) {
    for (v = 0; v < width; v++) {
      int uu = (u+height/2)%height;
      int vv = (v+width/2)%width;
      if (raio == 0) h =0.0;
      else
	h = 1/(1 + 0.001*pow(((uu-height/2)*(uu-height/2) + (vv-width/2)*(vv-width/2))/raio,R));
      fft_aux[u][v].real = h*fft_aux[u][v].real;
      fft_aux[u][v].imag = h*fft_aux[u][v].imag;
      fourier[u*width+v] = sqrt(fft_aux[u][v].real*fft_aux[u][v].real + fft_aux[u][v].imag*fft_aux[u][v].imag);
      if (fourier[u*width+v] > max) max = fourier[u*width+v];
      if (fourier[u*width+v] < min) min = fourier[u*width+v];

      }
   }

  // Gerar a imagem do espectro filtrado
  if (c == 0) {
    for (u=0; u<height; u++) {
      for (v=0; v<width; v++) {
	f_fft_imagem[((u+height/2)%height)*width+(v+width/2)%width] = (GLubyte)(255.*fourier[u*width+v]/max);
      }
    }
  } else {
    max = log(1+max);
    for (u=0; u<height; u++) {
      for (v=0; v<width; v++) {
	f_fft_imagem[((u+height/2)%height)*width+(v+width/2)%width] = (GLubyte)(255.*log(1+fourier[u*width+v])/max);
      }
    }
  }

  // Transformada rapida inversa
  FFT2D(fft_aux, width, height, -1); // FFT inversa

  // Gerar a imagem no dominio espacial
  for (u=0; u<height; u++) {
    for (v=0; v<width; v++) {
      ifft_imagem[u*width+v] = sqrt(fft_aux[u][v].real*fft_aux[u][v].real + fft_aux[u][v].imag*fft_aux[u][v].imag);
      if (sqrt(fft_aux[u][v].real*fft_aux[u][v].real + fft_aux[u][v].imag*fft_aux[u][v].imag) > 255) ifft_imagem[u*width+v] = 255;
    }
  }

}

// Filtro Passa-Altas de Butterworth (FPAB)
void FPAB (int raio)
{
  int u,v,R;
  double h;

  for (u = 0; u < height; u++) {
    for (v = 0; v < width; v++) {
	fft_aux[u][v].real = ifft_fourier[u][v].real;
	fft_aux[u][v].imag = ifft_fourier[u][v].imag;
      }
   }

  // Filtragem
  R = 2*n;
  float max=-1.E12, min=1E12;
  for (u = 0; u < height; u++) {
    for (v = 0; v < width; v++) {
      int uu = (u+height/2)%height;
      int vv = (v+width/2)%width;
      if (u==0 && v==0) h = 0.0;
      else
	h = 1/(1 + 0.414*pow(raio/((uu-height/2)*(uu-height/2) + (vv-width/2)*(vv-width/2)),R));
      fft_aux[u][v].real = h*fft_aux[u][v].real;
      fft_aux[u][v].imag = h*fft_aux[u][v].imag;
      fourier[u*width+v] = sqrt(fft_aux[u][v].real*fft_aux[u][v].real + fft_aux[u][v].imag*fft_aux[u][v].imag);
      if (fourier[u*width+v] > max) max = fourier[u*width+v];
      if (fourier[u*width+v] < min) min = fourier[u*width+v];

      }
   }

  // Gerar a imagem do espectro filtrado
  if (c == 0) {
    for (u=0; u<height; u++) {
      for (v=0; v<width; v++) {
	f_fft_imagem[((u+height/2)%height)*width+(v+width/2)%width] = (GLubyte)(255.*fourier[u*width+v]/max);
      }
    }
  } else {
    max = log(1+max);
    for (u=0; u<height; u++) {
      for (v=0; v<width; v++) {
	f_fft_imagem[((u+height/2)%height)*width+(v+width/2)%width] = (GLubyte)(255.*log(1+fourier[u*width+v])/max);
      }
    }
  }

  // Transformada rapida inversa

  FFT2D(fft_aux, width, height, -1); // FFT inversa

  // Gerar a imagem no dominio espacial
  for (u=0; u<height; u++) {
    for (v=0; v<width; v++) {
      ifft_imagem[u*width+v] = sqrt(fft_aux[u][v].real*fft_aux[u][v].real + fft_aux[u][v].imag*fft_aux[u][v].imag);
      if (sqrt(fft_aux[u][v].real*fft_aux[u][v].real + fft_aux[u][v].imag*fft_aux[u][v].imag) > 255) ifft_imagem[u*width+v] = 255;
    }
  }

}

/***********************************************************/
void init(void)
{
  int i;

   /* definir a cor em preto (0,0,0) como cor de "CLEAR" */
   glClearColor (0.0, 0.0, 0.0, 0.0);
  /* Space for complex image */
  fft_aux = (COMPLEX **) malloc(sizeof(COMPLEX*)*height);
  if (fft_aux == NULL) {
	printf("Error allocating memory for complex image data");
        exit(0);
	return;
  }
  for (i=0; i< height; i++) {
    fft_aux[i] = (COMPLEX *) malloc(sizeof(COMPLEX)*width);
    if (fft_aux[i] == NULL) {
	printf("Error allocating memory for complex image data");
        exit(0);
	return;
    }
  }

  ifft_fourier = (COMPLEX **) malloc(sizeof(COMPLEX*)*height);
  if (ifft_fourier == NULL) {
	printf("Error allocating memory for complex image data");
        exit(0);
	return;
  }

  for (i=0; i< height; i++) {
    ifft_fourier[i] = (COMPLEX *) malloc(sizeof(COMPLEX)*width);
    if (ifft_fourier[i] == NULL) {
	printf("Error allocating memory for complex image data");
        exit(0);
	return;
    }
  }

  /* Space for images*/
  fft_imagem = (GLubyte *) malloc(sizeof(GLubyte)*width*height);
  if (fft_imagem == NULL) {
	printf("Error allocating memory for fourier transform image data");
        exit(0);
	return;
  }
  f_fft_imagem = (GLubyte *) malloc(sizeof(GLubyte)*width*height);
  if (fft_imagem == NULL) {
	printf("Error allocating memory for fourier transform image data");
        exit(0);
	return;
  }
  ifft_imagem = (GLubyte *) malloc(sizeof(GLubyte)*width*height);
  if (ifft_imagem == NULL) {
	printf("Error allocating memory for fourier transform image data");
        exit(0);
	return;
  }
  /* Fourier Transform */
  fourier = (float *) malloc(sizeof(float)*width*height);
  if (fourier == NULL) {
	printf("Error allocating memory for fourier transform data");
        exit(0);
	return;
  }

  // Transformar para dominio da frequencia
  fft();

  // Filtragem e transformar de volta para dominio espacial
  FPBB(radius);
}

void
main_display(void)
{
    glClearColor(0.8, 0.8, 0.8, 0.0);
    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
    glColor3f(0.f, 0.f, 0.f);
}

void display_swin1()
{
  // Desenhar os pixels da image1
  glutSetWindow(swin1);
  // definir a forma como os dados sao desempacotados da memoria
  glPixelStorei(GL_UNPACK_ALIGNMENT, 1);
  // Desabilita dithering
  glDisable(GL_DITHER);
  // resetar as cores nos buffers de cor
  glClear(GL_COLOR_BUFFER_BIT);

  // definir a origem da imagem
  glRasterPos2i(0,0);
  // desenhar os pixels da imagem original
  glDrawPixels(width,height,GL_LUMINANCE, GL_UNSIGNED_BYTE, imagem);

  // forcar a execucao dos comandos enviados a OpenGL
  glFlush();
}

void display_swin2()
{
  // Desenhar o espectro de frequencia da image1
  glutSetWindow(swin2);
  glPixelStorei(GL_UNPACK_ALIGNMENT, 1);
  // Desabilita dithering
  glDisable(GL_DITHER);

  // resetar as cores nos buffers de cor
  glClear(GL_COLOR_BUFFER_BIT);

  glRasterPos2i(0,0);
  // desenhar os pixels da imagem no dominio da frequencia
  glDrawPixels(width,height,GL_LUMINANCE, GL_UNSIGNED_BYTE, fft_imagem);

  // forcar a execucao dos comandos enviados a OpenGL
  glFlush();
}

void display_swin3()
{
  // Desenhar o espectro de frequencia apos filtragem
  glutSetWindow(swin3);
  glPixelStorei(GL_UNPACK_ALIGNMENT, 1);
  // Desabilita dithering
  glDisable(GL_DITHER);

  // resetar as cores nos buffers de cor
  glClear(GL_COLOR_BUFFER_BIT);

  glRasterPos2i(0,0);
  // desenhar os pixels da imagem no dominio da frequencia
  glDrawPixels(width,height,GL_LUMINANCE, GL_UNSIGNED_BYTE, f_fft_imagem);

  // forcar a execucao dos comandos enviados a OpenGL
  glFlush();
}

void display_swin4()
{
  // Desenhar os pixels da image2
  glutSetWindow(swin4);
  glPixelStorei(GL_UNPACK_ALIGNMENT, 1);
  // Desabilita dithering
  glDisable(GL_DITHER);

  // resetar as cores nos buffers de cor
  glClear(GL_COLOR_BUFFER_BIT);

  glRasterPos2i(0,0);
  // desenhar os pixels da imagem no dominio da frequencia
  glDrawPixels(width,height,GL_LUMINANCE, GL_UNSIGNED_BYTE, ifft_imagem);

  // forcar a execucao dos comandos enviados a OpenGL
  glFlush();
}

void main_reshape(int w, int h)
{
  glViewport(0,0,(GLsizei)w,(GLsizei)h);
  glMatrixMode(GL_PROJECTION);
  glLoadIdentity();
  gluOrtho2D(0.0, (GLfloat) w, 0.0, (GLfloat) h);
  glMatrixMode(GL_MODELVIEW);
  glLoadIdentity();

  /* redefinir a posição e as dimensões da sub-janela 1 */
  glutSetWindow(swin1);
  glutPositionWindow(0, 0);
  glutReshapeWindow(width, height);
  /* redefinir a posição e as dimensões da sub-janela 2 */
  glutSetWindow(swin2);
  glutPositionWindow(width,0 );
  glutReshapeWindow(width, height);
  /* redefinir a posição e as dimensões da sub-janela 2 */
  glutSetWindow(swin3);
  glutPositionWindow(width*2,0 );
  glutReshapeWindow(width, height);
  glutSetWindow(swin4);
  glutPositionWindow(width*3,0 );
  glutReshapeWindow(width, height);
}

void reshape_swin1(int w, int h)
{
  glMatrixMode(GL_PROJECTION);
  glLoadIdentity();
  gluOrtho2D(0.0, (GLfloat) width, 0.0, (GLfloat) height);
  glMatrixMode(GL_MODELVIEW);
  glLoadIdentity();
}

void reshape_swin2(int w, int h)
{
  glMatrixMode(GL_PROJECTION);
  glLoadIdentity();
  gluOrtho2D(0.0, (GLfloat) width, 0.0, (GLfloat) height);
  glMatrixMode(GL_MODELVIEW);
  glLoadIdentity();
}

void reshape_swin3(int w, int h)
{
  glMatrixMode(GL_PROJECTION);
  glLoadIdentity();
  gluOrtho2D(0.0, (GLfloat) width, 0.0, (GLfloat) height);
  glMatrixMode(GL_MODELVIEW);
  glLoadIdentity();
}

void reshape_swin4(int w, int h)
{
  glMatrixMode(GL_PROJECTION);
  glLoadIdentity();
  gluOrtho2D(0.0, (GLfloat) width, 0.0, (GLfloat) height);
  glMatrixMode(GL_MODELVIEW);
  glLoadIdentity();
}

/* ARGSUSED1 */
void keyboard(unsigned char key, int x, int y)
{
  int i;

   switch (key) {
   case 'c':
     // Alterna compressao da escala dinamica
     if (c == 1) c=0;
     else c=1;
     fft();
     break;
   case 'a':
     // Filtro Passa-Altas de Butterworth
     opt = 'a';
     break;
   case 'b':
     // Filtro Passa-Baixas de Butterworth
     opt = 'b';
     break;
   case 'r':
     // Reseta
     n = 1;
     radius = 25;
     break;
   case '+':
     // Incrementa raio de frequencia
     radius += 5;
     if (width > height) {
       if (radius > height/2) radius = height/2;
     } else {
       if (radius > width/2) radius = width/2;
     }
     break;
   case '-':
     // Decrementa raio de frequencia
     radius -= 5;
     if (radius < 0) radius = 0;
     break;
   case 'i':
     // Incremanta a ordem
     n += 1;
     if (n > 6) n=6;
     break;
   case 'd':
     // Decremanta a ordem
     n -= 1;
     if (n < 1) n=1;
     break;
   case 27:
     for (i=0; i<height; i++) {
       free((char *) fft_aux[i]);
       free((char *) ifft_fourier[i]);
       free((char *) aux_imagem[i]);
     }
     free((char *)fft_aux);
     free((char *)ifft_fourier);
     free((char *)fft_imagem);
     free((char *)f_fft_imagem);
     free((char *)ifft_imagem);
     free((char *)fourier);
     free((char *)imagem);
     free((char *)aux_imagem);
     exit(0);
     return;
   }
   if (opt == 'a')
       FPAB(radius);
   else
       FPBB(radius);
   display_swin2();
   display_swin3();
   display_swin4();
}

int main(int argc, char **argv)
{
   int i, j;

  /* Alocar memoria para imagens */
  imagem = (GLubyte *) malloc (sizeof(GLubyte)*height*width);
  aux_imagem = (GLubyte **) malloc (sizeof(GLubyte *)*height);
  for (i=0; i<height; i++) {
    aux_imagem[i] = (GLubyte *)malloc(sizeof(GLubyte)*width);
  }

  /* Extrair imagem */
  data = header_data;
  for (i=(height-1); i >= 0; i--)
    for (j=(width-1); j >= 0; j--) {
      HEADER_PIXEL(data,pixel)
      imagem[i*width+j] = pixel[0];
      aux_imagem[i][j]=pixel[0];
    }

  glutInit(&argc, argv);
  glutInitDisplayMode(GLUT_SINGLE | GLUT_RGB );
  glutInitWindowSize((int)width*4, (int)height);
  glutInitWindowPosition(0,0);
  win = glutCreateWindow(argv[0]);
    glutReshapeFunc(main_reshape);
    glutDisplayFunc(main_display);
    glutKeyboardFunc(keyboard);
    swin1 = glutCreateSubWindow(win,0,0,(int)width,(int)height);
    glutReshapeFunc(reshape_swin1);
    glutDisplayFunc(display_swin1);
    glutKeyboardFunc(keyboard);
    swin2 = glutCreateSubWindow(win,(int)width,0,(int)width,(int)height);
    glutReshapeFunc(reshape_swin2);
    glutDisplayFunc(display_swin2);
    glutKeyboardFunc(keyboard);
    swin3 = glutCreateSubWindow(win,(int)width*2,0,(int)width,(int)height);
    glutReshapeFunc(reshape_swin3);
    glutDisplayFunc(display_swin3);
    glutKeyboardFunc(keyboard);
    swin4 = glutCreateSubWindow(win,(int)width*3,0,(int)width,(int)height);
    glutReshapeFunc(reshape_swin4);
    glutDisplayFunc(display_swin4);
    glutKeyboardFunc(keyboard);

  init();
  glutMainLoop();
  return(0);
}