/*+-------------------------------------------------------------------------+
  | MORFOLOGIA MATEMATICA -                                                 |
  | ---------------------                                                   |
  | Implementacao de um programa eficiente de RECONSTRUCAO, conforme sugeri-|
  | no artigo:                                                              |
  |                                                                         |
  | MORPHOLOGICAL GRAYSCALE RECONSTRUCTION IN IMAGE ANALYSIS: APPLICATIONS  |
  | AND EFFICIENT ALGORITHMS. Luc Vincent, IEEE on Image Processing, vol.2, |
  | nro. 2, April, 1993.        
  |  
     Le duas imagens, mascara e marcadora, e calcula a reconstrucao         |
  |                                                                         |
  | Por Luiz Eduardo da Silva.                                AGOSTO/1997   |
  +-------------------------------------------------------------------------+*/
#include <stdio.h>
 
#include "imagem.h"
#include "utils.h"

int max (int a, int b) {
  a = (a > b)? a : b;
  return a;
}

int min (int a, int b) {
  a = (a < b)? a : b;
  return a;
}

/*+------------------------------------------------------------------------+
  | RECONSTRUCAO.                                                          |
  | Este algoritmo tem dois passos: a Inicializacao e a Propagacao         |
  | A imagem F e' a imagem na qual sera' feita a reconstrucao e a imagem G |
  | e' usada como imagem MARCADORA, e' onde sera' devolvido o resultado    |
  | processamento.                                                         |
  +------------------------------------------------------------------------+*/
void reconstrucao (int *F, int *G, int *M, int nl, int nc, int mn) {
    int i, j, x, y, v1, v2, v3, v4, maior, fim, vz_existe, n;
    struct no *COM = NULL, *FIM = NULL;

/*+---------------------------- Inicializacao -----------------------------+
  |   F = Imagem Mascara.                                                  |
  |   G = Imagem Marcadora, ambas definidas no dominio D e tal que G <= F  |
  | A reconstrucao e' determinada diretamente na imagem G:                 |
  | 1) Percorrer em sentido RASTER:                                        |
  |    Seja p um pixel corrente, G(p)<-(max{G(p),q in N+(p) U {p}})^F(p)   |
  | 2) Percorrer em sentido ANTIRASTER:                                    |
  |    Seja p um pixel corrente, G(p)<-(max{G(p),q in N-(p) U {p}})^F(p)   |
  |    Se existe q in N-(p) tal que G(q) < G(p) and G(q) < F(q)            |
  |       Entao INSERE_FILA (p)                                            |
  | OBS: N+(p) = e' o conjunto de vizinhos de p que sao alcancados antes   |
  |      de p numa varredura RASTER.                                       |
  |      N-(p) = e' o conjunto de vizinhos de p que sao alcancados depois  |
  |      de p numa varredura RASTER.                                       |
  |      ^ = denota a operacao de minimo pontual, correspondente da inter- |
  |          secao binaria.                                                |
  +------------------------------------------------------------------------+*/

    for (i = 0; i < nl; i++) 
      for (j = 0; j < nc; j++) {
	v1 = *(G+i*nc+j);
	v2 = (j > 0)? *(G+i*nc+j-1):0;
	v3 = (i > 0)? *(G+(i-1)*nc+j):0;
	v4 = *(F+i*nc+j);
	maior = max (v1,v2);
	maior = max (maior,v3);
	*(M+i*nc+j) = min (maior,v4);
      }
    for (i = nl-1; i > 0; i--) 
      for (j = nc-1; j > 0; j--) {
	v1 = G[i*nc+j];
	v2 = (j < nc-1)? *(G+i*nc+j+1):0;
	v3 = (i < nl-1)? *(G+(i+1)*nc+j):0;
	v4 = *(F+i*nc+j);
	maior = max (v1,v2);
	maior = max (maior,v3);
	v1 = *(G+i*nc+j) = min (maior,v4);
	if ((v2 < v1 && v2 < v4) ||
	    (v3 < v1 && v3 < v4)) {
	  insere_fila (&COM, &FIM, i, j);
	}
      }

/*+------------------------------ Propagacao ----------------------------+
  | Enquanto a fila nao esta' vazia faca:                                |
  |   p <- RETIRA_FILA()                                                 |
  |   Para todo pixel q in N(p) faca:                                    |
  |      Se G(q) < G(p) and F(q) != G(q)                                 |
  |         Entao G(q) <- min {G(p), F(q)}                               |
  |         INSERE_FILA (q)                                              |
  +----------------------------------------------------------------------+*/

    while (!fila_vazia (COM)) {
      retira_fila (&COM, &FIM, &i, &j);
      v1 = *(G+i*nc+j);
      for (n = 0; n < 4; n++) {
	vz_existe = Viz (4, n, i, j, nl, nc, &y, &x);
	if (vz_existe) {
	  v2 = *(G+y*nc+x);
	  v3 = *(F+y*nc+x);
	  if (v2 < v1 && v3 != v2) {
	    *(M+y*nc+x) = min (v1, v3);
	    insere_fila (&COM, &FIM, y, x);
	  }
	}	  
      }
    }
}

void msg (char *s) {
    printf ("\nRECONSTRUCAO por Luc Vincent");
    printf ("\n----------------------------");
    printf ("\nUSO.: %s <IMG.PGM> <RAIO>", s);
    printf ("\nONDE:\n");
    printf ("    <IMG.PGM>  nome do arquivo .pgm para processar\n");
    printf ("    <RAIO>     raio da bola digital (elemento estruturante)\n");
}

/*+------------------------------------------------------------------------+
  |                  P R O G R A M A    P R I N C I P A L                  |
  +------------------------------------------------------------------------+*/
void main (int argc, char *argv[]) {
  int OK, OK1, nc, nl, mn, i, j, k, di, dj, min, raio;
  int *F, *G, *M;
  OK = OK1 = FALSE;
  if (argc >= 2){
     OK =  le_imagem_pgm (argv[1], &F, &nl, &nc, &mn);
     OK1 = le_imagem_pgm (argv[2], &G, &nl, &nc, &mn);
   //  if (argc == 3) raio = atoi (argv[2]);
   //  else
   //    raio = 1;
  }
  else msg (argv[0]);
  if (OK && OK1) {
      aloca_memo (&M, nl, nc);
      printf ("\nRECONSTRUCAO por Luc Vincent");
      info_imagem (argv[1], nl, nc, mn);
      info_imagem (argv[2], nl, nc, mn);
   
 




      /*+-------------------------------------------------------------+
        | Inicializacao da imagem marcadora G.                        |
        | Para encontrar maximos regionais.                           |
        +-------------------------------------------------------------+
      for (i = 0; i < nl; i++)
	for (j = 0; j < nc; j++)
	  *(G+i*nc+j) = *(F+i*nc+j) > 0? *(F+i*nc+j)-1:0;
      -----------------------------------------------------------------*/

      reconstrucao (F, G, M, nl, nc, mn);
      grava_imagem_pgm (M, "reconstroi.pgm", nl, nc, mn);
    /*  system ("xv reconstroi.pgm &"); */

      for (i = 0; i < nl; i++)
	for (j = 0; j < nc; j++)
	  *(G+i*nc+j) = *(F+i*nc+j) - *(M+i*nc+j) == 0? 1:0;
      grava_imagem_pbm (G, "diferenca.pbm", nl, nc);
   /*   system ("xv diferenca.pbm &"); */

      desaloca_memo (&F);
      desaloca_memo (&M);
      desaloca_memo (&G);
  }
}