/* Mathematiques de l'informatique graphique    */
/* Ligne polygonale en 3D                       */
/*                                              */
/* Auteur: Nicolas JANEY                        */
/* nicolas.janey@univ-fcomte.fr                 */
/* Novembre 2012                                */

#include <stdlib.h>
#include <stdio.h>
#include <math.h>

#include <GL/glut.h>
#include <GL/gl.h>
#include <GL/glu.h>

#include "LignePolygonale3D.h"
#include "Position3D.h"

/* Constructeurs                                */

LignePolygonale3D::LignePolygonale3D(void) {
  t = NULL;
  l = 0;
  n = 0;
}

LignePolygonale3D::LignePolygonale3D(int l) {
  init(l);
}

LignePolygonale3D::LignePolygonale3D(LignePolygonale3D *lp) {
  t = ( lp->n > 0 ) ? (Position3D **) calloc(lp->n,sizeof(Position3D *)) : NULL;
  n = 0;
  l = lp->n;
  for ( int i = 0 ; i < n ; i++ )
    add(lp->t[i]);
}

/* Destructeur                                  */

LignePolygonale3D::~LignePolygonale3D(void) {
  for ( int i = 0 ; i < l ; i++ ) {
    if ( t[i] ) {
      delete(t[i]); } }
  delete(t);
}

/* Allocation pour np Position3D disponibles    */

void LignePolygonale3D::init(int np) {
  t = ( np > 0 ) ? (Position3D **) calloc(np,sizeof(Position3D *)) : NULL;
  l = np;
  n = np;
  for ( int i = 0 ; i < np ; i++ )
    t[i] = new Position3D();
}

/* Ajout d'une Position3D                       */

void LignePolygonale3D::add(Position3D *p) {
  if ( l == 0 ) {
    t =(Position3D **) calloc(1,sizeof(Position3D *));
    l = 1; }
    else
    if ( n == l ) {
      t =(Position3D **) realloc(t,(l+1)*sizeof(Position3D *));
      l++; }
  t[n] = new Position3D(p);
  n++;
}

/* Dessin OpenGL                                */

void LignePolygonale3D::draw(void) {
  glBegin(GL_LINE_STRIP);
  for ( int i = 0 ; i < n ; i++ )
    glVertex3dv(t[i]->c);
  glEnd();
}

static void calculPositionSurBSpline(Position3D *p,Position3D **g,double t,const double m[4][4]) {
  double tt[4],ttt[4] ;
  tt[0] = t*t*t ;
  tt[1] = t*t ;
  tt[2] = t ;
  tt[3] = 1.0F ;
  for ( int j = 0 ; j < 4 ; j++ ) {
    ttt[j] = 0 ;
    for ( int k = 0 ; k < 4 ; k++ )
      ttt[j] += tt[k] * m[k][j] ; }
  p->c[0] = p->c[1] = p->c[2] = 0.0;
  for ( int j = 0 ; j < 4 ; j++ ) {
    p->c[0] += ttt[j] * g[j]->c[0] ;
    p->c[1] += ttt[j] * g[j]->c[1] ;
    p->c[2] += ttt[j] * g[j]->c[2] ; }
  p->c[3] = 1.0;
}

/* Construction d'une courbe de B-Spline        */
/* generique                                    */

static void construitBSpline(const double m[4][4],
                             Position3D **tp,int ntp,
                             Position3D **tcb,int ntcb) {
  for ( int i = 0 ; i < ntcb ; i++ ) {
    double t =(double) (ntp-3)*i/(ntcb-1);
    int crv = t;
    t -= crv;
    if ( crv == ntp-3 ) {
      crv = ntp-4;
      t = 1.0; }
    calculPositionSurBSpline(tcb[i],&tp[crv],t,m); }
}

/* Construction d'une courbe de B-Spline NRU    */

static void construitBSplineNRU(Position3D **tp,int ntp,Position3D **tcb,int ntcb) {
  static const double m[4][4] = { { -1.0/6.0, 3.0/6.0, -3.0/6.0, 1.0/6.0 },
                                  {  3.0/6.0,    -1.0,  3.0/6.0,     0.0 },
                                  { -3.0/6.0,     0.0,  3.0/6.0,     0.0 },
                                  {  1.0/6.0, 4.0/6.0,  1.0/6.0,     0.0 } };
  construitBSpline(m,tp,ntp,tcb,ntcb);
}

/* Construction d'une courbe de B-Spline        */
/* de Catmull-Rom                               */

static void construitBSplineCatmullRom(Position3D **tp,int ntp,Position3D **tcb,int ntcb) {
  static const double m[4][4] = { { -1.0/2.0, 3.0/2.0, -3.0/2.0, 1.0/2.0 },
                                  {      1.0,-5.0/2.0,  4.0/2.0,-1.0/2.0 },
                                  { -1.0/2.0,     0.0,  1.0/2.0,     0.0 },
                                  {      0.0,     1.0,      0.0,     0.0 } };
  construitBSpline(m,tp,ntp,tcb,ntcb);
}

//////////////////////////////////////////////////

/* Fonction factorielle                         */

static double factoriel(int n) {
  double fact = 1.0;
  for ( int i = 2 ; i <= n ; i++ ) {
    fact *= i; }
  return(fact);
}

/* Fonction de calcul des Cin                   */

static double *calculCinVersion1(int n) {
  double *cn =(double *) calloc(n,sizeof(double));
  cn[0] = cn[n-1] = 1.0;
  cn[1] = cn[n-2] = (double) (n-1);
  double factn = factoriel(n-1);
  for ( int i = 2 ; i < n-2 ; i++ ) {
    cn[i] = factn/(factoriel(i)*factoriel(n-i-1)); }
  return(cn);
}

/* Fonction de calcul des Cin                   */

static double *calculCinVersion2(int n) {
  double *cn =(double *) calloc(n,sizeof(double));
  cn[0] = cn[n-1] = 1.0;
  cn[1] = cn[n-2] = (double) (n-1);
  for ( int i = 2 ; i < n-2 ; i++ ) {
    cn[i] = (cn[i-1]*(n-i))/i; }
  return(cn);
}

/* Fonction power                               */

static double power(double v,int p) {
  return(( !v && !p ) ? 1.: pow(v,(double) p));
}

/* Construction d'une courbe de Bezier          */

static void construitBezier(Position3D **tp,int ntp,Position3D **tcb,int ntcb) {
  double *cn = calculCinVersion1(ntp);
  for ( int i = 0 ; i < ntcb ; i++ ) {
    double t =(double) i/(ntcb-1);
    double mt = 1.0-t;
    double x = 0.0,y = 0.0,z = 0.0;
    for ( int j = 0 ; j < ntp ; j++ ) {
      double fac = cn[j]*power(t,j)*power(mt,ntp-1-j);
      x += (fac*tp[j]->c[0]);
      y += (fac*tp[j]->c[1]);
      z += (fac*tp[j]->c[2]); }
    tcb[i]->c[0] = x;
    tcb[i]->c[1] = y;
    tcb[i]->c[2] = z; }
  free(cn);
}

/* Constructeur de creation de courbes lissees  */

LignePolygonale3D::LignePolygonale3D(LignePolygonale3D *lp,int np) {
  init(np);
  construitBezier(lp->t,lp->n,t,np);
}

//////////////////////////////////////////////////

/* Constructeur de creation de courbes lissees  */

LignePolygonale3D::LignePolygonale3D(LignePolygonale3D *lp,int np,type_lissage type) {
  init(np);
  switch(type) {
    case BEZIER :
      construitBezier(lp->t,lp->n,t,np);
      break;
    case BSPLINE_NRU :
      construitBSplineNRU(lp->t,lp->n,t,np);
      break;
    case BSPLINE_CATMULL_ROM :
      construitBSplineCatmullRom(lp->t,lp->n,t,np);
      break; }
}