
#include "spline.h"


/******************************************************************************/
/* RACINE CUBIQUE                                                             */
/******************************************************************************/

double cbrt(double x) {
  return (x>=0.?(pow(x,1./3.)):(-pow(-x,1./3.)));
}


/******************************************************************************/
/* FONCTION SPLINE : CALCUL DES COEFFICIENTS DES POLYNOMES SUR CHAQUE MORCEAU */
/*                                                                            */
/* Calcul des coefficients de degre 1, 2 et 3 des polynomes compris entre deux*/
/* points de controle consecutifs. Les points de controles sont stockes dans  */
/* le vecteur y (il s'agit en fait des composantes selon x ou y des points de */
/* controle du plan).                                                         */
/******************************************************************************/

void coeffs_spline(double *y,polynome_t *polynomes,long n) {

  /* variables locales */

  long i;
  double mu;
  double *diag  =(double *)malloc((n-2)*sizeof(double)),
         *pentes=(double *)malloc( n   *sizeof(double));


  /* premiere diagonale de la matrice privee de ses extremites */

  diag[0]=2.;
  for (i=1;i<n-2;i++) diag[i]=4.;


  /* construction du second membre -> initialisation des pentes */

  pentes[0]=(-5.*y[0]+4.*y[1]+y[2])/2.;
  for (i=1;i<n-1;i++) pentes[i]=3.*(y[i+1]-y[i-1]);
  pentes[n-1]=(-y[n-3]-4.*y[n-2]+5.*y[n-1])/2.;


  /* resolution du systeme tridagonal */

  pentes[1]-=pentes[0];
  for (i=1;i<n-2;i++) {
    mu=1./diag[i-1];
    diag[i]-=mu;
    pentes[i+1]-=mu*pentes[i];
  }
  mu=2./diag[n-3];
  pentes[n-1]=(pentes[n-1]-mu*pentes[n-2])/(1.-mu);
  for (i=n-2;i>=1;i--) {
    pentes[i]=(pentes[i]-pentes[i+1])/diag[i-1];
  }
  pentes[0]-=2.*pentes[1];
  pentes[n-1]=pentes[0]=(pentes[n-1]+pentes[0])/2.; /* je l'ai rajoute */


  /* calcul des coefficients des polynomes */

  for (i=0;i<n-1;i++) {
    polynomes[i].c1=pentes[i];
    polynomes[i].c2=3.*(y[i+1]-y[i])-2.*pentes[i]-pentes[i+1];
    polynomes[i].c3=pentes[i]-2.*(y[i+1]-y[i])+pentes[i+1];
  }


  /* liberation de la memoire */

  free(diag); free(pentes);
}


/******************************************************************************/
/* RACINES D'UN POLYNOME DE DEGRE TROIS                                       */
/******************************************************************************/

short racines(double c0,polynome_t polynome,double *r) {

  /* variables locales */

  double a,b,c,d,p,q,delta,racine_carre,theta;
  short retour;


  a=(polynome.c3!=0.?polynome.c3:EPSILON);
  b=polynome.c2/(3.*a); c=polynome.c1/a; d=c0/a;
  p=-3.*CARRE(b)+c;
  q=(2.*CARRE(b)-c)*b+d;
  delta=CUBE(p)*C_4_27+CARRE(q);
  if (delta>EPSILON) {

    /* une solution */

    racine_carre=sqrt(delta);
    r[0]=cbrt(.5*(racine_carre-q))-cbrt(.5*(racine_carre+q))-b;
    retour=1;

  } else if (delta<-EPSILON) {

    /* trois solutions */

    racine_carre=sqrt(C_M4_3*p);
    theta=acos(1.5*q/p*sqrt(-3./p))/3.;
    r[0]=racine_carre*cos(theta        )-b;
    r[1]=racine_carre*cos(theta+C_2PI_3)-b;
    r[2]=racine_carre*cos(theta+C_4PI_3)-b;
    retour=3;

  } else {

    /* deux solutions */

    r[0]=-cbrt(4.*q)-b;
    r[1]= cbrt(.5*q)-b;
    retour=2;

  }

  return retour;
}


/******************************************************************************/
/* RACINES DE LA DERIVEE D'UN POLYNOME DE DEGRE TROIS                         */
/*                                                                            */
/* On part des coefficients de degre 1, 2 et 3 des polynomes compris entre    */
/* deux points de controle consecutifs et on calcule des coefficients de degre*/
/* 0, 1 et 2 de la derivee. Puis on en calcule les racines.                   */
/******************************************************************************/

short racines_derivee(polynome_t polynome,double *r) {

  /* variables locales */

  double a,b,c,delta,rc_delta;
  short retour;


  /* ax^2+bx+c et discriminant */

  a=3.*(polynome.c3!=0.?polynome.c3:EPSILON);
  b=2.*polynome.c2;
  c=polynome.c1;
  delta=CARRE(b)-4.*a*c;
  

  if (delta>EPSILON) {

    /* deux solutions */
    
    rc_delta=sqrt(delta);
    r[0]=(-b-rc_delta)/(2.*a);
    r[1]=(-b+rc_delta)/(2.*a);
    retour=2;

  } else if (delta<-EPSILON) {
  
    /* aucune solution */

    retour=0;
  
  } else {
  
    /* une solution */
    
    r[0]=-b/(2.*a);
    retour=1;
  
  }

  return retour;
}


/******************************************************************************/
/* APPARTENANCE D'UN POINT A L'INTERIEUR D'UNE FORME                          */
/*                                                                            */
/* On compte le nombre d'intersections entre une demi-droite issue de (x,y)   */
/* et la courbe spline. Si ce nombre est impair, (x,y) appartient a la zone   */
/* (retour : 1) sinon il en est en dehors (retour : 0).                       */
/******************************************************************************/

short appartient_forme(polynome_t *polys_x,polynome_t *polys_y,double *points_x,
                       double *points_y,long nb_points,double x,double y) {

  /* variables locales */

  long nb=0,i;
  short nb_racines,j;
  double r[3],y_intersection;
  polynome_t poly_y;


  for (i=0;i<nb_points-1;i++) {
    nb_racines=racines(points_x[i]-x,polys_x[i],r);
    poly_y=polys_y[i];
    for (j=0;j<nb_racines;j++) {
      if (r[j]>=0. && r[j]<=1.) {
        y_intersection=points_y[i]+r[j]*(poly_y.c1+r[j]*(poly_y.c2+r[j]*poly_y.c3));
        if (y<=y_intersection) nb++;
      }
    }
  }
  
  return (nb%2);
}


/******************************************************************************/
/* BORNES DE LA FORME SELON UNE DIRECTION (X OU Y) : RECHERCHE SUR UN SEGMENT */
/******************************************************************************/

void bornes(polynome_t p,double point,double *inf,double *sup) {

  /* variables locales */

  long nb,i;
  double r[2],valeur;


  nb=racines_derivee(p,r);
  for (i=0;i<nb;i++) {
    if (r[i]>=0. && r[i]<=1.) {
      valeur=point+r[i]*(p.c1+r[i]*(p.c2+r[i]*p.c3));
      if (valeur<*inf) *inf=valeur;
      if (valeur>*sup) *sup=valeur;
    }
  }
}

