
#include "tirages.h"
#include "calculQr.h"
#include "zones.h"


/******************************************************************************/
/* SYNTHESE DE CASCADES INFINIMENT DIVISIBLES EN 2D                           */
/*                                                                            */
/* Generation de la texture de depart.                                        */
/*                                                                            */
/* ENTREES :                                                                  */
/*   -> x_max, y_max     : dimensions de la matrice de texture a generer      */
/*   -> r_min            : resolution de la cascade                           */
/*   -> dx               : pas d'espace                                       */
/*   -> type_loi         : identifiant du type de loi                         */
/*   -> parametres_loi   : parametres et notamment densite de la loi          */
/*   -> motif            : identifiant du motif                               */
/*   -> angle            : motifs a angle fixe/variable                       */
/*   -> parametres_motif : points de controles, parametres d'ellipse, ...     */
/* SORTIES :                                                                  */
/*   -> Qr               : matrice de texture                                 */
/*   -> va               : variables aleatoires utilisees pour la synthese    */
/*   -> zones            : zones presentes dans la texture                    */
/*   -> x,y              : vecteurs d'indice de la matrice Qr (facultatif)    */
/******************************************************************************/

void mexFunction(int nb_sorties,      mxArray *sorties[],
                 int nb_entrees,const mxArray *entrees[]) {

  /* variables locales */

  long nx,ny,L,l,N,indice_debut_non_shaper,i;
  char nb_formes=0,nb_ellipses=0,nb_shapers_quelconques=0;
  double x_max,y_max,r_min,u_max,dx,lambda,H1,Ttheo2,*Qr,quotient;
  loi_t loi;
  zone_t zone;
  forme_t forme;
  ellipse_t *ellipse;
  shaper_t shaper;
  positions_t positions;
  tirage_t *tirages;
  en_tete_t *en_tete;


  /* lecture des entrees */

  if (nb_entrees!=8 && nb_entrees!=9)
    mexErrMsgTxt("Il faut huit ou neuf arguments en entree !");
  if (nb_sorties!=3 && nb_sorties!=5)
    mexErrMsgTxt("Il faut trois ou cinq arguments en sortie !");
  if (mxGetM(entrees[0])!=1 || mxGetN(entrees[0])!=1 || (x_max=*mxGetPr(entrees[0]))<=0.)
    mexErrMsgTxt("x_max doit etre un scalaire strictement positif !");
  if (mxGetM(entrees[1])!=1 || mxGetN(entrees[1])!=1 || (y_max=*mxGetPr(entrees[1]))<=0.)
    mexErrMsgTxt("y_max doit etre un scalaire strictement positif !");
  if (mxGetM(entrees[2])!=1 || mxGetN(entrees[2])!=1 || (r_min=*mxGetPr(entrees[2]))<=0.)
    mexErrMsgTxt("r_min doit etre un scalaire strictement positif !");
  if (r_min>=1.) mexErrMsgTxt("r_min doit etre strictement inferieur a 1 !");
  u_max=1./CARRE(r_min)-1.;
  if (mxGetM(entrees[3])!=1 || mxGetN(entrees[3])!=1 || (dx=*mxGetPr(entrees[3]))<=0.)
    mexErrMsgTxt("dx doit etre un scalaire strictement positif !");
  nx=(long)(x_max/dx)+1; ny=(long)(y_max/dx)+1;


  /* lecture et verification des parametres de la loi */

  if (mxGetM(entrees[4])!=1 || mxGetN(entrees[4])!=1) mexErrMsgTxt("type_loi doit etre un scalaire !");
  zone.loi.type=(type_loi_t)*mxGetPr(entrees[4]);
  if (zone.loi.type<1 || zone.loi.type>3) mexErrMsgTxt("Loi inconnue !");
  L=MAX(mxGetM(entrees[5]),mxGetN(entrees[5]));
  l=MIN(mxGetM(entrees[5]),mxGetN(entrees[5]));
  if (l==1 && L==sizeof(parametres_t)/sizeof(double)+1)
    memcpy(&(zone.loi.parametres),mxGetPr(entrees[5]),L*sizeof(double));
  else mexErrMsgTxt("Mauvais nombre de parametres de loi !");
  loi=zone.loi;
  if (loi.C<=0.) mexErrMsgTxt("La parametre de densite de la loi doit etre strictement positif !");
  switch (loi.type) {
    case log_normal:
      if (loi.parametres.ln.sigma2<0.)
        mexErrMsgTxt("Le parametre sigma2 doit etre positif !");
      H1=loi.C*(1.-exp(loi.parametres.ln.mu+loi.parametres.ln.sigma2/2.));
      Ttheo2=-loi.C*(exp(2.*(loi.parametres.ln.sigma2+loi.parametres.ln.mu))-1.);
      break;
    case log_exponential:
      if (loi.parametres.le.Temp<=0.)
        mexErrMsgTxt("Le parametre Temp doit etre strictement positif !");
      H1=0.;
      Ttheo2=-loi.C*(CARRE(1.+loi.parametres.le.Temp)/(1.+2.*loi.parametres.le.Temp)-1.);
      break;
    case log_poisson:
      if (loi.parametres.lp.W0<0. || loi.parametres.lp.W0>=1.)
        mexErrMsgTxt("Le parametre W0 doit etre positif et strictement inferieur a 1 !");
      H1=loi.C*(1.-loi.parametres.lp.W0);
      Ttheo2=loi.C*(-2.*(1.-loi.parametres.lp.W0)+1.-CARRE(loi.parametres.lp.W0));
      break;
  }
  if (Ttheo2<-1.) mexWarnMsgTxt("La convergence n'est pas guarantie !");


  /* lecture des parametres du motif */

  if (mxGetM(entrees[6])!=1 || mxGetN(entrees[6])!=1) mexErrMsgTxt("motif doit etre un scalaire !");
  zone.motif=(motif_t)*mxGetPr(entrees[6]);
  if (zone.motif<1 || zone.motif>8) mexErrMsgTxt("type de motif inconnu !");
  if (mxGetM(entrees[7])!=1 || mxGetN(entrees[7])!=1) mexErrMsgTxt("angle doit etre un scalaire !");
  zone.angle=(motif_t)*mxGetPr(entrees[7]);
  if (zone.angle!=fixe && zone.angle!=variable) mexErrMsgTxt("angle doit etre egal a 0 ou 1 !");
  if (zone.motif==m_polygone || zone.motif==m_spline) {
    if (nb_entrees!=9) mexErrMsgTxt("Il faut neuf arguments en entree !");
    if (mxGetN(entrees[8])!=2 || mxGetM(entrees[8])<3+(zone.motif==m_spline?1:0))
      mexErrMsgTxt("La matrice des points de la forme est invalide !");
    forme.nb=mxGetM(entrees[8]);
    forme.x=mxGetPr(entrees[8]);
    forme.y=mxGetPr(entrees[8])+forme.nb;
    if (zone.motif==m_spline) {
      if (forme.x[forme.nb-1]!=forme.x[0]) forme.x[forme.nb-1]=forme.x[0];
      if (forme.y[forme.nb-1]!=forme.y[0]) forme.y[forme.nb-1]=forme.y[0];
    }
    parametres_formes(zone.motif,zone.angle,&forme);
  } else if (zone.motif==m_ellipse) {
    if (nb_entrees!=9) mexErrMsgTxt("Il faut neuf arguments en entree !");
    ellipse=(ellipse_t *)mxGetPr(entrees[8]);
  } else if (zone.motif==m_shaper_quelconque) {
    if (nb_entrees!=9) mexErrMsgTxt("Il faut neuf arguments en entree !");
    L=MAX(mxGetM(entrees[8]),mxGetN(entrees[8]));
    l=MIN(mxGetM(entrees[8]),mxGetN(entrees[8]));
    if (l!=1 || L<2) mexErrMsgTxt("Les points du shaper doivent etre dans un vecteur de longueur au moins 2 !");
    shaper.nb=L;
    shaper.points=(double *)mxGetPr(entrees[8]);
  }
  switch (zone.motif) {
    case m_carre:
      if (zone.angle==fixe)   zone.numero_fonction=n_carre_fixe;
      else                    zone.numero_fonction=n_carre_variable;    break;
    case m_cercle:            zone.numero_fonction=n_cercle;            break;
    case m_shaper_cos_carre:  zone.numero_fonction=n_shaper_cos_carre;  break;
    case m_shaper_exp:        zone.numero_fonction=n_shaper_exp;        break;
    case m_polygone:          zone.numero_fonction=n_polygone;          break;
    case m_spline:            zone.numero_fonction=n_spline;            break;
    case m_ellipse:           zone.numero_fonction=n_ellipse;           break;
    case m_shaper_quelconque: zone.numero_fonction=n_shaper_quelconque; break;
  }


  /* tirage des positions des multiplieurs */

  rand_init(time(NULL));
  lambda=zone.loi.C*2./T_PI*u_max*(x_max+1.)*(y_max+1.);
  N=(long)(lambda+.5); /* round(lambda); */
  positions.tableau=(position_t *)malloc(N*sizeof(position_t));
  positions.n=0;
  tirer_positions(N,0.,x_max,0.,y_max,u_max,r_min,1.,dx,&positions,nx,ny);
  if (est_shaper(zone.motif)) indice_debut_non_shaper=positions.n;
  else                        indice_debut_non_shaper=0;


  /* allocation memoire pour les sorties */

  allocation_sorties(sorties,nb_sorties,&tirages,positions.n,&Qr,nx,ny,dx);


  /* tirage des multiplieurs (des wi ainsi que des thetai) */

  tirer_multiplieurs(positions,tirages,zone.loi,zone.motif,zone.angle,zone.numero_fonction,0,0);


  /* calcul de Qr */

  calcul_Qr(dx,nx,ny,Qr,&forme,ellipse,&shaper,indice_debut_non_shaper,positions,tirages);


  /* initialisation de la liste des zones, a retourner sous matlab */

  if (zone.motif==m_spline || zone.motif==m_polygone) nb_formes=1;
  else if (zone.motif==m_ellipse)                     nb_ellipses=1;
  else if (zone.motif==m_shaper_quelconque)           nb_shapers_quelconques=1;
  zone.nb_points_controle=0;
  zone.masque=-1;
  zone.mort=0;
  construire_liste_zones(1,nb_formes,nb_ellipses,nb_shapers_quelconques,indice_debut_non_shaper,
                         &zone,&forme,ellipse,&shaper,sorties+2,&en_tete);


  /* liberation de la memoire temporaire */

  liberer_memoire_zones(en_tete,NULL,&forme,NULL,NULL,1);
  free(positions.tableau);


  /* normalisation de Qr */

  quotient=pow(r_min,H1);
  for (i=0;i<nx*ny;i++) Qr[i]/=quotient;
}

