/*----------------------------MegaWave2 module----------------------------*/
/* mwcommand
name = {classif_gk};
version = {"1.0"};
author = {"Lionel Moisan"};
function = {"Calcule la carte de classification (noyau gaussien)"};
usage = {
 's':[s=1.]->s        "standart deviation (default 1.)",
 'n':[n=500]->n       "taille de l'image (carrée), par défaut 500",
 e1->e1               "échantillons de la classe 1 (2-Flist)",
 e2->e2               "échantillons de la classe 2 (2-Flist)",
 out<-classif_gk      "résultat: carte de classification (Cimage)"
        };
*/

/* convention de remplissage de la carte de classification:
     80: classe 1 prédite
    170: classe 2 prédite 
*/
 
#include <stdio.h>
#include <math.h>
#include "mw.h"

Cimage classif_gk(e1,e2,n,s)
     Flist e1,e2;
     int *n;
     float *s;
{
  Cimage out;
  int x,y,i;
  float xx,yy,*p;
  double dd,var,p1,p2;

  out = mw_change_cimage(NULL,*n,*n);
  var = (double)(*s*(*s));

  /* parcours des pixels de la carte de classification (out) */
  for (x=0;x<*n;x++)
    for (y=0;y<*n;y++) {
      xx = (float)x/(float)(*n);
      yy = (float)y/(float)(*n);
      /* détermination de la proba pour la classe 1 */
      for (i=0,p1=0.,p=e1->values;i<e1->size;i++,p+=2) {
	dd = (double)( (*p-xx)*(*p-xx)+(*(p+1)-yy)*(*(p+1)-yy) );
	p1 += exp(-0.5*dd/var);
      }
      /* détermination de la proba pour la classe 2 */
      for (i=0,p2=0.,p=e2->values;i<e2->size;i++,p+=2) {
	dd = (double)( (*p-xx)*(*p-xx)+(*(p+1)-yy)*(*(p+1)-yy) );
	p2 += exp(-0.5*dd/var);
      }
      /* remplissage de la carte de classification */
      out->gray[*n*y+x] = ((p1>p2)?80:170);
    }

  return(out);
}