/* Surface B-Spline                             */
/*                                              */
/* Auteur: Nicolas JANEY                        */
/* nicolas.janey@univ-fcomte.fr                 */
/* Mars 2014                                    */

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

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

#include "Pos3D.h"
#include "Dir3D.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 = 100;
static int nP = 0;
static int aff = 0;

static Pos3D *tPos[4][4] = { { new Pos3D(-6.0, 4.0, 1.0),new Pos3D(-2.5, 5.0, 2.0),new Pos3D( 2.5, 3.5, 3.0),new Pos3D( 5.5, 4.5, 4.0) },
                             { new Pos3D(-4.0, 1.0,-1.0),new Pos3D(-1.5, 2.0,-2.0),new Pos3D( 2.0, 1.5,-1.0),new Pos3D( 5.0, 2.0, 1.0) },
                             { new Pos3D(-5.0,-1.0,-1.0),new Pos3D(-2.0,-1.5,-1.0),new Pos3D( 1.0,-2.0, 1.0),new Pos3D( 4.5,-3.0, 2.0) },
                             { new Pos3D(-5.5,-4.0, 1.0),new Pos3D(-3.0,-3.0,-1.5),new Pos3D( 0.5,-5.0,-3.5),new Pos3D( 3.5,-5.5,-1.0) } };

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 BEZIER[4][4] = { { -1.0,  3.0, -3.0,  1.0 },
                               {  3.0, -6.0,  3.0,  0.0 },
                               { -3.0,  3.0,  0.0,  0.0 },
                               {  1.0,  0.0,  0.0,  0.0 } };

static void determinationTangenteSurBSpline(Pos3D **g,double t,double (*mb)[4],Dir3D *p) {
  double vt[4] = { 3*t*t,2*t,1.0,0.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 ; }
}

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

static void calculTangentesBSpline(Pos3D **tPos,int n,int nb,double (*mb)[4],Dir3D **tRes) {
  for ( int i = 0 ; i < nb ; i++ ) {
    double pos = i/(nb-1.0)*(n-3);
    determinationTangenteSurBSpline(tPos,n,pos,mb,tRes[i]); }
}

static void determinationPositionSurBSpline(Pos3D **g,double t,double (*mb)[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 ; }
}

static void determinationPositionSurBSpline(Pos3D **tPos,int n,double pos,double (*mb)[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 void calculBSpline(Pos3D **tPos,int n,int nb,double (*mb)[4],Pos3D **tRes) {
  for ( int i = 0 ; i < nb ; i++ ) {
    double pos = i/(nb-1.0)*(n-3);
    determinationPositionSurBSpline(tPos,n,pos,mb,tRes[i]); }
}

/* 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;
        case GLUT_ACTIVE_CTRL :
          nP++;
          if ( nP == nb )
            nP = nb-1;
          glutPostRedisplay();
          break;
        default :
          rz++;
          glutPostRedisplay(); }
      break;
    case GLUT_KEY_PAGE_DOWN :
      switch (glutGetModifiers()) {
        case GLUT_ACTIVE_ALT :
          nb--;
          if ( nb == 1 )
            nb = 2;
          if ( nP == nb )
            nP = nb-1;
          glutPostRedisplay(); 
          break;
        case GLUT_ACTIVE_CTRL :
          nP--;
          if ( nP < 0 )
            nP = 0;
          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);
}

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

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

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

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

static Dir3D ***allocationDir3D(int n1,int n2) {
  Dir3D ***tT =(Dir3D ***) calloc(n1,sizeof(Dir3D **));
  for ( int j = 0 ; j < n1 ; j++ ) {
    tT[j] =(Dir3D **) calloc(n2,sizeof(Dir3D *));
    for ( int i = 0 ; i < n2 ; i++ )
      tT[j][i] = new Dir3D(); }
  return tT;
}

static void desallocationDir3D(Dir3D ***tT,int n1,int n2) {
  for ( int j = 0 ; j < n1 ; j++ ) {
    for ( int i = 0 ; i < n2 ; i++ )
      delete(tT[j][i]);
    free(tT[j]); }
  free(tT);
}

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);
}

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

static void dessiner4CourbesBSplineOpenGL(Pos3D *tP[4][4],int n) {
  Pos3D ***tPC = allocationPos3D(4,n);
  for ( int i = 0 ; i < 4 ; i++ ) {
    calculBSpline(tPos[i],4,n,BEZIER,tPC[i]); }
  for ( int i = 0 ; i < 4 ; i++ ) {
    affichage(tPC[i],n); }
  desallocationPos3D(tPC,4,n);
}

void affichage(Pos3D ***tP,int n1,int n2) {
  glBegin(GL_POINTS);
  for ( int i = 0 ; i < n1 ; i++ ) {
    for ( int j = 0 ; j < n2 ; j++ ) {
      glVertex3d(tP[i][j]->x,tP[i][j]->y,tP[i][j]->z); } }
  glEnd();
}

static Pos3D ***calculNCourbesBSpline(Pos3D *tPos[4][4],int n) {
  Pos3D ***tPC = allocationPos3D(4,n);
  for ( int i = 0 ; i < 4 ; i++ ) {
    calculBSpline(tPos[i],4,n,BEZIER,tPC[i]); }
  Pos3D **tPCtrl = allocationPos3D(4);
  Pos3D ***tPS = allocationPos3D(n,n);
  for ( int i = 0 ; i < n ; i++ ) {
    for ( int j = 0 ; j < 4 ; j++) {
      tPCtrl[j]->x = tPC[j][i]->x;
      tPCtrl[j]->y = tPC[j][i]->y;
      tPCtrl[j]->z = tPC[j][i]->z;
      tPCtrl[j]->w = tPC[j][i]->w; }
    calculBSpline(tPCtrl,4,n,BEZIER,tPS[i]); }
  desallocationPos3D(tPCtrl,4);
  desallocationPos3D(tPC,4,n);
  return tPS;
}

static void dessinerNCourbesBSplineOpenGL(Pos3D *tP[4][4],int n) {
  Pos3D ***tPS = calculNCourbesBSpline(tP,n);
  affichage(tPS,n,n);
  desallocationPos3D(tPS,n,n);
}

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

static void dessinerMaillageBSplineOpenGL(Pos3D *tP[4][4],int n) {
  Pos3D ***tPS = calculNCourbesBSpline(tP,n);
  dessinerMaillageOpenGL(tPS,n,n);
  desallocationPos3D(tPS,n,n);
}

static void dessinerSurfaceBSplineOpenGL(Pos3D ***tP,Dir3D ***tT,int n1,int n2) {
  for ( int i = 0 ; i < n1 ; i++ ) {
    glBegin(GL_LINE_STRIP);
    for ( int j = 0 ; j < n2 ; j++ ) {
      glNormal3d(tT[i][j]->x,tT[i][j]->y,tT[i][j]->z);
      glVertex3d(tP[i][j]->x,tP[i][j]->y,tP[i][j]->z); }
    glEnd(); }
  for ( int i = 0 ; i < n2 ; i++ ) {
    glBegin(GL_LINE_STRIP);
    for ( int j = 0 ; j < n1 ; j++ ) {
      glNormal3d(tT[j][i]->x,tT[j][i]->y,tT[j][i]->z);
      glVertex3d(tP[j][i]->x,tP[j][i]->y,tP[j][i]->z); }
    glEnd(); }
}

static void calculSurfaceBSpline(Pos3D *tPos[4][4],int n,Pos3D ***tP,Dir3D ***tT) {
  Pos3D ***tPC = allocationPos3D(4,n);
  for ( int i = 0 ; i < 4 ; i++ ) {
    calculBSpline(tPos[i],4,n,BEZIER,tPC[i]); }
  Pos3D **tPCtrl = allocationPos3D(4);
  for ( int i = 0 ; i < n ; i++ ) {
    for ( int j = 0 ; j < 4 ; j++) {
      tPCtrl[j]->x = tPC[j][i]->x;
      tPCtrl[j]->y = tPC[j][i]->y;
      tPCtrl[j]->z = tPC[j][i]->z;
      tPCtrl[j]->w = tPC[j][i]->w; }
    calculBSpline(tPCtrl,4,n,BEZIER,tP[i]); }
  desallocationPos3D(tPCtrl,4);
  desallocationPos3D(tPC,4,n);
}

static void dessinerSurfaceBSplineOpenGL(Pos3D *tP[4][4],int n) {
  Pos3D ***tPS = allocationPos3D(n,n);
  Dir3D ***tTS = allocationDir3D(n,n);
  calculSurfaceBSpline(tP,n,tPS,tTS);
  dessinerSurfaceBSplineOpenGL(tPS,tTS,n,n);
  desallocationPos3D(tPS,n,n);
  desallocationDir3D(tTS,n,n);
}

/* Scene dessinee                               */

static void dessinerOpenGL(Pos3D *tP[4][4],int n,int type) {
  switch (type) {
    case 0 :
      dessinerNuageDePointsOpenGL(tP);
      break;
    case 1 :
      dessinerMaillageOpenGL(tP);
      break;
    case 2 :
      dessiner4CourbesBSplineOpenGL(tP,n);
      break;
    case 3 :
      dessinerNCourbesBSplineOpenGL(tP,n);
      break;
    case 4 :
      dessinerMaillageBSplineOpenGL(tP,n);
      break; }
}

static void scene(void) {
  switch (aff) {
    case 0 :
      glColor3f(0.75F,0.75F,0.75F);
      dessinerOpenGL(tPos,0,0);
      break;
    case 1 :
      glColor3f(0.75F,0.75F,0.75F);
      dessinerOpenGL(tPos,0,1);
      break;
    case 2 :
      glColor3f(0.75F,0.75F,0.75F);
      dessinerOpenGL(tPos,50,2);
      break;
    case 3 :
      glColor3f(0.75F,0.75F,0.75F);
      dessinerOpenGL(tPos,50,3);
      break;
    case 4 :
      glColor3f(0.75F,0.75F,0.75F);
      dessinerOpenGL(tPos,50,4);
      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);
  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 taille 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(-6.0*ratio,6.0*ratio,-6.0,6.0,-10.0,10.0);
    else
    glOrtho(-6.0,6.0,-6.0/ratio,6.0/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 0x0D :
      glutPostRedisplay(); 
      break;
    case 0x20 :
      aff = (aff+1)%5;
      glutPostRedisplay();
      break;
    case 0x1B :
      exit(0); 
      break; }
}

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

/* Fonction principale                          */

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

  atexit(clean);

  glutInit(&argc,argv);
  glutInitDisplayMode(GLUT_RGBA|GLUT_DEPTH|GLUT_DOUBLE);
  glutInitWindowSize(450,300); 
  glutInitWindowPosition(50,50); 
  glutCreateWindow("Une surface B-Spline"); 
  init();
  glutKeyboardFunc(keyboard);
  glutSpecialFunc(special);
  glutReshapeFunc(reshape);
  glutKeyboardFunc(keyboard);
  glutDisplayFunc(display);
  glutMainLoop();
  return(0);
}
