/* Courbes B-Spline et de Bezier                */
/*                                              */
/* Auteur: Nicolas JANEY                        */
/* nicolas.janey@univ-fcomte.fr                 */
/* Mars 2019                                    */

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

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

#include "Pos3D.h"

/* Variables et constantes globales             */

static const float blanc[] = { 1.0F,1.0F,1.0F,1.0F };
static const float gris[]  = { 0.7F,0.7F,0.7F,1.0F };
static float rx = 0.0F;
static float ry = 0.0F;
static float rz = 0.0F;
static int nb = 500;
static int aff = 0;

static float tCouleur[7][3] = { { 1.0F,1.0F,1.0F },
                                { 1.0F,0.0F,1.0F },
                                { 0.0F,1.0F,1.0F },
                                { 1.0F,1.0F,0.0F },
                                { 1.0F,0.0F,0.0F },
                                { 0.0F,1.0F,0.0F },
                                { 0.0F,0.0F,1.0F } };

static int nbPoints = 37;
static Pos3D *tPos[] = { new Pos3D( 6.5, 4.9,-3.00), new Pos3D( 7.0, 5.0,-3.00), new Pos3D( 7.5, 4.9,-2.00),
                         new Pos3D( 8.0, 4.6, 0.00), new Pos3D( 7.5, 4.4, 2.00), new Pos3D( 7.0, 4.2, 3.50),
                         new Pos3D( 6.0, 3.8, 4.00), new Pos3D( 5.0, 2.8, 4.00), new Pos3D( 4.0, 1.8, 4.00),
                         new Pos3D( 3.0, 0.8, 4.00), new Pos3D( 0.0, 0.0, 4.00), new Pos3D(-1.4, 0.7, 3.75),
                         new Pos3D(-2.0, 2.0, 3.50), new Pos3D(-1.4, 3.3, 3.25), new Pos3D( 0.0, 4.0, 3.00),
                         new Pos3D( 1.4, 3.3, 2.75), new Pos3D( 2.0, 2.0, 2.50), new Pos3D( 1.0, 0.7, 2.25),
                         new Pos3D(-1.0, 0.0, 2.00), new Pos3D(-2.4, 0.7, 1.75), new Pos3D(-3.0, 2.0, 1.50),
                         new Pos3D(-2.4, 3.3, 1.25), new Pos3D(-1.0, 4.0, 1.00), new Pos3D( 0.4, 3.3, 0.75),
                         new Pos3D( 1.0, 2.0, 0.50), new Pos3D( 0.0, 0.7, 0.25), new Pos3D(-2.0, 0.2, 0.00),
                         new Pos3D(-7.0, 0.0, 1.50), new Pos3D(-8.0, 0.0,-1.30), new Pos3D(-7.0, 0.3,-1.90),
                         new Pos3D(-6.0, 1.8,-2.40), new Pos3D(-4.5, 4.9,-2.40), new Pos3D(-3.0, 3.2,-2.40),
                         new Pos3D( 0.5, 2.4,-2.80), new Pos3D( 3.5, 4.4,-2.80), new Pos3D( 4.5, 4.7,-2.90),
                         new Pos3D( 5.5, 4.8,-3.00) };

static double NRUBS[4][4]       = { { -1.0/6.0,  3.0/6.0, -3.0/6.0,  1.0/6.0 },
                                    {  3.0/6.0, -6.0/6.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 } };

static double CATMULL_ROM[4][4] = { { -1.0/2.0,  3.0/2.0, -3.0/2.0,  1.0/2.0 },
                                    {  2.0/2.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,  2.0/2.0,      0.0,      0.0 } };

static void vertex(Pos3D *p,int couleur,double nbPoints) {
  glPushMatrix();
  glColor3f(tCouleur[couleur%7][0],tCouleur[couleur%7][1],tCouleur[couleur%7][2]);
  glTranslated(p->x,p->y,p->z);
  glutSolidSphere(nbPoints,36,18);
  glPopMatrix();
}

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

/* Calcul la position d'un point sur une courbe  */
/* B-Spline controlee par quatre sommets         */
/* g : le tableau des 4 sommets de controle      */
/* t : la valeur de calcul du point              */
/*     t a prendre dans l'intervalle [0.0,1.0]   */
/* mb : la matrice de base                       */
/* p : le point resultat                         */

static void determinationPositionSurBSpline(Pos3D **g,double t,double mb[4][4],Pos3D *p) {
  double vt[4] = { t*t*t,t*t,t,1.0 };
  double vtmb[4] = { 0.0,0.0,0.0,0.0 };
  for ( int j = 0 ; j < 4 ; j++ ) {
    for ( int k = 0 ; k < 4 ; k++ )
      vtmb[j] += vt[k] * mb[k][j] ; }
  p->x = p->y = p->z = 0.0;
  for ( int j = 0 ; j < 4 ; j++ ) {
    p->x += vtmb[j] * g[j]->x ;
    p->y += vtmb[j] * g[j]->y ;
    p->z += vtmb[j] * g[j]->z ; }
}

/* Calcul les points definissant une courbe      */
/* B-Spline par morceaux definie par un ensemble */
/* de sommets de controle                        */
/* tPos : le tableau des sommets de controle     */
/* n : le nombre de sommets de contrôle          */
/* mb : la matrice de base                       */
/* nb : le nombre de points a calculer           */
/* tRes : le tableau de points resultat          */

static void calculBSpline(Pos3D **tPos,int n,double mb[4][4],int nb,Pos3D **tRes) {
  for ( int i = 0 ; i < nb ; i++ ) {
    double pos = i/(nb-1.0)*(n-3);
    int nb =(int) pos;
    if ( nb == n-3 )
      nb = n-4;
    double t = pos-nb;
    determinationPositionSurBSpline(&tPos[nb],t,mb,tRes[i]); }
}

/* Calcul de v a la puissance p                  */

float power(double v,int p) {
  return((float) (( v || p ) ? pow(v,(double) p) : 1.0));
}

/* Calcul les points definissant une courbe      */
/* de Bezier definie par un ensemble             */
/* de sommets de controle                        */
/* tPos : le tableau des sommets de controle     */
/* n : le nombre de sommets de contrôle          */
/* nb : le nombre de points a calculer           */
/* tRes : le tableau de points resultat          */

static void calculBezier(Pos3D **tPos,int n,int nb,Pos3D **tRes) {
  float t,mt;
  long long *cn;
  float x,y,z,fac;
  cn =(long long *) calloc(n,sizeof(long long));
  cn[0] = 1;
  cn[1] = n-1;
  for ( int i = 2 ; i < n ; i++ )
    cn[i] = (cn[i-1]*(n-i))/i;
  for ( int i = 0 ; i < nb ; i++ ) {
    t =(float) i/(nb-1);
    mt = 1.0F-t;
    x = y = z = 0.0F;
    for ( int j = 0 ; j < n ; j++ ) {
      fac = cn[j]*power(t,j)*power(mt,n-1-j);
      x += fac * tPos[j]->x;
      y += fac * tPos[j]->y;
      z += fac * tPos[j]->z; }
    tRes[i]->x = x;
    tRes[i]->y = y;
    tRes[i]->z = z; }
  free(cn);
}

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

static void determinationPositionSurBSpline(Pos3D **tPos,int n,double pos,double mb[4][4],Pos3D *p) {
  int nb =(int) pos;
  if ( nb == n-3 )
    nb = n-4;
  double t = pos-nb;
  determinationPositionSurBSpline(&tPos[nb],t,mb,p);
}

static Pos3D **allocationPos3D(int n) {
  Pos3D **tP =(Pos3D **) calloc(n,sizeof(Pos3D *));
  for ( int i = 0 ; i < n ; i++ )
    tP[i] = new Pos3D();
  return tP;
}

static void desallocationPos3D(Pos3D **tP,int n) {
  for ( int i = 0 ; i < n ; i++ )
    delete(tP[i]);
  free(tP);
}

static void modelisationBSpline(Pos3D **tPos,int n,int nb,double mb[4][4]) {
  Pos3D **tP = allocationPos3D(nb);
  calculBSpline(tPos,n,mb,nb,tP);
  for ( int i = 0 ; i < nb ; i++ ) {
    double pos = i/(nb-1.0)*(n-3);
    vertex(tP[i],( pos < n-3 ) ? (int) pos : n-4,0.05); }
  desallocationPos3D(tP,nb);
}

static void modelisationBezier(Pos3D **tPos,int n,int nb) {
  Pos3D **tP = allocationPos3D(nb);
  calculBezier(tPos,n,nb,tP);
  for ( int i = 0 ; i < nb ; i++ ) {
    vertex(tP[i],0,0.05); }
  desallocationPos3D(tP,nb);
}

static void modelisationBSplineUnie(Pos3D **tPos,int n,int nb,double mb[4][4],int couleur,double taille) {
  Pos3D **tP = allocationPos3D(nb);
  calculBSpline(tPos,n,mb,nb,tP);
  for ( int i = 0 ; i < nb ; i++ ) {
    double pos = i/(nb-1.0)*(n-3);
    vertex(tP[i],couleur,taille); }
  desallocationPos3D(tP,nb);
}

static void modelisationBezierUnie(Pos3D **tPos,int n,int nb,int couleur,double taille) {
  Pos3D **tP = allocationPos3D(nb);
  calculBezier(tPos,n,nb,tP);
  for ( int i = 0 ; i < nb ; i++ ) {
    vertex(tP[i],couleur,taille); }
  desallocationPos3D(tP,nb);
}

/* Fonction executee lors de l'appui            */
/* d'une touche de curseur ou d'une touche      */
/* page up ou page down                         */

static void special(int key,int x,int y) {
  switch (key) {
    case GLUT_KEY_UP :
      rx++;
      glutPostRedisplay();
      break;
    case GLUT_KEY_DOWN :
      rx--;
      glutPostRedisplay();
      break;
    case GLUT_KEY_LEFT :
      ry++;
      glutPostRedisplay();
      break;
    case GLUT_KEY_RIGHT :
      ry--;
      glutPostRedisplay();
      break;
    case GLUT_KEY_PAGE_UP :
      switch (glutGetModifiers()) {
        case GLUT_ACTIVE_ALT :
          nb++;
          glutPostRedisplay();
          break;
        default :
          rz++;
          glutPostRedisplay(); }
      break;
    case GLUT_KEY_PAGE_DOWN :
      switch (glutGetModifiers()) {
        case GLUT_ACTIVE_ALT :
          nb--;
          if ( nb == 1 )
            nb = 2;
          glutPostRedisplay(); 
          break;
        default :
          rz--;
          glutPostRedisplay(); }
      break; }
}

/* Fonction d'initialisation des parametres     */
/* OpenGL ne changeant pas au cours de la vie   */
/* du programme                                 */

static void init(void) {
  glDepthFunc(GL_LESS);
  glEnable(GL_DEPTH_TEST);
  glEnable(GL_NORMALIZE);
}

/* Scene dessinee                               */

static void dessinerOpenGL(Pos3D **tP,int n) {
  glBegin(GL_LINE_STRIP);
  for ( int i = 0 ; i < n ; i++ ) {
    glVertex3d(tP[i]->x,tP[i]->y,tP[i]->z); }
  glEnd();
}

static void scene(void) {
  switch (aff) {
    case 0 :
      glColor3f(0.75F,0.75F,0.75F);
      dessinerOpenGL(tPos,nbPoints);
      break;
    case 1 : {
      glColor3f(0.75F,0.75F,0.75F);
      dessinerOpenGL(tPos,nbPoints);
      modelisationBSpline(tPos,nbPoints,nb,NRUBS); }
      break;
    case 2 : {
      glColor3f(0.75F,0.75F,0.75F);
      dessinerOpenGL(tPos,nbPoints);
      modelisationBSpline(tPos,nbPoints,nb,CATMULL_ROM); }
      break;
    case 3 : {
      glColor3f(0.75F,0.75F,0.75F);
      dessinerOpenGL(tPos,nbPoints);
      modelisationBezier(tPos,nbPoints,nb); }
      break;
    case 4 : {
      glColor3fv(tCouleur[0]);
      glLineWidth(2.0F);
      dessinerOpenGL(tPos,nbPoints);
      glLineWidth(1.0F);
      modelisationBSplineUnie(tPos,nbPoints,nb,NRUBS,1,0.03);
      modelisationBSplineUnie(tPos,nbPoints,nb,CATMULL_ROM,2,0.03);
      modelisationBezierUnie(tPos,nbPoints,nb,3,0.03); }
      break; }
}

/* Fonction executee lors d'un rafraichissement */
/* de la fenetre de dessin                      */

static void display(void) {
  glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
  glPushMatrix();
  glRotatef(rx,1.0F,0.0F,0.0F);
  glRotatef(ry,0.0F,1.0F,0.0F);
  glRotatef(rz,0.0F,0.0F,1.0F);
  glTranslatef(0.0F,-2.5F,0.0F);
  scene();
  glPopMatrix();
  glFlush();
  glutSwapBuffers();
  int error = glGetError();
  if ( error != GL_NO_ERROR )
    printf("Erreur OpenGL: %d\n",error);
}

/* Fonction executee lors d'un changement       */
/* de la nbPoints de la fenetre OpenGL            */
/* -> Ajustement de la camera de visualisation  */

static void reshape(int tx,int ty) {
  glViewport(0,0,tx,ty); 
  glMatrixMode(GL_PROJECTION);
  glLoadIdentity();
  double ratio =(double) tx/ty;
  if ( ratio >= 1.0 )
    glOrtho(-3.5*ratio,3.5*ratio,-3.5,3.5,-10.0,10.0);
    else
    glOrtho(-3.5,3.5,-3.5/ratio,3.5/ratio,-10.0,10.0);
  glMatrixMode(GL_MODELVIEW);
  glLoadIdentity();
}

/* Fonction executee lors de la frappe          */
/* d'une touche du clavier                      */

static void keyboard(unsigned char key,int x,int y) {
  switch (key) {
    case 0x20 :
      aff = (aff+1)%5;
      glutPostRedisplay();
      break;
    case 0x1B :
      exit(0); 
      break; }
}

static void clean(void) {
  for ( int i = 0 ; i < nbPoints ; i++ ) {
    delete(tPos[i]); }
}

/* Fonction principale                          */

int main(int argc,char **argv) {

  atexit(clean);

  glutInit(&argc,argv);
  glutInitDisplayMode(GLUT_RGBA|GLUT_DEPTH|GLUT_DOUBLE);
  glutInitWindowSize(600,250); 
  glutInitWindowPosition(50,50); 
  glutCreateWindow("Courbes B-Spline et de Bézier"); 
  init();
  glutKeyboardFunc(keyboard);
  glutSpecialFunc(special);
  glutReshapeFunc(reshape);
  glutKeyboardFunc(keyboard);
  glutDisplayFunc(display);
  glutMainLoop();
  return(0);
}