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


void tirage_multiplieurs(positions_t *,positions_t,tirage_t *,zone_t *,long,long *);
void dans_masque(zone_t *,double,long,positions_t *,long);


/******************************************************************************/
/* ZOOM DANS UNE CASCADE INFINIMENT DIVISIBLE EN 2D                           */
/*                                                                            */
/* Zoom dans une texture en diminuant la resolution de la cascade. Permet ega-*/
/* lement de recalculer la matrice Qr en passant une liste de variables alea- */
/* toires modifiee (pour cela on effectue un zoom sur l'integralite de la tex-*/
/* ture en ne modifiant pas la resolution).                                   */
/*                                                                            */
/* ENTREES :                                                                  */
/*   -> cadre          : cadre de la zone de la texture a agrandir            */
/*   -> dx             : pas d'espace                                         */
/*   -> r_min          : resolution de la cascade                             */
/*   -> r_min_prec     : resolution precedente de la cascade                  */
/*   -> va             : variables aleatoires utilisees pour la synthese      */
/*   -> zones          : zones presentes dans la texture                      */
/* 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 L,l,nx,ny,nb_tirages_prec,nb_zones,nb_formes,nb_ellipses,nb_shapers_quelconques,
       nb_points,nb_anciens=0,indice_debut_non_shaper,nb_multiplieurs,deja_multiplieur,i,j,k;
  double x_min,y_min,delta_x,delta_y,dx,r_min,r_min_prec,u_max,lambda,
         inf_x,sup_x,inf_y,sup_y,x,y,*Qr;
  char stop;
  tirage_t *tirages_prec,*tirages;
  en_tete_t *en_tete;
  zone_t *zones=NULL;
  forme_t *formes=NULL;
  ellipse_t *ellipses=NULL;
  shaper_t *shapers=NULL;
  positions_t *positions=NULL,anciennes_positions,toutes_positions;
  position_t *tableau;
  polynome_t *polynomes_x,*polynomes_y;
  double *points_x,*points_y;


  /* lecture des entrees */

  if (nb_entrees!=6) mexErrMsgTxt("Il faut six arguments en entree !");
  if (nb_sorties!=3 && nb_sorties!=5)
    mexErrMsgTxt("Il faut trois ou cinq arguments en sortie !");
  L=MAX(mxGetM(entrees[0]),mxGetN(entrees[0]));
  l=MIN(mxGetM(entrees[0]),mxGetN(entrees[0]));
  if (l==1 && L==4) {
    x_min  =*(mxGetPr(entrees[0])  );
    y_min  =*(mxGetPr(entrees[0])+1);
    delta_x=*(mxGetPr(entrees[0])+2);
    delta_y=*(mxGetPr(entrees[0])+3);
  } else mexErrMsgTxt("Le cadre du zoom doit avoir quatre parametres !");
  if (mxGetM(entrees[1])!=1 || mxGetN(entrees[1])!=1 || (dx=*mxGetPr(entrees[1]))<=0.)
    mexErrMsgTxt("dx doit etre un scalaire strictement positif !");
  nx=(long)(delta_x/dx)+1; ny=(long)(delta_y/dx)+1;
  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 (mxGetM(entrees[3])!=1 || mxGetN(entrees[3])!=1 || (r_min_prec=*mxGetPr(entrees[3]))<=0.)
    mexErrMsgTxt("r_min_prec doit etre un scalaire strictement positif !");
  if (r_min_prec>=1.) mexErrMsgTxt("r_min_prec doit etre strictement inferieur a 1 !");
  if (r_min>r_min_prec) mexErrMsgTxt("r_min doit etre inferieur a r_min_prec !");
  u_max=1./CARRE(r_min)-1./CARRE(r_min_prec);
  if (mxGetM(entrees[4])!=sizeof(tirage_t)/sizeof(double))
    mexErrMsgTxt("Mauvais nombre de lignes pour la matrice des tirages !");
  tirages_prec=(tirage_t *)mxGetPr(entrees[4]);
  nb_tirages_prec=mxGetN(entrees[4]);
  en_tete=(en_tete_t *)mxGetPr(entrees[5]);


  /* lecture des differentes zones */

  nb_zones=lire_liste_zones(en_tete,&zones,&formes,&ellipses,&shapers);
  nb_formes             =en_tete->nb_formes;
  nb_ellipses           =en_tete->nb_ellipses;
  nb_shapers_quelconques=en_tete->nb_shapers_quelconques;
  
  
  /* determination d'un cadre pour la zone : determination du plus petit rectangle */
  /* dans lequel est incrite la zone                                               */

  cadre_zones(zones,&nb_zones,x_min,y_min,delta_x,delta_y,formes,&nb_formes,ellipses,
              &nb_ellipses,shapers,&nb_shapers_quelconques,tirages_prec,nb_tirages_prec);


  /* tirage des positions des nouveaux multiplieurs (entre r_min et */
  /* r_min_prec) pour chaque zone                                   */

  rand_init(time(NULL));
  positions=(positions_t *)malloc(nb_zones*sizeof(positions_t));
  for (i=0;i<nb_zones;i++) {
  
    /* on tire les points dans le cadre (le rectangle) de la zone i */

    lambda=zones[i].loi.C*2./T_PI*u_max*((zones[i].sup_x-zones[i].inf_x)+(i==0?r_min_prec:0))
                                       *((zones[i].sup_y-zones[i].inf_y)+(i==0?r_min_prec:0));
    positions[i].tableau=(position_t *)malloc(((long)(lambda+.5))*sizeof(position_t));
    positions[i].n=0;
    if (!(zones[i].mort))
      tirer_positions((long)(lambda+.5),zones[i].inf_x,zones[i].sup_x,zones[i].inf_y,zones[i].sup_y,
                      u_max,r_min,r_min_prec,dx,positions+i,nx,ny);


    /* on supprime les points en dehors des limites (la spline) de la zone i */

    tableau=positions[i].tableau;
    if (i>0) {
      polynomes_x=zones[i].spline.x; points_x=zones[i].points_x;
      polynomes_y=zones[i].spline.y; points_y=zones[i].points_y;
      nb_points=zones[i].nb_points_controle;
      k=0;
      for (j=0;j<positions[i].n;j++) {
        if (appartient_forme(polynomes_x,polynomes_y,points_x,points_y,nb_points,
                             tableau[j].x,tableau[j].y))
          tableau[k++]=tableau[j];
      }
      positions[i].n=k;
    }


    /* on supprime les points situes en dehors du "masque" de la zone i       */
    /* -> masque : une zone j constitue un masque pour une zone k si tous les */
    /*             multiplieurs de k doivent etre situes dans la zone j       */

    if (i>0 && (j=zones[i].masque)!=-1) dans_masque(zones,r_min_prec,j,positions,i);

    
    /* on supprime les points situes dans une zone j qui couvre la zone i */

    for (j=i+1;j<nb_zones;j++) {
      if (zones[j].masque==-1 || zones[j].masque==i) {
        inf_x=zones[j].inf_x-.5*r_min_prec; sup_x=zones[j].sup_x+.5*r_min_prec;
        inf_y=zones[j].inf_y-.5*r_min_prec; sup_y=zones[j].sup_y+.5*r_min_prec;
        polynomes_x=zones[j].spline.x; points_x=zones[j].points_x;
        polynomes_y=zones[j].spline.y; points_y=zones[j].points_y;
        nb_points=zones[j].nb_points_controle;
        l=0;
        for (k=0;k<positions[i].n;k++) {
          x=tableau[k].x; y=tableau[k].y;
          if (x<inf_x || x>sup_x || y<inf_y || y>sup_y ||
              !appartient_forme(polynomes_x,polynomes_y,points_x,points_y,nb_points,x,y))
            tableau[l++]=tableau[k];
        }
        positions[i].n=l;
      }
    }
  }


  /* determination, parmi les points deja tires avant cet appel à zcid2d, de ceux qui   */
  /* ont une influence sur la matrice Qr que l'on va construire, elimination des autres */

  for (i=0;i<nb_zones;i++) zones[i].nb_elements=0;
  for (i=0;i<nb_tirages_prec;i++) {
    if (tirages_prec[i].x>=x_min-.5 && tirages_prec[i].x<=x_min+delta_x+.5 &&
        tirages_prec[i].y>=y_min-.5 && tirages_prec[i].y<=y_min+delta_y+.5)
      nb_anciens++;
  }
  anciennes_positions.tableau=(position_t *)malloc(nb_anciens*sizeof(position_t));
  anciennes_positions.n=0;
  indice_debut_non_shaper=en_tete->indice_debut_non_shaper;
  j=0;
  for (i=0;i<nb_tirages_prec;i++) {
    if (tirages_prec[i].x>=x_min-.5 && tirages_prec[i].x<=x_min+delta_x+.5 &&
        tirages_prec[i].y>=y_min-.5 && tirages_prec[i].y<=y_min+delta_y+.5 &&
        zone_influence(tirages_prec[i].x-x_min,tirages_prec[i].y-y_min,
                       tirages_prec[i].r,nx,ny,dx,&anciennes_positions)) {
      tirages_prec[j]=tirages_prec[i];


      /* on decale le multiplieur pour compenser la translation des axes du repere */
      /* de la texture, lorsque l'on zoome                                         */

      tirages_prec[j].x-=x_min;
      tirages_prec[j].y-=y_min;


      zones[tirages_prec[j].numero_zone].nb_elements++;
      j++;
    } else if (est_shaper(tirages_prec[i].type_motif)) indice_debut_non_shaper--;
  }


  /* suppression des zones non representees (il n'existe plus de multiplieurs dans cette zone) */

  i=0;
  while (i<nb_zones) {
    if (positions[i].n+zones[i].nb_elements==0) {
      stop=0;
      if (i>0) for (j=i+1;j<nb_zones && !stop;j++) {
        if (zones[j].masque==i) { /* si la zone n'a plus de representants mais que */
          zones[i].mort=1;        /* l'on a besoin de la garder comme masque       */
          stop=1;
        }
      }
      if (!stop) {
        if (zones[i].spline.x) free(zones[i].spline.x);
        if (zones[i].spline.y) free(zones[i].spline.y);
        free(positions[i].tableau);
        for (j=i+1;j<nb_zones;j++) positions[j-1]=positions[j];
        supprimer_zone(zones,&nb_zones,i,formes,&nb_formes,ellipses,&nb_ellipses,shapers,
                       &nb_shapers_quelconques,tirages_prec,anciennes_positions.n,0);
      } else i++;
    } else i++;
  }


  /* allocation memoire pour les sorties */

  nb_multiplieurs=anciennes_positions.n;
  for (i=0;i<nb_zones;i++) {
    nb_multiplieurs+=positions[i].n;
  }
  mexPrintf("nombre de points tires : %d\n",nb_multiplieurs-anciennes_positions.n);
  allocation_sorties(sorties,nb_sorties,&tirages,nb_multiplieurs,&Qr,nx,ny,dx);


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

  toutes_positions.tableau=(position_t *)malloc(nb_multiplieurs*sizeof(position_t));
  toutes_positions.n=nb_multiplieurs;
  deja_multiplieur=0;
  for (i=0;i<nb_zones;i++) {
    if (est_shaper(zones[i].motif))
      tirage_multiplieurs(positions,toutes_positions,tirages,zones,i,&deja_multiplieur);
  }
  indice_debut_non_shaper+=deja_multiplieur;
  memcpy(tirages+deja_multiplieur,tirages_prec,anciennes_positions.n*sizeof(tirage_t));
  memcpy(toutes_positions.tableau+deja_multiplieur,anciennes_positions.tableau,
         anciennes_positions.n*sizeof(position_t));
  deja_multiplieur+=anciennes_positions.n;
  for (i=0;i<nb_zones;i++) {
    if (!est_shaper(zones[i].motif))
      tirage_multiplieurs(positions,toutes_positions,tirages,zones,i,&deja_multiplieur);
  }
  for (i=0;i<nb_zones;i++) free(positions[i].tableau);
  free(positions);
  free(anciennes_positions.tableau);


  /* calcul de Qr */

  calcul_Qr(dx,nx,ny,Qr,formes,ellipses,shapers,indice_debut_non_shaper,toutes_positions,tirages);
  free(toutes_positions.tableau);


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

  construire_liste_zones(nb_zones,nb_formes,nb_ellipses,nb_shapers_quelconques,indice_debut_non_shaper,
                         zones,formes,ellipses,shapers,sorties+2,&en_tete);


  /* liberation de la memoire temporaire */

  liberer_memoire_zones(en_tete,zones,formes,ellipses,shapers,0);
  for (i=0;i<nb_zones;i++) {
    if (zones[i].spline.x) free(zones[i].spline.x);
    if (zones[i].spline.y) free(zones[i].spline.y);
  }
}


/******************************************************************************/
/* TIRAGE DES MULTIPLIEURS ASSOCIES A DE NOUVEAUX POINTS D'UNE ZONE DONNEE    */
/******************************************************************************/

void tirage_multiplieurs(positions_t *positions,positions_t toutes_positions,tirage_t *tirages,
                         zone_t *zones,long i,long *deja_multiplieur) {

  tirer_multiplieurs(positions[i],tirages+(*deja_multiplieur),zones[i].loi,zones[i].motif,
                     zones[i].angle,zones[i].numero_fonction,zones[i].numero_motif,i);
  memcpy(toutes_positions.tableau+(*deja_multiplieur),positions[i].tableau,
         positions[i].n*sizeof(position_t));
  (*deja_multiplieur)+=positions[i].n;
}


/******************************************************************************/
/* ON VERIFIE QUE LES POINTS TIRES DE LA ZONE i SONT DANS LE MASQUE DE i      */
/******************************************************************************/

void dans_masque(zone_t *zones,double r_min_prec,long masque,positions_t *positions,long i) {

  /* variables locales */
  
  double inf_x,sup_x,inf_y,sup_y,*points_x,*points_y,x,y;
  polynome_t *polynomes_x,*polynomes_y;
  long nb_points,j,k;
  position_t *tableau=positions[i].tableau;
  

  /* appel recursif sur les zones masquant le masque courant */

  if ((j=zones[masque].masque)!=-1) dans_masque(zones,r_min_prec,j,positions,i);


  inf_x=zones[masque].inf_x-.5*r_min_prec; sup_x=zones[masque].sup_x+.5*r_min_prec;
  inf_y=zones[masque].inf_y-.5*r_min_prec; sup_y=zones[masque].sup_y+.5*r_min_prec;
  polynomes_x=zones[masque].spline.x; points_x=zones[masque].points_x;
  polynomes_y=zones[masque].spline.y; points_y=zones[masque].points_y;
  nb_points=zones[masque].nb_points_controle;
  k=0;
  for (j=0;j<positions[i].n;j++) {
    x=tableau[j].x; y=tableau[j].y;
    if (x>=inf_x && x<=sup_x && y>=inf_y && y<=sup_y &&
        appartient_forme(polynomes_x,polynomes_y,points_x,points_y,nb_points,x,y))
      tableau[k++]=tableau[j];
  }
  positions[i].n=k;
}

