#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include  <time.h>
#include "hdf.h"
#include <dirent.h>
#include <limits.h>
#define MIN_SIG -32.0
#define MAX_SIG -5.0

#define MAX_LON 359.99

#define MIN_PIX 5
#define MAX_PIX 255

#define TERRE 1027.0

#define TEMP_FILE_NAME "appendfile"

#define MAX_COLS  810



struct save {
	float64 lat;
	float64 lon;
	float64 sigma0;
	float64 incidence;
	float64 surface;
	float64 sigma0_attn_map;
};
typedef struct save data;


struct param {
	int date_deb;
	int nbjour;
	double maille;
	double lon[2];
    double lat[2];
    char * dossier;
    char * sortie;
    char * appendfilename;
};
typedef struct param parameter;
 
 
     int getParameter(int argc,char *argv[],parameter *p){
             if(argc==11){
                 p->date_deb=atoi(argv[1]);
                 p->nbjour=atoi(argv[2]);
                 p->maille=atof(argv[3]);
                 p->lon[0]=atof(argv[4]);
                 p->lon[1]=atof(argv[5]);
                 p->lat[0]=atof(argv[6]);
                 p->lat[1]=atof(argv[7]);
                 /* BUG
                  * attention, la dans le code original on faisait un strlen()
                  * sur l'adresse de argv[x] + 1, donc on fait un acces memoire illegal.
                  * Confusion avec strlen(argv[x]) + 1 pour allouer suffisamment de
                  * place pour la chaine de caracteres et son caractere NUL terminal.
                  */
                 p->dossier=(char*)malloc(sizeof(char)*(strlen(argv[8])+1));
                 strcpy(p->dossier,argv[8]);
                 p->sortie=(char*)malloc(sizeof(char)*(strlen(argv[9])+1));
                 strcpy(p->sortie,argv[9]);
                 p->appendfilename=(char*)malloc(sizeof(char)*(strlen(argv[10])+1));
                 strcpy(p->appendfilename,argv[10]);
             }   else{
                 fprintf(stderr,"./syntheseW [date_deb][nb_jours][maille][long_min][long_max][lat_min][lat_max][dossier entre][dossier sortie][appendfile]\n");
             return -1;
             }         
             return 0;
             }  

 float64 correctS(float64 sigma, float64 sigma_attn, float64 incidence){
	return (sigma + 2 * (sigma_attn/cos((double)((2*3.141592654*incidence)/360.))) );
}


double getPixelFromNumber(double val){
       return val;              
}                       

double getPixelValue4(double val){
       /* pour obtenir des couleurs claires*/
    return ((val-MIN_SIG)*(250))/(MAX_SIG-MIN_SIG)+5;
}  
double INVgetPixelValue4(double val){
      return ((val-5.0)*(MAX_SIG-MIN_SIG)/250.0)+MIN_SIG;
}

double getPixelValue5(double val){
       /* pour obtenir des couleurs sombres*/
    return (260-(((val-MIN_SIG)*(250))/(MAX_SIG-MIN_SIG)+5));
}  


int histoStat(char *picturename,char *statfilename,int linecol[2]){
    
    int i;
    FILE *file=fopen(statfilename,"ab");
    FILE *img=fopen(picturename,"r");
    
  float tmp;
  int p5pourcent;
  int p5p=-1;
  int p95pourcent;
  int p95p=-1;
  int histo[256];
    int histoC[256];
    int max=INT_MIN;
    int min=INT_MAX;
    
    for(i=0;i<256;++i){
          histo[i]=0;     
          histoC[i]=0;                                       
   }
                                     
    unsigned char pic[linecol[0]*linecol[1]];
    fread(pic,sizeof(char),linecol[0]*linecol[1],img);
    
   int current;
   int nbvalue;
   for(i=0;i<linecol[0]*linecol[1];++i){
   current=(int)pic[i];
   if(current>4){
   nbvalue++;
   if(current>max)  max=current;
   else if(current<min)  min=current;   
   }                                  
}
    
  for(i=0;i<linecol[0]*linecol[1];++i){
                                       
  current=(int)pic[i];
  if(current>4){
  histo[current]++;            
}        
  } 
  
  
  histoC[5]=histo[5];
  
  for(i=6;i<256;++i){
  histoC[i]=histoC[i-1]+ histo[i];                   
  }             
  
  
  /*On commence a 5 pour ne pas compter les points de la mer */
   for(i=5;i<256;++i){   
   if(histo[i]>histo[max]) max=i; 
   tmp=histoC[i]/nbvalue;
   tmp*=100;
   
   if(tmp>=5 && p5p==-1){
       p5pourcent=i;
       p5p=0;                             
   }
   
    if(tmp>=95 && p95p==-1){
       p95pourcent=i;
       p95p=0;                             
   }
   
   }
printf("%d\n",nbvalue);
for(i=0;i<50;++i){printf("%d: %d      %d     %f\n",i,histo[i],histoC[i],(histoC[i]*1.0)/(nbvalue*1.0));}

fprintf(file,"*********************\n");
fprintf(file,"Image: %s \n",picturename);
fprintf(file,"min pix: %d \t max pix: %d \n",min,max);
//float nb=histo[p5pourcent]/nbvalue;
float nb=(histo[p5pourcent]*1.0)/(nbvalue*1.0);
fprintf(file,"à 5pourcent:%d mesures(%f pourcent)\n",histo[p5pourcent],nb);
nb=(histo[p95pourcent]*1.0)/(nbvalue*1.0);
fprintf(file,"à 95pourcent:%d mesures(%f pourcent)\n",histo[p95pourcent],nb);
nb=(histo[max]*1.0)/(nbvalue*1.0);
fprintf(file,"max pour %d: %d mesures(%f pourcent)\n",max,histo[max],nb);
fprintf(file,"*********************\n");
return 1;                             
}    

int Statistiques(char *picturename,char *statfilename,int linecol[2],double incidence[2]){

int i;

int current;

int max=INT_MIN;
int min=INT_MAX;
double moyenne=0;
double variance=0;
double ecart_type;
int nbvalue=0;

FILE *file=fopen(statfilename,"ab");
FILE *img=fopen(picturename,"r");

unsigned char pic[linecol[0]*linecol[1]];
fread(pic,sizeof(char),linecol[0]*linecol[1],img);

for(i=0;i<linecol[0]*linecol[1];++i){
   current=(int)pic[i];
   if(current>4){
   nbvalue++;
   if(current>max)  max=current;
   else if(current<min)  min=current;   
   moyenne+=current; 
   }                                  
}
moyenne=moyenne/(nbvalue);

for(i=0;i<linecol[0]*linecol[1];++i){
current=(int)pic[i];
if(current>4){
  variance+=(pow(current-moyenne,2));   
}                                
}
variance=variance/(nbvalue);
ecart_type=sqrt(variance);

fprintf(file,"*********************\n");
fprintf(file,"Image: %s \n",picturename);
fprintf(file,"min pix: %d \t max pix: %d \n",min,max);
fprintf(file,"min val: %f \t max val: %f \n",INVgetPixelValue4(min),INVgetPixelValue4(max));
fprintf(file,"moyenne pix: %f \n",moyenne);
fprintf(file,"moyenne val: %f \n",INVgetPixelValue4(moyenne));
fprintf(file,"ecart-type pix: %f \n",ecart_type);
fprintf(file,"incidence min : %f \t incidence max : %f \n",incidence[0],incidence[1]);


fclose(img);
fclose(file);

return 0;
}   

 
 /*
 *non alloué lors de l'appel de la fonction*
  double **somme: tableau de pixel contenant( apres execution de la methode) pour chaque pixel i la somme des sigma0 naturelles
  double **nbval: tableau indiquant (apres execution de la methode) le nombre de valeurs additionnées pour le pixel i
  ====> Il ne faut donc pas oublié de les liberer apres l'execution
 
 *directement utilisable*
 double *lat:  [0]= lat min,[1]=lat max
 double *lon:  [0]= lon min,[1]=lon max
 int *lignecol:  [0]= nbr de ligne,[1]= nbr de colonne
 double maille: taille de la maille
 FILE *file: fichier d'entré
 
 Return: nombre d'entré utilisable
 
 Cette fonction récupère les données dans le flux *file et remplit le tableau somme et nbval
 Ainsi pour un pixel i donné il suffit de faire somme[i]/nbval[i] pour obtenir la moyenne des sigma0 naturelles associées à ce pixel
 */    
 int readFile(double **sommeHH,double **nbvalHH,double **sommeVV,double **nbvalVV,parameter p,int *lignecol,char *filename,double **incidence){
     data mydata;
     /* FIXME
      * Ajout d'une variable i comme dans CreateImgEcart() pour simplifier
      * les expressions.
      */
     int x,y,i;
     int k=0;
     FILE *file;
     
     /*printf("coco: %s\n",filename);*/
     /*file=fopen(filename,"r");*/
     
     file=fopen(filename,"r");
     if(file==NULL) return 1;
     
     /*
     NE PAS OUBLIER DE LIBERER CES TABLEAUX
     */
     /* BUG
      * Cette fonction lit et ecrit un float64 trop loin dans ces tableaux.
      * Comme il s'agit peut etre d'un comportement correct, par exemple
      * si on part de x=1 au lieu de x=0, on alloue ici la place necessaire
      * a cette valeur supplementaire pour eviter un acces memoire illegal.
      */
     *nbvalHH=(double*)calloc(lignecol[0]*lignecol[1] +1,sizeof(double));
     *sommeHH=(double*)calloc(lignecol[0]*lignecol[1] +1,sizeof(double));
     *nbvalVV=(double*)calloc(lignecol[0]*lignecol[1] +1,sizeof(double));
     *sommeVV=(double*)calloc(lignecol[0]*lignecol[1] +1,sizeof(double));
  
   printf("******PROCESSING********\n");
   
   
     while(feof(file)==0)
     {

     /*Pour chaque entrée dans le fichier concatené*/
       fread(&mydata, sizeof(data),1,file);
       double longi; 
       /*cette manipulation est necessaire pour passer la longitude dans le meme systeme que la lat*/

       if(mydata.lon>p.lon[1]){longi=(mydata.lon-MAX_LON); }
         else {longi=mydata.lon;}
         
      /*On verifie que l'entrée est dans la zone souhaité*/ 
      
     if((mydata.lat<=p.lat[1] && mydata.lat>=p.lat[0])&&(longi<=p.lon[1] && longi>=p.lon[0]) )
         {
        /* printf("surface %f  \n",mydata.surface);  */ 
     /* calcule des coordonnées en pixel a partir des long/lat de l'entré*/            
         x=fabs((p.lon[0]-longi)/p.maille);
         y=fabs((p.lat[1]-mydata.lat)/p.maille);
         /* FIXME
          * Pour simplifier les expressions on utilise i pour les coordonnees
          * du pixel courant.
          */
         i = y*lignecol[1]+x;
              
         /*printf("********************\n");
         printf("lat= %f  lon=%f\n",mydata.lat,mydata.lon);
         printf("sigma0db= %f \n",mydata.sigma0);
         printf("sigma0= %f \n",pow(10.0,mydata.sigma0/10));
         printf("%d %d  \n",x,y);*/
          if (mydata.surface==TERRE){
          ++k;
          if(mydata.incidence<50){
         (*sommeHH)[i]+=pow(10.0,correctS(mydata.sigma0,mydata.sigma0_attn_map,mydata.incidence)/10);
         (*nbvalHH)[i]+=1;
         if(mydata.incidence<(*incidence)[0]) (*incidence)[0]=mydata.incidence;
         else if (mydata.incidence>(*incidence)[1]) (*incidence)[1]=mydata.incidence;
         }
         else{
         (*sommeVV)[i]+=pow(10.0,correctS(mydata.sigma0,mydata.sigma0_attn_map,mydata.incidence)/10);
         (*nbvalVV)[i]+=1;
         if(mydata.incidence<(*incidence)[2]) (*incidence)[2]=mydata.incidence;
         else if (mydata.incidence>(*incidence)[3]) (*incidence)[3]=mydata.incidence;
         }     
         
         }
         /*Si l'entré est dans la mer*/
          else{  
              if(mydata.incidence<50){
              (*sommeHH)[i]=-1;
              (*nbvalHH)[i]=-1;
              }  
              else{
              (*sommeVV)[i]=-1;
              (*nbvalVV)[i]=-1;
              }
         } 
         
          
         }
  
     } 
      fclose(file);
     return k;
 }    
 
     
 int writePicture(double *somme,double *nbval,int lignecol[2],char *filename, char *nbvaluefilename ){
     
     int tmmp=0;
     
      FILE *out=fopen(filename,"w");
      FILE *out2=fopen(nbvaluefilename,"w");
      int i;
      for(i=0;i<lignecol[0]*lignecol[1];++i){
                                             
      if(nbval[i]>0){            
          double tmp=getPixelValue4(10.0*log10(somme[i]/nbval[i]));
          unsigned char p=(unsigned char)tmp;
          fprintf(out,"%c",p); 
          fprintf(out2,"%c",(unsigned char)getPixelFromNumber(nbval[i])); 
      }
      
      /*Pas de valeur pour ce pixel*/
      else if(nbval[i]==0) {
          tmmp++;
          fprintf(out,"%c",4);
          fprintf(out2,"%c",0);
      } 
             
      /*Mer*/     
      else if(nbval[i]<0) {
     /* printf("dans l eau: %d\n",i);*/
         fprintf(out,"%c",0);
         fprintf(out2,"%c",0);
         }  
      }  
      fclose(out);
      fclose(out2);
      printf("il ya %d pixel sans valeur dans %s\n",tmmp,filename);
     return 1;
}               


int CreateImgEcart(double *sommeHH,double *nbvalHH,double *sommeVV,double *nbvalVV,parameter p,int *lignecol,char *filename,char *FecartypeVV,char *FecartypeHH,double **incidenceEcart)
{
double *ecartypeHH;
double *ecartypeVV;
data mydata;
int x,y;
/* BUG
 * i etait utilise sans etre initialise ici, d'ou la segmentation fault !
 * On a suppose que i correspondait aux coordonnees du pixel en cours de
 * traitement; si ce n'est pas le cas, il faut lui affecter la valeur
 * qu'il devait avoir pour ce calcul.
 */
int i;   
FILE *outHH=fopen(FecartypeHH,"w");
FILE *outVV=fopen(FecartypeVV,"w");
FILE *file=fopen(filename,"r");

     /* BUG
      * Meme probleme que dans ReadFile(), on accede a un float64 de trop.
      * Meme bidouille pour eviter de lire de la memoire qui ne nous appartient
      * pas.
      */
     ecartypeHH=(double*)calloc(lignecol[0]*lignecol[1] +1,sizeof(double));
     ecartypeVV=(double*)calloc(lignecol[0]*lignecol[1] +1,sizeof(double));
  
     while(feof(file)==0)
     {
     /*Pour chaque entrée dans le fichier concatené*/
       fread(&mydata, sizeof(data),1,file);

       double longi; 
       /*cette manipulation est necessaire pour passer la longitude dans le meme systeme que la lat*/

       if(mydata.lon>p.lon[1])longi=(mydata.lon-MAX_LON); 
         else longi=mydata.lon;
         
      /*On verifie que l'entrée est dans la zone souhaité*/ 
      
     if((mydata.lat<=p.lat[1] && mydata.lat>=p.lat[0])&&(longi<=p.lon[1] && longi>=p.lon[0]) )
         {
        /* printf("surface %f  \n",mydata.surface);  */ 
     /* calcule des coordonnées en pixel a partir des long/lat de l'entré*/            
         x=fabs((p.lon[0]-longi)/p.maille);
         y=fabs((p.lat[1]-mydata.lat)/p.maille);
         i = y*lignecol[1]+x;
              
         /*printf("********************\n");
         printf("lat= %f  lon=%f\n",mydata.lat,mydata.lon);
         printf("sigma0db= %f \n",mydata.sigma0);
         printf("sigma0= %f \n",pow(10.0,mydata.sigma0/10));
         printf("%d %d  \n",x,y);*/

         /* FIXME
          * Pour simplifier les expressions et toujours sur la supposition que
          * i correspondait aux coordonnees du pixel courant, on utilise i partout.
          */
          if (mydata.surface==TERRE){
          if(mydata.incidence<50){
         ecartypeHH[i]+=pow(pow(10.0,correctS(mydata.sigma0,mydata.sigma0_attn_map,mydata.incidence)/10)-sommeHH[i]/nbvalHH[i],2);
           if(mydata.incidence<(*incidenceEcart)[0]) (*incidenceEcart)[0]=mydata.incidence;
           else if (mydata.incidence>(*incidenceEcart)[1]) (*incidenceEcart)[1]=mydata.incidence;
         }
         else{
         ecartypeVV[i]+=pow(pow(10.0,correctS(mydata.sigma0,mydata.sigma0_attn_map,mydata.incidence)/10)-sommeVV[i]/nbvalVV[i],2);
         if(mydata.incidence<(*incidenceEcart)[2]) (*incidenceEcart)[2]=mydata.incidence;
         else if (mydata.incidence>(*incidenceEcart)[3]) (*incidenceEcart)[3]=mydata.incidence;
         }     
         
         }
         /*Si l'entré n'est pas dans la zone souhaité*/
          else{  
              if(mydata.incidence<50){
              ecartypeHH[i]=-1;
              }  
              else{
              ecartypeVV[i]=-1;    
              }
         } 
         }
  
     } 
     
     
     for(i=0;i<lignecol[0]*lignecol[1];++i){
                                             
     if(nbvalHH[i]>0){            
          double tmp=getPixelValue4(10.0*log10(ecartypeHH[i]/nbvalHH[i]));
          unsigned char p=(unsigned char)tmp;
          fprintf(outHH,"%c",p); 
        }  
        /*Pas de valeur pour ce pixel*/
        else if(nbvalHH[i]==0) {
          fprintf(outHH,"%c",4);
      }
          
      else if(nbvalHH[i]<0) {
   /* printf("dans l eau: %d\n",i);*/
         fprintf(outHH,"%c",0);
         }
          
     if(nbvalVV[i]>0){            
          double tmp=getPixelValue4(10.0*log10(ecartypeVV[i]/nbvalVV[i]));
          unsigned char p=(unsigned char)tmp;
          fprintf(outVV,"%c",p); 
        }  
        /*Pas de valeur pour ce pixel*/
        else if(nbvalVV[i]==0) {
          fprintf(outVV,"%c",4);
      }
       else if(nbvalVV[i]<0) {
     /* printf("dans l eau: %d\n",i);*/
         fprintf(outVV,"%c",0);
         }
         
      }
     
    
      fclose(file);
      fclose(outHH);
      fclose(outVV);
      free(ecartypeVV);
      free(ecartypeHH);
     return 1;   
}
 

int main(int argc,char *argv[])
{
    double *incidence;
    double *incidenceEcart;
    incidence=(double*)malloc(sizeof(double)*4);
    incidenceEcart=(double*)malloc(sizeof(double)*4);
    incidence[0]=INT_MAX;
    incidence[1]=INT_MIN;
    incidence[2]=INT_MAX;
    incidence[3]=INT_MIN;
    incidenceEcart[0]=INT_MAX;
    incidenceEcart[1]=INT_MIN;
    incidenceEcart[2]=INT_MAX;
    incidenceEcart[3]=INT_MIN;
    parameter p;
    int m=getParameter(argc,argv,&p);
    if(m==-1) return -1;
    int lignecol[2];
    
     lignecol[0]=fabs((p.lat[1]-p.lat[0])/p.maille);
     lignecol[1]=fabs((p.lon[1]-p.lon[0])/p.maille);
  
    double *nbvalHH;
    double *sommeHH;
    double *nbvalVV;
    double *sommeVV;
   
    int i,j;
    int tail;
  
    tail=strlen(p.appendfilename)+strlen(p.sortie)+1;
    char appendfile[tail];
    for(i=0;i<tail;++i)  appendfile[i]='\0';
    strcat(appendfile,p.sortie);
    strcat(appendfile,p.appendfilename);
    
    i=strlen(p.appendfilename)-10;
     char suffix[i];
    for(j=0;j<i;++j) suffix[j]='\0';
    strcpy(suffix,p.appendfilename+i);
   
    printf("********************\n");
    printf("Chargement des données\n");

    int k=-1;
    k=readFile(&sommeHH,&nbvalHH,&sommeVV,&nbvalVV,p,lignecol,appendfile,&incidence);
    if(k==1) return -1;
    printf("********************\n");
    printf("Creation de l'image\n");
    
  /* fprintf(stderr,"suffixe:%s\n",suffix);*/
   
    tail=strlen("_imageHH.raw")+strlen(p.sortie)+1+strlen(suffix); 
    char imageHH[tail];
     for(i=0;i<tail;++i)  imageHH[i]='\0';
    strcat(imageHH,p.sortie);
    strcat(imageHH,suffix);
    strcat(imageHH,"_imageHH.raw");


    tail=strlen("_imageVV.raw")+strlen(p.sortie)+1+strlen(suffix);
    char imageVV[tail];
    for(i=0;i<tail;++i)  imageVV[i]='\0';
    strcat(imageVV,p.sortie);
    strcat(imageVV,suffix);
    strcat(imageVV,"_imageVV.raw");

    
    tail=strlen("_nbvalueHH.raw")+strlen(p.sortie)+1+strlen(suffix);
    char nbvalueHH[tail];
    for(i=0;i<tail;++i)  nbvalueHH[i]='\0';
    strcat(nbvalueHH,p.sortie);
    strcat(nbvalueHH,suffix);
    strcat(nbvalueHH,"_nbvalueHH.raw");

    
    tail=strlen("_nbvalueVV.raw")+strlen(p.sortie)+1+strlen(suffix); 
    char nbvalueVV[tail];
    for(i=0;i<tail;++i)  nbvalueVV[i]='\0';
    strcat(nbvalueVV,p.sortie);
    strcat(nbvalueVV,suffix);
    strcat(nbvalueVV,"_nbvalueVV.raw");
     
   
    tail=strlen("_ecartypeHH.raw")+strlen(p.sortie)+1+strlen(suffix);    
    char ecartypeHH[tail];
    for(i=0;i<tail;++i)  ecartypeHH[i]='\0';
    strcat(ecartypeHH,p.sortie);
    strcat(ecartypeHH,suffix);
    strcat(ecartypeHH,"_ecartypeHH.raw");  
    
    char ecartypeVV[tail];
    for(i=0;i<tail;++i)  ecartypeVV[i]='\0';
    strcat(ecartypeVV,p.sortie);
    strcat(ecartypeVV,suffix);
    strcat(ecartypeVV,"_ecartypeVV.raw");
    
     fprintf(stderr,"%s\n%s\n%s\n%s\n",ecartypeHH,ecartypeVV,nbvalueVV,nbvalueHH);
    
    writePicture(sommeHH,nbvalHH,lignecol,imageHH,nbvalueHH);
    writePicture(sommeVV,nbvalVV,lignecol,imageVV,nbvalueVV);
    CreateImgEcart(sommeHH,nbvalHH,sommeVV,nbvalVV,p,lignecol,appendfile,ecartypeVV,ecartypeHH,&incidenceEcart);

    printf(" col=%d ligne=%d\n\n",lignecol[1],lignecol[0]);
    printf("Image cree à partir de:  %d valeurs\n",k);     
    /*printf("Incidence minimal: %f\nIncidence maximal:%f\n\n\n",incidence[0],incidence[1]);                      
     */
     tail=strlen("_stat.txt")+strlen(p.sortie)+1+strlen(suffix);
    char stat[tail];
    for(i=0;i<tail;++i)  stat[i]='\0';
    strcat(stat,p.sortie);
    strcat(stat,suffix);
    strcat(stat,"_stat.txt");
     
    Statistiques(imageHH,stat,lignecol,incidence);
    Statistiques(imageVV,stat,lignecol,incidence+2);
    Statistiques(ecartypeHH,stat,lignecol,incidenceEcart);
    Statistiques(ecartypeVV,stat,lignecol,incidenceEcart+2);
    //histoStat(nbvalueHH,stat,lignecol);
    //histoStat(nbvalueVV,stat,lignecol);
     
free(nbvalHH);
free(sommeHH);
free(nbvalVV);
free(sommeVV);
free(p.dossier); 
free(p.sortie); 
free(p.appendfilename); 
free(incidence);
free(incidenceEcart);
/*free(suffix);*/ /* BUG suffix n'a pas ete alloue ! */


return 1;    
}