/*------------------------------------------------------*/
/* Prog    : TpIFT6150-3-D.c                            */
/* Auteur  : Francois Destrempes                        */
/* Date    :                                            */
/* version :                                            */ 
/* langage : C                                          */
/* labo    : DIRO                                       */
/*------------------------------------------------------*/

/*------------------------------------------------*/
/* FICHIERS INCLUS -------------------------------*/
/*------------------------------------------------*/
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <string.h>

#include "FonctionDemo3.h"

/*------------------------------------------------*/
/* DEFINITIONS -----------------------------------*/  
/*------------------------------------------------*/
#define NAME_IMG_IN  "photograph"

#define NAME_IMG_OUT1 "Image-depart"
#define NAME_IMG_OUT2 "Image-bruitee" 
#define NAME_IMG_OUT3 "Buades-Image-debruitee" 
#define NAME_IMG_OUT4 "Wang-Image-debruitee" 

/*------------------------------------------------*/
/* FONCTIONS -------------------------------------*/                     
/*------------------------------------------------*/
void egal(float** matin, float** matout,int lgth,int wdth)
{int i,j;
     for(i=0;i<lgth;i++) for(j=0;j<wdth;j++) {matout[i][j]=matin[i][j];}}
  

void Filtre_Gaussien(float** Gauss,int foret, int ecartype)
{
  int i,j;
   for(i=0;i<foret;i++) for(j=0;j<foret;j++) //les deplacements dans chaque filtre est la meme//
			   {
			    		    
			    if(i!=4 && j!=4) //condition de ne pas prendre le centre du filtre// 
			      {
				Gauss[i][j]=exp(-1*(pow(i-4,2)+pow(j-4,2))/2*ecartype); //filtre gaussien pareille pour tous//
			      } else {}
			   }}

  
void Buades_debruitage(float** matout, float** matin, float cstebruit,int lgth,int wdth)
{
int i,j;
 int k,h;  //indices virtuelles des filtres voisinages du pixel dinteret//
 int g,f;  //indices des centres de filtres//
 int foret,nbfiltre; 
 int x,y;  //indices des voisins du centre du filtre dans un filtre//
 int w,z; //indices des voisins du centre du filtre dans limage globale//
float** Euclid; 
float** Gauss; 
float**Moy;

 float squarediff,sigma, sumproduit2; 

 float** Produit; 
 float** poids; 

 foret=7.0;
 sigma=1;
 nbfiltre=lgth/foret;//nombre filtre par ligne et par colonne//

    //Allocation memoire pour la matrice image//
 poids=fmatrix_allocate_2d(nbfiltre,nbfiltre);
  Euclid=fmatrix_allocate_2d(lgth,wdth);
  Gauss=fmatrix_allocate_2d(foret,foret);
  Moy=fmatrix_allocate_2d(lgth,wdth);
  Produit=fmatrix_allocate_2d(nbfiltre,nbfiltre);

 //Calcul de la gaussienne// 
 for(i=0;i<foret;i++) for(j=0;j<foret;j++)
			{
			  Gauss[i][j]=0.0;
			}   
for(i=0;i<foret;i++) for(j=0;j<foret;j++) //les deplacements dans chaque filtre est la meme//
			   {
			    		    
			    if(i!=4 && j!=4) //condition de ne pas prendre le centre du filtre// 
			      {
				Gauss[i][j]=exp(-1*(pow(i-4,2)+pow(j-4,2))/2*sigma); //filtre gaussien pareille pour tous//
			      } else {}}
 
 //Calcul de la ponderation//

  for(f=3;f<lgth-3;f++)for(g=3;g<wdth-3;g++) //les centres des filtres//
 {
   squarediff=0.0;
   sumproduit2=0.0;

   //Initialisation a zero//
   for(i=0;i<lgth;i++) for(j=0;j<wdth;j++)
   {
     Euclid[i][j]=0.0;
     Moy[i][j]=0.0;
   }  
   for(k=0;k<nbfiltre;k++) for(h=0;h<nbfiltre;h++)
   {
     poids[k][h]=0.0;
     Produit[k][h]=0.0;
     } 

for(k=0;k<nbfiltre;k++)  for(h=0;h<nbfiltre;h++)//Pour chaque voisinage du pixel//
				  {
	    for(i=0;i<foret;i++) for(j=0;j<foret;j++)
				  {
				    w=f-3+i; //coordonnees reelles du pixel dans limage//
				    z=g-3+j;
					    
	if(i+k*foret < lgth && j+h*foret < wdth) //Ne pas depasser les limites de l'image-Pour tous les filtres en dessous et a droit//
			  {  
			    Produit[k][h]=Produit[k][h]+Gauss[i][j]*(matin[w][z]-matin[i+k*foret][j+h*foret]);//Calcul de la difference Euclidienne//
			  }   
			    			
	if(i-k*foret > 0 && j-h*foret > 0 && i-k*foret<lgth && j-h*foret<wdth )  //les filtres en amont du filtre et au dessus //
			  {
			    Produit[k][h]=Produit[k][h]+Gauss[i][j]*(matin[w][z]-matin[i-k*foret][j-h*foret]); //Calcul de la difference Euclidienne//
			    } 

	if(i+k*foret < lgth && j-h*foret > 0 && j-h*foret<wdth) //Ne pas depasser les limites de l'image-les filtres en amont et en dessous//
			  {  
			    Produit[k][h]=Produit[k][h]+Gauss[i][j]*(matin[w][z]-matin[i+k*foret][j-h*foret];//Calcul de la difference Euclidienne//
			  } 

	if(abs(i-k*foret) > 0 && j+h*foret < wdth) //Ne pas depasser les limites de l'image-les filtres au dessus et a droite//
			  {  
			    Produit[k][h]=Produit[k][h]+Gauss[i][j]*(matin[w][z]-matin[i-k*foret][j+h*foret]);//Calcul de la difference Euclidienne//
			  } 
			  
 else
			       {}

				  }}     
//Calcul de la moyenne des difference euclidiennes//
for(k=0;k<nbfiltre;k++) for(h=0;h<nbfiltre;h++)
      {
	squarediff=squarediff+exp(-1*pow(Produit[k][h],2));
      }				 				 

//Calcul final du poids//
	 for(k=0;k<nbfiltre;k++) for(h=0;h<nbfiltre;h++)
				   { 							         
				     w=k*foret+3; 
				     z=h*foret+3;
			   	     poids[k][h]=1/squarediff*exp(-1*Produit[k][h]/0.15);
				     sumproduit2=sumproduit2+poids[k][h]*matin[w][z];
				   }
       matout[f][g]=sumproduit2;}

  //On bouche les trous//
  for(f=0;f<3;f++) for (g=0;g<3;g++)
		     {matout[f][g]=matin[f][g];}
  for (f=lgth-3;f<lgth;f++) for(g=wdth-3;g<wdth;g++)
			      {matout[f][g]=matin[f][g];}
}




  

int main(int argc,char** argv)
 {
  int i,j,k,l;
  int length,width;
  float var,Som1,Som2,ISNR;
  int size_filter;
  float** MatriceImage;  /* image de depart */
  float** MatB;  /*image bruitee*/
  float** IF1; /*Image debruitee par Buades*/
  float** w; /*le poid moyen de tous les pixels dans image*/
   
  
  /*Allocation memoire pour la matrice image*/
  MatriceImage=LoadImagePgm(NAME_IMG_IN,&length,&width);
  MatB=fmatrix_allocate_2d(length,width);
  IF1=fmatrix_allocate_2d(length,width);
  w=fmatrix_allocate_2d(length,width);

  
/*Initialisation des matrices*/
for(i=0;i<length;i++) for(j=0;j<width;j++)
{MatB[i][j]=0.0; IF1[i][j]=0.0;w[i][j]=0.0;}
  
  printf("Entrez la variance du bruit : ");
  scanf("%f",&var);

  /* ouvrir l'image d'entree */
  //system("gimp photograph.pgm&");
 
  /* ajouter du bruit a l'image d'entree (add_gaussian_noise) */
 //Ajout du bruit a l'image//
 egal(MatriceImage,MatB,length,width);
 add_gaussian_noise(MatB, length, width, var);
 

/*Debruitage Buades*/
Buades_debruitage(IF1,MatB,var,length,width);



 

  /* Sauvegarde des matrices sous forme d'image pgm */
SaveImagePgm(NAME_IMG_OUT1,MatriceImage,length,width);  
SaveImagePgm(NAME_IMG_OUT2,MatB,length,width);
SaveImagePgm(NAME_IMG_OUT3,IF1,length,width);

  
  /* Liberation memoire pour les matrices */
  free_fmatrix_2d(MatriceImage); 
    free_fmatrix_2d(MatB); 
      free_fmatrix_2d(IF1); 
        free_fmatrix_2d(w); 
 
  /*Commande systeme: visualisation de Ingout.pgm*/
	/* system("gimp Image-bruitee.pgm&"); 
 system("gimp Buades-Image-debruitee.pgm&");
 
  /*retour sans probleme*/ 
  printf("\n C'est fini ... \n\n\n");
  return 0; 	 
}