3.8 Subdivision Surfaces

The last shape representation that we’ll define in this chapter implements subdivision surfaces, a representation that is particularly well suited to describing complex smooth shapes. The subdivision surface for a particular mesh is defined by repeatedly subdividing the faces of the mesh into smaller faces and then finding the new vertex locations using weighted combinations of the old vertex positions.

For appropriately chosen subdivision rules, this process converges to give a smooth limit surface as the number of subdivision steps goes to infinity. In practice, just a few levels of subdivision typically suffice to give a good approximation of the limit surface. Figure 3.24 shows a simple example of a subdivision, where a tetrahedron has been subdivided zero, one, two, and six times.

Figure 3.25 shows the effect of applying subdivision to the Killeroo model; on the top is the original control mesh, and below is the subdivision surface that the control mesh represents.

Although originally developed in the 1970s, subdivision surfaces have seen widespread use in recent years thanks to some important advantages over polygonal and spline-based representations of surfaces. The advantages of subdivision include the following:

• Subdivision surfaces are smooth, as opposed to polygon meshes, which appear faceted when viewed close up, regardless of how finely they are modeled.
• Much of the existing infrastructure in modeling systems can be retargeted to subdivision. The classic toolbox of techniques for modeling polygon meshes can be applied to modeling subdivision control meshes.
• Subdivision surfaces are well suited to describing objects with complex topology, since they start with a control mesh of arbitrary (manifold) topology. Parametric surface models generally don’t handle complex topology well.
• Subdivision methods are often generalizations of spline-based surface representations, so spline surfaces can often just be run through general subdivision surface renderers.
• It is easy to add detail to a localized region of a subdivision surface simply by adding faces to appropriate parts of the control mesh. This is much harder with spline representations.

Here, we will describe an implementation of Loop subdivision surfaces. ^† The Loop subdivision rules are based on triangular faces in the control mesh; faces with more than three vertices are triangulated at the start. At each subdivision step, all faces split into four child faces (Figure 3.26). New vertices are added along all of the edges of the original mesh, with positions computed using weighted averages of nearby vertices. Furthermore, the position of each original vertex is updated with a weighted average of its position and its new neighbors’ positions. The implementation here uses weights based on improvements to Loop’s method developed by Hoppe et al. (1994). We will not include discussion here about how these weights are derived. They must be chosen carefully to ensure that the limit surface actually has particular desirable smoothness properties, although subtle mathematics is necessary to prove that they indeed do this.

Rather than being implemented as a Shape in pbrt, subdivision surfaces are generated by a function, LoopSubdivide(), that applies the subdivision rules to a mesh represented by a collection of vertices and vertex indices and returns a vector of Triangles that represent the final subdivided mesh.

3.8.1 Mesh Representation

The parameters to LoopSubdivide() specify a triangle mesh in exactly the same format used in the TriangleMesh constructor (Section 3.6): each face is described by three integer vertex indices, giving offsets into the vertex array p for the face’s three vertices. We will need to process this data to determine which faces are adjacent to each other, which faces are adjacent to which vertices, and so on, in order to implement the subdivision algorithm.

We will shortly define SDVertex and SDFace structures, which hold data for vertices and faces in the subdivision mesh. LoopSubdivide() starts by allocating one instance of the SDVertex class for each vertex in the mesh and an SDFace for each face. For now, these are mostly uninitialized.

The Loop subdivision scheme, like most other subdivision schemes, assumes that the control mesh is manifold—no more than two faces share any given edge. Such a mesh may be closed or open: a closed mesh has no boundary, and all faces have adjacent faces across each of their edges. An open mesh has some faces that do not have all three neighbors. The implementation here supports both closed and open meshes.

In the interior of a triangle mesh, most vertices are adjacent to six faces and have six neighbor vertices directly connected to them with edges. On the boundaries of an open mesh, most vertices are adjacent to three faces and four vertices. The number of vertices directly adjacent to a vertex is called the vertex’s valence. Interior vertices with valence other than six, or boundary vertices with valence other than four, are called extraordinary vertices; otherwise, they are called regular. ^† Loop subdivision surfaces are smooth everywhere except at their extraordinary vertices.

Each SDVertex stores its position p, a Boolean that indicates whether it is a regular or extraordinary vertex, and a Boolean that records if it lies on the boundary of the mesh. It also holds a pointer to an arbitrary face adjacent to it; this pointer gives a starting point for finding all of the adjacent faces. Finally, there is a pointer to store the corresponding SDVertex for the next level of subdivision, if any.

The SDFace structure is where most of the topological information about the mesh is maintained. Because all faces are triangular, faces always store pointers to their three vertices and pointers to the adjacent faces across its three edges. The corresponding face neighbor pointers will be nullptr if the face is on the boundary of an open mesh.

The face neighbor pointers are indexed such that if we label the edge from v[i] to v[(i+1)%3] as the ith edge, then the neighbor face across that edge is stored in f[i] (Figure 3.27). This labeling convention is important to keep in mind. Later when we are updating the topology of a newly subdivided mesh, we will make extensive use of it to navigate around the mesh. Similarly to the SDVertex class, the SDFace also stores pointers to child faces at the next level of subdivision.

The SDFace constructor is straightforward—it simply sets these various pointers to nullptr—so it is not shown here.

In order to simplify navigation of the SDFace data structure, we’ll provide macros that make it easy to determine the vertex and face indices before or after a particular index. These macros add appropriate offsets and compute the result modulo three to handle cycling around.

In addition to requiring a manifold mesh, the subdivision code expects that the control mesh specified by the user will be consistently ordered—each directed edge in the mesh can be present only once. An edge that is shared by two faces should be specified in a different direction by each face. Consider two vertices, and , with an edge between them. We expect that one of the triangular faces that has this edge will specify its three vertices so that is before , and that the other face will specify its vertices so that is before (Figure 3.28). A Möbius strip is one example of a surface that cannot be consistently ordered, but such surfaces come up rarely in rendering so in practice this restriction is not a problem. Poorly formed mesh data from other programs that don’t create consistently ordered meshes can be troublesome, however.

Given this assumption about the input data, the LoopSubdivide() can now initialize the mesh’s topological data structures. It first loops over all of the faces and sets their v pointers to point to their three vertices. It also sets each vertex’s SDVertex::startFace pointer to point to one of the vertex’s neighboring faces. It doesn’t matter which of its adjacent faces is used, so the implementation just keeps resetting it each time it comes across another face that the vertex is incident to, thus ensuring that all vertices have a non-nullptr face pointer by the time the loop is complete.

Now it is necessary to set each face’s f pointer to point to its neighboring faces. This is a bit trickier, since face adjacency information isn’t directly specified in the data passed to LoopSubdivide(). The implementation here loops over the faces and creates an SDEdge object for each of their three edges. When it comes to another face that shares the same edge, it can update both faces’ neighbor pointers.

The SDEdge constructor takes pointers to the two vertices at each end of the edge. It orders them so that v[0] holds the one that is first in memory. This code may seem strange, but it is simply relying on the fact that pointers in C++ are effectively numbers that can be manipulated like integers ^† and that the ordering of vertices on an edge is arbitrary. Sorting the two vertices based on the addresses of their pointers guarantees that the edge (v, v) is correctly recognized as the same as the edge (v, v), regardless of what order the vertices are provided in.

The class also defines an ordering operation for SDEdge objects so that they can be stored in other data structures that rely on ordering being well defined.

Now the LoopSubdivide() function can get to work, looping over the edges in all of the faces and updating the neighbor pointers as it goes. It uses a set to store the edges that have only one adjacent face so far. The set makes it possible to search for a particular edge in time.

For each edge in each face, the loop body creates an edge object and sees if the same edge has been seen previously. If so, it initializes both faces’ neighbor pointers across the edge. If not, it adds the edge to the set of edges. The indices of the two vertices at the ends of the edge, v0 and v1, are equal to the edge index and the edge index plus one.

Given an edge that hasn’t been encountered before, the current face’s pointer is stored in the edge object’s f[0] member. Because the input mesh is assumed to be manifold, there can be at most one other face that shares this edge. When such a face is discovered, it can be used to initialize the neighboring face field. Storing the edge number of this edge in the current face allows the neighboring face to initialize its corresponding edge neighbor pointer.

When the second face on an edge is found, the neighbor pointers for each of the two faces are set. The edge is then removed from the edge set, since no edge can be shared by more than two faces.

Now that all faces have proper neighbor pointers, the boundary and regular flags in each of the vertices can be set. In order to determine if a vertex is a boundary vertex, we’ll define an ordering of faces around a vertex (Figure 3.29). For a vertex v[i] on a face f, we define the vertex’s next face as the face across the edge from v[i] to v[NEXT(i)] and the previous face as the face across the edge from v[PREV(i)] to v[i].

By successively going to the next face around v, we can iterate over the faces adjacent to it. If we eventually return to the face we started at, then we are at an interior vertex; if we come to an edge with a nullptr neighbor pointer, then we’re at a boundary vertex (Figure 3.30). Once the initialization routine has determined if this is a boundary vertex, it computes the valence of the vertex and sets the regular flag if the valence is 6 for an interior vertex or 4 for a boundary vertex; otherwise, it is an extraordinary vertex.

Because the valence of a vertex is frequently needed, we provide the method SDVertex::valence().

To compute the valence of a nonboundary vertex, this method counts the number of the adjacent faces starting by following each face’s neighbor pointers around the vertex until it reaches the starting face. The valence is equal to the number of faces visited.

For boundary vertices we can use the same approach, although in this case, the valence is one more than the number of adjacent faces. The loop over adjacent faces is slightly more complicated here: it follows pointers to the next face around the vertex until it reaches the boundary, counting the number of faces seen. It then starts again at startFace and follows previous face pointers until it encounters the boundary in the other direction.

SDFace::vnum() is a utility function that finds the index of a given vertex pointer. It is a fatal error to pass a pointer to a vertex that isn’t part of the current face—this case would represent a bug elsewhere in the subdivision code.

Since the next face for a vertex v[i] on a face f is over the ith edge (recall the mapping of edge neighbor pointers from Figure 3.27), we can find the appropriate face neighbor pointer easily given the index i for the vertex, which the vnum() utility function provides. The previous face is across the edge from PREV(i) to i, so the method returns f[PREV(i)] for the previous face.

It is also useful to be able to get the next and previous vertices around a face starting at any vertex. The SDFace::nextVert() and SDFace::prevVert() methods do just that (Figure 3.31).

3.8.2 Subdivision

Now we can show how subdivision proceeds with the modified Loop rules. The implementation here applies subdivision a fixed number of times to generate a triangle mesh for rendering; Exercise 3.9.7 at the end of the chapter discusses adaptive subdivision, where each original face is subdivided enough times so that the result looks smooth from a particular viewpoint rather than just using a fixed number of levels of subdivision, which may over-subdivide some areas while simultaneously under-subdividing others.

The <<Refine subdivision mesh into triangles>> fragment repeatedly applies the subdivision rules to the mesh, each time generating a new mesh to be used as the input to the next step. After each subdivision step, the f and v arrays are updated to point to the faces and vertices from the level of subdivision just computed. When it’s done subdividing, a triangle mesh representation of the surface is returned.

An instance of the MemoryArena class is used to allocate temporary storage through this process. This class, defined in Section A.4.3, provides a custom memory allocation method that quickly allocates memory, automatically freeing the memory when it goes out of scope.

The main loop of a subdivision step proceeds as follows: it creates vectors to store the vertices and faces at the current level of subdivision and then proceeds to compute new vertex positions and update the topological representation for the refined mesh. Figure 3.32 shows the basic refinement rules for faces in the mesh. Each face is split into four child faces, such that the ith child face is next to the ith vertex of the input face and the final face is in the center. Three new vertices are then computed along the split edges of the original face.

First, storage is allocated for the updated values of the vertices already present in the input mesh. The method also allocates storage for the child faces. It doesn’t yet do any initialization of the new vertices and faces other than setting the regular and boundary flags for the vertices since subdivision leaves boundary vertices on the boundary and interior vertices in the interior and it doesn’t change the valence of vertices in the mesh.

Computing New Vertex Positions

Before worrying about initializing the topology of the subdivided mesh, the refinement method computes positions for all of the vertices in the mesh. First, it considers the problem of computing updated positions for all of the vertices that were already present in the mesh; these vertices are called even vertices. It then computes the new vertices on the split edges. These are called odd vertices.

Different techniques are used to compute the updated positions for each of the different types of even vertices—regular and extraordinary, boundary and interior. This gives four cases to handle.

For both types of interior vertices, we take the set of vertices adjacent to each vertex (called the one-ring around it, reflecting the fact that it’s a ring of neighbors) and weight each of the neighbor vertices by a weight (Figure 3.33). The vertex we are updating, in the center, is weighted by , where is the valence of the vertex. Thus, the new position v for a vertex v is

This formulation ensures that the sum of weights is one, which guarantees the convex hull property of Loop subdivision surfaces, which ensures that the final mesh is in the convex hull of the control mesh. The position of the vertex being updated is only affected by vertices that are nearby; this is known as local support. Loop subdivision is particularly efficient because its subdivision rules all have this property.

The specific weight used for this step is a key component of the subdivision method and must be chosen carefully in order to ensure smoothness of the limit surface, among other desirable properties. ^† The beta() function that follows computes a  value based on the vertex’s valence that ensures smoothness. For regular interior vertices, beta() returns . Since this is a common case, the implementation uses directly instead of calling beta() every time.

The weightOneRing() function loops over the one-ring of adjacent vertices and applies the given weight to compute a new vertex position. It uses the SDVertex::oneRing() method, defined in the following, which returns the positions of the vertices around the vertex vert.

Because a variable number of vertices are in the one-rings, we use the ALLOCA() macro to efficiently allocate space to store their positions.

The oneRing() method assumes that the pointer passed in points to an area of memory large enough to hold the one-ring around the vertex.

It’s relatively easy to get the one-ring around an interior vertex by looping over the faces adjacent to the vertex and for each face retaining the vertex after the center vertex. (Brief sketching with pencil and paper should convince you that this process returns all of the vertices in the one-ring.)

The one-ring around a boundary vertex is a bit trickier. The implementation here carefully stores the one-ring in the given Point3f array so that the first and last entries in the array are the two adjacent vertices along the boundary. This ordering is important because the adjacent boundary vertices will often be weighted differently from the adjacent vertices that are in the interior of the mesh. Doing so requires that we first loop around neighbor faces until we reach a face on the boundary and then loop around the other way, storing vertices one by one.

For vertices on the boundary, the new vertex’s position is based only on the two neighboring boundary vertices (Figure 3.34). Not depending on interior vertices ensures that two abutting surfaces that share the same vertices on the boundary will have abutting limit surfaces. The weightBoundary() utility function applies the given weighting on the two neighbor vertices v and v to compute the new position v as

The same weight of is used for both regular and extraordinary vertices.

The weightBoundary() utility function applies the given weights at a boundary vertex. Because the SDVertex::oneRing() function orders the boundary vertex’s one-ring such that the first and last entries are the boundary neighbors, the implementation here is particularly straightforward.

Now the refinement method computes the positions of the odd vertices—the new vertices along the split edges of the mesh. It loops over each edge of each face in the mesh, computing the new vertex that splits the edge (Figure 3.35). For interior edges, the new vertex is found by weighting the two vertices at the ends of the edge and the two vertices across from the edge on the adjacent faces. It loops through all three edges of each face, and each time it comes to an edge that hasn’t been seen before it computes and stores the new odd vertex for the edge in the edgeVerts associative array.

As was done when setting the face neighbor pointers in the original mesh, an SDEdge object is created for the edge and checked to see if it is in the set of edges that have already been visited. If it isn’t, the new vertex on this edge is computed and added to the map, which is an associative array structure that performs efficient lookups.

In Loop subdivision, the new vertices added by subdivision are always regular. (This means that the proportion of extraordinary vertices with respect to regular vertices will decrease with each level of subdivision.) Therefore, the regular member of the new vertex can immediately be set to true. The boundary member can also be easily initialized, by checking to see if there is a neighbor face across the edge that is being split. Finally, the new vertex’s startFace pointer can also be set here. For any odd vertex on the edge of a face, the center child (child face number three) is guaranteed to be adjacent to the new vertex.

For odd boundary vertices, the new vertex is just the average of the two adjacent vertices. For odd interior vertices, the two vertices at the ends of the edge are given weight , and the two vertices opposite the edge are given weight (Figure 3.35). These last two vertices can be found using the SDFace::otherVert() utility function, which returns the vertex opposite a given edge of a face.

The SDFace::otherVert() method is self-explanatory:

Updating Mesh Topology

In order to keep the details of the topology update as straightforward as possible, the numbering scheme for the subdivided faces and their vertices has been chosen carefully (Figure 3.36). Review the figure carefully; the conventions shown are key to the next few pages.

There are four main tasks required to update the topological pointers of the refined mesh:

1. The odd vertices’ SDVertex::startFace pointers need to store a pointer to one of their adjacent faces.
2. Similarly, the even vertices’ SDVertex::startFace pointers must be set.
3. The new faces’ neighbor f[i] pointers need to be set to point to the neighboring faces.
4. The new faces’ v[i] pointers need to point to the appropriate vertices.

The startFace pointers of the odd vertices were already initialized when they were first created. We’ll handle the other three tasks in order here.

If a vertex is the ith vertex of its startFace, then it is guaranteed that it will be adjacent to the ith child face of startFace. Therefore, it is just necessary to loop through all the parent vertices in the mesh, and for each one find its vertex index in its startFace. This index can then be used to find the child face adjacent to the new even vertex.

Next, the face neighbor pointers for the newly created faces are updated. We break this into two steps: one to update neighbors among children of the same parent, and one to do neighbors across children of different parents. This involves some tricky pointer manipulation.

For the first step, recall that the interior child face is always stored in children[3]. Furthermore, the st child face (for ) is across the kth edge of the interior face, and the interior face is across the st edge of the kth face.

We’ll now update the children’s face neighbor pointers that point to children of other parents. Only the first three children need to be addressed here; the interior child’s neighbor pointers have already been fully initialized. Inspection of Figure 3.36 reveals that the kth and PREV(k)th edges of the ith child need to be set. To set the kth edge of the kth child, we first find the kth edge of the parent face, then the neighbor parent f2 across that edge. If f2 exists (meaning we aren’t on a boundary), the neighbor parent index for the vertex v[k] is found. That index is equal to the index of the neighbor child we are searching for. This process is then repeated to find the child across the PREV(k)th edge.

Finally, we handle the fourth step in the topological updates: setting the children faces’ vertex pointers.

For the kth child face (for ), the kth vertex corresponds to the even vertex that is adjacent to the child face. For the noninterior child faces, there is one even vertex and two odd vertices; for the interior child face, there are three odd vertices. This vertex can be found by following the child pointer of the parent vertex, available from the parent face.

To update the rest of the vertex pointers, the edgeVerts associative array is reused to find the odd vertex for each split edge of the parent face. Three child faces have that vertex as an incident vertex. The vertex indices for the three faces are easily found, again based on the numbering scheme established in Figure 3.36.

After the geometric and topological work has been done for a subdivision step, the newly created vertices and faces are moved into the v and f arrays:

To the Limit Surface and Output

One of the remarkable properties of subdivision surfaces is that there are special subdivision rules that give the positions that the vertices of the mesh would have if we continued subdividing forever. We apply these rules here to initialize an array of limit surface positions, pLimit. Note that it’s important to temporarily store the limit surface positions somewhere other than in the vertices while the computation is taking place. Because the limit surface position of each vertex depends on the original positions of its surrounding vertices, the original positions of all vertices must remain unchanged until the computation is complete.

The limit rule for a boundary vertex weights the two neighbor vertices by and the center vertex by . The rule for interior vertices is based on a function loopGamma(), which computes appropriate vertex weights based on the valence of the vertex.

In order to generate a smooth-looking triangle mesh with per-vertex surface normals, a pair of nonparallel tangent vectors to the limit surface is computed at each vertex. As with the limit rule for positions, this is an analytic computation that gives the precise tangents on the actual limit surface.

Figure 3.37 shows the setting for computing tangents in the mesh interior. The center vertex is given a weight of zero, and the neighbors are given weights .

To compute the first tangent vector , the weights are

where is the valence of the vertex. The second tangent is computed with weights

Tangents on boundary vertices are a bit trickier. Figure 3.38 shows the ordering of vertices in the one-ring expected in the following discussion.

The first tangent, known as the across tangent, is given by the vector between the two neighboring boundary vertices:

The second tangent, known as the transverse tangent, is computed based on the vertex’s valence. The center vertex is given a weight and the one-ring vertices are given weights specified by a vector . The transverse tangent rules we will use are

For valences of 5 and higher, and

where

Although we will not prove it here, these weights sum to zero for all values of . This guarantees that the weighted sum is in fact a tangent vector.

Finally, the fragment <<Create triangle mesh from subdivision mesh>> initializes a vector of Triangles corresponding to the triangulation of the limit surface. We won’t include it here, since it’s just a straightforward transformation of the subdivided mesh into an indexed triangle mesh.