/* NURB curve (1-dimension)

Testing whether we really get a quarter circle with the NURBS given in
Rogers and also in Hill, and whether the curve goes outside the convex
hull.

Written by Caroline Geiersbach and Scott D. Anderson
scott.anderson@acm.org
Summer 2003
*/

#include <stdio.h>
#include <stdlib.h>
#include <GL/glut.h>
#include <math.h>

GLUnurbsObj *myNurb;
int uOrder=3; 
const int numPoints = 3;

GLfloat W = sqrt(2)/2;

GLfloat points[numPoints][4];
GLfloat pointsWeights[numPoints][4] = {
   {1,0,0,1},
   {1,1,0,W},
   {0,1,0,1}
};


// evaluate the ith b-spline function using Cox-de Boor recursion.
// The index "k" is from 0 to L, where there are L control points;
// that is, it chooses the control point that corresponds to this
// basis spline.  This code is Hill, page 631.

bool trace = false;

GLfloat bSpline(int k, int order, float t, float knots[]) {
  char* indent[4] = {"","        ","      ","    "};
  float term1=0, term2=0, denom1, denom2, bk, bk1;
  if(trace) printf("%s bSpline(%d,%d,%f)\n",indent[order],k,order,t);
  if(1==order) {
    float result = (t >= knots[k] && t < knots[k+1]);
    if(trace) printf("%s bSpline(%d,%d,%f)=>%f\n",indent[order],
                     k,order,t,result);
    return result;
  }
  denom1 = knots[k+order-1] - knots[k];
  if(denom1 == 0.0) {
    if(trace) printf("%s denom1 is zero\n",indent[order]);
  } else {
    bk=bSpline(k,order-1,t,knots);
    term1 = (t-knots[k])*bk/denom1;
    if(trace) printf("%s t-knots[k]=%f bk=%f denom1=%f term1=%f\n",indent[order],
                     t-knots[k],bk,denom1,term1);
  }
  denom2 = knots[k+order] - knots[k+1];
  if(denom2 == 0.0) {
    if(trace) printf("%s denom2 is zero\n",indent[order]);
  } else {
    float bk1=bSpline(k+1,order-1,t,knots);
    term2 = (knots[k+order]-t)*bk1/denom2;
    if(trace) printf("%s knots[k+order]-t=%f bk1=%f denom2=%f term2=%f\n",indent[order],
                     knots[k+order]-t,bk1,denom2,term2);
  }
  if(trace) printf("%s bSpline(%d,%d,%f)=>%f\n",indent[order],k,order,t,term1+term2);
  return term1+term2;
}

// This is the rational version of the bSpline function for use with
// non-homogenous coordinates.

GLfloat nurbBasis(int k, int order, float t, int numPts, GLfloat weights[], GLfloat knots[]) {
  float denom = 0.0;
  for(int i=0; i<numPts; i++) {
    denom += weights[i]*bSpline(i,order,t,knots);
  }
  if(trace) printf("nurbBasis: t=%f k=%d weights[k]=%f denom = %f\n",t,k,weights[k],denom);
  float result = weights[k]*bSpline(k,order,t,knots)/denom;
  if(trace) printf("nurbBasis: result = %f\n",result);
  return result;
}

// Evaluate a nurb curve using the nurbBasis function
void nurbEval(int order, float t, int numPts, GLfloat points[][4], GLfloat weights[], GLfloat knots[], GLfloat P[]) {
  P[0] = P[1] = P[2] = 0;
  float sum = 0;
  for(int cp=0; cp<=numPts; cp++) {
    if(trace) printf("nurbEval: t=%f cp=%d\n",t,cp);
    float coef = nurbBasis(cp,order,t,numPts,weights,knots);
    sum += coef;
    if(coef < 0) 
      printf("coef = %f < 0 when k=%d\n",coef,cp);
    if(coef > 1) 
      printf("coef = %f > 1 when k=%d\n",coef,cp);
    for(int i=0; i<3; i++) 
      P[i] += coef*points[cp][i];
  }
  if(sum > 1.01 || sum < 0.99) 
    printf("sum isn't unity:  %f\n",sum);
}

// eval a nurb curve using the homogeneous coordinate idea of Hill
void nurbEval2(int order, float t, int numPts, GLfloat points[][4], GLfloat weights[], GLfloat knots[], GLfloat P[]) {
  P[0] = P[1] = P[2] = P[3] = 0;
  for(int cp=0; cp<numPts; cp++) {
    for(int i=0; i<4; i++) {
      if(trace) 
        printf("cp=%d i=%d weights[cp]=%f points[cp][i]=%f bspline=%f\n",
               cp,i,weights[cp],points[cp][i],bSpline(cp,order,t,knots));
      P[i] += weights[cp]*points[cp][i]* bSpline(cp,order,t,knots);
    }
    if(trace) printf("cp=%d P=(%f,%f,%f,%f)\n",P[0],P[1],P[2],P[3]);
  }
}

int whatToDisplay = 15;

void display(void) {
    int i,j;
    glClearColor(0.7,0.7,0.7,1);
    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
    glMatrixMode(GL_PROJECTION);
    glLoadIdentity();
    // since the figure is in the unit square, this should work.
    float margin = 0.2;
    glOrtho(0-margin,1+margin,
            0-margin,1+margin,
            0,1);
    glMatrixMode(GL_MODELVIEW);
    glLoadIdentity();
   
    pointsWeights[1][3] = W;

    int numKnots = uOrder+numPoints;

    GLfloat uKnot[12] = {0,0,0,1,1,1}; // for square knot
    
    for (i = 0; i < numPoints; i++) {
      points[i][3] = pointsWeights[i][3];
      printf("weight = %f\n",points[i][3]);
      for(j = 0; j<3; j++) {
        points[i][j] = pointsWeights[i][j]*points[i][3];
      }
      printf("%f %f %f %f\n",points[i][0],points[i][1],points[i][2],points[i][3]);
    }
    glLineWidth(3);

    if(whatToDisplay & 1) {
      glColor3f(1,0.5,0);                // orange
      gluBeginCurve(myNurb);
      /*
        gluNurbsCurve(GLUnurbsObj *myNurb, GLint uKnotCount, GLfloat *uKnot,
        GLint uStride, GLfloat *points, GLint uOrder(degree+1), GLenum type)
      */
      // this definitely works.
      gluNurbsCurve(myNurb,numKnots,uKnot,4,&points[0][0],uOrder,
                    GL_MAP1_VERTEX_4);
      gluEndCurve(myNurb);
    }
 
    if(whatToDisplay & 2) {
      glColor3f(1,0,0);                // red
      gluBeginCurve(myNurb);
      // this definitely works.
      gluNurbsCurve(myNurb,numKnots,uKnot,4,&pointsWeights[0][0],uOrder,
                    GL_MAP1_VERTEX_3);
      gluEndCurve(myNurb);
    }
 
    if(whatToDisplay & 4) {
      glColor3f(0,0,0);                // black
      glLineWidth(1);
      glBegin(GL_LINE_STRIP); //lines between points
      for(i=0;i<numPoints;i++) {
        glVertex3fv(pointsWeights[i]);
      }
      glEnd();
    }

    if(whatToDisplay & 8) {
      glColor3f(0,0,1);                // blue
      glBegin(GL_LINE_STRIP);
      for(i=0;i<numPoints;i++) {
        glVertex3fv(pointsWeights[i]);
      }
      glEnd();
      glBegin(GL_LINE_STRIP);
      for(float theta=0.0; theta<M_PI_2; theta += 0.1) 
        glVertex2f(cos(theta),sin(theta));
      glEnd();
    }
    
    if(whatToDisplay & 16) {
      glColor3f(0,1,0);                // green
      glBegin(GL_LINE_STRIP);
      for(i=0;i<numPoints;i++) {
        glVertex4fv(points[i]);
      }
      glEnd();
    }
    
    // ================================================================
    // by hand

    GLfloat points[3][4] = {
      {1,0,0,1},
      {1,1,0,1},
      {0,1,0,1}};

    int numPts=3;
    int order=3;
    glPointSize(4);
    GLfloat P[4];
    GLfloat weights[] = {1,W,1};
    float knots[6] = {0,0,0,1,1,1};

    if(false){
      for(int i=0;i<10;i++) {
        float c = 0;
        for(int j=0; j<numPts; j++) 
          c += bSpline(j,order,i/10.0,knots);
        if( c>1.001 || c < 0.999 ) {
          printf("bSpline sum for t=%f is %f\n",i/10.0,c);
          trace=true;
          for(int j=0; j<numPts; j++) 
            c += bSpline(j,2,i/10.0,knots);
          trace=false;
        }
      }
    }

    float result;
    if(false) {
      trace=true;
      result=0;
      result+=bSpline(0,3,0,knots);
      result+=bSpline(1,3,0,knots);
      result+=bSpline(2,3,0,knots);
      printf("result = %f\n",result);
      exit(1);   
    } 

    trace=false;
    printf("call nurbBasis\n");
    result=nurbBasis(1,order,0,numPts,weights,knots);
    printf("result = %f\n",result);

    glColor3f(1,1,0);                // yellow
    printf("yellow\n");
    glBegin(GL_POINTS);
    trace=false;
    for(i=0;i<=10;i++) {
      // printf("t=%f\n",i/10.0);
      nurbEval(order,i/10.0,numPts,pointsWeights,weights,knots,P);
      printf("t=%f P=(%f,%f)\n",i/10.0,P[0],P[1]);
      float theta=M_PI/2*i/10.0;
      printf("T=%f P=(%f,%f)\n",theta,cos(theta),sin(theta));
      glVertex3fv(P);
    }
    glEnd();

    glColor3f(0,1,1);                // cyan
    glBegin(GL_POINTS);
    pointsWeights[0][3] = pointsWeights[1][3] = pointsWeights[2][3] = 1;
    for(i=0;i<10;i++) {
      nurbEval2(order,i/10.0,numPts,pointsWeights,weights,knots,P);
      printf("P=(%f,%f,0,%f)=(%f,%f)\n",P[0],P[1],P[3],P[0]/P[3],P[1]/P[3]);
      glVertex4fv(P);
    }
    glEnd();

    glFlush();
    glutSwapBuffers(); 
}

int tmp;

void keys(unsigned char k, int x, int y) {
  if(k=='q') exit(0);
  if(k=='?') printf("W = %f\n",W);
  if(k=='w') {
    char buf[80];
    printf("W = %f; what new value? ",W);
    fgets(buf,sizeof(buf),stdin);
    sscanf(buf,"%f",&W);
    glutPostRedisplay();
  }
  if(k >= '0' && k <= '9') tmp = tmp*10 + k-'0';
  if(k == '\r') { whatToDisplay=tmp; tmp=0; printf("what to display = %d\n",whatToDisplay); glutPostRedisplay(); }
}

int main(int argc, char **argv) {
  glutInit(&argc, argv);
  glutInitDisplayMode(GLUT_DOUBLE | GLUT_RGB | GLUT_DEPTH);
  glutInitWindowSize(650, 650);
  glutCreateWindow("NURBS!");
  glutDisplayFunc(display);
  glutKeyboardFunc(keys);
  myNurb=gluNewNurbsRenderer();
  gluNurbsProperty(myNurb,GLU_DISPLAY_MODE,GLU_FILL);
  glEnable(GL_MAP1_VERTEX_4);
  glutMainLoop();
  return 0;
}