
#include "spline.h"
#include "calculQr.h"


/******************************************************************************/
/* SELECTION D'UNE PARTIE DE LA TEXTURE                                       */
/*                                                                            */
/* Selection des points de la matrice Qr inclus dans une zone delimitee par   */
/* une courbe spline fermee. On determine le plus petit rectangle dans lequel */
/* s'inscrit la zone et la matrice retournee a la taille de ce rectangle.     */
/*                                                                            */
/* ENTREES :                                                                  */
/*   -> points : points de controle de la spline                              */
/*   -> Qr     : matrice de texture                                           */
/* SORTIES :                                                                  */
/*   -> Vr     : matrice de texture extraite :                                */
/*               Vr(i-rec_i,j-rec_j)=|Qr(i,j) si (i,j) est dans la zone selec-*/
/*                                   |   tionnee ou (rec_i,rec_j) est le coin */
/*                                   |   superieur gauche du rectangle        */
/*                                   |0 sinon                                 */
/******************************************************************************/

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

  /* variables locales */

  long nb_points,m,n,deb_i,fin_i,deb_j,fin_j,m2,n2,i,j;
  short appartient;
  double *sx,*sy,*Qr,sup_x,sup_y,inf_x,inf_y,*Vr;
  polynome_t *polynomes_x,*polynomes_y;


  /* lecture des entrees */

  if (nb_entrees!=2) mexErrMsgTxt("Il faut deux entrees !");
  if (nb_sorties>1) mexErrMsgTxt("Il faut au maximum une sortie !");
  if (mxGetN(entrees[0])!=2 || mxGetM(entrees[0])<3)
    mexErrMsgTxt("La matrice des points de la forme est invalide !");
  nb_points=mxGetM(entrees[0]);
  sx=mxGetPr(entrees[0]);
  sy=mxGetPr(entrees[0])+nb_points;
  polynomes_x=(polynome_t *)malloc((nb_points-1)*sizeof(polynome_t));
  polynomes_y=(polynome_t *)malloc((nb_points-1)*sizeof(polynome_t));
  coeffs_spline(sx,polynomes_x,nb_points);
  coeffs_spline(sy,polynomes_y,nb_points);
  Qr=mxGetPr(entrees[1]);
  m=mxGetM(entrees[1]);
  n=mxGetN(entrees[1]);


  /* cadre de la selection */

  sup_x=sup_y=0.;
  inf_x=(double)m;
  inf_y=(double)n;
  for (i=0;i<nb_points-1;i++) {
    bornes(polynomes_x[i],sx[i],&inf_x,&sup_x);
    bornes(polynomes_y[i],sy[i],&inf_y,&sup_y);
  }
  if (inf_x<0.) inf_x=0.; if (inf_y<0.) inf_y=0.;
  if (sup_x>(double)m) sup_x=(double)m;
  if (sup_y>(double)n) sup_y=(double)n;
  deb_i=(long)inf_x; fin_i=MIN((long)sup_x+1,m);
  deb_j=(long)inf_y; fin_j=MIN((long)sup_y+1,n);
  m2=fin_i-deb_i; n2=fin_j-deb_j;


  /* allocation memoire pour les sorties */

  sorties[0]=mxCreateDoubleMatrix(m2,n2,mxREAL);
  Vr=mxGetPr(sorties[0]);


  for (i=deb_i;i<fin_i;i++) {
    for (j=deb_j;j<fin_j;j++) {

      /* appartenance de (i,j) a l'interieur de la courbe */

      appartient=appartient_forme(polynomes_x,polynomes_y,sx,sy,nb_points,(double)i,(double)j);


      /* si le point est dans la spline -> Vr(i-deb_i,j-deb_j):=Qr(i,j) */

      if (appartient) Vr[(j-deb_j)*m2+(i-deb_i)]=Qr[j*m+i];
      else            Vr[(j-deb_j)*m2+(i-deb_i)]=0.;
    }
  }
  

  /* liberation de la memoire */

  free(polynomes_x);
  free(polynomes_y);
}

