/* gcc -Wall -lm snake_couleur.c -o snake_couleur
*/
#include<stdio.h>
#include<stdlib.h>
#include<string.h>
#include<math.h>

#define Alpha 2.      /* Poids accorde a la continuite du snake */
#define Beta .8      /* Poids accorde a la "courbure" du snake */
#define W_grad 4.    /* Poids accorde au gradient de l'image   */
#define nb_points 22  /* Nombre de points du snake */
#define nb_iter  50   /* Nombre d'iterations pour le deplacemt du snake */
#define INPUT_DIM "fruits.dim"
#define INPUT_IMA "fruits.rvb"
#define INPUT_FLOU_DIM "fruits_flou.dim"
#define INPUT_FLOU_IMA "fruits_flou.rvb"
#define OUTPUT_DIM "snake.dim"
#define OUTPUT_IMA "snake.ima"
#define OUTPUT_RVB "snake.rvb"

int snake[nb_points][2];
int Decision[81][nb_points];

int ecrit_image ( unsigned char *image, int ncol, int nlig );
int ecrit_image_couleur ( unsigned char *image_r, unsigned char *image_v,
			  unsigned char *image_b, int ncol, int nlig );
int segment ( int nlig, int ncol, unsigned char *image, 
	      int i0, int j0, int i1, int j1 );
 /* i -> indice de ligne et j -> indice de colonne */
int nouvelle_ligne ( int t, int code );
int nouvelle_colonne ( int t, int code );
double derivee ( int t, int code_precedent, int code );
double derivee2 ( int t, int code_precedent, int code, int code_suivant );
double Eext ( int t, int code, unsigned char *image_r, unsigned char *image_v, 
	      unsigned char *image_b, int ncol );
double Etot ( int t, int code_precedent, int code, int code_suivant,
	      double Energie_prec, unsigned char* image_r,  unsigned char *image_v,
	      unsigned char *image_b, int ncol );
int distance ( int t, int code_precedent, int code );
int get_V0 ( int avant_dernier, int dernier );

/* valeur de l'etat : etat = 9 * code_suivant + code_courant */


void main(int argc, char *argv[])
{

  int debug1,debug2;
  double debug;

  double Energie[81][nb_points];
  double Energie_courante, Energie_min, Energie_min_finale, Energie_prec;
  /*  int Decision[81][nb_points];*/
  /*  int snake[nb_points][2];*/
  int i,j, iter, t, etat, code_prec, deplacement;
  int decision_finale, decision, V0;
  unsigned char *image_r, *image_v, *image_b;
  unsigned char *imageTravail_r, *imageTravail_v, *imageTravail_b;
  unsigned char *buffer;
  int nlig, ncol;
  FILE *fp;
  
  /* Initialisation du snake */
  /* Rq : snake [i][0] = indice de ligne, snake [i][1] = indice de colonne */

  /*  snake[28][0] = 217; snake[28][1] = 74;
  snake[29][0] = 206; snake[29][1] = 95;
  snake[0][0] = 220; snake[0][1] = 112;
  snake[1][0] = 206; snake[1][1] = 133; 
  snake[2][0] = 222; snake[2][1] = 150;
  snake[3][0] = 208; snake[3][1] = 173; 
  snake[4][0] = 221; snake[4][1] = 194;
  snake[5][0] = 199; snake[5][1] = 204;
  snake[6][0] = 180; snake[6][1] = 222;
  snake[7][0] = 166; snake[7][1] = 206;
  snake[8][0] = 137; snake[8][1] = 218;
  snake[9][0] = 118; snake[9][1] = 204;
  snake[10][0] = 96; snake[10][1] = 224;
  snake[11][0] = 84; snake[11][1] = 206;
  snake[12][0] = 60; snake[12][1] = 216;
  snake[13][0] = 50; snake[13][1] = 202;
  snake[14][0] = 48; snake[14][1] = 184;
  snake[15][0] = 34; snake[15][1] = 168;
  snake[16][0] = 47; snake[16][1] = 146;
  snake[17][0] = 32; snake[17][1] = 126;
  snake[18][0] = 46; snake[18][1] = 110;
  snake[19][0] = 32; snake[19][1] = 86 ;
  snake[20][0] = 46; snake[20][1] = 77;
  snake[21][0] = 42; snake[21][1] = 56;
  snake[22][0] = 57; snake[22][1] = 50;
  snake[23][0] = 73; snake[23][1] = 40; 
  snake[24][0] = 94; snake[24][1] = 46;
  snake[25][0] = 108; snake[25][1] = 36;
  snake[26][0] = 126; snake[26][1] = 42;
  snake[27][0] = 206; snake[27][1] = 54;*/


  snake[0][0] = 108; snake[0][1] = 16;
  snake[1][0] = 128; snake[1][1] = 27;
  snake[2][0] = 132; snake[2][1] =  47;
  snake[3][0] = 143; snake[3][1] = 60;
  snake[4][0] = 140; snake[4][1] = 84;
  snake[5][0] = 138; snake[5][1] = 108;
  snake[6][0] = 136; snake[6][1] = 131;
  snake[7][0] = 118; snake[7][1] = 156;
  snake[8][0] = 100; snake[8][1] = 169;
  snake[9][0] = 77 ; snake[9][1] = 175;
  snake[10][0] = 55; snake[10][1] = 162;
  snake[11][0] = 39; snake[11][1] = 154; 
  snake[12][0] = 29; snake[12][1] = 132;
  snake[13][0] = 14; snake[13][1] = 115;
  snake[14][0] = 7 ; snake[14][1] = 89;
  snake[15][0] = 12; snake[15][1] = 72;
  snake[16][0] = 8 ; snake[16][1] = 48;
  snake[17][0] = 30; snake[17][1] = 30;
  snake[18][0] = 41; snake[18][1] = 11;
  snake[19][0] = 55; snake[19][1] = 4;
  snake[20][0] = 78; snake[20][1] = 8;
  snake[21][0] = 93; snake[21][1] = 5;

  /* snake[0][0] = 273 ; snake[0][1] = 111;
  snake[1][0] = 243; snake[1][1] = 141;
  snake[2][0] = 224; snake[2][1] = 157;
  snake[3][0] = 216; snake[3][1] = 193;
  snake[4][0] = 216; snake[4][1] = 230;
  snake[5][0] = 232; snake[5][1] = 267;
  snake[6][0] = 258; snake[6][1] = 284;
  snake[7][0] = 281; snake[7][1] = 295;
  snake[8][0] = 310; snake[8][1] = 296;
  snake[9][0] = 325; snake[9][1] = 281;
  snake[10][0] = 355; snake[10][1] = 236;
  snake[11][0] = 355; snake[11][1] = 188;
  snake[12][0] = 353; snake[12][1] = 168;
  snake[13][0] = 349; snake[13][1] = 150;
  snake[14][0] = 319; snake[14][1] = 131;
  snake[15][0] = 294; snake[15][1] = 114;*/

  /* points pour le singe  */
  /*  snake[0][0]  = 26 ; snake[0][1]  = 28 ;
  snake[1][0]  = 43 ; snake[1][1]  = 13 ;
  snake[2][0]  = 69 ; snake[2][1]  = 23 ;
  snake[3][0]  = 77 ; snake[3][1]  = 40 ;
  snake[4][0]  = 89 ; snake[4][1]  = 48 ;
  snake[5][0]  = 91 ; snake[5][1]  = 82 ;
  snake[6][0]  = 95 ; snake[6][1]  = 136;
  snake[7][0]  = 96; snake[7][1]  = 165;
  snake[8][0]  = 106; snake[8][1]  = 192; 
  snake[9][0]  = 111; snake[9][1]  = 210;
  snake[10][0]  = 99 ; snake[10][1]  = 219; 
  snake[11][0] = 81 ; snake[11][1] = 205;
  snake[12][0] = 72 ; snake[12][1] = 177;
  snake[13][0] = 67 ; snake[13][1] = 156;
  snake[14][0] = 63 ; snake[14][1] = 122;
  snake[15][0] = 49 ; snake[15][1] = 98 ;
  snake[16][0] = 44 ; snake[16][1] = 62 ;*/

  /* Lecture de l'image */
  /* Ouverture du .dim */
  fp = fopen( INPUT_DIM,"r");
  if(fp == 0)
    {
      printf("Le fichier %s n'existe pas.\nArret du programme.\n", INPUT_DIM);
      exit(0);
    }

  /* lecture des dimensions de l'image */
  fscanf( fp, "%d %d", &ncol, &nlig);
  fclose(fp);
  printf(" Format de l'image :\n");
  printf("            largeur : %d, hauteur : %d\n", ncol, nlig);
  
  /* Allocation memoire du buffer de lecture */
  buffer = (unsigned char *) malloc(3*ncol*sizeof(char));
  printf(" Lecture...\n");
  
  /* allocation de l'image*/ 
  image_r = (unsigned char *) malloc(nlig*ncol*sizeof(char));
  imageTravail_r = (unsigned char *) malloc(nlig*ncol*sizeof(char));
  if (image_r == NULL)
    {
      printf("Pas assez de memoire pour le chargement de l'image source.\n");
      exit(0);
    }
  if (imageTravail_r == NULL)
    {
      printf("Pas assez de memoire pour le chargement de l'image source.\n");
      exit(0);
    }
   
  image_v = (unsigned char *) malloc(nlig*ncol*sizeof(char));
  imageTravail_v = (unsigned char *) malloc(nlig*ncol*sizeof(char));
  if (image_v == NULL)
    {
      printf("Pas assez de memoire pour le chargement de l'image source.\n");
      exit(0);
    }
  if (imageTravail_v == NULL)
    {
      printf("Pas assez de memoire pour le chargement de l'image source.\n");
      exit(0);
    }
   
  image_b = (unsigned char *) malloc(nlig*ncol*sizeof(char));
  imageTravail_b = (unsigned char *) malloc(nlig*ncol*sizeof(char));
  if (image_b == NULL)
    {
      printf("Pas assez de memoire pour le chargement de l'image source.\n");
      exit(0);
    }
  if (imageTravail_b == NULL)
    {
      printf("Pas assez de memoire pour le chargement de l'image source.\n");
      exit(0);
    }
   
  /* lecture de l'image source */
  fp = fopen(INPUT_FLOU_IMA,"r");
  for(i=0; i<nlig; i++) 
    {
      fread(buffer, sizeof(char),3*ncol, fp); 
      for(j=0;j<ncol;j++)
	{
	  image_r[i*ncol + j] = buffer[3*j+0];
	  image_v[i*ncol + j] = buffer[3*j+1];
	  image_b[i*ncol + j] = buffer[3*j+2];
	}
    }

  fclose(fp); 

  fp = fopen(INPUT_IMA,"r");
  for(i=0; i<nlig; i++) 
    {
      fread(buffer, sizeof(char),3*ncol, fp); 
      for(j=0;j<ncol;j++)
	{
	  imageTravail_r[i*ncol + j] = buffer[3*j+0];
	  imageTravail_v[i*ncol + j] = buffer[3*j+1];
	  imageTravail_b[i*ncol + j] = buffer[3*j+2];
	}
    }
  fclose(fp);
 
  ecrit_image_couleur( imageTravail_r, imageTravail_v, imageTravail_b, ncol, nlig);

  /* on trace le snake */
  /* Composante rouge */
  /*  for(i=0; i< nb_points-1; i++)
    {
      segment(nlig, ncol, imageTravail_r, snake[i][1], snake[i][0],
	      snake[i+1][1], snake[i+1][0]);
    }
  
  segment(nlig, ncol, imageTravail_r, snake[0][1], snake[0][0],
	  snake[nb_points-1][1], snake[nb_points-1][0]);*/
  
  /* Composante verte */
  /*  for(i=0; i< nb_points-1; i++)
    {
      segment(nlig, ncol, imageTravail_v, snake[i][1], snake[i][0],
	      snake[i+1][1], snake[i+1][0]);
    }
  
   segment(nlig, ncol, imageTravail_v, snake[0][1], snake[0][0],
	  snake[nb_points-1][1], snake[nb_points-1][0]);*/
  
  /* Composante bleue */
  /*  for(i=0; i< nb_points-1; i++)
    {
      segment(nlig, ncol, imageTravail_b, snake[i][1], snake[i][0],
	      snake[i+1][1], snake[i+1][0]);
    }
  
   segment(nlig, ncol, imageTravail_b, snake[0][1], snake[0][0],
	  snake[nb_points-1][1], snake[nb_points-1][0]);
  
  ecrit_image_couleur(imageTravail_r, imageTravail_v, imageTravail_b, ncol, nlig);*/


  /* On commence par une boucle histoire de voir si ca marche un peu... */
  /* A changer apres par un test sur les variations de l'energie du snake */
 
 /* A l'attaque !!! */
  /*ecrit_image(image_v, ncol, nlig);*/
  for (iter=0; iter < nb_iter; iter++)
    {
      
      /* Initialisation de la matrice d'energie et de decision */
      for (i=0 ; i<81; i++)
	for (j=0; j<nb_points; j++)
	  {
	    Energie[i][j] = 0.0;
	    Decision[i][j] = 0;
	  }
      
      /* Premier point du snake !!!! */
      /* On ne calcule donc que son energie !!! */
            for (etat = 0; etat < 81; etat ++)
	{
	  Energie[etat][0] = Etot(0, 4, etat%9, (int)(etat/9), 0, image_r, image_v, image_b, ncol );
	  debug=Energie[etat][0];
	}
	
      /* On commence notre voyage le long du snake; t = temps ! */
      for (t=1; t<nb_points-1; t++)
	{
	  for (etat=0; etat < 81; etat ++)
	    {
	      /* On calcule les energies possibles et on choisit */
	      /* la plus petite -> Decision                      */
	      Energie_prec = Energie[9*(etat%9)+0][t-1];
	      Energie_courante = Etot (t, 0, etat%9, (int)(etat/9), Energie_prec, image_r, image_v, image_b, ncol);

	      /* on initialise Energie_min et decision */
	      Energie_min = Energie_courante;
	      decision = 0;
	      
	      for (code_prec=1; code_prec<9; code_prec++)
		{
		  Energie_prec =  Energie[9*(etat%9)+code_prec][t-1];
		  Energie_courante = Etot(t, code_prec, etat%9, (int)(etat/9), Energie_prec, image_r, image_v, image_b, ncol);
				  
		  if ( Energie_courante <= Energie_min )
		    if ( distance (t, code_prec, etat%9) > 8)
		    {
		      Energie_min = Energie_courante;
		      decision = code_prec;
		    }
		}
	      Energie[etat][t] = Energie_min;
	      /*Decision[etat][t] = decision;*/
	      Decision[etat][t]= 9 * (etat%9) + decision;
	      debug1=Decision[etat][t];
	    }
	}


      /*****************************************************/
      /* pour t = nb_points-1 (le dernier point  du snake) */
      /*****************************************************/

      for (etat=0; etat < 9 ; etat ++)
	{
	  /* On calcule les energies possibles et on choisit */
	  /* la plus petite -> Decision                      */
	  /* (code du deplacement du point precedent (t-1)   */ 
	  

	  /* On va chercher le premier point pour calculer la derivee2 */
	  V0 = get_V0 (0, etat);
	  Energie_prec =  Energie[9*(etat%9)+0][nb_points-2];
	  Energie_courante = Etot(nb_points-1, 0, etat%9, V0  , Energie_prec, image_r, image_v, image_b, ncol);
	 	  
	  /* on initialise Energie_min et decision */ 
 	  
	  Energie_min = Energie_courante;
	  decision = 0;
	  for (code_prec=1; code_prec<9; code_prec++)
	    { 
	      V0 = get_V0 (code_prec, etat);
	      Energie_prec =  Energie[9*(etat%9)+code_prec][nb_points-2];
	      Energie_courante = Etot( nb_points-1, code_prec, etat%9, V0, Energie_prec, image_r, image_v, image_b, ncol);
	      
	      if ( Energie_courante <= Energie_min )
		if ( distance (nb_points-1, code_prec, etat%9) > 8)
		{
		  Energie_min = Energie_courante;
		  decision = code_prec;
		  
		}
	    }
	  Energie[etat][nb_points-1] = Energie_min;
	  Decision[etat][nb_points-1] = 9 * (etat%9) + decision;
	  debug1=Decision[etat][nb_points-1];
	}
      
      /* on cherche l'energie totale minimum */
      /*Energie_min_finale = Energie[0][nb_points-1]; 
      decision_finale = 0;
      for (etat=1; etat<9; etat++)
	{
	  if (Energie[etat][nb_points-1] <= Energie_min_finale)
	    {
	      Energie_min_finale = Energie[etat][nb_points-1];
	      decision_finale = etat;
	    }
	}*/

      /* ecriture de l'image (.ima) */
      /*      fp = fopen("decision.txt","w");
      for(etat=0;etat<81;etat++)
	    fprintf(fp,"%d  %d  %d  %d\n", Decision[etat][0], Decision[etat][1], Decision[etat][2], Decision[etat][3]);
      fclose(fp);
      */

      Energie_min_finale = Energie[0][nb_points-1];
      decision = 0;
      for (etat=1; etat<9; etat++)
	{
	  Energie_courante= Energie [etat][nb_points-1];	      
	  if ( Energie_courante <= Energie_min_finale )
	    {
	      Energie_min_finale = Energie_courante;
	      decision = etat;
	    }
	}
      /*decision_finale = decision%9;*/
      
      if((iter%5)==0)
	{
	  printf("%f  ",(float)Energie_min_finale);
	  printf("%d  \n", decision);
	}

      /* on remonte le tableau de decision */
      /* et on recalcule la position du snake */  
     
      deplacement = decision;
      etat = deplacement;
      for (t=nb_points-1; t>=0; t--)
	{
	  snake[t][0] = nouvelle_colonne(t, deplacement); 
	  snake[t][1] = nouvelle_ligne(t, deplacement); 
	  etat = Decision[etat][t];
	  deplacement = etat%9;
	}

      
      /* deplacement = decision_finale;
      for (t=nb_points-1; t>=0; t--)
	{
	  snake[t][0] = nouvelle_colonne(t, deplacement); 
	  snake[t][1] = nouvelle_ligne(t, deplacement);
	  deplacement = Decision[deplacement][t];
	  if (t==1) printf("%d\n", deplacement);
	}*/

      /* trace du snake apres une iteration */
      /*      if ( (iter%5) == 0)
	{
	  for(i=0; i< nb_points-1; i++)
	    {
	      segment(nlig, ncol, imageTravail_r, snake[i][1], snake[i][0],
		      snake[i+1][1], snake[i+1][0]);
	    }
	    
	  segment(nlig, ncol, imageTravail_r, snake[0][1], snake[0][0],
		  snake[nb_points-1][1], snake[nb_points-1][0]);

      	  for(i=0; i< nb_points-1; i++)
	    {
	      segment(nlig, ncol, imageTravail_v, snake[i][1], snake[i][0],
		      snake[i+1][1], snake[i+1][0]);
	    }
	    
	  segment(nlig, ncol, imageTravail_v, snake[0][1], snake[0][0],
		  snake[nb_points-1][1], snake[nb_points-1][0]);

      	  for(i=0; i< nb_points-1; i++)
	    {
	      segment(nlig, ncol, imageTravail_b, snake[i][1], snake[i][0],
		      snake[i+1][1], snake[i+1][0]);
	    }
	    
	  segment(nlig, ncol, imageTravail_b, snake[0][1], snake[0][0],
		  snake[nb_points-1][1], snake[nb_points-1][0]);

      	  ecrit_image_couleur(imageTravail_r, imageTravail_v, imageTravail_b, ncol, nlig);
	}*/

      if ( iter == ( nb_iter-1 ) )
	{
	  for(i=0; i< nb_points-1; i++)
	    {
	      segment(nlig, ncol, imageTravail_r, snake[i][1], snake[i][0],
		      snake[i+1][1], snake[i+1][0]);
	    }
	  
	  segment(nlig, ncol, imageTravail_r, snake[0][1], snake[0][0],
		  snake[nb_points-1][1], snake[nb_points-1][0]);

      	  for(i=0; i< nb_points-1; i++)
	    {
	      segment(nlig, ncol, imageTravail_v, snake[i][1], snake[i][0],
		      snake[i+1][1], snake[i+1][0]);
	    }
	  
	  segment(nlig, ncol, imageTravail_v, snake[0][1], snake[0][0],
		  snake[nb_points-1][1], snake[nb_points-1][0]);
	  
      	  for(i=0; i< nb_points-1; i++)
	    {
	      segment(nlig, ncol, imageTravail_b, snake[i][1], snake[i][0],
		      snake[i+1][1], snake[i+1][0]);
	    }
	  
	  segment(nlig, ncol, imageTravail_b, snake[0][1], snake[0][0],
		  snake[nb_points-1][1], snake[nb_points-1][0]);
	  
	  ecrit_image_couleur(imageTravail_r, imageTravail_v, imageTravail_b, ncol, nlig);
	}
	
    }
 

  free(buffer);
  free(image_r);
  free(imageTravail_r);
  free(image_v);
  free(imageTravail_v);
  free(image_b);
  free(imageTravail_b);
 
}


/* Fonctions de calcul du deplacement du point du snake */
/* Matrice de codage du deplacement : */
/*  0 1 2  */  
/*  3 4 5  */  
/*  6 7 8  */  

int nouvelle_colonne (int t, int code)
{
  int i;

  i = snake[t][0] + code%3 - 1;
  return(i);
}

int nouvelle_ligne (int t, int code)
{
  int j;
 
  j = snake[t][1] + (int)code/3 - 1;
  return(j);

}

/* trouve le V0 pour le calcul de l'energie du dernier point */
int get_V0 (int avant_dernier, int dernier)
{
  int t, etat;

  etat = 9 * dernier + avant_dernier;
  for (t=nb_points-2; t>0; t--)
    {
      etat = Decision[etat][t];
    }
  return (etat%9);
}

/* Fonction de calcul de la norme de la derivee du snake */

double derivee (int t, int code_precedent, int code)
{
  int i0, j0, i1, j1;
  double res;

  if (t==0)
    {
      i0 = nouvelle_colonne(nb_points-1, code_precedent);
      j0 = nouvelle_ligne(nb_points-1, code_precedent);
    }
  else
    {
      i0 = nouvelle_colonne(t-1, code_precedent);
      j0 = nouvelle_ligne(t-1, code_precedent); 
    }
  i1 = nouvelle_colonne(t, code);
  j1 = nouvelle_ligne(t, code); 
  res = (i0-i1)*(i0-i1) + (j0-j1)*(j0-j1);
 
  return(res);
}

int distance (int t, int code_precedent, int code)
{
    int i0, j0, i1, j1;
    int res;

    if (t==0)
      {
	i0 = nouvelle_colonne(nb_points-1, code_precedent);
	j0 = nouvelle_ligne(nb_points-1, code_precedent);
      }
    else
      {
	i0 = nouvelle_colonne(t-1, code_precedent);
	j0 = nouvelle_ligne(t-1, code_precedent); 
      }
    i1 = nouvelle_colonne(t, code);
    j1 = nouvelle_ligne(t, code); 
    res = (int) ( sqrt ( (double)( (i0-i1)*(i0-i1) + (j0-j1)*(j0-j1) ) ) );
    
    return(res);
}

/* Fonction de calcul de la norme de la derivee seconde du snake */

double derivee2 (int t, int code_precedent, int code, int code_suivant)
{
  int i0, j0, i1, j1, i2, j2;
  double res;

  if (t==0)
    {
      i0 = nouvelle_colonne(nb_points-1, code_precedent);
      j0 = nouvelle_ligne(nb_points-1, code_precedent);
    }
  else 
    { 
      i0 = nouvelle_colonne(t-1, code_precedent);
      j0 = nouvelle_ligne(t-1, code_precedent); 
    }
  i1 = nouvelle_colonne(t, code);
  j1 = nouvelle_ligne(t, code);
 
  if (t==(nb_points-1))
    {
      i2 = nouvelle_colonne(0, code_suivant);
      j2 = nouvelle_ligne(0, code_suivant); 
    }
  else
    {
      i2 = nouvelle_colonne(t+1, code_suivant);
      j2 = nouvelle_ligne(t+1, code_suivant);
    }
  
  res = (i2-2*i1+i0)*(i2-2*i1+i0) + (j2-2*j1+j0)* (j2-2*j1+j0);
  return(res);
}

/**********************************************/
/* Fonction de calcul de l'energie exterieure */
/**********************************************/

double Eext (int t, int code, unsigned char *image_r, unsigned char *image_v, unsigned char *image_b, int ncol)
{
  int i0, j0, i1, j1;
  float Y0, Y1, Y2;
  double gradient;
  double intensite;

  i0 = snake[t][0];
  j0 = snake[t][1];
  i1 = nouvelle_colonne(t, code);
  j1 = nouvelle_ligne(t, code);

  /* Pour l'instant : Energie exterieure = - W_grad * |gradient| */
  /* Ici : on prend le gradient de la luminance de TV */
  
  /* Y0 = 0.3*image_r[j1*ncol + i1] + 0.6*image_v[j1*ncol + i1] + 0.1*image_b[j1*ncol + i1];
  Y1 = 0.3*image_r[j1*ncol + (i1-1)] + 0.6*image_v[j1*ncol + (i1-1)] + 0.1*image_b[j1*ncol + (i1-1)];
  Y2 = 0.3*image_r[(j1-1)*ncol + i1] + 0.6*image_v[(j1-1)*ncol + i1] + 0.1*image_b[(j1-1)*ncol + i1];
  
  gradient = - (double) ( (Y0-Y1)*(Y0-Y1) + (Y0-Y2)*(Y0-Y2) );*/
  
  
  /* Ici : on prend le max de gradient selon les 3 couleurs */
 
  /* Y0 = - (double)
    ( (image_r[j1*ncol + i1]-image_r[j1*ncol + (i1-1)])
      *(image_r[j1*ncol + i1]-image_r[j1*ncol + (i1-1)]) +
      (image_r[j1*ncol + i1]-image_r[(j1-1)*ncol + i1])
      *(image_r[j1*ncol + i1]-image_r[(j1-1)*ncol + i1]) );

  Y1 = - (double)
    ( (image_v[j1*ncol + i1]-image_v[j1*ncol + (i1-1)])
      *(image_v[j1*ncol + i1]-image_v[j1*ncol + (i1-1)]) +
      (image_v[j1*ncol + i1]-image_v[(j1-1)*ncol + i1])
      *(image_v[j1*ncol + i1]-image_v[(j1-1)*ncol + i1]) );

  Y2 = - (double)
    ( (image_b[j1*ncol + i1]-image_b[j1*ncol + (i1-1)])
      *(image_b[j1*ncol + i1]-image_b[j1*ncol + (i1-1)]) +
      (image_b[j1*ncol + i1]-image_b[(j1-1)*ncol + i1])
      *(image_b[j1*ncol + i1]-image_b[(j1-1)*ncol + i1]) );

  gradient =Y0;
  if (Y1 > gradient)
    gradient=Y1;
  if (Y2 > gradient)
    gradient=Y2;*/
    

  /* Ici : on reste dans un seul canal : le vert car c'est en general */
  /*   le plus contraste et representatif */
 
  gradient = - (double)
    ( (image_v[j1*ncol + i1]-image_v[j1*ncol + (i1-1)])
      *(image_v[j1*ncol + i1]-image_v[j1*ncol + (i1-1)]) +
      (image_v[j1*ncol + i1]-image_v[(j1-1)*ncol + i1])
      *(image_v[j1*ncol + i1]-image_v[(j1-1)*ncol + i1]) );

  gradient *= W_grad;
 

 /* intensite = (double) Y0;*/

  return (gradient);

}


/******************/
/* Energie totale */
/******************/
double Etot ( int t, int code_precedent, int code, int code_suivant,
	      double Energie_prec, unsigned char* image_r, unsigned char *image_v, unsigned char *image_b, int ncol)
{
  double NRJ;
  NRJ = Energie_prec
    + Alpha * derivee(t, code_precedent, code)
    + Beta * derivee2(t, code_precedent, code, code_suivant)
    + Eext(t, code, image_r, image_v, image_b, ncol);
  return (NRJ);
}

/*******************************************/
/* Fonction d'ecriture de l'image resultat */
/*******************************************/

int ecrit_image (unsigned char *image, int ncol, int nlig)
{
  FILE *fp;
  int i, j;


 /* ecriture de l'image (.ima) */
  fp=fopen(OUTPUT_IMA,"w");
  for (i=0;i<nlig;i++)
    for(j=0;j<ncol;j++)
      fwrite(image+i*ncol+j,sizeof(char),1,fp);
  fclose(fp);
  
  /* ecriture du fichier .dim associe */

  fp=fopen(OUTPUT_DIM,"w");
  fprintf(fp,"%d %d",ncol,nlig);
  fclose(fp);

  return(0);
}

/***********************************************************/
/* Fonction d'ecriture de l'image resultat en couleur !!!! */
/***********************************************************/

int ecrit_image_couleur (unsigned char *image_r, unsigned char *image_v, unsigned char *image_b, int ncol, int nlig)
{
  FILE *fp;
  int i, j;


  /* ecriture du fichier .dim associe */

  fp=fopen(OUTPUT_DIM,"w");
  fprintf(fp,"%d %d",ncol,nlig);
  fclose(fp);

 /* ecriture de l'image (.rvb) */
  fp=fopen(OUTPUT_RVB,"w");
  for (i=0;i<nlig;i++)
    for(j=0;j<ncol;j++)
      {
	fwrite(image_r+i*ncol+j,sizeof(char),1,fp);
	fwrite(image_v+i*ncol+j,sizeof(char),1,fp);
	fwrite(image_b+i*ncol+j,sizeof(char),1,fp);
      }
  fclose(fp);
  

  return(0);
}


/********************************/
/* Fonction de trace de segment */
/********************************/

int segment(int nlig, int ncol,unsigned char *image, int i0, int j0, int i1, int j1)
{
  int i, j;
  float a1, b1, a2, b2;
  
  if (i0 > i1)
    {
      i = i1; i1 = i0; i0 = i;
      i = j1; j1 = j0; j0 = i;
    }

  /* On fait des points pour voir les extremites du segment */
  image [i0*ncol + j0] = 255;
  image [i0*ncol + j0-1] = 255;
  image [i0*ncol + j0+1] = 255;
  image [(i0-1)*ncol + j0] = 255;
  image [(i0-1)*ncol + j0-1] = 255;
  image [(i0-1)*ncol + j0+1] = 255;
  image [(i0+1)*ncol + j0] = 255;
  image [(i0+1)*ncol + j0-1] = 255;
  image [(i0+1)*ncol + j0+1] = 255;

  image [i1*ncol + j1] = 255;
  image [i1*ncol + j1-1] = 255;
  image [i1*ncol + j1+1] = 255;
  image [(i1-1)*ncol + j1] = 255;
  image [(i1-1)*ncol + j1-1] = 255;
  image [(i1-1)*ncol + j1+1] = 255;
  image [(i1+1)*ncol + j1] = 255;
  image [(i1+1)*ncol + j1-1] = 255;
  image [(i1+1)*ncol + j1+1] = 255;

  /* Si les 2 points sont confondus, on s'arrete !!! */
  if( (i0==i1) && (j0==j1))
      return(0);
  else
    {

      /* Calcul des parametres de la droite : J = a1*I + b1 */
      a1 = ((float)j1 - (float)j0) / ( (float)i1 - (float)i0);
      b1 = ((float)i1 * (float)j0 - (float)j1 * (float)i0)/((float)i1 - (float)i0);
      
      /* Modification de l'image */
      if (a1 >= 1)
	{
	  
	  /* j decrira l'espace et non i sinon il y a des trous !! */
	  
	  a2 = ((float)i1 - (float)i0) / ( (float)j1 - (float)j0);
	  b2 = ((float)j1 * (float)i0 - (float)i1 * (float)j0)/((float)j1 - (float)j0); 
	  for (j = j0; j<= j1; j++)
	    {
	      i = (int)(a2* (float)j + b2);
	      image [i*ncol + j] = 255;
	    }
	}
      else
	{
	  if (a1 <= -1)
	    
	    /* j decrira l'espace et non x sinon il y a des trous !! */
	    
	    {
	      a2 = ((float)i1 - (float)i0) / ( (float)j1 - (float)j0);
	      b2 = ((float)j1 * (float)i0 - (float)i1 * (float)j0)/((float)j1 - (float)j0); 
	      for (j = j0; j>= j1; j--)
		{
		  i = (int)(a2* (float)j + b2);
		  image [i*ncol + j] = 255;
		}
	    }
	  else
	    {
	      for (i = i0; i<= i1; i++)
		{
		  j = (int)(a1* (float)i + b1);
		  image [i*ncol + j] = 255;
		}
	    }
	}
      return(0);

    }

}