OpenGL Programming/Scientific OpenGL Tutorial 04

A three-dimensional plot


Now that we have a grip on drawing two-dimensional plots, it is time to tackle three-dimensional plots. It is not too different from drawing two-dimensional plots, it is just a matter of adding the third dimension and choosing a suitable model-view-projection (MVP) matrix to present the three-dimensional data in a clear way.

One significant difference however is that we now have much more data to plot, as the number of data points is now squared. The strategy from the second graph tutorial to separate the number of vertices drawn from the number of data points will now really pay off, so we will use it in this tutorial as well.

We will also see that when drawing grid lines, vertices are used more than once. To ensure we reuse vertices, and to work around some other problems, we will use an Index Buffer Object (IBO).

Putting our function in a textureEdit

We will use a 3D version of the Mexican hat function. Basically, this is a function of only one variable, but we will make a rotation symmetrical version of it by using the distance to the origin as that variable:


We will evaluate this function in a grid of N by N points. We will make it so it is easy to change the exact number of points at compile time:

#define N 256
GLbyte graph[N][N];

for(int i = 0; i < N; i++) {
  for(int j = 0; j < N; j++) {
    float x = (i - N / 2) / (N / 2.0);
    float y = (j - N / 2) / (N / 2.0);
    float t = hypotf(x, y) * 4.0;
    float z = (1 - t * t) * expf(t * t / -2.0);
    graph[i][j] = roundf(z * 127 + 128);

The hypot*() functions are not well known but they are part of the C99 standard, and are very convenient. We scaled the function somewhat to ensure the hat fits nicely within the range -1..1. We can now tell OpenGL to use this data as a two dimensional texture:

glGenTextures(1, &texture_id);
glBindTexture(GL_TEXTURE_2D, texture_id);


  • Try different values of N. What is the maximum size supported by your video card? How much memory does the texture consume?
  • We don't need to evaluate all the N x N points in one go. Make it so you evaluate only N points at a time, and put the partial results in the texture with glTexSubImage2D().
  • Try using an image instead of a function (you can use the texture from Modern OpenGL Tutorial 06 for example.
  • Make it so you can switch texture interpolation and wrapping on and off with the F1 and F2 keys.


The vertex shader we will use works more or less the same as the one from the second graph tutorial. However, there we had the luxury of plotting a 2D function on a 2D screen, so we didn't have to transform our vertex coordinates, we just moved the texture around a little. When plotting a 3D function, we really have to project our vertices somehow to make sense of all the three dimensions. Also, we can move the texture around in two dimensions now. So it makes sense two just have two generic transformation matrices; one for the texture coordinates and one for the vertex coordinates. This is what our shader will look like:

attribute vec2 coord2d;
uniform mat4 texture_transform;
uniform mat4 vertex_transform;
uniform sampler2D mytexture;
varying vec4 graph_coord;

void main(void) {
  graph_coord = texture_transform * vec4(coord2d, 0, 1);
  graph_coord.z = texture2D(mytexture, graph_coord.xy / 2.0 + 0.5).r;
  gl_Position = vertex_transform * vec4(coord2d, graph_coord.z, 1);

The attribute coord2d has the same function as the attribute coord1d from the second graph tutorial. The uniform matrix texture_transform takes over the role of the uniforms offset_x and scale_x. The uniform matrix vertex_transform is new, and will be used to change our view on the graph. In the main function, we recover the graph coordinates by applying the texture_transform matrix to the 2D coordinates that we feed it. Once we know that, we can recover the z coordinate by doing a texture lookup with those coordinates. We keep the graph coordinates in a varying vec4 so they can be used by the fragment shader to give the graph nice colors. The gl_Position variable is calculated the same was as it was done in the second tutorial, except for the new transformation matrix that is applied.

We will use the following fragment shader:

varying vec4 graph_coord;

void main(void) {
  gl_FragColor = graph_coord / 2.0 + 0.5;

Calculating the texture and vertex transformation matricesEdit

If you have followed the previous tutorials, you have already seen how we created a transformation matrix from offset and scale variables. This time is no different, except we have two offset variables, for both the x and y axis:

glm::mat4 texture_transform = glm::translate(glm::scale(glm::mat4(1.0f), glm::vec3(scale, scale, 1)), glm::vec3(offset_x, offset_y, 0));
glUniformMatrix4fv(uniform_texture_transform, 1, GL_FALSE, glm::value_ptr(texture_transform));

Also, from Modern OpenGL Tutorial 05 we have seen how to create model, view and projection matrices. Let's keep our vertices in a box from (-1, -1, -1) to (1, 1, 1), so our model transformation matrix is just the identity matrix. Then we can position the camera at say, (0, -2, 2), where the vector (0, 0, 1) is the up direction. So we are a bit in front of and above the graph, looking down on it. Finally, we will use the same perspective projection as used in the other tutorials. The resulting MVP matrix is calculated as follows:

glm::mat4 model = glm::mat4(1.0f);
glm::mat4 view = glm::lookAt(glm::vec3(0.0, -2.0, 2.0), glm::vec3(0.0, 0.0, 0.0), glm::vec3(0.0, 0.0, 1.0));
glm::mat4 projection = glm::perspective(45.0f, 1.0f * 640 / 480, 0.1f, 10.0f);

glm::mat4 vertex_transform = projection * view * model;
glUniformMatrix4fv(uniform_vertex_transform, 1, GL_FALSE, glm::value_ptr(vertex_transform));


  • Make it so you can change offset_y using the up and down keys, and change the scale with page-up and page-down.
  • Change the model matrix such that the graph is slowly rotating around the z axis in time.
  • Make it so you can toggle rotation by pressing the F3 key.

Drawing a gridEdit

One way to plot a three-dimensional function is to draw grid lines. This is exactly what gnuplot does by default when you use the splot command, so we will do that as well. We used 101 points before for the two-dimensional plot, so we will use a grid of 101 by 101 points now, and put that in our VBO:

struct point {
  GLfloat x;
  GLfloat y;

point vertices[101][101];

for(int i = 0; i < 101; i++) {
  for(int j = 0; j < 101; j++) {
    vertices[i][j].x = (j - 50) / 50.0;
    vertices[i][j].y = (i - 50) / 50.0;

glBindBuffer(GL_ARRAY_BUFFER, vbo);
glBufferData(GL_ARRAY_BUFFER, sizeof vertices, vertices, GL_STATIC_DRAW);

The vertices are in the right order to draw horizontal lines (i.e., those with constant y coordinates). The only problem is that we cannot just make a single call to glDrawArrays(GL_LINE_STRIP), as it would not know when it has reached the right edge of the graph and moves back to the left edge. Instead, it would create a zig-zag line pattern. The most simple solution is to just manually draw 101 lines:

glVertexAttribPointer(attribute_coord2d, 2, GL_FLOAT, GL_FALSE, 0, 0);
for(int i = 0; i < 101; i++)
  glDrawArrays(GL_LINE_STRIP, 101 * i, 101);

That works, although we did make 101 OpenGL calls, which is not so much, but one would rather avoid doing that. We also need to draw vertical lines, but the vertices are not in the right order! Although, in this case, we can cheat by using the stride and pointer parameters:

for(int i = 0; i < 101; i++) {
  glVertexAttribPointer(attribute_coord2d, 2, GL_FLOAT, GL_FALSE, 101 * sizeof(point), (void *)(i * sizeof(point)));
  glDrawArrays(GL_LINE_STRIP, 0, 101);


  • Why did we have to use an offset pointer in glVertexAttribPointer()? Couldn't we have left it 0 and used glDrawArrays(GL_LINE_STRIP, i, 101) instead?
  • Think of a way to create a vertex array with duplicate vertices such that you can draw all horizontal and vertical lines with just one call to glDrawArrays(GL_LINE_STRIP).
  • If you are not limited by OpenGL ES, have a look at the glMultiDrawArrays() command.
  • Suppose you want to draw triangles to fill the whole 101 x 101 square. Can you order the vertices such that you can draw everything with one call to glDrawArrays(GL_TRIANGLE_STRIP) without wasting vertices? What about when using multiple calls to glDrawArrays()?

Using an IBO to prevent wasting vertices and OpenGL callsEdit

If we wanted to draw the graph with triangles, to get a filled surface instead of a grid, we could not reuse our VBO anymore, since the order of the vertices for grid lines is completely different than for triangles. If we wanted to draw both a filled surface and grid lines on top, we would need multiple copies for all of our vertices, just because of ordering problems, and we would need many glDrawArrays() commands.

Luckily there is a way to decouple a set of vertices from the order in which to draw them. With the glDrawElements() function, we can have a second array which contains indices to the vertex array (or any other attribute array for that matter). Unfortunately, there is still no way to draw with GL_LINE_STRIP, because the index array also cannot tell OpenGL where to start and end the line strips. But we can draw with GL_LINES! Now, you may think that then we have a lot of duplication again, because we would have to draw a line from vertex index 0 to 1, from 1 to 2, and so on. However, indices are small numbers, usually one or two bytes big, while attributes are usually much bigger. Even in our simple case, our 2D vertex attributes are already 8 bytes big. So the overhead is much smaller. The advantages are that we can draw all the line segments for both horizontal and vertical grid lines with one call to glDrawElements(), without any unnecessary pixels drawn. Of course, we can also store the indices in the memory of the GPU, by using Index Buffer Objects. Here goes:

GLushort indices[2 * 100 * 101 * 2];
int i = 0;

// Horizontal grid lines
for(int y = 0; y < 101; y++) {
  for(int x = 0; x < 100; x++) {
    indices[i++] = y * 101 + x;
    indices[i++] = y * 101 + x + 1;

// Vertical grid lines
for(int x = 0; x < 101; x++) {
  for(int y = 0; y < 100; y++) {
    indices[i++] = y * 101 + x;
    indices[i++] = (y + 1) * 101 + x;

GLint ibo;
glGenBuffers(1, &ibo);
glBufferData(GL_ELEMENT_ARRAY_BUFFER, sizeof indices, indices, GL_STATIC_DRAW);

The amount of indices is two per line segment, and we have 100 segments per grid line. Then, we have two times 101 grid lines. Here is how we finally draw the grid using vertices from our VBO and indices from our IBO:

glBindBuffer(GL_ARRAY_BUFFER, vbo);
glVertexAttribPointer(attribute_coord2d, 2, GL_FLOAT, GL_FALSE, 0, 0);
glDrawElements(GL_LINES, 2 * 100 * 101 * 2, GL_UNSIGNED_SHORT, 0);


  • Find out how you can draw only the vertical grid lines with a single call to glDrawElements() without changing anything else.
  • Is there a limit to the number of vertices that can be drawn with glDrawElements()?
  • Create an array of indices to draw a filled surface with GL_TRIANGLES.
  • Can you reuse the same IBO with another VBO? Or the same VBO with another IBO?
  • Suppose you have two attribute arrays, one for vertices and one for colors. Find out how glDrawElements() works in that case.

< OpenGL Programming

Browse & download complete code