Dmitry Kurtaev/ C++: Grid surface

Problem: drawing surface using heights map.
Example: written on C++, freeglut.
Languages: EN, RU

Let we have function z = f(x, y) which defines surface in (x, y, z) space. Nodes (x, y) from uniform grid. z is height for each node.

Table of contents

Coordinate axises

Initially we write program with coordinate axises drawing.
Plane of (x, y) points is horizontal base, by z axis values of surface function.

Surface definition

     ^ y
  y3 |   *-----*-----*-----*-----*  min_x = x0; max_x = x4
     |   |     |     |     |     |  min_y = y0; max_y = y3
     |   |     |     |     |     |  n_nodes_by_x = 5
     |   |     |     |     |     |  n_nodes_by_y = 4
  y2 |   *-----*-----*-----*-----*
     |   |     |     |     |     |
     |   |     |     |     |     |
     |   |     |     |     |     |
  y1 |   *-----*-----*-----*-----*
     |   |     |     |     |     |
     |   |     |     |     |     |
     |   |     |     |     |     |
  y0 |   *-----*-----*-----*-----*
     +-------------------------------> x
         x0    x1    x2    x3    x4
  

                                         ^ z
                                         |
                                         |
                 f(x3, y1)               |
                @              f(x1, y0) |
                |              @         |
                |              |         |
             x4 |  x3    x2    | x1  x0  |
       x <------|--------------|---------+
             *--|--*-----*-----*-----*  / y0
            /   | /     /     /     /  /
           /    |/     / @ f(x1, y2)  /
          *-----*-----*--|--*-----*  / y1
         /     /     /   | /     /  /
        /     /     /    |/     /  /
       *-----*-----*-----*-----*  / y2
      /     /     /     /     /  /
     /     /     /     /     /  /
    *-----*-----*-----*-----*  / y3
                            y v
    

Define surface as set of points (x, y, z) where (x, y) are nodes of uniform grid. Grid called as uniform when distance between neighboring nodes along each axis is constant.
In every node we knows z as function of (x, y): z = f(x, y). We call it value as height. Set of heights we call as heights map.

Surface via spheres

Let's build some surface using f(x, y) function. At the start our surface is a set of small spheres.

Surface via triangles

For filling empty space among nodes we need to draw OpenGL's primitives with vertices on nodes. Split cells of grid by triangles.
Now we need to change drawing function only due all vertices data stored at corresponding array.

     ^ y
  y3 |   *-----*-----*-----*-----*
     |   |   / |   / |   / |   / |
     |   |  /  |  /  |  /  |  /  |
     |   | /   | /   | /   | /   |
  y2 |   *-----*-----*-----*-----*
     |   |   / |   / |   / |   / |
     |   |  /  |  /  |  /  |  /  |
     |   | /   | /   | /   | /   |
  y1 |   *-----*-----*-----*-----*
     |   |   / |   / |   / |   / |
     |   |  /  |  /  |  /  |  /  |
     |   | /   | /   | /   | /   |
  y0 |   *-----*-----*-----*-----*
     +-------------------------------> x
         x0    x1    x2    x3    x4
  

  (x,y+dy) ____top____(x+dx,y+dy)
          |  <--     /|
          | |   ^   / |
          | v   |  /  |
         l|  -->  /   |r
         e|      /    |i
         f|     /     |g
         t|    / <--  |h
          |   / |   ^ |t
          |  /  V   | |
          | /    -->  |
          |/__________|
      (x,y)  bottom (x+dx,y)
    

Normals computation

For robust shading our surface set normals to it for every vertex.

                  ^ normal
(x2, y2, z2) _____|____
             \    |   / (x1, y1, z1)
              \   |  /
               \    /
                \  /
                 \/ (x3, y3, z3)

For triangles with vertices (x1, y1, z1), (x2, y2, z2), (x3, y3, z3), in counter-clockwise, normal is a vector (nx, ny, nz):

nx = y1 * (z2 - z3) + y2 * (z3 - z1) + y3 * (z1 - z2)
ny = x1 * (z3 - z2) + x2 * (z1 - z3) + x3 * (z2 - z1)
nz = x1 * (y2 - y3) + x2 * (y3 - y1) + x3 * (y1 - y2)

For each internal vertex (node of which not on top, left, bottom or right border of grid) computing normal as mean of normals for triangles includes this vertex.
Boundary normals may be copyied from neighboring layer of internal vertices.

^                               1) nx = dy*(z1-z0); ny = dx*(z2-z1); nz = dx*dy;
|                               2) nx = dy*(z2-z3); ny = dx*(z3-z0); nz = dx*dy;
y+dy   ______z6_____z5          3) nx = dy*(z0-z4); ny = dx*(z3-z0); nz = dx*dy;
|     |****/ |    / |           4) nx = dy*(z0-z4); ny = dx*(z4-z5); nz = dx*dy;
|     |***/  | 5 /  |           5) nx = dy*(z6-z5); ny = dx*(z0-z6); nz = dx*dy;
|     |**/   |  /   |           6) nx = dy*(z1-z0); ny = dx*(z0-z6); nz = dx*dy;
|     |*/  6 | /  4 |
y     z1-----z0-----z4
|     |  1 / |  3 /*|
|     |   /  |   /**|     nx = 1/6 * dy * (2 * z1 + z2 - z3 - 2 * z4 - z5 + z6)
|     |  / 2 |  /***|     ny = 1/6 * dx * (-z1 + z2 + 2 * z3 + z4 - z5 - 2 * z6)
|     | /    | /****|     nz = dx * dy
y-dy  z2-----z3-----*
|
+--- x-dx -- x --- x+dx ---->

Drawing by indices

Useful to use OpenGL usage with drawing primitives not by explit calls of glVertex and glNormal but via indices of corresponding vertices at our arrays.

Drawing by triangle strip

Using certain indices order we can more efficiently draw triangles which forms strip. It is directly our case because all grid could be divided by strips.

0_____2_____4_____6         GL_TRIANGLES: 1, 2, 0, 1, 3, 2, 3, 4, 2, 3, 5, 4,
|    /|    /|    /|                       5, 6, 4, 5, 7, 6
|  /  |  /  |  /  |         GL_TRIANGLE_STRIP: 0, 1, 2, 3, 4, 5, 6, 7
|/____|/____|/____|
1     3     5     7

Least it gives us less number of stored indices. But more profit achieves due OpenGL use our tip about generic vertices. For us it decreasing time of frame drawing.

During jump between strips we changing indexing direction (from right to left and from bottom to top or vice versa). We can understand it via imaging grid as single folded strip:

8_____9_____10____11   indices: 4, 0, 5, 1, 6, 2, 7, 3, 7, 11, 6, 10, 5, 9, 4, 8
|    /|    /|    /|
|  /  |  /  |  /  |
|/____|/____|/____|
4    /5    /6    /7
|  /  |  /  |  /  |
|/____|/____|/____|
0     1     2     3
                           /\
                          / |\
                         /  | \
           .------      /\  | / 8
          /            / |\ |/
         /            /  | \/
        \/           /\  | / 9
                    / |\ |/
4_____5_____6_____7/  | \/
|    /|    /|    /|\  | / 10
|  /  |  /  |  /  | \ |/
|/____|/____|/____|  \/
0     1     2     3   11

If comparing number of frames per second using code snippets from this page, we can demonstrating advantages of each approach.
(For estimation drawing 100 surfaces simultaneously, grids is 100x100, surfaces from similar function, with different scales along (x, y). Code from examples not optimized thus values from table below are approximated. Each number is accumulated during 5 seconds).

+---------------+-------------------+
| Method        | Frames per second |
+---------------+-------------------+
| via triangles | 9                 |
| by indices    | 38-40             |
| by strip      | 59                |
+----------------+------------------+