Motivation: memory and connectedness.

So far, we have been representing our model geometries using a "soup" of triangles. In other words, we have just been given an array of triangles without any information about how they are connected to each other. This is fine for the ray tracing we have been doing, as well as most rendering algorithms, but when we want to actually create this set of triangles to represent a model, this idea of "connectedness" is important.

There are also some advantages to representing a surface as a connected set of triangles in a rendering pass. For example, consider the mesh below with 5 triangles. How much memory does it take to store this mesh if we represent each triangle by three 3d vectors? Assume floating-point values are stored as Float32 (4 bytes).

The issue with our previous representation is that the vertices shared by multiple triangles are duplicated in memory. We can solve this problem by storing a single vertex in memory and have each triangle reference that vertex. For the example above, there are 6 vertices and 5 triangles. The 6 vertices consume 6 x 3 x 4 = 72 bytes of memory. The memory used to store the triangles depends on how we reference the vertices to define the triangles. Assuming a 32-bit integer is used, then we have 3 x 4 = 12 bytes per triangle, thus 12 x 5 = 60 bytes for all 5 triangles. The total memory used to store the mesh is then 72 + 60 = 132 bytes. The difference with the 180 bytes we had before might not seem like a big difference, but this difference is much greater for larger meshes (more triangles).

Note that for this tiny mesh, we could even have gotten away with using a single byte for the triangle vertex indices (as an unsigned integer with a range from 0 - 255), which would have consumed 5 x 3 = 15 bytes for the triangles and 87 bytes total for the vertices and triangles.

In addition to memory consumption, another advantage of the new proposed representation is that triangles are connected. If the shapes in a spider web weren't connected, it would just fall to the ground. For us, this idea of connection is important to devise algorithms to search through and modify meshes.

Mesh: geometry and topology.

We will represent our models using a mesh, which contains two core ingredients: the geometry and topology. The geometry of the mesh refers to the points of the model in 3d space, whereas the topology refers to how those points are connected to represent the model. In our course, we will always use triangle meshes, but it's good to know that other types of meshes exist, for example, quadrilateral or more general polygons for surface meshes. It's also possible to have volumetric elements like tetrahedra, hexahedra, prisms, pyramids or general polyhedra. The fundamental shape of a mesh is usually called an element.

We will describe the geometry using an array of vertices (sometimes called nodes). These will form the vertices of each triangle in the mesh. We will also assign a unique number to each vertex which corresponds to its location in our array of vertices. Here is an example of a two-triangle, four-vertex mesh.

Storing a mesh with 1d arrays using a stride.

Instead of storing vertices in a 2d array (maybe, where each row corresponds to a single vertex), we will store the coordinates in a 1d array. Since we are working with 3d coordinates, each chunk of 3 values in this array represents the coordinates of a particular vertex. In the example mesh above, the coordinates of the vertices can be represented in a 1d array called vertices like this:

vertices = [x0, y0, z0, x1, y1, z1, x2, y2, z2, x3, y3, z3];

Similarly, the two triangles in the mesh can be represented in a 1d array called triangles like this:

triangles = [3, 2, 1, 2, 3, 0];

To access the y-coordinate of the third vertex (y2), we would use:

y2 = vertices[3 * 2 + 1];

where the 3 represents the stride to take when accessing each vertex. We use a stride of 3 since there are 3 coordinates for each vertex. The + 1 accesses the y coordinate here, which would be + 0 for the x-component and + 2 for the z-component. We could also access all three components using the slice method:

p2 = vertices.slice(3 * 2, 3 * 2 + 3);

The variable p2 will be an array of length 3, which can also be interpreted as a vec3 in glMatrix (passing this variable to glMatrix vec3 functions will work).

To more generally access the d'th component of the i'th vertex, we should use:

vertices[3 * i + d]

Triangles should always be oriented with a counterclockwise (CCW) ordering, which is shown in the blue arrows in the image above. This is important when we shade our models and need the normal vector (which depends on this ordering) to be pointing out of the model. If you trace the order of the triangle vertices with your right hand, then your thumb will point in the direction of the normal. The first triangle references vertices 3, 2 and 1 and the second triangle references vertices 2, 3 and 0.

Similar to how we access the coordinates of a vertex, we can access the vertex indices of each triangle using a stride of 3 (since there are three vertices per triangle). For example, the third vertex of the second triangle (which is 0) would be accessed using:

triangles[3 * 1 + 2]

where the 3 again represents the stride, the 1 represents the triangle index and the + 2 is used to access the third vertex index. Similary to what we did before, we can access the j'th vertex of the k'th triangle using:

triangles[3 * k + j]

To then access the vertex coordinates in a particular triangle, we can combine what we just derived for vertices and triangles. We can access the d'th coordinate of the jth vertex in triangle k using:

vertices[3 * triangles[3 * k + j] + d]

Traversing the vertex indices of a triangle.

To go through all the vertices of the k'th triangle, we can do something like this:

for (let j = 0; j < 3; j++) {
  const v = triangles[3 * k + j];
}

But what if we want to go through the edges (triangle borders)? We know that when we are at vertex j, the one next to it (next in the CCW ordering) is at location j + 1, so what about this?

for (let j = 0; j < 3; j++) {
  const p = triangles[3 * k + j];
  const q = triangles[3 * k + j + 1];
}

Bug! This will work for j = 0 and j = 1, but the "next" vertex around the triangle when we are at index j = 2 is actually 0. There are different ways to fix this, and here's one option using modulo 3, since there are 3 vertices per triangle:

for (let j = 0; j < 3; j++) {
  const p = triangles[3 * k + j];
  const q = triangles[3 * k + (j + 1) % 3];
}

Similarly, if we want to access a vertex "opposite" an edge (opp), we can use an offset of 2 before the modulus operator:

for (let j = 0; j < 3; j++) {
  const p = triangles[3 * k + j];
  const q = triangles[3 * k + (j + 1) % 3];
  const opp = triangles[3 * k + (j + 2) % 3];
}

If you really want to optimize your code, you could use a look-up table instead of doing the modulus over and over again. For example, we could have:

const edgeLookup = [[0, 1], [1, 2], [2, 0]];
const oppositeLookup = [2, 1, 0];

for (let j = 0; j < 3; j++) {
  const p = triangles[3 * k + edgeLookup[j][0]];
  const q = triangles[3 * k + edgeLookup[j][1]];
  const opp = triangles[3 * k + oppositeLookup[j]];
}

You can also use a 1d array for edgeLookup with a stride of 2 since there are two vertices per edge.

In the image above, $q_0$ is the index of the vertex opposite edge $(e_0, e_1)$ when traversing the left triangle and $q_1$ is the index of the vertex opposite edge $(e_1, e_0)$ when traversing the right triangle (remember to go in CCW order).

Extracting the edges of a mesh.

Some algorithms may require us to extract the edges in the mesh. It's also useful to extract edges so we can directly visualize the borders of all the triangle elements. We know how to loop through the edges of a single triangle, which we'll use to build a global list of edges. When encountering an edge of a triangle, we'll check if we have added this edge yet, and if we haven't, then add it to our global list.

Internal edges always appear twice (for the types of meshes we will see today). Since we are using a consistent CCW ordering of the triangles, these edges will appear once in one direction and another in the opposite direction. For example, in the two-triangle mesh above, we will have edge 3-2 from triangle 0 and edge 2-3 from triangle 1. Since each vertex is associated with a unique index, we can uniquely identify edges by labelling the first vertex of the edge as the min of these two values, and the second vertex of the edge as the max of these two values. This provides us with a sort of "key" that we can use to keep track of whether we have visited an edge or not.

Here is one way to extract the edges in a mesh. Note that this operation is entirely topological - we do not need any information about the vertex coordinates to extract the edges. The integers in the 1d array of edges reference the same vertices array as the triangles.

const extractEdges(triangles) => {
  let edgeMap = {}; // initialize dictionary to track edges
  let edges = []; // the list of edges we will extract and return
  const nTriangles = triangles.length / 3;

  // loop through all the triangles
  for (let k = 0; k < nTriangles; k++) {
    // loop through all 3 edges of triangle k
    for (let j = 0; j < 3; j++) {
      // extract the two vertex indices of edge j in triangle k
      const p = triangles[3 * k + j];
      const q = triangles[3 * k + (j + 1) % 3];

      // uniquely represent this edge 
      const edge = [Math.min(p, q), Math.max(p, q)];
      const edgeKey = JSON.stringify(edge); // unique edge key
      
      // does this edge exist in our map yet?
      if (edgeKey in edgeMap) {
        // this edge already exists, nothing to do!
      } else {
        // this edge does not exist yet, add it to the map and to the list
        const edgeIndex = edges.length / 2;
        edgeMap[edgeKey] = edgeIndex;
        edges.push(p, q); // not necessarily sorted, but it doesn't matter
      }
    }
  }
  return edges;
}

Note that the push method for Array accepts multiple arguments, so you can add multiple items to our 1d edges array in one function call.

Subdividing a mesh.

A lot of the models we have worked so far have had a low-poly look to them, primarily because we have only used at most a few hundred triangles to represent them. Subdivision is the process of dividing triangles to create a finer mesh. There are different schemes you can use to perform the subdivision and, again, we need to make the distinction between geometry and topology.

We will focus on subdividing a mesh such that every triangle is split up into four new triangles (discarding the original triangles). This involves introducing a new vertex along every (unique) edge of the mesh. This topological modification will always be the same, regardless of the scheme you are implementing. However, the geometry of the subdivision differs based on the scheme you choose. We will implement a simple subdivision scheme together today, and will leave an important subdivision scheme for this week's lab.

The topology of subdivision.

As mentioned earlier, we will add a new vertex along every edge of the mesh. Every existing vertex will still be there, so we will just append new vertex coordinates to our vertices array. Let nVertices be the initial number of vertices before the subdivision and let edgeIndex represent the index of the edge we are processing (same variable as in the code above). After appending the coordinates to the vertices array, the index of the vertex we inserted is nVertices + edgeIndex.

For every triangle, we can then identify the indices of the three vertices that were inserted along the edge:

let edgeIndices = new Array(3);

// loop through all 3 edges of triangle k
for (let j = 0; j < 3; j++) {
  // extract the two vertex indices of edge j in triangle k
  const p = triangles[3 * k + j];
  const q = triangles[3 * k + (j + 1) % 3];

  // uniquely represent this edge 
  const edge = [Math.min(p, q), Math.max(p, q)];
  const edgeKey = JSON.stringify(edge); // unique edge key

  // does this edge exist in our map yet?
  if (edgeKey in edgeMap) {
    // this edge already exists, store the index we had assigned to it
    edgeIndices[j] = nVertices + edgeMap[edgeKey];
  } else {
    // this edge does not exist yet, add it to the map and save the index
    const edgeIndex = edges.length / 2;
    edgeMap[edgeKey] = edgeIndex;
    edgeIndices[j] = nVertices + edgeIndex;
  }
}

Now the last thing we need to do is create four new triangles. All of these triangles need to maintain the same consistent CCW orientation.

What are the vertex indices of the four triangles created in this figure?

Exercise: subdividing an icosahedron to create a better approximation of a sphere.

Let's practice with the concepts introduced above by subdividing the mesh of an icosahedron (initially 12 vertices and 20 triangles) to produce a better approximation of a unit sphere (radius = 1) centered at the origin (0, 0, 0). We will use a loop very similar to the loop used to extract edges and create vertices and triangles on the fly.

To define the geometry of the subdivision, we will leave existing vertices untouched, since we will assume that they are already exactly on the sphere. When creating vertices along edges, we need to make sure they are on the unit sphere or, stated another way, that they have a distance of 1 from the origin. We can achieve this by first calculating the midpoint of every edge and then normalizing the position vector of this midpoint.

Please see the repl for this exercise here.

Loop subdivision scheme.

The example we just did was restricted to approximating a sphere. For more general shapes, we can use well-established schemes for creating Subdivision Surfaces. For quadrilateral meshes, we can use the Catmull-Clark subdivision scheme, and for triangle meshes, we can use the Loop Subdivision scheme, named after Charles Loop who invented the scheme as part of a Master's thesis.

The technology of subdivision surfaces first made an appearance in the 1997 Pixar short called Geri's Game (see the image at the top-left of this lecture). Here, subdivision surfaces were used to obtain much smoother surfaces from very coarse meshes, like we see in Geri's hand and the chess pieces.

We will restrict our attention to closed manifold meshes, which means that every edge in the mesh is shared by exactly two triangles. A mesh with boundary will have some edges adjacent to a single triangle and some non-manifold meshes might have some edges adjacent to three or more triangles.

The topological modifications in the Loop subdivision scheme are almost exactly the same as what we just did in the sphere exercise. However, the geometry of the new vertices is different and requires two modifications:

1. New vertices along edges: The blue dot in the leftmost image below is the weighted sum of the coordinates of the black dots. In other words, this new edge vertex receives 3/8 of the edge endpoint vertices and 1/8 of the opposite vertices (use the the same "opposite" method defined above).
2. Existing vertices: All of the black dots in the rightmost image below exist previously in the mesh. Consider treating the one in the center with the outer red ring, and denote this as $\vec{p}_i$ which is the vertex at index i. This vertex will move to a location which is the weighted sum: $\vec{p}_i^* = (1 - k\beta) \vec{p}_i + \sum_j\beta \vec{p}_j$.

where $k$ is the number of vertices surrounding $\vec{p}_i$ and $\beta$ is equal to 0.1875 if $k$ is equal to 3 and $\frac{0.375}{k}$ otherwise.

The image in the middle describes how boundary edge vertices are treated, which you will not need for this week's lab (it's only here in case you want to extend your algorithm to treat more general meshes with boundary).

Note that step 2 requires knowing all the vertices that are adjacent to the red-ring vertex in the center. We'll call this data structure v2v (for vertex-to-vertex), which can also be calculated by extending our previous loop to extract edges:

// initialize an array for each vertex to place the v2v information
let v2v = new Array(nVertices);
for (let i = 0; i < nVertices; i++) v2v[i] = [];

for (let i = 0; i < nTriangles; i++) {
  for (let j = 0; j < 3; j++) {
    let p = this.triangles[3 * i + j];
    let q = this.triangles[3 * i + ((j + 1) % 3)];

    // add q as a neighbor to p
    // (p will be added as a neighbor to q when looping through the adjacent triangle)
    v2v[p].push(q);
  }
}

Note that $k = \texttt{v2v[i].length}$.

Here is what the Loop subdivision scheme produces on an initial mesh of a toroidal tetrahedron.

© Philip Claude Caplan, 2023 (Last updated: 2023-10-18)