Syllabus   Blank Homework   Quizzes  
Notes   Labs   Scores   Blank

Lecture Notes
Dr. Tong Lai Yu, 2010
    1. Introduction
    2. OpenGL Shading Language ( GLSL ) I
    3. GLSL II
    4. Curve and Surface Design
    5. Modeling Shapes with Polygonal Meshes
    6. Texture Mapping
    7. Casting Shadows
    8. Tools for Raster Display
    9. Parsing External Objects

    Curves and Surfaces

    1. Representation of Curves and Surfaces

      Explicit Representation

    2. curve
        y = f ( x )

        z = g ( x )

        Examples:

        2-D line: y = m.x + b

        2-D circle: y = ( r2 - x2 ) 1/2     0 ≤ |x| < r

    3. surface ( needs two variables )
        z = f ( x, y )

        e.g. a spherical surface: z = ( r2 - x2 - y2 ) 1/2

      a curve or surface may not have an explicit representation
    4. Implicit Representation

    5. curve
        f ( x, y ) = 0

        Examples:

        2-D line: a.x + b.y + c = 0

        2-D circle: x2 + y2 - r2 = 0

    6. surface
        f ( x, y, z ) = 0

        e.g. a plane: a.x + b.y + c.z + d = 0

        e.g. a spherical surface: x2 + y2 + z2 - r2 = 0

      We can represent a curve as the intersection of two surfaces:

        f( x, y, z ) = 0

        g( x, y, z ) = 0

      Algebraic surfaces are those for which f( x, y, z ) is the sum of polynomials in the three variables. Of particular importance are the quadratic surface where each term in f can have degree up to 2 ( e.g. x, xy, or z2 but not xy2 ).

      See

      and
    7. Parametric form

    8. curve

      expresses the value of each spatial variable for points in terms of an independent variable, u, the parameter.

        x = x(u),

        y = y(u),

        z = z(u).

      It is in the same form in two and three dimensions.

      A useful interpretation of the parametric form is to visualize the locus of points

        p(u) = [x(u), y(u), z(u)]T
      being drawn as u varies.

      Parametric Continuity

      A curve is

      • 0-smooth in an interval if it is continuous.
      • 1-smooth in an interval if its first derivative exists and is continuous in the interval.
      • 2-smooth in an interval if its first and second derivatives exist and are continuous in the interval.
      • 3-smooth in an interval if its first, second and third derivatives exist and are continuous in the interval.

    9. surface

      requires two parameters

        x = x ( u, v )

        y = y ( u, v )

        z = z ( u, v )

    10. Parametric Polynomial Curves

      Parametric Polynomial Surfaces

      Why

    11. uses much less memory; just needs to save a few control points rather than the whole image

    12. local control of shape

    13. smoothness and continuity

    14. ability to evaluate derivatives

    15. stability

    16. ease of rendering
    17. Parametric Cubic Polynomial Curves

      p(u) =
      3

      k = 0
      ck uk = c0 + c1 u + c2 u2 + c3 u3 = uTc

      where

      c = c0
      c1
      c2
      c3
         
      u = 1
      u
      u2
      u3
         
      ck = ( ckx cky ckz )

      The design of a particular type of cubic will be based on data given at some values of parameter.

      Need

    18. control points

    19. join points
    20. Interpolation

      Interpolation polynomial:

      p(x) = anxn + an-1xn-1 + ⋅ ⋅ ⋅ + a1x1 + a0.

      If p(x) interpolates the data points ( xi, yi ), then

      p(xi) = yi     for all i ∈ {0, 1, ..., n}

      i.e. the polynomial curve passes through the data points.
      The coefficients ai can be found using Vandermonde_matrix

       

      Cubic Interpolating Polynomial

    21. Suppose we have four control points, P0, P1, P2, P3 where
        Pk   =   xk
        yk
        zk

    22. Find coefficients so that polynomail p(u) passes through ( interpolates ) the four control points.

    23. For simplicity, consider when u = 0, 1/3, 2/3, 1 p(u) passes through the four control points.

    24. That is

      Or in matrix form

        P = AC
      where
        P   =   P0
        P1
        P2
        P3
         
        C   =   C0
        C1
        C2
        C3

      and

        and
        A-1 =
         

      The desired coefficients are found by

        C = A-1P

      ( Note in our derivation, Pi can be substituted by xi or yi or zi ). For example,

    25. Lagrange's Method

      A simplest form of numerical solution for polynomial interpolation.

      /*  polyint.cpp			
       *  This program demonstrates using polynomial interpretation to draw a curve
       *  using Lagrange's method.
       */
      #include <GL/glut.h>
      #include <stdlib.h>
      #include <stdio.h>
      
      using namespace std;
      
      GLfloat ctrlpoints[4][3] = {
      	{ -4.0, -4.0, 0.0}, { -2.0, 4.0, 0.0}, 
      	{2.0, -4.0, 0.0}, {4.0, 4.0, 0.0}};
      
      void init(void)
      {
         glClearColor(0.0, 0.0, 0.0, 0.0);
         glShadeModel(GL_FLAT);
      }
      
      
      //polynomial interpretation for N points
      float polyint ( float  points[][3], float x, int N )
      {
        float y;
      
        float num = 1.0, den = 1.0;
        float sum = 0.0;
      
        for ( int i = 0; i < N; ++i ) {
          num = den = 1.0;
          for ( int j = 0; j < N; ++j ) {
            if ( j == i ) continue;
      
            num = num * ( x - points[j][0] );		 	//x - xj
          }
          for ( int j = 0; j < N; ++j ) {
            if ( j == i ) continue;
            den = den * ( points[i][0] - points[j][0] );	//xi - xj
          }
          sum += num / den * points[i][1];
        }
        y = sum;
      
        return y;
      }
      
      void display(void)
      {
         int i;
         float x, y;
      
         glClear(GL_COLOR_BUFFER_BIT);
         glColor3f(1.0, 1.0, 1.0);
         glBegin(GL_LINE_STRIP);
            for (i = -40; i <= 40; i++) {
      	 x = (float) i /10.0;
               y = polyint( ctrlpoints, x, 4);
      	 glVertex2f ( x, y );
            }
         glEnd();
         /* The following code displays the control points as dots. */
         glPointSize(5.0);
         glColor3f(1.0, 1.0, 0.0);
         glBegin(GL_POINTS);
            for (i = 0; i < 4; i++) 
               glVertex3fv(&ctrlpoints[i][0]);
         glEnd();
         glFlush();
      }
      
      void reshape(int w, int h)
      {
         glViewport(0, 0, (GLsizei) w, (GLsizei) h);
         glMatrixMode(GL_PROJECTION);
         glLoadIdentity();
         if (w <= h)
            glOrtho(-5.0, 5.0, -5.0*(GLfloat)h/(GLfloat)w, 
                     5.0*(GLfloat)h/(GLfloat)w, -5.0, 5.0);
         else
            glOrtho(-5.0*(GLfloat)w/(GLfloat)h, 
                     5.0*(GLfloat)w/(GLfloat)h, -5.0, 5.0, -5.0, 5.0);
         glMatrixMode(GL_MODELVIEW);
         glLoadIdentity();
      }
      
      void keyboard(unsigned char key, int x, int y)
      {
         switch (key) {
            case 27:
               exit(0);
               break;
         }
      }
      
      int main(int argc, char** argv)
      {
         glutInit(&argc, argv);
         glutInitDisplayMode (GLUT_SINGLE | GLUT_RGB);
         glutInitWindowSize (500, 500);
         glutInitWindowPosition (100, 100);
         glutCreateWindow (argv[0]);
         init ();
         glutDisplayFunc(display);
         glutReshapeFunc(reshape);
         glutKeyboardFunc (keyboard);
         glutMainLoop();
         return 0;
      }
        

      Neville's Algorithm and Barycentric Formula

      More practical, calculate values recursively, can be easily extended to calculate derivative values.

      Barycentric Formula extends the idea further. Let

      Equation (1) can be rewritten as:

      Hermite Interpolation

      1. Consider both derivatives at data points as well as data points themselves; closely related to "Newton divided difference" method.


        Curves defined by two end points and the tangent vectors.

      2. The tangents may be determined by other control points.

      3. Hermite method may not be very effective as sometimes it requires a large variation of the magnitude of the tangent vector in order to perturb a small change of the curve.

    26. Bezier Curve

    27. A spline curve is a smooth curve specified succinctly in terms of a few points.

    28. Bezier curves and B-spline curves are two main classes of splines

    29. Bezier curves were first developed by automobile designers to describe the shape of exterior car panels in the 1960s and 70s.

    30. Degree Three Bezier Curves

      • most commonly used

      • degree three polynomial curves specified by 4 control points with a curve passing through first and last point

      • Defined parametrically by a function p(u): as u varies from 0 to 1, the values of p(u) sweep out the curve

          p(u) = B0 ( u )P0 + B1 ( u ) P1 + B2 ( u ) P2 + B3 ( u ) P3

        where

        So
          B0(u) = (1 -u)3   B2(u) = 3u2(1 -u)
          B1(u) = 3u(1 -u)2   B3(u) = u3
        Note

        Therefore, p(u) is indeed a valid point.

        Also,

          p(0) = P0,   p(1) = P3

        We can also express p(u) as

          p(u) = (1 - 3u + 3u2 - u3)P0 + (3u - 6u2 + 3u3)P1 + (3u2 - 3u3)P2 + u3P3

        Or in matrix form

          a)
          p(u) = ( 1, u, u2, u3 ) 1
          -3
          3
          -1
          0
          3
          -6
          3
          0
          0
          3
          -3
          0
          0
          0
          1
          P0
          P1
          P2
          P3

          b)

    31. Derivative of p(u) above is given by

        p'(u) = B'0(u)P0 + B'1(u)P1 + B'2(u)P2 + B'3(u)P3

      Thus

        p'(0) = 3(P1 - P0 )       p'(1) = 3(P3 - P2 )

      This means that the curve p(u) starts at u = 0 traveling int the direction of the vector from P0 to P1 and at the end, it travels in the direction of P2 to P3

    32. May join together multiple cubic curves to approximate a curve fit to an arbitrary number of control points
    33. Evaluators

      Evaluators

    34. provide a way to specify points on a curve or surface using the control points

    35. allow you to describe any polynomial splines or surfaces of any degree, including B-splines, NURBS ( Non-Uniform Rational B-Spline ) surfaces, Bezier curves and surfaces and Hermite splines
    36. One-dimensional Evaluators

      Example: A simple Bezier curve

      /*  bezcurve.c
       *  This program uses evaluators to draw a Bezier curve.
       */
      #include <GL/glut.h>
      #include <stdlib.h>
                                                                                      
      GLfloat ctrlpoints[4][3] = {
              { -4.0, -4.0, 0.0}, { -2.0, 4.0, 0.0},
              {2.0, -4.0, 0.0}, {4.0, 4.0, 0.0}};
                                                                                      
      void init(void)
      {
         glClearColor(0.0, 0.0, 0.0, 0.0);
         glShadeModel(GL_FLAT);
         /*   
          GL_MAP1_VERTEX_3 -- specifies that 3-dimensional control points are 
      			provided and 3-D vertices should be produced
          0.0 -- low value of parameter u
          1.0 -- high value of parameter u
          3   -- number of floating-point values to advance in the data between two
      	   consecutive control points
          4   -- order of the spline ( = degree + 1 )
         */
         glMap1f(GL_MAP1_VERTEX_3, 0.0, 1.0, 3, 4, &ctrlpoints[0][0]);
         glEnable(GL_MAP1_VERTEX_3);
      }
                                                                                      
      void display(void)
      {
         int i;
                                                                                      
         glClear(GL_COLOR_BUFFER_BIT);
         glColor3f(1.0, 1.0, 1.0);
         glBegin(GL_LINE_STRIP);
            for (i = 0; i <= 30; i++)
               glEvalCoord1f((GLfloat) i/30.0);
         glEnd();
         /* The following code displays the control points as dots. */
         glPointSize(5.0);
         glColor3f(1.0, 1.0, 0.0);
         glBegin(GL_POINTS);
            for (i = 0; i < 4; i++)
               glVertex3fv(&ctrlpoints[i][0]);
         glEnd();
         glFlush();
      }
       
      void reshape(int w, int h)
      {
         glViewport(0, 0, (GLsizei) w, (GLsizei) h);
         glMatrixMode(GL_PROJECTION);
         glLoadIdentity();
         if (w <= h)
            glOrtho(-5.0, 5.0, -5.0*(GLfloat)h/(GLfloat)w,
                     5.0*(GLfloat)h/(GLfloat)w, -5.0, 5.0);
         else
            glOrtho(-5.0*(GLfloat)w/(GLfloat)h,
                     5.0*(GLfloat)w/(GLfloat)h, -5.0, 5.0, -5.0, 5.0);
         glMatrixMode(GL_MODELVIEW);
         glLoadIdentity();
      }
                                                                                      
      void keyboard(unsigned char key, int x, int y)
      {
         switch (key) {
            case 27:
               exit(0);
               break;
         }
      }
                                                                                      
      int main(int argc, char** argv)
      {
         glutInit(&argc, argv);
         glutInitDisplayMode (GLUT_SINGLE | GLUT_RGB);
         glutInitWindowSize (500, 500);
         glutInitWindowPosition (100, 100);
         glutCreateWindow (argv[0]);
         init ();
         glutDisplayFunc(display);
         glutReshapeFunc(reshape);
         glutKeyboardFunc (keyboard);
         glutMainLoop();
         return 0;
      }
      

    37. Bezier Surface Patches

    38. Bezier patches of degree 3 are defined by 4 x 4 array of control points Pi,j
        p(u, v ) =
        3

        i = 0
        3

        j = 0
        Bi(u) Bj(v) Pi,j

    39. Bezier Surface

      Lit, shaded Bezier Surface Draw with mesh

      /*  bezsurf.c
       *  This program renders a wireframe Bezier surface,
       *  using two-dimensional evaluators.
       */
      #include <GL/glut.h>
      #include <stdlib.h>
      
      GLfloat ctrlpoints[4][4][3] = {
         {{-1.5, -1.5, 4.0}, {-0.5, -1.5, 2.0}, 
          {0.5, -1.5, -1.0}, {1.5, -1.5, 2.0}}, 
         {{-1.5, -0.5, 1.0}, {-0.5, -0.5, 3.0}, 
          {0.5, -0.5, 0.0}, {1.5, -0.5, -1.0}}, 
         {{-1.5, 0.5, 4.0}, {-0.5, 0.5, 0.0}, 
          {0.5, 0.5, 3.0}, {1.5, 0.5, 4.0}}, 
         {{-1.5, 1.5, -2.0}, {-0.5, 1.5, -2.0}, 
          {0.5, 1.5, 0.0}, {1.5, 1.5, -1.0}}
      };
      
      void display(void)
      {
         int i, j;
      
         glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
         glColor3f(1.0, 1.0, 1.0);
         glPushMatrix ();
         glRotatef(85.0, 1.0, 1.0, 1.0);
         for (j = 0; j <= 8; j++) {
            glBegin(GL_LINE_STRIP);
            for (i = 0; i <= 30; i++)
               glEvalCoord2f((GLfloat)i/30.0, (GLfloat)j/8.0);
            glEnd();
            glBegin(GL_LINE_STRIP);
            for (i = 0; i <= 30; i++)
               glEvalCoord2f((GLfloat)j/8.0, (GLfloat)i/30.0);
            glEnd();
         }
         //glEvalMesh2( GL_LINE, 0, 20, 0, 20 );
         glPopMatrix ();
         glFlush();
      }
      
      void init(void)
      {
         glClearColor (0.0, 0.0, 0.0, 0.0);
         /*   
          GL_MAP1_VERTEX_3 -- specifies that 3-dimensional control points are 
      			provided and 3-D vertices should be produced
          0   -- min u
          1   -- max u
          3   -- ustride, number of values to advance in the data between two
      	   consecutive control points ( see diagram below )
          4   -- order of u ( = degree + 1 )
          0   -- min v
          1   -- max v
          12  -- vstride ( see diagram below )
          4   -- order of v ( = degree + 1 )
         */
         glMap2f(GL_MAP2_VERTEX_3, 0, 1, 3, 4,
                 0, 1, 12, 4, &ctrlpoints[0][0][0]);
         glEnable(GL_MAP2_VERTEX_3);
         /*   
          glMapGrid2f() -- defines a grid going from u1 to u2; v1 to v2 in
      		     evenly-spaced steps
          20  --  number of u steps
          0.0 --  u1 ( starting u )
          1.0 --  u2 ( ending u )
          20  --  number of v steps
          0.0 --  v1 ( starting v )
          1.0 --  v2 ( ending v )
       
          use with glEvalMesh2()
         */
         //glMapGrid2f(20, 0.0, 1.0, 20, 0.0, 1.0); 
         glEnable(GL_DEPTH_TEST);
         glShadeModel(GL_FLAT);
      }
      
      void reshape(int w, int h)
      {
         glViewport(0, 0, (GLsizei) w, (GLsizei) h);
         glMatrixMode(GL_PROJECTION);
         glLoadIdentity();
         if (w <= h)
            glOrtho(-4.0, 4.0, -4.0*(GLfloat)h/(GLfloat)w, 
                    4.0*(GLfloat)h/(GLfloat)w, -4.0, 4.0);
         else
            glOrtho(-4.0*(GLfloat)w/(GLfloat)h, 
                    4.0*(GLfloat)w/(GLfloat)h, -4.0, 4.0, -4.0, 4.0);
         glMatrixMode(GL_MODELVIEW);
         glLoadIdentity();
      }
      
      void keyboard(unsigned char key, int x, int y)
      {
         switch (key) {
            case 27:
               exit(0);
               break;
         }
      }
      
      int main(int argc, char** argv)
      {
         glutInit(&argc, argv);
         glutInitDisplayMode (GLUT_SINGLE | GLUT_RGB | GLUT_DEPTH);
         glutInitWindowSize (500, 500);
         glutInitWindowPosition (100, 100);
         glutCreateWindow (argv[0]);
         init ();
         glutDisplayFunc(display);
         glutReshapeFunc(reshape);
         glutKeyboardFunc(keyboard);
         glutMainLoop();
         return 0;
      }
      

      Each control point is 3-D ( has 3 values ). So ustride = 3, vstride = 3 x 4 = 12
       u →
      v ↓
      . . . .
      . . . .
      . . . .
      . . . .

      
      /*  bezmesh.c
       *  This program renders a lighted, filled Bezier surface,
       *  using two-dimensional evaluators.
       */
      #include <stdlib.h>
      #include <GL/glut.h>
      
      GLfloat ctrlpoints[4][4][3] = {
         { {-1.5, -1.5, 4.0},
           {-0.5, -1.5, 2.0},
           {0.5, -1.5, -1.0},
           {1.5, -1.5, 2.0}},
         { {-1.5, -0.5, 1.0},
           {-0.5, -0.5, 3.0},
           {0.5, -0.5, 0.0},
           {1.5, -0.5, -1.0}},
         { {-1.5, 0.5, 4.0},
           {-0.5, 0.5, 0.0},
           {0.5, 0.5, 3.0},
           {1.5, 0.5, 4.0}},
         { {-1.5, 1.5, -2.0},
           {-0.5, 1.5, -2.0},
           {0.5, 1.5, 0.0},
           {1.5, 1.5, -1.0}}
      };
      
      void initlights(void)
      {
         GLfloat ambient[] = {0.2, 0.2, 0.2, 1.0};
         GLfloat position[] = {0.0, 0.0, 2.0, 1.0};
         GLfloat mat_diffuse[] = {0.6, 0.6, 0.6, 1.0};
         GLfloat mat_specular[] = {1.0, 1.0, 1.0, 1.0};
         GLfloat mat_shininess[] = {50.0};
      
         glEnable(GL_LIGHTING);
         glEnable(GL_LIGHT0);
      
         glLightfv(GL_LIGHT0, GL_AMBIENT, ambient);
         glLightfv(GL_LIGHT0, GL_POSITION, position);
      
         glMaterialfv(GL_FRONT, GL_DIFFUSE, mat_diffuse);
         glMaterialfv(GL_FRONT, GL_SPECULAR, mat_specular);
         glMaterialfv(GL_FRONT, GL_SHININESS, mat_shininess);
      }
      
      void display(void)
      {
         glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
         glPushMatrix();
         glRotatef(85.0, 1.0, 1.0, 1.0);
         glEvalMesh2(GL_FILL, 0, 20, 0, 20);
         glPopMatrix();
         glFlush();
      }
      
      void init(void)
      {
         glClearColor(0.0, 0.0, 0.0, 0.0);
         glEnable(GL_DEPTH_TEST);
         glMap2f(GL_MAP2_VERTEX_3, 0, 1, 3, 4,
                 0, 1, 12, 4, &ctrlpoints[0][0][0]);
         glEnable(GL_MAP2_VERTEX_3);
         glEnable(GL_AUTO_NORMAL);
         glMapGrid2f(20, 0.0, 1.0, 20, 0.0, 1.0);
         initlights();       /* for lighted version only */
      }
      
      void reshape(int w, int h)
      {
         glViewport(0, 0, (GLsizei) w, (GLsizei) h);
         glMatrixMode(GL_PROJECTION);
         glLoadIdentity();
         if (w <= h)
            glOrtho(-4.0, 4.0, -4.0*(GLfloat)h/(GLfloat)w,
                    4.0*(GLfloat)h/(GLfloat)w, -4.0, 4.0);
         else
            glOrtho(-4.0*(GLfloat)w/(GLfloat)h,
                    4.0*(GLfloat)w/(GLfloat)h, -4.0, 4.0, -4.0, 4.0);
         glMatrixMode(GL_MODELVIEW);
         glLoadIdentity();
      }
      
      void keyboard(unsigned char key, int x, int y)
      {
         switch (key) {
            case 27:
               exit(0);
               break;
         }
      }
      
      int main(int argc, char **argv)
      {
         glutInit(&argc, argv);
         glutInitDisplayMode(GLUT_SINGLE | GLUT_RGB | GLUT_DEPTH);
         glutInitWindowSize (500, 500);
         glutInitWindowPosition (100, 100);
         glutCreateWindow(argv[0]);
         init();
         glutReshapeFunc(reshape);
         glutDisplayFunc(display);
         glutKeyboardFunc(keyboard);
         glutMainLoop();
         return 0;
      }
      

    40. B-Splines

      Powerful tool for generating curves with many control points. B stands for "Basis" Advantages

    41. A single B-spline can specify a long complicated curve

    42. B-splines can be designed with sharp bends and even "corners"

    43. Can translate piecewise Bezier curves into B-splines

    44. Act more flexibly and intuitively with a large number of control points

    45. B-splines are a lot more sensitive to the placement of control points than Bezier cuves
    46. Uniform B-Splines of Degree Three

    47. one of the simplest and most useful cases of B-splines

    48. Degree 3 B-Spline with n + 1 control points:

        q(u) = n

        i=0
        Ni(u).Pi   3 ≤ u ≤ n + 1

      where P0, ..., Pn are the control points, N0(u), ..., Nn(u) are the blending ( or basis ) functions

    49. For degree 3,
        Ni(u) = 0   if u ≤ i or i + 4 ≤ u

      So
        q(u) = j

        i=j-3
        Ni(u).Pi   u ∈ [j, j+1],   3 ≤ j ≤ n

      When a single control point pi is moved, only the portion of the curve q(u) with i < u < i + 4 is changed --> local control

      The blending functions have the following properties:

      • a) they are translates of each other, i.e.
          Ni(u) = N0(u-i)

      • b) they are piecewise degree three polyonmials

      • c) C2-continuous, i.e. Ni(u) have continuous second derivatives

      • d) partition of unity
          Ni(u) = 1
        for 3 ≤ u ≤ n + 1. This is necessary as q(u) is a valid point.

      • e) Ni ≥ 0, ( thus Ni ≤ 1, ) for all u

      • f)local properties: Ni = 0 for u ≤ i and i + 4 ≤ u
         
    50. Non Uniform Rational B-Splines ( NURB )

      Degree m B-Splines

    51. B-spline curve with n control points and n blending functions Nk,m(t) of order m:

      where P0, ..., Pn-1 are the control points. ( for uniform B-splines, u0 = 0, u1 = 1, ... )

    52. U = ( u0, u1, ...., un+m-1 ) is a knot vector.

      The B-spline blending function can be calculated recursively:

      Note that the functions sum up to 1 at any u and thus its legitimate to use them to form combination of points.

      Uniform knot vector: e.g. ( 0, 1, 2, 3, ...., n+m-1)

      Standard knot vector for n control points and order-m B-splines ( m ≤ n ):

      1. Totally n + m knots, denoted as u0, ..., un+m-1
      2. First m knots, u0, ..., um-1 all have value 0
      3. Knots um, ..., un-1 increases in increments of value 1, from 1 to n - m.
      4. The final m knots, un, ..., un+m-1 are all equal to n - m + 1.

      Example

      Knot vector of 8 control points ( n = 8 ), order 4 ( m = 4 ):

      u0 = u1 = u2 =u3 = 0

      u4 = 1
      u5 = 2
      u6 = 3
      u7 = 4

      u8 = u9 = u10 =u11 = 5

      The number of times a knot value is duplicated is called the knot's multiplicity ( ≤ order of the spline ).

      One can show that

      Nk,n+1(u) = Bnk(u) ( Bernstein coefficients ) i.e. Bezier curves are special cases of B-splines.

      /*
       * Build standard knot vector for n control points
       * and B-splines of order m
       */
      void buildKnots ( int m, int n, float knot[] )
      {
        if ( n < m  ) return;         //not enough control points
        for ( int i = 0; i < n + m; ++i ){
          if (i < m) knot[i] = 0.0;
          else if (i < n) knot[i] = i - m + 1;        //i is at least m here
          else knot[i] = n - m + 1;
        }
      }
      
      //evaluate blending functions recurvsively
      float bSpline ( int k, int m, float u, float knot[] )
      {
        float d1, d2, sum = 0.0;
      
        if ( m == 1 )
          return ( knot[k] < u &&  u <= knot[k+1] );   //1 or 0
      
        //m larger than 1, so recurse
         d1 = knot[k+m-1] - knot[k];
        if ( d1 != 0 )
          sum = (u - knot[k]) * bSpline(k,m-1,u, knot) / d1;
        d2 = knot[k+m] - knot[k+1];
        if ( d2 != 0 )
          sum += (knot[k+m] - u) * bSpline(k+1, m-1, u, knot) / d2;
      
        return sum;
      }
      
      //non uniform rational B-splines, n control points, order m, p[] is the output point 
      void nurb ( float control_points[][3], float u, float p[] )
      {
        const int n = 4, m = 4;               //this reduces to cubic B-spline
        float knot[n+m];
        buildKnots ( m,  n, knot );           //construct the knob
      
        // sum the control points mulitplied by their respective blending functions
        for ( int i = 0; i < 3; ++i ){        //x, y, z components
          p[i] = 0;
          for ( int k = 0; k < n; ++k )
            p[i] += bSpline(k, m, u, knot ) * control_points[k][i];
        }
      }

    53. A simple NURBS Example
      /*
       *  surface.c
       *  This program draws a NURBS surface in the shape of a 
       *  symmetrical hill.  The 'c' keyboard key allows you to 
       *  toggle the visibility of the control points themselves.  
       *  Note that some of the control points are hidden by the  
       *  surface itself.
       */
      #include <GL/glut.h>
      #include <stdlib.h>
      #include <stdio.h>
      
      #ifndef CALLBACK
      #define CALLBACK
      #endif
      
      GLfloat ctlpoints[4][4][3];
      int showPoints = 0;
      
      GLUnurbsObj *theNurb;
      GLenum errorCode;
      
      /*
       *  Initializes the control points of the surface to a small hill.
       *  The control points range from -3 to +3 in x, y, and z
       */
      void init_surface(void)
      {
         int u, v;
         for (u = 0; u < 4; u++) {
            for (v = 0; v < 4; v++) {
               ctlpoints[u][v][0] = 2.0*((GLfloat)u - 1.5);
               ctlpoints[u][v][1] = 2.0*((GLfloat)v - 1.5);
      
               if ( (u == 1 || u == 2) && (v == 1 || v == 2))
                  ctlpoints[u][v][2] = 3.0;
               else
                  ctlpoints[u][v][2] = -3.0;
            }
         }				
      }				
      
      void CALLBACK nurbsError(GLenum errorCode)
      {
         const GLubyte *estring;
      
         estring = gluErrorString(errorCode);
         fprintf (stderr, "Nurbs Error: %s\n", estring);
         exit (0);
      }
      			
      /*  Initialize material property and depth buffer.
       */
      void init(void)
      {
         GLfloat mat_diffuse[] = { 0.7, 0.7, 0.7, 1.0 };
         GLfloat mat_specular[] = { 1.0, 1.0, 1.0, 1.0 };
         GLfloat mat_shininess[] = { 100.0 };
      
         glClearColor (0.0, 0.0, 0.0, 0.0);
         glMaterialfv(GL_FRONT, GL_DIFFUSE, mat_diffuse);
         glMaterialfv(GL_FRONT, GL_SPECULAR, mat_specular);
         glMaterialfv(GL_FRONT, GL_SHININESS, mat_shininess);
      
         glEnable(GL_LIGHTING);
         glEnable(GL_LIGHT0);
         glEnable(GL_DEPTH_TEST);
         glEnable(GL_AUTO_NORMAL);
         glEnable(GL_NORMALIZE);
      
         init_surface();
      
         theNurb = gluNewNurbsRenderer();
         gluNurbsProperty(theNurb, ( GLenum ) GLU_SAMPLING_TOLERANCE, 25.0);
         gluNurbsProperty(theNurb, ( GLenum ) GLU_DISPLAY_MODE, GLU_FILL);
         gluNurbsCallback(theNurb, ( GLenum ) GLU_ERROR,  nurbsError);
      
      }
      
      void display(void)
      {
         GLfloat knots[8] = {0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0};
         int i, j;
      
         glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
      
         glPushMatrix();
         glRotatef(330.0, 1.,0.,0.);
         glScalef (0.5, 0.5, 0.5);
      
         gluBeginSurface(theNurb);
         gluNurbsSurface(theNurb, 
                         8, knots, 8, knots,
                         4 * 3, 3, &ctlpoints[0][0][0], 
                         4, 4, GL_MAP2_VERTEX_3);
         gluEndSurface(theNurb);
      
         if (showPoints) {
            glPointSize(5.0);
            glDisable(GL_LIGHTING);
            glColor3f(1.0, 1.0, 0.0);
            glBegin(GL_POINTS);
            for (i = 0; i < 4; i++) {
               for (j = 0; j < 4; j++) {
      	    glVertex3f(ctlpoints[i][j][0], 
                     ctlpoints[i][j][1], ctlpoints[i][j][2]);
               }
            }
            glEnd();
            glEnable(GL_LIGHTING);
         }
         glPopMatrix();
         glFlush();
      }
      
      void reshape(int w, int h)
      {
         glViewport(0, 0, (GLsizei) w, (GLsizei) h);
         glMatrixMode(GL_PROJECTION);
         glLoadIdentity();
         gluPerspective (45.0, (GLdouble)w/(GLdouble)h, 3.0, 8.0);
         glMatrixMode(GL_MODELVIEW);
         glLoadIdentity();
         glTranslatef (0.0, 0.0, -5.0);
      }
      
      void keyboard(unsigned char key, int x, int y)
      {
         switch (key) {
            case 'c':
            case 'C':
               showPoints = !showPoints;
               glutPostRedisplay();
               break;
            case 27:
               exit(0);
               break;
            default:
               break;
         }
      }
      
      int main(int argc, char** argv)
      {
         glutInit(&argc, argv);
         glutInitDisplayMode(GLUT_SINGLE | GLUT_RGB | GLUT_DEPTH);
         glutInitWindowSize (500, 500);
         glutInitWindowPosition (100, 100);
         glutCreateWindow(argv[0]);
         init();
         glutReshapeFunc(reshape);
         glutDisplayFunc(display);
         glutKeyboardFunc (keyboard);
         glutMainLoop();
         return 0; 
      }
        

    54. Creating Surfaces by Recurvise Subdivision -- the Utah Teapot

      //bteapot.c
      /* Shaded teapot using Bezier surfaces */
      
      #include <stdlib.h>
      #include <GL/glut.h>
      
      typedef GLfloat point[3];
      
      point data[32][4][4];
      
      /* 306 vertices */
      
      #include "vertices.h"
      
      /*
       32 patches each defined by 16 vertices, arranged in a 4 x 4 array 
       NOTE: teapot vertices labeled from 1 to 306 
       remnent of the days of FORTRAN 
      */
      
      #include "patches.h"
      
      void display(void)
      {
         int i, j, k;
      
         glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
         glColor3f(1.0, 1.0, 1.0);
      
         glLoadIdentity();
         glTranslatef(0.0, 0.0, -10.0); 
         glRotatef(-35.26, 1.0, 0.0, 0.0);
         glRotatef(-45.0, 0.0, 1.0, 0.0);  
      
         /* gluLookAt(-1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0); */
      
         /* data aligned along z axis, rotate to align with y axis */
      
         glRotatef(-90.0, 1.0,0.0, 0.0); 
         for(k=0;k<32;k++) 
         {
           glMap2f(GL_MAP2_VERTEX_3, 0, 1, 3, 4,
                 0, 1, 12, 4, &data[k][0][0][0]);
           for (j = 0; j <= 8; j++) {
      	glBegin(GL_LINE_STRIP);
           	  for (i = 0; i <= 30; i++)
                  glEvalCoord2f((GLfloat)i/30.0, (GLfloat)j/8.0);
              glEnd();
              glBegin(GL_LINE_STRIP);
                for (i = 0; i <= 30; i++)
                   glEvalCoord2f((GLfloat)j/8.0, (GLfloat)i/30.0);
                glEnd();
           }
         }
         glFlush();
      }
      
      void myReshape(int w, int h)
      {
        glViewport(0, 0, w, h);
        glMatrixMode(GL_PROJECTION);
        glLoadIdentity();
        if (w <= h)
          glOrtho(-4.0, 4.0, -4.0 * (GLfloat) h / (GLfloat) w,
               4.0 * (GLfloat) h / (GLfloat) w, -20.0, 20.0);
        else
          glOrtho(-4.0 * (GLfloat) w / (GLfloat) h,
             	4.0 * (GLfloat) w / (GLfloat) h, -4.0, 4.0, -20.0, 20.0);
          glMatrixMode(GL_MODELVIEW);
        display();
      }
      
      void myinit()
      {
        glEnable(GL_MAP2_VERTEX_3);
        glMapGrid2f(20, 0.0, 1.0, 20, 0.0, 1.0);
        glEnable(GL_DEPTH_TEST);
      }
      
      main(int argc, char *argv[])
      {
        int i,j, k, m, n;
      
        for(i=0;i<32;i++) 
          for(j=0;j<4;j++) 
            for(k=0;k<4;k++) 
      	for(n=0;n<3;n++){
      
      	/* put teapot data into single array for subdivision */
      
                m=indices[i][j][k];
                for(n=0;n<3;n++) data[i][j][k][n]=vertices[m-1][n];
           	}
          glutInit(&argc, argv);
          glutInitDisplayMode(GLUT_SINGLE | GLUT_RGB | GLUT_DEPTH);
          glutInitWindowSize(500, 500);
          glutCreateWindow("teapot");
          myinit();
          glutReshapeFunc(myReshape);
          glutDisplayFunc(display);
      
          glutMainLoop();
      }
      

      //teapot.c
      /* Shaded teapot using recursive subdivision of Bezier surfaces */
      /* Depth of recursion entered as command line argument */
      
      #include <stdlib.h>
      #include <GL/glut.h>
      
      typedef GLfloat point[3];
      
      void draw_patch(point p[4][4]);
      void divide_curve(point c[4], point *r, point *l);
      void divide_patch(point p[4][4], int n);
      void transpose(point p[4][4]);
      void normal(point n, point p, point q, point r);
      
      point data[32][4][4];
      
      /* 306 vertices */
      
      #include "vertices.h"
      
      /* 32 patches each defined by 16 vertices, arranged in a 4 x 4 array */
      /* NOTE: numbering scheme for teapot has vertices labeled from 1 to 306 */
      /* remnent of the days of FORTRAN */
      
      #include "patches.h"
      
      int level; 
      
      void normal(point n, point p, point q, point r)
      {
        n[0]=(q[1]-p[1])*(r[2]-p[2])-(q[2]-p[2])*(r[1]-p[1]);
        n[1]=(q[2]-p[2])*(r[0]-p[0])-(q[0]-p[0])*(r[2]-p[2]);
        n[2]=(q[0]-p[0])*(r[2]-p[2])-(q[2]-p[2])*(r[0]-p[0]);
      }
      
      void transpose(point a[4][4])
      {
      
        /* transpose wastes time but makes program more readable */
      
        int i,j, k;
        GLfloat tt;
        for(i=0;i<4;i++) for(j=i;j<4; j++) for(k=0;k<3;k++)
        {
              tt=a[i][j][k];
              a[i][j][k]=a[j][i][k];
              a[j][i][k]=tt;
        }
      }
      
      void divide_patch(point p[4][4], int n)
      {
         point q[4][4], r[4][4], s[4][4], t[4][4];
         point a[4][4], b[4][4];
         int i,j, k;
         if(n==0) draw_patch(p); /* draw patch if recursion done */
      
      /* subdivide curves in u direction, transpose results, divide
      in u direction again (equivalent to subdivision in v) */
      
         else {
             for(k=0; k<4; k++) divide_curve(p[k], a[k], b[k]);
             transpose(a);
             transpose(b); 
             for(k=0; k<4; k++) {
                 divide_curve(a[k], q[k], r[k]);
                 divide_curve(b[k], s[k], t[k]);
             }
      
      /* recursive division of 4 resulting patches */
      
             divide_patch(q, n-1);
             divide_patch(r, n-1);
             divide_patch(s, n-1);
             divide_patch(t, n-1);
         }
      }
      
      void divide_curve(point c[4], point r[4], point l[4])
      {
      
      /* division of convex hull of Bezier curve */
      
         int i;
         point t;
         for(i=0;i<3;i++)
         {
             l[0][i]=c[0][i];
             r[3][i]=c[3][i];
             l[1][i]=(c[1][i]+c[0][i])/2;
             r[2][i]=(c[2][i]+c[3][i])/2;
             t[i]=(l[1][i]+r[2][i])/2;
             l[2][i]=(t[i]+l[1][i])/2;
             r[1][i]=(t[i]+r[2][i])/2;
             l[3][i]=r[0][i]=(l[2][i]+r[1][i])/2;
         }
      }
      
      void draw_patch(point p[4][4])
      {
          point n;
          normal(n, p[0][0], p[3][0], p[3][3]);
          glBegin(GL_QUADS);
          glNormal3fv(n);
          glVertex3fv(p[0][0]);
          glVertex3fv(p[3][0]);
          glVertex3fv(p[3][3]);
          glVertex3fv(p[0][3]);
          glEnd();
      }
      
      void display(void)
      {
          int i;
      
          glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
          glLoadIdentity();
      
      /* data aligned along z axis, rotate to align with y axis */
      
          glRotatef(-90.0, 1.0,0.0, 0.0);
          for(i=0;i<32;i++) divide_patch(data[i], level); /* divide all 32 patches */
      
          glFlush();
      
      }
      
      
      void myReshape(int w, int h)
      {
          glViewport(0, 0, w, h);
          glMatrixMode(GL_PROJECTION);
          glLoadIdentity();
          if (w <= h)
              glOrtho(-4.0, 4.0, -4.0 * (GLfloat) h / (GLfloat) w,
                  4.0 * (GLfloat) h / (GLfloat) w, -10.0, 10.0);
          else
              glOrtho(-4.0 * (GLfloat) w / (GLfloat) h,
                  4.0 * (GLfloat) w / (GLfloat) h, -4.0, 4.0, -10.0, 10.0);
          glMatrixMode(GL_MODELVIEW);
          display();
      }
      
      void myinit()
      {
          GLfloat mat_specular[]={1.0, 1.0, 1.0, 1.0};
          GLfloat mat_diffuse[]={1.0, 1.0, 1.0, 1.0};
          GLfloat mat_ambient[]={1.0, 1.0, 1.0, 1.0};
          GLfloat mat_shininess={100.0};
          GLfloat light_ambient[]={0.0, 0.0, 0.0, 1.0};
          GLfloat light_diffuse[]={1.0, 1.0, 1.0, 1.0};
          GLfloat light_specular[]={1.0, 1.0, 1.0, 1.0};
          GLfloat light_position[]={10.0, 10.0, 10.0, 0.0};
      
          glLightfv(GL_LIGHT0, GL_POSITION, light_position);
          glLightfv(GL_LIGHT0, GL_AMBIENT, light_ambient);
          glLightfv(GL_LIGHT0, GL_DIFFUSE, light_diffuse);
          glLightfv(GL_LIGHT0, GL_SPECULAR, light_specular);
      
          glMaterialfv(GL_FRONT, GL_SPECULAR, mat_specular);
          glMaterialfv(GL_FRONT, GL_AMBIENT, mat_ambient);
          glMaterialfv(GL_FRONT, GL_DIFFUSE, mat_diffuse);
          glMaterialf(GL_FRONT, GL_SHININESS, mat_shininess);
      
          glShadeModel(GL_SMOOTH);
          glEnable(GL_LIGHTING);
          glEnable(GL_LIGHT0);
      
          glEnable(GL_DEPTH_TEST);
          glClearColor (0.0, 0.0, 0.0, 1.0);
          glColor3f (1.0, 1.0, 1.0);
          glEnable(GL_NORMALIZE); /* automatic normaization of normals */
          glEnable(GL_CULL_FACE); /* eliminate backfacing polygons */
          glCullFace(GL_BACK); 
      
      }
      
      main(int argc, char *argv[])
      {
        int i,j, k, m, n;
      
      
          level=4;
      
          for(i=0;i<32;i++) for(j=0;j<4;j++) for(k=0;k<4;k++) for(n=0;n<3;n++)
          {
      
      /* put teapot data into single array for subdivision */
      
              m=indices[i][j][k];
              for(n=0;n<3;n++) data[i][j][k][n]=vertices[m-1][n];
          }
          glutInit(&argc, argv);
          glutInitDisplayMode(GLUT_SINGLE | GLUT_RGB | GLUT_DEPTH);
          glutInitWindowSize(500, 500);
          glutCreateWindow("teapot");
          myinit();
          glutReshapeFunc(myReshape);
          glutDisplayFunc(display);
      
          glutMainLoop();
      }