/*==============================================================
mesure du rapport S/B entre une image originale et une image bruitée
--------------------------------------------------------------*/

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

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

int src_mx,src_my; /* dimension de l'image originale et de l'image filtree */

gris *src;         /* image originale */
gris *srcb;        /* image dilatee (pour eviter les effets de bords */

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

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


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

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

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

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

  if (argc!=5) {
    printf("   Mesure du rapport S/B\n");
	printf("___________________________\n");
    printf("ssb <image originale> <image_X> <image_Y> <image bruitée> \n\n");
    printf("<image originale> <image bruitée>  : fichiers au format TIF 24bits\n");
    printf("<image_X> <image_Y> : taille de l'image\n");
    exit(0);
  }

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

  /* ouverture du fichier - calcul taille du fichier -> taille de l'entete */
  f=fopen(fin,"rb");
  if (f==NULL) {
	  printf("Impossible d'ouvrir l'image originale\n");
	  exit(0);
  }
  fb=fopen(fbruit,"rb");
  if (fb==NULL) {
	  printf("Impossible d'ouvrir l'image bruitée\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);
  }
  fclose(f);f=fopen(fin,"rb");

  /* allocations memoire values, src, im, dest, entete, et precalc */
  entete=(char *)malloc(lentete);
  src=(gris *)malloc(src_mx*src_my*3);
  srcb=(gris *)malloc(src_mx*src_my*3);
  if (!src) exit(1);
  if (!srcb) exit(1);

  /* lecture de l'entete puis de l'image */
  fread(entete,lentete,1,f);
  fread(src,src_mx,src_my*3,f);
  fclose(f);
  fread(entete,lentete,1,fb);
  fread(srcb,src_mx,src_my*3,fb);
  fclose(fb);

  printf("%s : %d octets = %d (entete) + %d*%d*3 (image)\n",fin,taillef,lentete,src_mx,src_my);

  /* mesure */
  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++) {
	  int s= src[(x+y*src_mx)*3+plan];
	  int d=srcb[(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) * ((double)src[(x+y*src_mx)*3+plan]-moyA);
	  double d =(srcb[(x+y*src_mx)*3+plan]-(double)src[(x+y*src_mx)*3+plan]+moyA-moyB) \
	  * (double)(srcb[(x+y*src_mx)*3+plan]-(double)src[(x+y*src_mx)*3+plan]+moyA-moyB);
	  varA+=s;varB+=d;
  }
  varA/=(src_mx*src_my*3);  varB/=(src_mx*src_my*3);

  printf("référence : S/B pour une image 8bits : 6*8=48dB\n");
  printf("Variance de l'image originale = %d = %ddB\n",(int)varA,(int)(10*log10(varA)));
  printf("Variance du bruit = %d = %ddB\n",(int)varB,(int)(10*log10(varB)));

  return 0;
}