The teapot is a famous model among 3D developers.

You can find various versions of the model as vertices, but did you know that the original version was actually made of (3,3) Bézier surfaces? Bézier surface are described by control points, and we can evaluate the surface with any level of precision to produce a set of vertices.

See:

The Origins of the Teapot (PDF) - Computer Graphics & Applications - January 1987 (vol. 7 no. 1) - a venerable article, with Bézier data ("patches")!
The History of The Teapot - more historical details

So how about we create a high-definition version of the teapot? :)

Math to code

Wikipedia has related information at Bézier surface#Equation

Wikipedia has related information at Bézier_curve

From the Wikipedia article (see below for explanations):

A given Bézier surface of order (n, m) is defined by a set of (n + 1)(m + 1) control points k_i,j. [...]

A two-dimensional Bézier surface can be defined as a parametric surface where the position of a point p as a function of the parametric coordinates u, v is given by:

\mathbf {p} (u,v)=\sum _{i=0}^{n}\sum _{j=0}^{m}B_{i}^{n}(u)\;B_{j}^{m}(v)\;\mathbf {k} _{i,j}

evaluated over the unit square, where

B_{i}^{n}(u)={n \choose i}\;u^{i}(1-u)^{n-i}

is a Bernstein polynomial, and

{n \choose i}={\frac {n!}{i!(n-i)!}}

is the binomial coefficient.

OK, well actually this is pretty simple.

the big "E" means "sum" from value a to value b, step 1 (i.e.: it's a for loop that does additions)
our data is formed of 4x4 points, so our Bézier surface order is (3,3) - that's a 3D surface
we create a grid of vertices, and we index them by (u,v); u and v are between 0 and 1 (equivalent to t in Bézier curves, we can say it's a percentage of completion along an axis)

What's complicated is all the code we'll need around the maths ;)

Getting the control points

Teapot with its control points

The PDF article presents the data as a big set of control point vertices, and then several patches of 4x4 vertex indices.

We want our control points somehow like this (we'll use a vertex structure for clarity):

  struct vertex { GLfloat x, y, z; };
  ...
  #define ORDER 3
  struct vertex control_points_k[ORDER+1][ORDER+1] = {
    { { 1, 2, 3}, { 4, 5, 6}, { 7, 8, 9}, {10,11,12} },
    { {13,14,15}, {16,17,18}, {19,20,21}, {22,23,24} },
    ...
  };

Moreover, we want this array for each of the 28 patches in the teapot.

We note that the data provided in the article is not directly usable:

we don't have directly the vertices, but instead we have indices to the vertices
the indices start at 1, instead of 0 as in C/C++

We'll store the data as a C array, then we'll write a function to convert it to the format we want. We do it in a first, separate step, because doing everything as once would make our code look complicated.

The vertices:

struct vertex teapot_cp_vertices[] = {
  {  1.4   ,   0.0   ,  2.4     },
  {  1.4   ,  -0.784 ,  2.4     },
  {  0.784 ,  -1.4   ,  2.4     },
  ...

The indices:

#define TEAPOT_NB_PATCHES 28
GLushort teapot_patches[][ORDER+1][ORDER+1] = {
  // rim
  { {   1,   2,   3,   4 }, {   5,   6,   7,   8 }, {   9,  10,  11,  12 }, {  13,  14,  15,  16, } },
  { {   4,  17,  18,  19 }, {   8,  20,  21,  22 }, {  12,  23,  24,  25 }, {  16,  26,  27,  28, } },
  { {  19,  29,  30,  31 }, {  22,  32,  33,  34 }, {  25,  35,  36,  37 }, {  28,  38,  39,  40, } },
  { {  31,  41,  42,   1 }, {  34,  43,  44,   5 }, {  37,  45,  46,   9 }, {  40,  47,  48,  13, } },
  // body
  ...

Now our function to get the set of control points:

void build_control_points_k(int p, struct vertex control_points_k[][ORDER+1]) {
  for (int i = 0; i <= ORDER; i++)
    for (int j = 0; j <= ORDER; j++)
      control_points_k[i][j] = teapot_cp_vertices[teapot_patches[p][i][j] - 1];
}

Computing the vertices

Now we'll evaluate the Bézier surface, with a resolution of 10x10 (or any other precision you want).

For each 4x4 patch, we compute each point in our 10x10 grid (so with u and v progressing with 1/10 steps):

#define RESU 10
#define RESV 10
struct vertex teapot_vertices[TEAPOT_NB_PATCHES * RESU*RESV];
...
void build_teapot() {
  // Vertices
  for (int p = 0; p < TEAPOT_NB_PATCHES; p++) {
    struct vertex control_points_k[ORDER+1][ORDER+1];
    build_control_points_k(p, control_points_k);
    for (int ru = 0; ru <= RESU-1; ru++) {
      float u = 1.0 * ru / (RESU-1);
      for (int rv = 0; rv <= RESV-1; rv++) {
	float v = 1.0 * rv / (RESV-1);
	teapot_vertices[p*RESU*RESV + ru*RESV + rv] = compute_position(control_points_k, u, v);
      }
    }
  }

  // Elements
  ...
}

For a vertex at (u,v) we compute the "EE" sum from the equation ("formula") above:

struct vertex compute_position(struct vertex control_points_k[][ORDER+1], float u, float v) {
  struct vertex result = { 0.0, 0.0, 0.0 };
  for (int i = 0; i <= ORDER; i++) {
    for (int j = 0; j <= ORDER; j++) {
      float poly_i = bernstein_polynomial(i, ORDER, u);
      float poly_j = bernstein_polynomial(j, ORDER, v);
      result.x += poly_i * poly_j * control_points_k[i][j].x;
      result.y += poly_i * poly_j * control_points_k[i][j].y;
      result.z += poly_i * poly_j * control_points_k[i][j].z;
    }
  }
  return result;
}

Note: we can optimize the code by computing poly_i only once per i loop: move it between the two for lines at the beginning.

The bernstein_polynomial and binomial_coefficient functions are tedious, but straighforward:

float bernstein_polynomial(int i, int n, float u) {
  return binomial_coefficient(i, n) * powf(u, i) * powf(1-u, n-i);
}

float binomial_coefficient(int i, int n) {
  assert(i >= 0); assert(n >= 0);
  return 1.0f * factorial(n) / (factorial(i) * factorial(n-i));
}
int factorial(int n) {
  assert(n >= 0);
  int result = 1;
  for (int i = n; i > 1; i--)
    result *= i;
  return result;
}

Note : the article presents Pascal code to do the job. You may have noticed that the authors hard-coded the equation for n=m=3. We didn't use that method, because it makes the code actually less easy to compare with the equation, and does not really make the code clearer nor shorter.

To be able to declare our functions in any order, we need to pre-declare them at the top of our file:

void build_control_points_k(int p, struct vertex control_points_k[][ORDER+1]);
struct vertex compute_position(struct vertex control_points_k[][ORDER+1], float u, float v);
float bernstein_polynomial(int i, int n, float u);
float binomial_coefficient(int i, int n);
int factorial(int n);

Drawing the vertices

Now that we have our grid of vertices, we can draw each of its squares, using the elements technique :

GLushort teapot_elements[TEAPOT_NB_PATCHES * (RESU-1)*(RESV-1) * 2*3];
...
void build_teapot() {
  // Vertices
  ...

  // Elements
  int n = 0;
  for (int p = 0; p < TEAPOT_NB_PATCHES; p++)
    for (int ru = 0; ru < RESU-1; ru++)
      for (int rv = 0; rv < RESV-1; rv++) {
	// 1 square ABCD = 2 triangles ABC + CDA
	// ABC
	teapot_elements[n] = p*RESU*RESV +  ru   *RESV +  rv   ; n++;
	teapot_elements[n] = p*RESU*RESV +  ru   *RESV + (rv+1); n++;
	teapot_elements[n] = p*RESU*RESV + (ru+1)*RESV + (rv+1); n++;
	// CDA
	teapot_elements[n] = p*RESU*RESV + (ru+1)*RESV + (rv+1); n++;
	teapot_elements[n] = p*RESU*RESV + (ru+1)*RESV +  rv   ; n++;
	teapot_elements[n] = p*RESU*RESV +  ru   *RESV +  rv   ; n++;
      }
}

We can draw it as usual:

  glBindBuffer(GL_ARRAY_BUFFER, vbo_teapot_vertices);
  glVertexAttribPointer(
    attribute_coord3d, // attribute
    3,                 // number of elements per vertex, here (x,y,z)
    GL_FLOAT,          // the type of each element
    GL_FALSE,          // take our values as-is
    0,                 // no extra data between each position
    0                  // offset of first element
  );
  glBindBuffer(GL_ELEMENT_ARRAY_BUFFER, ibo_teapot_elements);
  glDrawElements(GL_TRIANGLES, sizeof(teapot_elements)/sizeof(teapot_elements[0]), GL_UNSIGNED_SHORT, 0);

We've got our flying teapot!

Making precision arbitrary

Teapot with 4x4 resolution

Currently, we can modify the resolution and the Bernstein degree by modifying the #defines.

If we want to change these values dynamically (for instance, change the resolution when the user click a +/- button), we won't be able to use static arrays anymore, due to a limitation in C/C++. We chose to use static arrays, because it makes the code easier to understand. To make the change, you'll need to either:

Create an array of pointers to arrays of floats (multi-dimensional)
Use a single-dimension array and use maths to compute the right index. For instance, in a 4x4 array, array[2][3] is equivalent to array[2*4+3] - that's just what we did for the teapot_elements[] array.

This is left as an exercise to the reader ;)

Limit

When we make the precision very high, for instance 49x49 (67228 vertices and 774144 elements), some vertices seem to merge. Remember that we use GL_UNSIGNED_SHORT to index vertices? This means we can only address up to 65536 vertices. If we want to draw more, we need to split the teapot into several arrays of vertices.

Going further

In The History of The Teapot page, you'll find a link to an archive with Bézier patches for other teaset elements, notably a spoon and a cup. Display them around the teapot!

Note: there is a built-in teapot shipped in GLUT, with a static resolution (glutSolidTeapot() - made of vertices instead of Bézier surfaces). We didn't use it in this tutorial, because it's not fun, because it's old-style/1.x OpenGL, and also because we want to use GLUT as little as possible: you may not have GLUT available for mobile development, for instance.

This tutorial does not discuss normals. They are necessary to compute lighting (feel free to contribute a new section).

- Comment on this page

- Recent stats

< OpenGL Programming

Browse & download complete code