/*+-------------------------------------------------------------------------+
  | MORFOLOGIA MATEMATICA -                                                 |
  | ---------------------                                                   |
  | Implementacao de um algoritmo otimo para Linha de Divisor de Agua (LDA),|
  | usando a ideia de crescimento de regioes a partir de marcadores conforme|
  | sugerido no artigo de MEYER:                                            |
  |                                                                         |
  | F.MEYER. Color image segmentation. Technical report, E'cole des Mines   |
  | de Paris, Centre de Morpholigie Mathe'matique, France, 1992.            |
  |                                                                         |
  |    A vantagem deste algoritmo e' que somente os pixeis que serao modifi-|
  | cados sao avaliados em cada passo do algoritmo.                         |
  |    A desvantagem e' que e' necessario um passo adicional para identifi- |
  | car os minimos regionais que sao utilizados como marcadores para obten- |
  | cao da LDA completa. Para este passo estamos utilizando uma modificacao |
  | do algoritmo de LDA sugerido num artigo de Luc Vincent.                 |
  |                                                                         |
  | Versao 3.0                                                              |
  |                                                                         |
  | Por Luiz Eduardo da Silva.                                   ABRIL/1997 |
  +-------------------------------------------------------------------------+*/

#include <stdio.h>
 

#include "imagem.h"
#include "utils.h"

#define TAM_FILA 256
#define WSHED   -4
#define MARCA   -3
#define MASK    -2
#define INIT    -1
#define TAM_ORD 500
#define NRO_VIZ 4

/*+-------------------------------------------------------------------------+
  | Estrutura de Fila = usada para reproduzir a ordem de preenchimento das  |
  | bacias no relevo da imagem.                                             |
  +-------------------------------------------------------------------------+*/

struct {
  struct no *com;
  struct no *fim;
} Filas[TAM_FILA];

int Max_Prior;

/*+-------------------------------------------------------------------------+
  | Estruturas globais utilizadas no programa.                              |
  +-------------------------------------------------------------------------+*/

struct tord {
  int freq;
  struct no *list;
};


/*+-------------------------------------------------------------------------+
  | Rotina que ordena os pixel's da imagem.                                 |
  |    Depois da execucao desta rotina e atraves da estrutura  montada T e' |
  | possivel ter acesso direto aos pixel's de cada nivel (Threshold=Limiar) |
  | da imagem.                                                              |
  | Parametros:                                                             |
  |     1. I = Imagem de entrada.                                           |
  |     2. nc = Nro de colunas na imagem.                                   |
  |     3. nl = Nro de linhas na imagem.                                    |
  |     4. mn = Maximo nivel de cinza.                                      |
  |     5. T = Estrutura montada, um vetor para acesso direto aos pixel's de|
  |        cada nivel da imagem.                                            |
  |     6. Hmin = O minimo nivel na imagem.                                 |
  |     7. Hmax = O maximo nivel na imagem.                                 |
  +-------------------------------------------------------------------------+*/
void ordena_pixels (int *I, int nc, int nl, int mn,
                    struct tord *T, int *Hmin, int *Hmax) {
    int i, j;
    struct no *aux;
    *Hmin = mn;
    *Hmax = -3;
    for (i = 0; i < TAM_ORD; i++) {
	T[i].freq = 0;
	T[i].list = NULL;
    }
    for (i = nl-1; i >= 0; i--)
	for (j = nc-1; j >= 0; j--) {
	    if (*Hmin > *(I+i*nc+j))
	       *Hmin = *(I+i*nc+j);
	    if (*Hmax < *(I+i*nc+j))
	       *Hmax = *(I+i*nc+j);
	    T[*(I+i*nc+j)].freq++;
	    aux = malloc (sizeof (struct no));
	    aux->i = i;
	    aux->j = j;
	    aux->next = T[*(I+i*nc+j)].list;
	    T[*(I+i*nc+j)].list = aux;
	}
}

/*+-------------------------------------------------------------------------+
  | Destroi as listas de pixel's montadas pela rotina anterior              |
  +-------------------------------------------------------------------------+*/
void destroi_listas (struct tord *T, int mn) {
    int i;
    struct no *aux1, *aux2;
    for (i = 0; i < mn; i++)
	if (T[i].freq != 0) {
	    aux1 = T[i].list;
	    while (aux1) {
		aux2 = aux1;
		aux1 = aux1->next;
		free (aux2);
	    }
	}
}


/*+-------------------------------------------------------------------------+
  |   Implementacao de um algoritmo para encontra MINIMOS REGIONAIS baseada |
  | na simulacao de imersao sugerida num paper de Luc Vincent.              |
  | 1. Os pixel's sao ordenadas pelo seu nivel de cinza para permitir o     |
  |    acesso direto dos pixel's em cada nivel.                             |
  |                                                                         |
  |    Suponha que hmin e' o menor nivel de cinza na imagem.                |
  |                hmax e' o maior nivel de cinza.                          |
  |                T_h(I) sao os ptos da imagem I no nivel h.               |
  |                                                                         |
  | 2. Para todo h em [hmin, hmax-1].                                       |
  |  A) Todos os pontos T_h(I) que sao vizinhos de regioes marcadas tambem  |
  |     sao marcados. Estes NAO sao minimos.                                |
  |  B) Os ponto T_h(I) que nao foram marcados sao minimos regionais da img |
  |                                                                         |
  | USA...: L = img de label's                                              |
  | PRODUZ: M = img de MINIMOS, B = img de BACIAS, W = incializa img WSHED. |
  +-------------------------------------------------------------------------+*/
void minimos (int *I, int nl, int nc, int mn, 
	      int *M, int *B, int *W, int *NRO_MIN) {
    int h, i, j, i1, j1, 
	existe, vz_existe, n,
	Hmin, Hmax, label_corrente,
	com, fim, tam, nro_ptos;
    struct pto { int i, j;} *fila;
    struct tord T[TAM_ORD];
    struct no *aux;
    int *L;

    ordena_pixels (I, nc, nl, mn, T, &Hmin, &Hmax);
    aloca_memo (&L, nl, nc);
    for (i = 0; i < nl; i++)
	for (j = 0; j < nc; j++) {
	    *(L+i*nc+j) = INIT;  /*-- Imagem de label's   --*/
	    *(M+i*nc+j) = 0;     /*-- Imagem dos MINIMOS  --*/
	    *(B+i*nc+j) = 0;     /*-- Imagem de BACIAS    --*/
	    *(W+i*nc+j) = 0;
	}

    tam = nl * nc;
    fila = malloc (nl * nc * sizeof (struct pto));
    nro_ptos = 0;
    com = fim = 0;

    *NRO_MIN = 0;

    for (h = Hmin; h < Hmax + 1; h++) {
	if (T [h].freq != 0) {
	    /*+-------------------------------------------------------+
              | Atribuir o valor MARK para os pontos do nivel h, e    |
              | identificar os pontos que sao vizinhos de regioes mar-|
              | cadas.                                                |
              +-------------------------------------------------------+*/
	    aux = T [h].list;
	    while (aux) {
		i = aux->i;
		j = aux->j;
		*(L+i*nc+j) = MASK;
		existe = FALSE;
		for (n = 0; n < NRO_VIZ; n++) {
		    vz_existe = Viz (NRO_VIZ, n, i, j, nl, nc, &i1, &j1);
		    if (vz_existe && *(L+i1*nc+j1) == MARCA) {
			existe = TRUE;
		    }
		}
		if (existe) {
		    /*--- Insere um pixel vizinho de regiao marcada ---*/
		    fila[fim].i = i;
		    fila[fim].j = j;
		    fim = (fim + 1) % tam;
		    nro_ptos++;
		}
		aux = aux->next;
	    }

	    /*+-------------------------------------------------------+
              | Marcar com o valor MARCA os ptos do nivel h que sao   |
              | vizinhos de regioes marcadas, estes nao sao MINIMOS.  |
              +-------------------------------------------------------+*/
	    while (nro_ptos != 0 && nro_ptos <= tam) {
		
		/*--- Retira um pixel da fila para marcacao ---*/
		i = fila[com].i;
		j = fila[com].j;
		nro_ptos--;
		com = (com + 1) % tam;

		*(L+i*nc+j) = MARCA;
		for (n = 0; n < NRO_VIZ; n++) {
		    vz_existe = Viz(NRO_VIZ, n, i, j, nl, nc, &i1, &j1);
		    if (vz_existe && (*(L+i1*nc+j1) == MASK)) {
			
			/*--- Insere um pixel vizinho que nao foi marcado ---*/
			fila[fim].i = i1;
			fila[fim].j = j1;
			fim = (fim + 1) % tam;
			nro_ptos++;
			
			*(L+i1*nc+j1) = MARCA;
		    }
		}
	    }

	    
	    /*+--------------------------------------------------------+
	      | Os ptos que nao foram marcados no passo anterior sao   |
	      | MINIMOS REGIONAIS DA IMAGEM.                           |
	      +--------------------------------------------------------+*/
	    aux = T[h].list;
	    while (aux) {
		i = aux->i;
		j = aux->j;
		if (*(L+i*nc+j) == MASK) {

		    (*NRO_MIN)++;

		    fila[fim].i = i;
		    fila[fim].j = j;
		    fim = (fim + 1) % tam;
		    nro_ptos++;

		    while (nro_ptos != 0 && nro_ptos <= tam) {
			
			i = fila[com].i;
			j = fila[com].j;
			nro_ptos--;
			com = (com + 1) % tam;
			
			*(L+i*nc+j) = MARCA;
			*(M+i*nc+j) = *NRO_MIN;
			*(B+i*nc+j) = h;

			for (n = 0; n < NRO_VIZ; n++) {
			    vz_existe=Viz(NRO_VIZ, n, i, j, nl, nc, &i1, &j1);
			    if (vz_existe && (*(L+i1*nc+j1) == MASK)) {
				
				fila[fim].i = i1;
				fila[fim].j = j1;
				fim = (fim + 1) % tam;
				nro_ptos++;
				
				*(L+i1*nc+j1) = MARCA;
				*(M+i1*nc+j1) = *NRO_MIN;
				*(B+i1*nc+j1) = h;
			    }
			}
		    }
		}
		aux = aux->next;
	    }
	} /*-- end of ... if (T.freq... --*/
    } /*-- end of ... for (h = Hmin;... --*/

    destroi_listas (T, mn);
    desaloca_memo (&L);
}

/*+-------------------------------------------------------------------------+ 
  | Primitivas de manipulacao da fila de prioridade                         |
  +-------------------------------------------------------------------------+*/
void inicia_filas (void) {
  int i;
  for (i = 0; i < TAM_FILA; i++) {
      Filas[i].com = NULL;
      Filas[i].fim = NULL;
  }
}

void mostra_filas (void) {
  int i;
  struct no *aux;
  printf ("\n\nFILAS:");
  printf ("\nMax Prioridade: %d", Max_Prior);
  for (i = 0; i < TAM_FILA; i++) {
      printf ("\n%d ==>", i);
      aux = Filas[i].com;
      while (aux) {
         printf ("(%2d,%2d)", aux->i, aux->j);
         aux = aux->next;
      }
  }
  getchar ();
}

void insere (int nro_fila, int i, int j, int inicio) {
   struct no *aux;
   aux = (struct no *) malloc (sizeof (struct no));
   aux->i = i;
   aux->j = j;
   aux->next = NULL;
   if ((nro_fila < Max_Prior) && inicio)
       nro_fila = Max_Prior;
   if (Filas[nro_fila].fim)
       Filas[nro_fila].fim->next = aux;
   Filas[nro_fila].fim = aux;
   if (!Filas[nro_fila].com)
      Filas[nro_fila].com = aux;
}

void retira (int *i, int *j) {
   struct no *aux;
   *i = Filas[Max_Prior].com->i;
   *j = Filas[Max_Prior].com->j;
   aux = Filas[Max_Prior].com;
   Filas[Max_Prior].com = aux->next;
   free (aux);
   if (!Filas[Max_Prior].com) {
      Filas[Max_Prior].fim = NULL;
   }
}

int f_vazia (int nro_fila) {
   return Filas[nro_fila].com == NULL;
}

/*+-------------------------------------------------------------------------+
  | LDA com MARCADORES.                                                     |
  |    Tem dois passo: a Inicializacao e o Crescimento de regioes, que esta'|
  | explicado no seu corpo.                                                 |
  +-------------------------------------------------------------------------+*/
void lda_marcador (int *I, int nl, int nc, int mn,
		   int *M, int *B, int *W) {
   int i, j, y, x, n, MARK, BACIA, fim, vz_existe, eh_lda, achou;

/*+-------------------- Inicializacao -------------------+
  | A fila ordenada de prioridade e' criada com tantos   |
  | niveis de prioridade quantos niveis de cinza da ima- |
  | gem. Os marcadores sao identificados por label's. To-|
  | dos os ptos vizinhos dos marcadores sao colocados na |
  | na fila com a "prioridade apropriada". O valor de ca-|
  | da pto na imagem determina o nivel de prioridade da- |
  | quele ponto.                                         |
  +------------------------------------------------------+*/
   inicia_filas ();

   Max_Prior = TAM_FILA;
   for (i = 0; i < nl; i++)
       for (j = 0; j < nc; j++) {
	   if (*(M+i*nc+j) == 0) {
	       achou = FALSE;
	       for (n = 0; !achou && n < 4; n++) {
		   vz_existe = Viz (4, n, i, j, nl, nc, &y, &x);
		   if (vz_existe && *(M+y*nc+x) > 0)
		       achou = TRUE;
	       }
	       if (achou) {
		   /*--- Esta marcacao e' para evitar que ---*/
		   /*--- o pto seja colocado 2 vezes na fila */
		   *(M+i*nc+j) = MARCA; 
		   insere (*(I+i*nc+j), i, j, 0);
		   if (*(I+i*nc+j) < Max_Prior)
		       Max_Prior = *(I+i*nc+j);
	       }
	   }
       }

/*+------------- Crescimento dos marcadores -------------+
  | Para um ponto (i,j) extraido da fila de prioridade:  |
  | 1) Se este ponto tem em sua vizinhanca um e somente  |
  |    uma regiao rotulada, ele e' agregado a esta re-   |
  |    giao. Seus vizinhos sem label sao colocados na fi-|
  |    la ordenada, novamente "com a prioridade apropri- |
  |    ada" como no passo de inicializacao.              |
  | 2) Se este ponto e' vizinho de duas regioes com ro-  |
  |    tulos diferentes, este ponto e' ponto de borda e  |
  |    obtem um rotulo especial para caracterizar a fron-|
  |    teira.                                            | 
  +------------------------------------------------------+*/
   fim = FALSE;
   while (!fim) {
       retira(&i, &j);
       MARK = -1;
       eh_lda = FALSE;
       for (n = 0; n < 4; n++) {
	   vz_existe = Viz (4, n, i, j, nl, nc, &y, &x);
	   /*-- Se o vizinho existe e tem rotulo --*/
	   if (vz_existe && *(M+y*nc+x) > 0) {
	       if (MARK == -1) {
		   MARK = *(M+y*nc+x);
		   BACIA = *(B+y*nc+x);
	       }
	       else if (MARK != *(M+y*nc+x))
		   eh_lda = TRUE;
	   }
       }
       /*--- O pto (i,j) tem algum vizinho rotulado? ---*/
       if (eh_lda) {
	   *(M+i*nc+j) = WSHED;
	   *(B+i*nc+j) = *(I+i*nc+j);
	   *(W+i*nc+j) = 1;
       }
       else {
	   *(M+i*nc+j) = MARK;
	   *(B+i*nc+j) = BACIA;
	   *(W+i*nc+j) = 0;
	   for (n = 0; n < 4; n++) {
	       vz_existe = Viz (4, n, i, j, nl, nc, &y, &x);
	       /*--- Se o vizinho existe e nao tem marca ---*/
	       if (vz_existe && *(M+y*nc+x) == 0) {
		   *(M+y*nc+x) = MARCA;
		   insere (*(I+y*nc+x), y, x, 1);
	       }
	   }
       }
	   
       while (Max_Prior < TAM_FILA  && f_vazia (Max_Prior))
	   Max_Prior++;
       if (Max_Prior == TAM_FILA)
	   fim = TRUE;
   }
   for (i = 0; i < nl; i++)
       for (j = 0; j < nc; j++) {
	   if (*(M+i*nc+j) == 0) {
              *(M+i*nc+j) = WSHED;
              *(B+i*nc+j) = *(I+i*nc+j);
	      *(W+i*nc+j) = 1;
	   }
       }
}

void msg (char *s) {
    printf ("\nLDA de imagem numerica com MARCADORES");
    printf ("\n-------------------------------------");
    printf ("\nUSO.: %s <IMG.PGM>", s);
    printf ("\nONDE:\n");
    printf ("    <IMG.PGM>  nome do arquivo .pgm para processar\n\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, nc, nl, mn, nc1, nl1, nro_minimos, i, j;
  int *I,  /*-- Imagem original --*/
      *M,  /*-- Marcadores = Minimos regionais da imagem para LDA completa --*/
      *W,  /*-- Imagem Watershed resultante --*/
      *B;  /*-- Imagem de Bacias topograficas=minimo propagado delimitado  --*/
           /*--                               por caminho de ptos maximais --*/
  OK = FALSE;
  if (argc == 2){
     OK = le_imagem_pgm (argv[1], &I, &nl, &nc, &mn);
  }
  else msg (argv[0]);
  if (OK) {
      aloca_memo (&W, nl, nc);
      aloca_memo (&M, nl, nc);
      aloca_memo (&B, nl, nc);
      printf ("\n LDA completa com MARCADORES");
      info_imagem (argv[1], nl, nc, mn);
      minimos (I, nl, nc, mn, M, B, W, &nro_minimos);
      lda_marcador (I, nl, nc, mn, M, B, W);
      grava_imagem_pbm (W, "lda.pbm", nl, nc);
      /* system ("xv lda.pbm &");  */
      desaloca_memo (&B);
      desaloca_memo (&W);
      desaloca_memo (&M);
      desaloca_memo (&I);
  }
}