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

#define MAXPTS 20
#define EPS 0.000001
#define EQUIDISTANT_WIREFRAME 1
#define INTERPOLATE_WIREFRAME 2
#define EQUIDISTANT_ANALYTIC 3
#define INTERPOLATE_ANALYTIC 4
#define EQUIDISTANT_DEBOER 5
#define INTERPOLATE_DEBOER 6

int mousex, mousey;
int sirka=800, vyska=600;
int N,M,D1,D2;
double V[MAXPTS][MAXPTS][3],w[MAXPTS][MAXPTS];
int sel_pt_i=0,sel_pt_j=0,color=0,mode=0;
int is_ctrl,is_alt,is_shift;
double ZOOM=2.0,XOFFSET=-1.0,YOFFSET=0.0,ROTATEX=15.0,ROTATEY=0.0,STEP=0.05;
GLuint wlist,slist1,slist2;
int SHOW_PTS=1;

double s[100],t[100],N1[50][50],N2[50][50];

void PisVektorovyText(GLfloat x, GLfloat y, char *format,...) {
  va_list args;
  char buffer[200], *p;

  va_start(args, format); vsprintf(buffer, format, args); va_end(args); 
  
  glPushMatrix(); glLoadIdentity(); glTranslatef(-1.0,0.0,-6.0f); // glRotatef(10,1.0,0.0,0.0);
  glRotatef(0,0.0,1.0,0.0); glScalef(2,2,-2); glTranslatef(x, y, 0); glScalef(0.0007,0.0007,0.0007);
  for (p = buffer; *p; p++) glutStrokeCharacter(GLUT_STROKE_ROMAN, *p);
  glPopMatrix();
}
void SetVertexColor(double u, double v) { glColor3d(u,v,1.0-u-v); }

int ComputeAnalyticPoint(double u, double v, double *x, double *y, double *z){
  int i,j;
  double ss,tt,cit,men;

  if (u>=1.0) u=0.9999;
  if (v>=1.0) v=0.9999;
  ss=s[D1]+u*(s[N+1]-s[D1]); tt=t[D2]+v*(t[M+1]-t[D2]); 
  
  for (i=0;i<=N+D1;i++) if ((ss>=s[i])&&(ss<s[i+1])) N1[i][0]=1.0; else N1[i][0]=0.0;
  for (j=1;j<=D1;j++) for (i=0;i<=N+D1-j;i++) {
    N1[i][j]=0.0;
    if (s[i+j]  -s[i]  >EPS) N1[i][j] += N1[i][j-1]   * ( (ss-s[i])     / (s[i+j]-s[i])     );
    if (s[i+j+1]-s[i+1]>EPS) N1[i][j] += N1[i+1][j-1] * ( (s[i+j+1]-ss) / (s[i+j+1]-s[i+1]) );
  }
  
  for (i=0;i<=M+D2;i++) if ((tt>=t[i])&&(tt<t[i+1])) N2[i][0]=1.0; else N2[i][0]=0.0;
  for (j=1;j<=D2;j++) for (i=0;i<=M+D2-j;i++) {
    N2[i][j]=0.0;
    if (t[i+j]  -t[i]  >EPS) N2[i][j] += N2[i][j-1]   * ( (tt-t[i])     / (t[i+j]-t[i])     );
    if (t[i+j+1]-t[i+1]>EPS) N2[i][j] += N2[i+1][j-1] * ( (t[i+j+1]-tt) / (t[i+j+1]-t[i+1]) );
  }
  
  cit=men=0.0;
  for (i=0;i<=N;i++) for (j=0;j<=M;j++) men+=w[i][j]*N1[i][D1]*N2[j][D2];
  if (fabs(men)<EPS) return 0;
  for (i=0;i<=N;i++) for (j=0;j<=M;j++) cit+=w[i][j]*N1[i][D1]*N2[j][D2]*V[i][j][0]; *x = cit / men; cit=0.0;
  for (i=0;i<=N;i++) for (j=0;j<=M;j++) cit+=w[i][j]*N1[i][D1]*N2[j][D2]*V[i][j][1]; *y = cit / men; cit=0.0;
  for (i=0;i<=N;i++) for (j=0;j<=M;j++) cit+=w[i][j]*N1[i][D1]*N2[j][D2]*V[i][j][2]; *z = cit / men; cit=0.0;
  return 1;
}

void GenerateNURBWireframe(void){
  double u,v;
  double x[4],y[4],z[4];
  int p,q,i,ok;

  p=N+D1+1; q=M+D2+1;
  if (mode==1) {
    for (i=0;i<=p;i++) s[i]=i; for (i=0;i<=q;i++) t[i]=i; 
  } else {
    for (i=0;i<D1;i++) s[i]=0; for (i=0;i<D2;i++) t[i]=0; 
    for (i=D1;i<=N+1;i++) s[i]=i-D1; for (i=D2;i<=M+1;i++) t[i]=i-D2; 
    for (i=N+2;i<=p;i++) s[i]=N+1-D1; for (i=M+2;i<=q;i++) t[i]=M+1-D2; 
  }

  glNewList(wlist,GL_COMPILE);
  glColor3d(0.0,0.0,0.7); glLineWidth(3.0);
  for (i=0;i<=N;i++) { 
    u=(double)i/(double)N;
    glBegin(GL_LINE_STRIP);
    for (v=0.0;v<1.0+STEP-EPS;v+=STEP) {
      ok =ComputeAnalyticPoint(u,v,&x[0],&y[0],&z[0]);
      if (ok) glVertex3d(x[0],y[0],z[0]);
    }
    glEnd();
  }
  for (i=0;i<=M;i++) { 
    v=(double)i/(double)N;
    glBegin(GL_LINE_STRIP);
    for (u=0.0;u<1.0+STEP-EPS;u+=STEP) {
      ok =ComputeAnalyticPoint(u,v,&x[0],&y[0],&z[0]);
      if (ok) glVertex3d(x[0],y[0],z[0]);
    }
    glEnd();
  }
  glEndList();
}

void GenerateNURBSolidAnalytic(void){
  double u,v;
  double x[4],y[4],z[4];
  int p,q,i,ok;

  p=N+D1+1; q=M+D2+1;
  if (mode==3) {
    for (i=0;i<=p;i++) s[i]=i; for (i=0;i<=q;i++) t[i]=i; 
  } else {
    for (i=0;i<D1;i++) s[i]=0; for (i=0;i<D2;i++) t[i]=0; 
    for (i=D1;i<=N+1;i++) s[i]=i-D1; for (i=D2;i<=M+1;i++) t[i]=i-D2; 
    for (i=N+2;i<=p;i++) s[i]=N+1-D1; for (i=M+2;i<=q;i++) t[i]=M+1-D2; 
  }

  glNewList(slist1,GL_COMPILE);
  glBegin(GL_QUADS);
  for (u=0.0;u<1.0-EPS;u+=STEP) 
  for (v=0.0;v<1.0-EPS;v+=STEP) {
    ok =ComputeAnalyticPoint(u,v,&x[0],&y[0],&z[0]);
    ok&=ComputeAnalyticPoint(u,v+STEP,&x[1],&y[1],&z[1]);
    ok&=ComputeAnalyticPoint(u+STEP,v+STEP,&x[2],&y[2],&z[2]);
    ok&=ComputeAnalyticPoint(u+STEP,v,&x[3],&y[3],&z[3]);
    if (ok) {
      SetVertexColor(u,v);           glVertex3d(x[0],y[0],z[0]);
      SetVertexColor(u,v+STEP);      glVertex3d(x[1],y[1],z[1]);
      SetVertexColor(u+STEP,v+STEP); glVertex3d(x[2],y[2],z[2]);
      SetVertexColor(u+STEP,v);      glVertex3d(x[3],y[3],z[3]);
    } // else printf("Fault at %.2f %.2f\n",u,v);
  }
  glEnd();
  glEndList();
}

void GenerateNURBSolidDeBoer(void){
}

void display(void) {
  int i,j;

  glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT );
  glMatrixMode(GL_MODELVIEW);
  glLoadIdentity();
  glTranslatef(XOFFSET,YOFFSET,-6.0f); // move into the screen 6.0
  glRotatef(ROTATEX,1.0,0.0,0.0);
  glRotatef(ROTATEY,0.0,1.0,0.0);
  glScalef(ZOOM,ZOOM,-ZOOM);

  glColor3d(1.0,1.0,1.0); glLineWidth(2.0);
  PisVektorovyText(-0.9,1.0,"N=%d M=%d d1=%d d2=%d point=[%d,%d] wt=%.1f step=%.2f",
      N,M,D1,D2,sel_pt_i,sel_pt_j,w[sel_pt_i][sel_pt_j],STEP);

  if (SHOW_PTS) {
    glColor3d(0.5,0.5,0.5); glLineWidth(1.0);
    for (i=0;i<=N;i++) { glBegin(GL_LINE_STRIP); for (j=0;j<=M;j++) glVertex3dv(V[i][j]); glEnd(); }
    for (j=0;j<=M;j++) { glBegin(GL_LINE_STRIP); for (i=0;i<=N;i++) glVertex3dv(V[i][j]); glEnd(); }
  }

  if ((mode==1)||(mode==2)) glCallList(wlist);
  if ((mode==3)||(mode==4)) glCallList(slist1);
  if ((mode==5)||(mode==6)) glCallList(slist2);
  
  if (SHOW_PTS) {
    glColor3d(0.5,0.5,0.5);
    glPointSize(5.0);
    glBegin(GL_POINTS); for (i=0;i<=N;i++) for (j=0;j<=M;j++) glVertex3dv(V[i][j]); glEnd();
    glColor3d(0.0,1.0,0.0);
    glPointSize(7.0);
    glBegin(GL_POINTS); glVertex3dv(V[sel_pt_i][sel_pt_j]); glEnd();
  }
  
  glutSwapBuffers();
}

void reshape(int w, int h) {
  sirka=w; vyska=h;
  glViewport(0, 0, w, h);
  glMatrixMode(GL_PROJECTION);
  glLoadIdentity();
  gluPerspective(45.0f,(GLfloat)sirka/(GLfloat)vyska,0.1f,100.0f);
  glMatrixMode(GL_MODELVIEW);
}

void keypressed(unsigned char key, int x, int y) { 
//  int i,j,k;

  is_ctrl=glutGetModifiers(); 
  is_alt=is_ctrl&GLUT_ACTIVE_ALT; 
  is_shift=is_ctrl&GLUT_ACTIVE_SHIFT; 
  is_ctrl=is_ctrl&GLUT_ACTIVE_CTRL;
  switch (key) {
    case 27:  exit(0); break;
    case '1': D1--; 
      if ((mode==1)||(mode==2)) GenerateNURBWireframe();
      if ((mode==3)||(mode==4)) GenerateNURBSolidAnalytic();
      if ((mode==5)||(mode==6)) GenerateNURBSolidDeBoer();
      break;
    case '!': D1++; 
      if ((mode==1)||(mode==2)) GenerateNURBWireframe();
      if ((mode==3)||(mode==4)) GenerateNURBSolidAnalytic();
      if ((mode==5)||(mode==6)) GenerateNURBSolidDeBoer();
      break;
    case '2': D2--; 
      if ((mode==1)||(mode==2)) GenerateNURBWireframe();
      if ((mode==3)||(mode==4)) GenerateNURBSolidAnalytic();
      if ((mode==5)||(mode==6)) GenerateNURBSolidDeBoer();
      break;
    case '@': D2++;
      if ((mode==1)||(mode==2)) GenerateNURBWireframe();
      if ((mode==3)||(mode==4)) GenerateNURBSolidAnalytic();
      if ((mode==5)||(mode==6)) GenerateNURBSolidDeBoer();
      break;
    case 'm': if (M-2) { M--; sel_pt_i=sel_pt_j=0; }
      if ((mode==1)||(mode==2)) GenerateNURBWireframe();
      if ((mode==3)||(mode==4)) GenerateNURBSolidAnalytic();
      if ((mode==5)||(mode==6)) GenerateNURBSolidDeBoer();
      break;
    case 'M': if (!M-MAXPTS+5) M++;
      if ((mode==1)||(mode==2)) GenerateNURBWireframe();
      if ((mode==3)||(mode==4)) GenerateNURBSolidAnalytic();
      if ((mode==5)||(mode==6)) GenerateNURBSolidDeBoer();
      break;
    case 'n': if (N-2) { N--; sel_pt_i=sel_pt_j=0; }
      if ((mode==1)||(mode==2)) GenerateNURBWireframe();
      if ((mode==3)||(mode==4)) GenerateNURBSolidAnalytic();
      if ((mode==5)||(mode==6)) GenerateNURBSolidDeBoer();
      break;
    case 'N': if (!N-MAXPTS+5) N++;
      if ((mode==1)||(mode==2)) GenerateNURBWireframe();
      if ((mode==3)||(mode==4)) GenerateNURBSolidAnalytic();
      if ((mode==5)||(mode==6)) GenerateNURBSolidDeBoer();
      break;
    case '+': ZOOM+=STEP; break;
    case '-': ZOOM-=STEP; break;
    case 'S': 
      STEP+=0.01; 
      if ((mode==1)||(mode==2)) GenerateNURBWireframe();
      if ((mode==3)||(mode==4)) GenerateNURBSolidAnalytic();
      if ((mode==5)||(mode==6)) GenerateNURBSolidDeBoer();
    break;
    case 's': 
      if (STEP>0.01) STEP-=0.01; 
      if ((mode==1)||(mode==2)) GenerateNURBWireframe();
      if ((mode==3)||(mode==4)) GenerateNURBSolidAnalytic();
      if ((mode==5)||(mode==6)) GenerateNURBSolidDeBoer();
    break;
    case 'p': SHOW_PTS^=1; break;
  }
  glutPostRedisplay();
}

void keypressed_special(int key, int x, int y) {
  is_ctrl=glutGetModifiers(); 
  is_alt=is_ctrl&GLUT_ACTIVE_ALT; 
  is_shift=is_ctrl&GLUT_ACTIVE_SHIFT; 
  is_ctrl=is_ctrl&GLUT_ACTIVE_CTRL;
 
  if (is_alt) {
    switch (key) {
      case GLUT_KEY_UP:    ROTATEX+=2; break;
      case GLUT_KEY_DOWN:  ROTATEX-=2; break;
      case GLUT_KEY_LEFT:  ROTATEY+=2; break;
      case GLUT_KEY_RIGHT: ROTATEY-=2; break;
    }
  }
  if (is_ctrl) {
    switch (key) {
      case GLUT_KEY_UP:    YOFFSET-=STEP; break;
      case GLUT_KEY_DOWN:  YOFFSET+=STEP; break;
      case GLUT_KEY_LEFT:  XOFFSET+=STEP; break;
      case GLUT_KEY_RIGHT: XOFFSET-=STEP; break;
    }
  }
  if (is_shift) {
    switch (key) {
      case GLUT_KEY_PAGE_UP:    V[sel_pt_i][sel_pt_j][1]+=0.025; break;
      case GLUT_KEY_PAGE_DOWN:  V[sel_pt_i][sel_pt_j][1]-=0.025; break;
      case GLUT_KEY_UP:         V[sel_pt_i][sel_pt_j][2]+=0.025; break;
      case GLUT_KEY_DOWN:       V[sel_pt_i][sel_pt_j][2]-=0.025; break;
      case GLUT_KEY_LEFT:       V[sel_pt_i][sel_pt_j][0]-=0.025; break;
      case GLUT_KEY_RIGHT:      V[sel_pt_i][sel_pt_j][0]+=0.025; break;
    }
    if ((mode==1)||(mode==2)) GenerateNURBWireframe();
    if ((mode==3)||(mode==4)) GenerateNURBSolidAnalytic();
    if ((mode==5)||(mode==6)) GenerateNURBSolidDeBoer();
  }
  if ((!is_alt)&&(!is_ctrl)&&(!is_shift)) {
    switch (key) {
      case GLUT_KEY_UP:    if (sel_pt_j<M) sel_pt_j++; break;
      case GLUT_KEY_DOWN:  if (sel_pt_j>0) sel_pt_j--; break;
      case GLUT_KEY_LEFT:  if (sel_pt_i>0) sel_pt_i--; break;
      case GLUT_KEY_RIGHT: if (sel_pt_i<N) sel_pt_i++; break;
    }
  }
  glutPostRedisplay();
}

void motion(int x, int y) { mousex = x; mousey = y; }
void mouse(int b, int s, int x, int y) { 
  mousex = x; mousey = y; 
  if ((b==3)||(b==4)) {
    w[sel_pt_i][sel_pt_j]+=-0.05+0.1*(4-b);
    if ((mode==1)||(mode==2)) GenerateNURBWireframe();
    if ((mode==3)||(mode==4)) GenerateNURBSolidAnalytic();
    if ((mode==5)||(mode==6)) GenerateNURBSolidDeBoer();
    glutPostRedisplay();
  }
  motion(x,y); 
}
void menu(int choice) { 
  mode=choice; 
  if ((mode==1)||(mode==2)) GenerateNURBWireframe();
  if ((mode==3)||(mode==4)) GenerateNURBSolidAnalytic();
  if ((mode==5)||(mode==6)) GenerateNURBSolidDeBoer();
  glutPostRedisplay(); 
}

int main(int argc, char **argv) {
  int i,j; 
  
  glutInit(&argc, argv);
  glutInitDisplayMode(GLUT_DOUBLE | GLUT_RGB | GLUT_DEPTH );
  
  glutInitWindowSize(800, 600); glutInitWindowPosition(0, 0);
  glutCreateWindow("Miso Forisek -- Krivky a plochy v PG");
  glClearColor(0.0, 0.0, 0.0, 0.0); glColor3d(1.0, 0.0, 0.0); glLineWidth(1.0);
 
  glClearDepth(10.0);                           // Enables Clearing Of The Depth Buffer
  glDepthFunc(GL_LESS);                         // The Type Of Depth Test To Do
  glEnable(GL_DEPTH_TEST);                      // Enables Depth Testing
  glShadeModel(GL_SMOOTH);                      // Enables Smooth Color Shading

  glMatrixMode(GL_PROJECTION);
  glLoadIdentity();                             // Reset The Projection Matrix
  gluPerspective(45.0f,(GLfloat)sirka/(GLfloat)vyska,0.1f,100.0f);     // Calculate The Aspect Ratio Of The Window
  glMatrixMode(GL_MODELVIEW);
 
  glutDisplayFunc(display);
  glutReshapeFunc(reshape);
  glutKeyboardFunc(keypressed);
  glutSpecialFunc(keypressed_special);
  glutMotionFunc(motion);
  glutMouseFunc(mouse);
  
  glutCreateMenu(menu); 
  glutAddMenuEntry("Grid only", 0); 
  glutAddMenuEntry("Wire-frame model (equidistant)", 1); 
  glutAddMenuEntry("Wire-frame model (interpolates corners)", 2); 
  glutAddMenuEntry("Solid model -- analytic (equidistant)", 3); 
  glutAddMenuEntry("Solid model -- analytic (interpolates corners)", 4); 
//  glutAddMenuEntry("Solid model -- deBoer (equidistant)", 5); 
//  glutAddMenuEntry("Solid model -- deBoer (interpolates corners)", 6); 
  glutAttachMenu(GLUT_RIGHT_BUTTON);

  M=N=5; D1=D2=3;
  for (i=0;i<MAXPTS;i++) for (j=0;j<MAXPTS;j++) 
    { V[i][j][0]=i/5.0; V[i][j][2]=j/5.0; V[i][j][1]=0.0; w[i][j]=1.0; }
  wlist=glGenLists(1); 
  slist1=glGenLists(1); 
  slist2=glGenLists(1); 

  glutMainLoop();
  return 0;
}