/* 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.

GLfloat bSpline(int k, int order, float t, float knots[]) {
  float denom1, denom2, sum = 0.0;
  if(1==order) 
    return (t >= knots[k] && t < knots[k+1]);
  denom1 = knots[k+order-1] - knots[k];
  if(denom1 != 0.0)
    sum += (t-knots[k])*bSpline(k,order-1,t,knots)/denom1;
  denom2 = knots[k+order] - knots[k+1];
  if(denom2 != 0.0)
    sum += (knots[k+order]-t)*bSpline(k+1,order-1,t,knots)/denom2;
  return sum;
}

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

GLfloat nurbBasis(int k, int order, float t, int L, GLfloat weights[], GLfloat knots[]) {
  float denom = 0.0;
  for(int i=0; i<=L; i++) 
    denom += weights[i]*bSpline(k,order,t,knots);
  return weights[k]*bSpline(k,order,t,knots);
}

// Evaluate a nurb curve using the nurbBasis function
void nurbEval(int order, float t, int L, GLfloat points[][4], GLfloat weights[], GLfloat knots[], GLfloat P[]) {
  P[0] = P[1] = P[2] = 0;
  float sum = 0;
  for(int cp=0; cp<=L; cp++) {
    float coef = nurbBasis(cp,order,t,L,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 L, GLfloat points[][4], GLfloat weights[], GLfloat knots[], GLfloat P[]) {
  P[0] = P[1] = P[2] = P[3] = 0;
  for(int cp=0; cp<=L; cp++) 
    for(int i=0; i<4; i++) 
      P[i] = weights[cp]*points[cp][i]* bSpline(cp,order,t,knots);
}

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,2,2,3,3,4,4,4}; //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);

    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)
    */
    gluNurbsCurve(myNurb,numKnots,uKnot,4,&points[0][0],uOrder,
                  GL_MAP1_VERTEX_4);
    gluEndCurve(myNurb);
 
    glColor3f(1,0,0);                // red
    gluBeginCurve(myNurb);
    gluNurbsCurve(myNurb,numKnots,uKnot,4,&pointsWeights[0][0],uOrder,
                  GL_MAP1_VERTEX_3);
    gluEndCurve(myNurb);
 
    glColor3f(0,0,0);                // black
    glLineWidth(1);
    glBegin(GL_LINE_LOOP); //lines between points
    for(i=0;i<numPoints;i++) {
      glVertex3fv(pointsWeights[i]);
    }
    glEnd();

    glColor3f(0,0,1);                // blue
    glBegin(GL_LINE_LOOP);
    for(i=0;i<numPoints;i++) {
      glVertex3fv(points[i]);
    }
    glEnd();
    
    glColor3f(0,1,0);                // green
    glBegin(GL_LINE_LOOP);
    for(i=0;i<numPoints;i++) {
      glVertex4fv(points[i]);
    }
    glEnd();
    
    // ================================================================
    // by hand

    glPointSize(4);
    GLfloat P[4];
    GLfloat weights[] = {1,W,1};
    GLfloat knots[12] = {0,0,0,0.25,0.25,0.5,0.5,0.75,0.75,1,1,1};

    for(int i=0;i<10;i++) {
      float c = 0;
      for(int j=0; j<12; j++) 
        c += bSpline(j,2,i/10.0,knots);
      if( c != 1 ) 
        printf("bSpline sum for t=%f is %f\n",i/10.0,c);
    }

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

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

    glFlush();
    glutSwapBuffers(); 
}

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

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