#include <string.h>
#include <stdio.h>
#include <stdarg.h>
#include <math.h>
#include <GL/glut.h>

#define MAXPTS 50
#define MAXM 100
#define EPS 0.000001
#define STEP 0.001
#define ANALYTIC 1
#define DEBOOR 2

int mousex, mousey;
int sirka=640, vyska=480;
int pts=0;
double P[MAXPTS][3];
int sel_pt=-1;

double W[MAXPTS],T[MAXM],Tdisp[MAXM][2];
int col,i,D=3,POPISKY=0,MODE;
double N[MAXM][MAXM], alpha[MAXM][MAXM], V[MAXM][MAXM][2], w[MAXM][MAXM];

void PisBitmapovyText(int x, int y, char *string, void *font) {
  int len, i;
  glRasterPos2f(x, y);
  len = (int) strlen(string);
  for (i = 0; i < len; i++) glutBitmapCharacter(font, string[i]);
}

int RatajAnalytic(double t, double *x, double *y){ 
  int i,k;
  double cit,men;

  for (i=0;i<pts+D;i++) 
    if ((t>=T[i]) && (t<T[i+1])) { N[i][0]=1.0; col=i; } else N[i][0]=0.0;
  for (k=1;k<=D;k++) for (i=0;i<pts+D-k;i++) {
    N[i][k] = 0;
    if (T[i+k]  -T[i]  >EPS) N[i][k] += N[i][k-1]   * ( (t-T[i])     / (T[i+k]-T[i])     );
    if (T[i+k+1]-T[i+1]>EPS) N[i][k] += N[i+1][k-1] * ( (T[i+k+1]-t) / (T[i+k+1]-T[i+1]) );
  }
  cit=men=0.0;
  for (i=0;i<pts;i++) men+=W[i]*N[i][D]; if (men<EPS) return 0;
  for (i=0;i<pts;i++) cit+=W[i]*N[i][D]*P[i][0]; *x = cit / men; cit=0.0;
  for (i=0;i<pts;i++) cit+=W[i]*N[i][D]*P[i][1]; *y = cit / men;
  return 1;
}

int RatajDeBoor(double t, double *x, double *y){ 
  int r,i,j;
  
  r=0; while (r<=pts+D-1) if ((T[r]<=t) && (t<T[r+1])) break; else r++;
  for (i=0;i<pts;i++) { V[i][0][0]=P[i][0]; V[i][0][1]=P[i][1]; w[i][0]=W[i]; }
  for (j=1;j<=D;j++) for (i=r-D+j;i<=r;i++) {
    alpha[i][j]=T[i+D-j+1]-T[i];
    if (fabs(alpha[i][j])<EPS) alpha[i][j]=0; else {
      alpha[i][j]=1/alpha[i][j]; alpha[i][j]*=t-T[i];
    }
    w[i][j]=(1-alpha[i][j])*w[i-1][j-1] + alpha[i][j]*w[i][j-1];
    if (fabs(w[i][j])<EPS) { V[i][j][0]=0.0; V[i][j][1]=0.0; } else {
      V[i][j][0] = (1-alpha[i][j])*w[i-1][j-1]*V[i-1][j-1][0] + alpha[i][j]*w[i][j-1]*V[i][j-1][0];
      V[i][j][1] = (1-alpha[i][j])*w[i-1][j-1]*V[i-1][j-1][1] + alpha[i][j]*w[i][j-1]*V[i][j-1][1];
      V[i][j][0] /= w[i][j];
      V[i][j][1] /= w[i][j];
    }
  }
  *x = V[r][D][0]; *y = V[r][D][1]; col=r;
  return 1;
}

void display(void) {
  double t,x,y;
  int i,ok;
  char pom[100];
  
  //
  // HORNY POKEC
  // 
  glClear(GL_COLOR_BUFFER_BIT); glColor3d(0.7, 0.7, 0.7);
  strcpy(pom,"Left button: add/remove points. Middle button (drag&drop): move points.");
  strcat(pom," Wheel: adjust point weights. Right button: select method."); 
  PisBitmapovyText(5,15,pom,GLUT_BITMAP_TIMES_ROMAN_10);
  sprintf(pom,"Curve dimension: %d (adjust using +-).",D);
  PisBitmapovyText(5,27,pom,GLUT_BITMAP_TIMES_ROMAN_10);
  PisBitmapovyText(5,39,"Adjust the knot vector by moving the points below. (S toggles showing values.)"
    ,GLUT_BITMAP_TIMES_ROMAN_10);
  
  //
  // UZLOVY VEKTOR
  // 
  glColor3d(0.4, 0.4, 0.0); glPointSize(5); 
  glBegin(GL_LINE_STRIP); glVertex2d(10.0,50.0); glVertex2d(sirka-60.0,50.0); glEnd();
  glBegin(GL_POINTS);
    for (i=0;i<=pts+D;i++) {
      Tdisp[i][0]=10.0+(sirka-70.0)*T[i]/(T[pts+D]-T[0]); Tdisp[i][1]=50.0;
      glVertex2d(Tdisp[i][0],Tdisp[i][1]);
    }
  glEnd();
 
  //
  // POPIS UZLOVEHO VEKTORU
  // 
  glColor3d(0.7, 0.7, 0.7); 
  if (POPISKY) for (i=0;i<=pts+D;i++) {      
    sprintf(pom,"t_%d = %.1f",i,T[i]);
    PisBitmapovyText(10.0+(sirka-70.0)*T[i]/(T[pts+D]-T[0]),62.0,pom,GLUT_BITMAP_TIMES_ROMAN_10);
  }
  
  // 
  // VODOROVNA CIARA
  // 
  glBegin(GL_LINE_STRIP); glVertex2d(0.0,67.0); glVertex2d(sirka,67.0); glEnd();

  //
  // SIVE CIARY MEDZI RIADIACIMI BODMI
  // 
  glColor3d(0.3, 0.3, 0.3); glBegin(GL_LINE_STRIP); for (i=0;i<pts;i++) { glVertex3dv(P[i]); } glEnd();

  //
  // SAMOTNA KRIVKA
  // 
  if ((pts>1) && (T[pts]>T[D])) {
  glColor3d(1.0, 1.0, 0.0);
  glBegin(GL_LINE_STRIP);
    for (t=T[D];t<=T[pts];t+=STEP*(T[pts]-T[D])) {
      switch (MODE) {
        case ANALYTIC: ok=RatajAnalytic(t,&x,&y); break;
        case DEBOOR:   ok=RatajDeBoor(t,&x,&y); break;
      }
      if (ok) {
        if ((col+MODE)%2) glColor3d(1.0,1.0,0.0); else glColor3d(1.0,0.33,0.0); glVertex3d(x,y,0.0);
      }
    }
  glEnd();
  }

  // 
  // RIADIACE BODY
  // 
  glPointSize(10); glColor3d(0.6, 0.0, 0.0); glBegin(GL_POINTS); for (i=0;i<pts;i++) glVertex3dv(P[i]); glEnd(); 
  
  //
  // VAHY RIADIACICH BODOV
  // 
  glColor3d(0.5, 0.5, 0.5);
  for (i=0;i<pts;i++) {
    sprintf(pom,"%.1f",W[i]);
    PisBitmapovyText((int)(P[i][0]+5),(int)(P[i][1]+5),pom,GLUT_BITMAP_TIMES_ROMAN_10);
  }

  //
  // KONIEC PREDSTAVENIA
  // 
  glutSwapBuffers();
}

void reshape(int w, int h) {
  sirka=w; vyska=h;
  glViewport(0, 0, w, h);
  glMatrixMode(GL_PROJECTION);
  glLoadIdentity();
  gluOrtho2D(0, w, 0, h);
  glScalef(1, -1, 1);
  glTranslatef(0, -h, 0);
  glMatrixMode(GL_MODELVIEW);
}

void keyPressed(unsigned char key, int x, int y) { 
  if (key == 27) exit(0); 
  if ((key == '-') && (D>1)) { D--; glutPostRedisplay(); }
  if ((key == '+') && (D<MAXM-MAXPTS)) { D++; T[D+pts]=T[D+pts-1]+1.0; glutPostRedisplay(); }
  if ((key == 's') || (key == 'S')) { POPISKY=!POPISKY; glutPostRedisplay(); }
}

void motion(int x, int y) {
  int p;
  double newt;
    
  if (sel_pt>-1) { P[sel_pt][0]=x; P[sel_pt][1]=y; }
  if (sel_pt<-1) { 
    p=-(sel_pt+2); // hybeme knotbodom p
    newt=T[0] + (T[pts+D]-T[0])*(x-10.0)/(sirka-70.0); 
    if (newt<T[p-1]) newt=T[p-1]; if (newt>T[p+1]) newt=T[p+1]; 
    T[p]=newt;
  }
  mousex = x; mousey = y;
  glutPostRedisplay();
}

void mouse(int b, int s, int x, int y) {
  int i,j,ok;
  double add;

  mousex = x; mousey = y;
  if (b==GLUT_LEFT_BUTTON){
    if (s==GLUT_DOWN) {
      ok=1;
      for (i=0;i<pts;i++) if ((abs(P[i][0]-x)<=4) && (abs(P[i][1]-y)<=4)) {
        ok=0;  
        for (j=i;j<pts-1;j++) { P[j][0]=P[j+1][0]; P[j][1]=P[j+1][1]; P[j][2]=P[j+1][2]; }
        pts--; 
        break;
      }
      if ((ok) && (pts<MAXPTS)) { P[pts][0]=x; P[pts][1]=y; P[pts][2]=0; pts++; }
    }
  }
  
  if (b==GLUT_MIDDLE_BUTTON) {
    if (s==GLUT_DOWN) {
      sel_pt=-1;
      for (i=0;i<pts;i++) if ((abs(P[i][0]-x)<=4) && (abs(P[i][1]-y)<=4)) { sel_pt=i; break; }
      for (i=1;i<pts+D;i++) if ((abs(Tdisp[i][0]-x)<=4) && (abs(Tdisp[i][1]-y)<=4)) { sel_pt=-i-2; break; }
    } else { sel_pt=-1; }
  }

  if ((b==3) || (b==4)) {
    // kolecko
    if (b==3) add=0.05; else add=-0.05;
    for (i=0;i<pts;i++) if ((abs(P[i][0]-x)<=4) && (abs(P[i][1]-y)<=4)) {
      W[i]+=add; if (W[i]<0.0) W[i]=0.0; break;
    }
  }
  
  glutPostRedisplay();
  motion(x,y);
}

void menu(int choice) {
  switch(choice) {
    case 1: MODE=ANALYTIC; break;
    case 2: MODE=DEBOOR; break;
  }
  glutPostRedisplay();
}

int main(int argc, char **argv) {
  glutInit(&argc, argv);
  glutInitDisplayMode(GLUT_DOUBLE | GLUT_RGB);
  
  glutInitWindowSize(640, 480); glutInitWindowPosition(0, 0);
  glutCreateWindow("Miso Forisek -- DU #3");
  glClearColor(0.0, 0.0, 0.0, 0.0); glColor3d(1.0, 0.0, 0.0); glLineWidth(1.0);
  
  glutDisplayFunc(display);
  glutReshapeFunc(reshape);
  glutKeyboardFunc(keyPressed);
  glutMotionFunc(motion);
  glutMouseFunc(mouse);
  
  glutCreateMenu(menu); glutAddMenuEntry("Analytic", 1); glutAddMenuEntry("Deboor", 2);
  glutAttachMenu(GLUT_RIGHT_BUTTON);

  for (i=0;i<200;i++) { W[i]=1.0; T[i]=i; } MODE=ANALYTIC;

  glutMainLoop();
  return 0;
}