///////////////////////////////////////////////////////
//DETECTABILITE D'UNE ANOMALIE DANS UNE IMAGE DE BRUIT
//Auteur: F. Richard
//Date: 28/04/2006
///////////////////////////////////////////////////////

//Ce programme permet d'evaluer l'algorithme EM  
//pour la detection d'une structure anormale présente dans une image de bruit

//Le probleme est le suivant:
//On simule une image de bruit contenant une anomalie. La region occupee par l'anomalie est un cercle dont le centre et le rayon sont tires de maniere aleatoire. Cette region est decrite par l'image anomalie_carte. L'image de bruit contenant l'anomalie est definie par 
//image(s)= {b0(s)    (s'il n'y a pas d'anomalie a la position s)
//          {b1(s)    (s'il y'a une anomalie a la position s),
//ou b0(s) est la realisation d'une v.a. gaussienne de variance sig0 et de moyenne nulle, 
//et b1(s) est la realisation d'une v.a. gaussienne de variance sig1 et de moyenne ng_anomalie. 

//L'objectif est de detecter l'anomalie  dans l'image, tout en estimant correctement les differents parametres du modele!

//On compare les resultats de detection obtenus  
//- par test du rappport de vraisemblance, lorsque les parametres du modele sont connus.
//- par algorithme EM, lorsque les parametres sont inconnus.
//On compare egalement les resultats d'estimation des parametres obtenus 
//- par maximum de vraisemblance, lorsque l'anomalie est localisee
//- par algorithme EM, lorsque l'anomalie n'est pas localisee






//////////////////////////////////////////////
//PARAMETRES GENERAUX D'IMAGES
//////////////////////////////////////////////
//taille de l'image (nb de lignes=nb de colonnes=taille_image)
taille_image=128;
//valeur maximale des niveaux de gris
ng_max=256;



///////////////////////////////////////////////
//PARAMETRES DE SIMULATION 
///////////////////////////////////////////////

//variance du bruit sur l'image
sig0=80;
//parametre anomalie
//niveau de gris (additif) d'une anomalie
ng_anomalie_max=50;
//variance du bruit sur l'anomalie 
sig1=80;
//tailles minimale et maximale de l'anomalie (rayon max)
taille_anomalie_min=5;
taille_anomalie_max=20;



//nb d'experiences
nb_expe=100;

//nb d'iterations maximum pour l'algorithme EM
nb_iter_max=50;



getf("test_rapport_vraisemblance.sci");
getf("max_vraisemblance.sci");
getf("algo_em_mixture.sci");

xset("window",0);
xbasc();
xset("window",1);
xbasc();
xset("window",2);
xbasc();



//initialisation des variables contenant les resultats
//vraies/fausses detections/erreurs de classement (VD/FD/EC) avec le test du rapport de vraisemblance (RV) et l'algo EM (version 1 et 2)
VD_RV=[]; VD_EM=[]; VD_EM2=[]; VD_EM3=[]; FD_RV=[]; FD_EM=[]; FD_EM2=[]; FD_EM3=[]; EC_RV=[]; EC_EM=[]; EC_EM2=[]; EC_EM3=[];
//erreurs d'estimations de chq parametre du modele avec le max. de vraisemblance (MV) et l'algo EM
err_moy0_MV=[]; err_moy1_MV=[]; err_sig0_MV=[]; err_sig1_MV=[]; err_p_anom_MV=[];
err_moy0_EM=[]; err_moy1_EM=[]; err_sig0_EM=[]; err_sig1_EM=[]; err_p_anom_EM=[];
err_moy1_EM2=[];err_moy1_EM3=[];
 
//I. SIMULATION
//I.1. Simulation de la position d'une anomalie

Mi=[1:taille_image]'*ones(1,taille_image);
Mj=Mi';

anomalie_centre=grand(2,1,'unf',1,taille_image);
anomalie_rayon=grand(1,1,'unf',taille_anomalie_min,taille_anomalie_max);

anomalie_carte=zeros(taille_image,taille_image);
anomalie_rayon=anomalie_rayon^2;
anomalie_carte((Mi-anomalie_centre(1)).^2+(Mj-anomalie_centre(2)).^2<anomalie_rayon)=1;

Mi=[]; Mj=[];

//I.2. Creation de l'image bruitee contenant l'anomalie

image0=grand(taille_image,taille_image,'nor',0,1);
image0(anomalie_carte==0)=sqrt(sig0)*image0(anomalie_carte==0);
image0(anomalie_carte==1)=sqrt(sig1)*image0(anomalie_carte==1);

for ng_anomalie=ng_anomalie_max:-1:1, //experiences de detection en reduisant la valeur de ng_anomalie a chq iteration

image=image0;
image(anomalie_carte==1)=ng_anomalie+image0(anomalie_carte==1);

//II. ESTIMATIONS
//II.1 TEST DU RAPPORT DE VRAISEMBLANCE : estimation de la carte de l'anomalie (parametres du modele connus).

//II.1.1 Taille de l'anomalie et probabilité a priori d'une lésion
taille_anomalie=sum(anomalie_carte);
p_anomalie=taille_anomalie/taille_image^2;

//II.1.2 Estimation
anomalie_carte_estimee=test_rapport_vraisemblance(image,sig0,ng_anomalie,sig1,p_anomalie);

//II.1.3 Evaluation de la detection
//pourcentage de points detectes
VD_RV=[VD_RV sum(anomalie_carte_estimee.*anomalie_carte)/taille_anomalie*100];
//nb de fausses alarmes
FD_RV=[FD_RV sum(anomalie_carte_estimee.*(1-anomalie_carte))];
//pourcentage d'erreurs de classement
EC_RV=[EC_RV (sum(anomalie_carte_estimee.*(1-anomalie_carte))+sum((1-anomalie_carte_estimee).*anomalie_carte))/taille_image^2*100];

//II.2. MAXIMUM DE VRAISEMBLANCE : estimation des paramètres du modele (carte d'anomalie connue)
 
//II.2.1 Estimation des parametres
[moy0_MV,sig0_MV,moy1_MV,sig1_MV,p_anom_MV]=max_vraisemblance(image,anomalie_carte);

//II.2.2 Calcul de l'erreur d'estimation
err_moy0_MV=[err_moy0_MV abs(moy0_MV)];
err_moy1_MV=[err_moy1_MV abs(moy1_MV-ng_anomalie)];
err_sig0_MV=[err_sig0_MV abs(sig0_MV-sig0)];
err_sig1_MV=[err_sig1_MV abs(sig1_MV-sig1)];
err_p_anom_MV=[err_p_anom_MV 0];



//II.3. ALGORITHME EM : estimation conjointe de tous les paramètres du modele et de la carte d'anomalie

//II.3.1 Estimation 
[moy0_EM,sig0_EM,moy1_EM,sig1_EM,p_anom_EM,proba_anomalie]=algo_em(image,nb_iter_max);




anomalie_carte_estimee=zeros(proba_anomalie);
anomalie_carte_estimee(proba_anomalie>0.5)=1;


//II.3.2 Evaluation 
//pourcentage de points detectes
VD_EM=[VD_EM sum(anomalie_carte_estimee.*anomalie_carte)/taille_anomalie*100];
//nb de fausses alarmes
FD_EM=[FD_EM sum(anomalie_carte_estimee.*(1-anomalie_carte))];
//pourcentage d'erreurs de classement
EC_EM=[EC_EM (sum(anomalie_carte_estimee.*(1-anomalie_carte))+sum((1-anomalie_carte_estimee).*anomalie_carte))/taille_image^2*100];

//erreur d'estimation des parametres
err_moy0_EM=[err_moy1_EM abs(moy0_EM)];
err_sig0_EM=[err_sig1_EM abs(sig0_EM-sig0)];
err_moy1_EM=[err_moy1_EM abs(moy1_EM-ng_anomalie)];
err_sig1_EM=[err_sig1_EM abs(sig1_EM-sig1)];
err_p_anom_EM=[err_p_anom_EM abs(p_anom_EM-p_anomalie)];


//II.4. ALGORITHME EM version 2 : estimation de la carte d'anomalie et des parametres lies a l'anomalie (ng_anomalie,sig1 et p_anomalie)

//II.4.1 Estimation 


[moy1_EM,sig1_EM,p_anom_EM,proba_anomalie]=algo_em_v2(image,0,sig0,nb_iter_max); 
moy0_EM=0; sig0_EM=sig0;

anomalie_carte_estimee=zeros(proba_anomalie);
anomalie_carte_estimee(proba_anomalie>0.5)=1;


//II.4.2 Evaluation 
//pourcentage de points detectes
VD_EM2=[VD_EM2 sum(anomalie_carte_estimee.*anomalie_carte)/taille_anomalie*100];
//nb de fausses alarmes
FD_EM2=[FD_EM2 sum(anomalie_carte_estimee.*(1-anomalie_carte))];
//pourcentage d'erreurs de classement
EC_EM2=[EC_EM2 (sum(anomalie_carte_estimee.*(1-anomalie_carte))+sum((1-anomalie_carte_estimee).*anomalie_carte))/taille_image^2*100];

//erreur d'estimation du parametre
err_moy1_EM2=[err_moy1_EM2 abs(moy1_EM-ng_anomalie)];



//II.4. ALGORITHME EM version 3 : estimation de la carte d'anomalie et d'un seul parametre (ng_anomalie)

//II.4.1 Estimation 


[moy1_EM,proba_anomalie]=algo_em_v3(image,0,sig0,sig1,p_anomalie,nb_iter_max); 
moy0_EM=0; sig0_EM=sig0; sig1_EM=sig1; p_anom_EM=p_anomalie;


anomalie_carte_estimee=zeros(proba_anomalie);
anomalie_carte_estimee(proba_anomalie>0.5)=1;


//II.4.2 Evaluation 
//pourcentage de points detectes
VD_EM3=[VD_EM3 sum(anomalie_carte_estimee.*anomalie_carte)/taille_anomalie*100];
//nb de fausses alarmes
FD_EM3=[FD_EM3 sum(anomalie_carte_estimee.*(1-anomalie_carte))];
//pourcentage d'erreurs de classement
EC_EM3=[EC_EM3 (sum(anomalie_carte_estimee.*(1-anomalie_carte))+sum((1-anomalie_carte_estimee).*anomalie_carte))/taille_image^2*100];

//erreur d'estimation du parametre
err_moy1_EM3=[err_moy1_EM3 abs(moy1_EM-ng_anomalie)];




//III. Affichage des resultats 

T=ng_anomalie_max:-1:ng_anomalie;

xset("window",0);
xbasc();
xset("colormap",graycolormap(255));
grayplot(1:taille_image,1:taille_image,image,rect=[1,1,taille_image,taille_image]);

xset("window",1);
xbasc();
plot2d(T,[EC_RV',EC_EM',EC_EM2',EC_EM3']);
xtitle("pourcentage d erreurs de classement");
legends(['RV'; 'EM';'EM2'; 'EM3'],[1,2,3,4],1);


err_moy1_MV=err_moy1_MV./T;
err_moy1_EM=err_moy1_EM./T;
err_moy1_EM2=err_moy1_EM2./T;
err_moy1_EM3=err_moy1_EM3./T;
xset("window",2);
xbasc();
plot2d(T,[err_moy1_MV',err_moy1_EM',err_moy1_EM2',err_moy1_EM3']);
xtitle("ordre de l erreur d estimation du parametre ng anomalie");
legends(['MV'; 'EM'; 'EM2'; 'EM3'],[1,2,3,4],1);


end;
