/*==============================================================
      filtre d-alpha (pretraitement)
--------------------------------------------------------------*/

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>

typedef unsigned char gris;		/* niveau de gris */

float alpha;		/* valeur de alpha pour le filtre d-alpha */
int lf,lf2;			/* lf=taille de la fenetre (impair) et lf2=lf/2 */

int src_mx,src_my; /* dimension de l'image originale et de l'image filtree */
int mx,my;         /* dimension de l'image originale dilatee (pour eviter les effets de bords) */

gris *src;         /* image originale */
gris *im;          /* image dilatee (pour eviter les effets de bords */
gris *dest;        /* image filtree */

char *entete;      /* entete du fichier image non compresse */
int lentete;       /*    et sa taille */

gris *values;      /* valeurs des pixels dans la fenetre */
float *precalc;   /* precalculs de la fonction "mise a la puissance alpha" */

double moyA,varA,moyB,varB;  /* moments d'ordre 1 et 2 en entree et en sortie */

/*===================================
    Calcul de la fonction f-alpha
-----------------------------------*/
float falpha(gris p) {
  int k;
  float somme;

  somme=0;
  for (k=0;k<lf*lf;k++)
    somme+=precalc[abs((int)(p)-values[k])];

  return somme;
}

/*============================================
    Filtre en utilisant la fonction f-alpha
--------------------------------------------*/
gris filtre() {
  gris fmin,fmax;		/* intervalle de valeurs a tester par fonction f-alpha */
  float falpha_min;		/* minimum trouve */
  gris falpha_kmin;		/* valeur ou ce minimum est atteint */

  int c;int k;float f;

  /* calcul du minimum et du maximum pour determiner l'intervalle de test */
  fmin=values[0];fmax=values[0];
  for (k=1;k<lf*lf;k++)
    if (values[k]<fmin) fmin=values[k];
    else if (values[k]>fmax) fmax=values[k];

  /* calcul du minimum de falpha sur cet intervalle */
  falpha_min=falpha(fmin); falpha_kmin=fmin;
  for (c=fmin+1;c<=fmax;c++)
	{
      f=falpha((gris)c);
      if (f<falpha_min) { falpha_min=f; falpha_kmin=c; }
	}

  return falpha_kmin;
}

/*===================================
    Programme principal
-----------------------------------*/
int main(int argc, char *argv[])
{
  char *fin,fout[300]; /* noms des fichiers source et destination */

  FILE *f;         /* fichier */
  int taillef;     /* taille de ce fichier */

  fpos_t posd,posf; /* utilises pour calculer la taille du fichier */

  int x,y,xx,yy;  /* pour balayer l'image */
  int i;
  int plan;

  if (argc!=6) {
    printf("     Filtrage d-alpha\n");
	printf("____________________________\n");
    printf("dalpha <image> <image_X> <image_Y> <fenetre> <alpha>\n\n");
    printf("<image>   : fichier a filtrer au format TIF 24bits\n");
    printf("<image_X> <image_Y> : taille de l'image\n");
    printf("<fenetre> : taille de la fenetre pour le filtre (entier impair)\n");
    printf("<alpha>   : parametre du filtre d-alpha (reel positif)\n");
    exit(0);
  }

  /* lecture des parametres sur la ligne de commande -> */
  fin=argv[1];
  sscanf(argv[2],"%d",&src_mx);
  sscanf(argv[3],"%d",&src_my);
  sscanf(argv[4],"%d",&lf); 
  sscanf(argv[5],"%f",&alpha);

  if (strchr(fin,'.')) {
	char st[100]; char sta[100];
	strcpy(st,fin);
	*strchr(st,'.')=0;
    sprintf(sta,"%.1f",alpha);
	if (strchr(sta,'.')) *strchr(sta,'.')=',';
	sprintf(fout,"%s_%d_%s.%s",st,lf,sta,st+strlen(st)+1);
  } else
	sprintf(fout,"%s_%d_%.1f",fin,lf,alpha);
  

  /* ouverture du fichier - calcul taille du fichier -> taille de l'entete */
  f=fopen(fin,"rb");
  if (f==NULL) {
	  printf("Impossible d'ouvrir le fichier\n");
	  exit(0);
  }
    fgetpos(f,&posd);
    fseek(f,0,SEEK_END);
    fgetpos(f,&posf);
    fsetpos(f,&posd);
  taillef=(int)(posf-posd);
  lentete=taillef-src_mx*src_my*3;
  if (lentete<0) {
    printf("Taille de l'image incompatible avec la taille du fichier !\n");
    exit(0);
  }

  /* allocations memoire values, src, im, dest, entete, et precalc */
  lf2=(lf-1)/2;
  if (lf!=(2*lf2+1)) {
    printf("La taille de la fenetre doit etre impaire !\n");
    exit(0);
  }
  mx=src_mx+2*lf2;
  my=src_my+2*lf2;
  values=(gris *)malloc(lf*lf);
  src=(gris *)malloc(src_mx*src_my*3);
  dest=(gris *)malloc(src_mx*src_my*3);
  im=(gris *)malloc(mx*my*3);
  entete=(char *)malloc(lentete);
  precalc=(float *)malloc(520*sizeof(float));
  for (x=0;x<520;x++) precalc[x]=(float)pow(x,alpha);

  /* lecture de l'entete puis de l'image */
  fread(entete,lentete,1,f);
  fread(src,src_mx,src_my*3,f);
  fclose(f);
  printf("%s : %d octets = %d (entete) + %d*%d*3 (image)\n",fin,taillef,lentete,src_mx,src_my);
  printf("fenetre de taille %d avec alpha = %.2f\n",lf,alpha);

  /* recopie de src dans im */
  for (y=0;y<my;y++)
	for (x=0;x<mx;x++)
	{
      if (x<=lf2) xx=0;
      else if (x>=src_mx+lf2) xx=src_mx;
      else xx=x-lf2;
      if (y<=lf2) yy=0;
      else if (y>=src_my+lf2) yy=src_my;
      else yy=y-lf2;
      for (plan=0;plan<3;plan++)
		im[(x+mx*y)*3+plan]=src[(xx+src_mx*yy)*3+plan];
    }
  
  /* boucle de filtrage */
  for (y=lf2;y<my-lf2;y++) {
	printf("\rtraitement en cours ... %02d%%",((y-lf2)*100)/src_my);fflush(stdout);
    for (x=lf2;x<mx-lf2;x++)
	  for (plan=0;plan<3;plan++)
    {
      /* calcul de values */
      i=0;
      for (xx=x-lf2;xx<=x+lf2;xx++)
       for (yy=y-lf2;yy<=y+lf2;yy++)
         values[i++]=im[(xx+mx*yy)*3+plan];
      dest[(x-lf2+src_mx*(y-lf2))*3+plan]=filtre();
	}
  }
  printf("\renregistrement dans %s ...        \n",fout);

  /* ecriture ligne par ligne du fichier de sortie */
  f=fopen(fout,"wb");
  fwrite(entete,lentete,1,f);
  fwrite(dest,src_mx,src_my*3,f);
  fclose(f);

  moyA=0;varA=0;moyB=0;varB=0;
  for (y=0;y<src_my;y++)
  for (x=0;x<src_mx;x++)
for (plan=0;plan<3;plan++)
 {
	  unsigned int s=src[(x+y*src_mx)*3+plan];
	  unsigned int d=dest[(x+y*src_mx)*3+plan];
	  moyA+=s;moyB+=d;
  }
  moyA/=(src_mx*src_my*3);  moyB/=(src_mx*src_my*3);
  for (y=0;y<src_my;y++)
  for (x=0;x<src_mx;x++)
for (plan=0;plan<3;plan++)
 {
	  double s=(src[(x+y*src_mx)*3+plan]-moyA)*(src[(x+y*src_mx)*3+plan]-moyA);
	  double d=(dest[(x+y*src_mx)*3+plan]-moyB)*(dest[(x+y*src_mx)*3+plan]-moyB);
	  varA+=s;varB+=d;
  }
  varA/=(src_mx*src_my*3);  varB/=(src_mx*src_my*3);
  printf("\rVariance en entree et en sortie : %.2e et %.2e\n",varA,varB);



  return 0;
}