/*
 *  quantiza.c
 *
 *  Este programa ilustra o uso da funcao 
 *  glHistogram para obter a distribuicao de frequencia
 *  das intensidades de uma imagem de luminancia 
 *  e diferentes tecnicas de quantizacao
 *  (subdivisao uniforme (u), algoritmo de populosidade (p),
 *   algoritmo de corte mediano (m), dithering bayer (b),
 *   dithering floyd-steinberg (f)) 
 *
 *  Importante: deve-se inicializar glew com glewInit 
 *    para utilizar as extensoes (GL_EXTENSIONS). Certifique-se
 *    ainda se suporta a extensao GL_ARB_imaging
 *
 *            Ting (09/09/08)
 */

//
#define _USE_MATH_DEFINES

#ifdef _WIN32
#include <GL/glew.h>
#endif
#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.

#ifndef _WIN32
extern void glGetHistogram (GLenum, GLboolean, GLenum, GLenum, GLvoid *);
extern void glHistogram (GLenum, GLsizei, GLenum, GLboolean);
#endif

/* ascii code for the escape key */
#define ESCAPE 27

/* dimensoes da imagem */
#define ALTURA 64
#define LARGURA 64
#define HISTOGRAM_SIZE 256
#define QTDE 4

GLubyte imagem[2*ALTURA][2*LARGURA];
GLubyte q_imagem[2*ALTURA][2*LARGURA];
GLuint values[HISTOGRAM_SIZE], accum[HISTOGRAM_SIZE];
int nq[QTDE+1];
GLint dither[2][2];  //QTDE = 2*2

typedef struct AMOSTRA_struct {   /* pontos */
  int valor;
  int freq;
} AMOSTRA;

GLuint opcao;
int win;  // id da janela principal
int swin1, swin2; // id das subjanelas do glut

AMOSTRA samples[HISTOGRAM_SIZE];

void  mediana(double, int, int *, int *);
//*********************************************************************/
// Imagem com faixa vertical em tons de cinza
void makeHorSenoidalImage(void)
{
   int i, j, c;
    
   for (j = 0; j < (2*LARGURA-1); j=j+2) {
      c = 255 * sin((j*M_PI)/40);
      if (c < 0) c = 0;
      else if (c > 255) c = 255;
      for (i = 0; i < (2*ALTURA-1); i=i+2) {
         imagem[i][j] = imagem[i+1][j] = imagem[i+1][j+1] = imagem[i+1][j] = (GLubyte) c;
      }
   }
}

static int cmpfreq(const void *p1, const void *p2) 
{
  return (((AMOSTRA *)p1 )-> freq - ((AMOSTRA *)p2 )-> freq); 
}

static int cmpvalor(const void *p1, const void *p2) 
{
  return (((AMOSTRA *)p1 )-> valor - ((AMOSTRA *)p2 )-> valor); 
}

/****************************************************************
 *   Quantizacao
 ****************************************************************/
// Quantizacao Uniforme
void quantizacao_uniforme ()
{
  int i, j, k, delta;
 
  // obter os extremos das celulas de quantizacao
  delta = HISTOGRAM_SIZE/QTDE;
  for (i=0; i< QTDE+1; i++) {
    nq[i] = i*delta;
  }

  // transformar valores
  for (i = 0; i < 2*ALTURA; i++) {
      for (j = 0; j < 2*LARGURA; j++) {
	q_imagem[i][j] = 0;
	for (k=1; k< QTDE+1; k++) {
	  if (imagem[i][j] <= nq[k] && imagem[i][j] > nq[k-1])
	    q_imagem[i][j] = nq[k]-1;
       }
      }
  }
}

// Quantizacao por Populosidade
void quantizacao_populosidade ()
{
  int i,j,k;

  // ordenar as amostras em relacao a frequencia de ocorrencia
  qsort((void *)samples, HISTOGRAM_SIZE, sizeof(AMOSTRA), cmpfreq);

  // tomar os valores de maior ocorrencia como niveis de quantizacao
  for (i=0; i< QTDE; i++) {
       nq[i] = samples[HISTOGRAM_SIZE-1-i].valor;
  }

  // obter limites dos intervalos/celulas de quantizacao
  for (i=QTDE-1; i> 0 ; i--) {
    nq[i] = (nq[i]+nq[i-1])/2;
  }

  nq[QTDE] = 255;
  nq[0] = 0;

  // transformar valores   
  for (i = 0; i < 2*ALTURA; i++) {
      for (j = 0; j < 2*LARGURA; j++) {
	q_imagem[i][j] = 0;
	for (k=1; k< QTDE+1; k++) {
	  if (imagem[i][j] <= nq[k] && imagem[i][j] > nq[k-1])
	    q_imagem[i][j] = nq[k]-1;
       }
      }
  }
}

// Quantizacao por Corte Mediano
void  mediana(double tamanho, int max, int *qtde, int *limiares)
{
  int i, max_ind, maior, accum, valor_mediano;

  if (*qtde >= max) return;

  // determinar o maior sub-intervalo
  maior = 0;
  for (i=0; i<(*qtde-1); i++) {
      if ((samples[limiares[i+1]].valor - samples[limiares[i]].valor) > maior) {
	maior = samples[limiares[i+1]].valor - samples[limiares[i]].valor;
        max_ind = i;
      }
  }

  // determinar mediana para amostras de size <tamanho>
  accum = 0;

  for (i = limiares[max_ind]; accum < tamanho; i++) {
    accum += samples[i].freq;
    valor_mediano = i+1;
  }

  // insere a nova mediana
  for (i=*qtde-1; i>0 && samples[limiares[i]].valor > samples[valor_mediano].valor; i--)
    limiares[i+1] = limiares[i];

  limiares[i+1] = valor_mediano;
  (*qtde)++;

  // computar recorrentemente mediana do maior sub-intervalo 
  // (maior variacao de valores)
  mediana(tamanho/2, max, qtde, limiares);
  
  return;
}

void quantizacao_mediano ()
{
  int i,j, k, aux, ind[QTDE+1];
  float cor;

  // ordenar as amostras em relacao aos valores de intensidade
  qsort((void *)samples, HISTOGRAM_SIZE, sizeof(AMOSTRA), cmpvalor);

  // Determinar medianas/limiares das celulas de quantizacao
  ind[0] = 0;
  ind[1] = 255;
  aux = 2;  
  mediana((double)(2*LARGURA*ALTURA), 5, &aux, ind);

  // Formar celulas de quantizacao
  for (i=0; i<QTDE+1; i++) nq[i] = ind[i];
  
  // Computar os niveis de quantizacao (centroide)
  for (i=0; i< QTDE; i++) {
    for (aux=0,j=nq[i]; j<nq[i+1]; j++) aux += samples[j].freq;
    for (cor=0,j=nq[i]; j<nq[i+1]; j++)
      cor += ((double)(samples[j].freq)/aux)*samples[j].valor;
    ind[i] = cor;
  }

  // Substituir os valores 
  for (i = 0; i < 2*ALTURA; i++) {
    for (j = 0; j < 2*LARGURA; j++){
            for (k=0; k< QTDE+1; k++) {
	      if (imagem[i][j] == nq[0]) 
		q_imagem[i][j] = ind[0];
	      else if (imagem[i][j] <= nq[k] && imagem[i][j] > nq[k-1])
	          q_imagem[i][j] = ind[k-1];
	    }
    }
  }
}

/*********************************************************************
 * Dithering
 *********************************************************************/
// Dithering Ordenado: D_2 
void dithering_bayer(void)
{
  int i, j, ii, jj, aux, ind[QTDE+1];
  int dither[2][2];

  // ordenar as amostras em relacao aos valores de intensidade
  qsort((void *)samples, HISTOGRAM_SIZE, sizeof(AMOSTRA), cmpvalor);

  // Determinar medianas/limiares das celulas de quantizacao
  ind[0] = 0;
  ind[1] = 255;
  aux = 2;  
  mediana((double)(2*LARGURA*ALTURA), 5, &aux, ind);

  dither[0][0] = ind[1];
  dither[1][1] = ind[2];
  dither[0][1] = ind[3];
  dither[1][0] = ind[4];

  // Substituir os valores 
  for (i = 0; i < 2*ALTURA; i++) {
    for (j = 0; j < 2*LARGURA; j++){
      ii = i%2; jj = j%2;
      if (imagem[i][j] <= dither[ii][jj]) 
	q_imagem[i][j] = 0;
      else 
	q_imagem[i][j] = 255;
    }
  }
}

// Dithering Floyd-Steinberg 
void dithering_floyd_steinberg(void)
{
  int i, j, limiar, erro;
  int aux_imagem[2*ALTURA][2*LARGURA];

  // Utilizar o ponto médio como limiar
  limiar = 128;

  // Copiar a imagem
  for (i = 0; i < 2*ALTURA; i++) {
    for (j = 0; j < 2*LARGURA; j++) 
      aux_imagem[i][j] = imagem[i][j];
  }

  // Comparar com o limiar (threshold)
  for (i = 0; i < 2*ALTURA; i++) { 
    for (j = 0; j < 2*LARGURA; j++) {
      // assegurar que os valores estejam na faixa valida
      if (aux_imagem[i][j] < 0) aux_imagem[i][j] = 0;
      else if (aux_imagem[i][j] > 255) aux_imagem[i][j] = 255;

      if (aux_imagem[i][j] <= limiar) {
	    erro = aux_imagem[i][j];
	    q_imagem[i][j] = 0;
      } else {
        erro = aux_imagem[i][j] - 255;
        q_imagem[i][j] = 255;
      }
      // e difundir os erros entre os pixels vizinhos
      if (j < (2*LARGURA)-1)
		aux_imagem[i][j+1] += 7./16*erro;
	  if (i < (2*ALTURA)-1) {
		aux_imagem[i+1][j-1] += 3./16*erro;  
		aux_imagem[i+1][j] += 5./16*erro;  
		if (j< (2*LARGURA)-1) aux_imagem[i+1][j+1] += 1./16*erro;  
	  }
    }
  }
}

/********************************************************************/
// OpenGL functions
//
void init(void)
{
  glClearColor (0.0, 0.0, 0.0, 0.0);
  glShadeModel(GL_FLAT);
  makeHorSenoidalImage();

  glPixelStorei(GL_UNPACK_ALIGNMENT, 1);
  opcao = 0;
}

void display_swin1()
{
  int i;

  // Desenhar os pixels da imagem
   glutSetWindow(swin1);
   glDisable(GL_DITHER);

   /* Define o formato da imagem da qual será feita a distribuicao
     de frequencias  */
   glHistogram(GL_HISTOGRAM, HISTOGRAM_SIZE, GL_LUMINANCE, GL_FALSE);

   /* Habilita histograma */
   glEnable(GL_HISTOGRAM);

   glClear(GL_COLOR_BUFFER_BIT);
   glRasterPos2i(0, 0);
   glDrawPixels(LARGURA*2, ALTURA*2, GL_LUMINANCE, GL_UNSIGNED_BYTE, imagem);
   glGetHistogram(GL_HISTOGRAM, GL_TRUE, GL_LUMINANCE, 
		  GL_UNSIGNED_INT, values);

  glDisable(GL_HISTOGRAM);

  for (i=0; i< HISTOGRAM_SIZE; i++) { 
    samples[i].valor = i;
    samples[i].freq = values[i];
  }

  glFlush();

}

void display_swin2()
{
  glutSetWindow(swin2);
  glDisable(GL_DITHER);

  glClear(GL_COLOR_BUFFER_BIT);
  glRasterPos2i(0, 0);

   switch (opcao) {
   case 0:
     quantizacao_uniforme();
     break;
   case 1:
     quantizacao_populosidade();
     break;
   case 2:
     quantizacao_mediano();
     break;
   case 3:
     dithering_bayer();
     break;
   case 4:
     dithering_floyd_steinberg();
     break;
   }

  glDrawPixels(LARGURA*2, ALTURA*2, GL_LUMINANCE, GL_UNSIGNED_BYTE, q_imagem);

  glFlush();
}

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

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

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

void reshape_swin1(int w, int h)
{
  glMatrixMode(GL_PROJECTION);
  glLoadIdentity();
  gluOrtho2D(0.0, (GLfloat) LARGURA*2, 0.0, (GLfloat) ALTURA*2);
  glMatrixMode(GL_MODELVIEW);
  glLoadIdentity();
}

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

void keyboard(unsigned char key, int x, int y)
{
   switch (key) {
   case 'u':
     opcao = 0;
     break;
   case 'p':
     opcao = 1;
     break;
   case 'm':
     opcao = 2;
     break;
   case 'b':
     opcao = 3;
     break;
   case 'f':
     opcao = 4;
     break;
   case 27:
     exit(0);
   }
   glutPostRedisplay();
}

/************************************************************************/
// Main configuration
//
int main(int argc, char** argv)
{
#ifdef _WIN32 
	GLenum err;
#endif

   glutInit(&argc, argv);
   glutInitDisplayMode(GLUT_SINGLE | GLUT_RGB );
   glutInitWindowSize(LARGURA*4, ALTURA*2);
   glutInitWindowPosition(0,0);
   win = glutCreateWindow(argv[0]);

#ifdef _WIN32 
   // Inicializar glew
  err = glewInit();
  if (GLEW_OK != err)
  {
    /* Problem: glewInit failed, something is seriously wrong. */
    printf ("glewInit failed.\n");
  }
  
  if (glewGetExtension("GL_ARB_imaging"))
  {
    printf ("GL_ARB_Imaging is supported.\n");
  }

#endif

    glutDisplayFunc(main_display);
    glutReshapeFunc(main_reshape);
    glutKeyboardFunc(keyboard);
   swin1 = glutCreateSubWindow(win,0,0,LARGURA*2,ALTURA*2);
    glutReshapeFunc(reshape_swin1);
    glutDisplayFunc(display_swin1);
    glutKeyboardFunc(keyboard);
   swin2 = glutCreateSubWindow(win,LARGURA*2,0,LARGURA*2,ALTURA*2);
    glutReshapeFunc(reshape_swin2);
    glutDisplayFunc(display_swin2);
    glutKeyboardFunc(keyboard);

   init();

   glutMainLoop();
   return 0; 
}