## 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.

*(Model courtesy of headus/Rezard.)*

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
`Triangle`s that represent the final subdivided mesh.

`LoopSubdiv`vertices and faces>>

`faces`>>

`edgeNum`>>

`f`and

`v`for next level of subdivision>>

`k`th edge>>

`f`pointers for siblings>>

`f`pointers for neighbor children>> } }

### 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.

`LoopSubdiv`vertices and faces>>=

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 `i`th 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 modulus 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.

`faces`>>=

`edgeNum`>>

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.

`edgeNum`>>=

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
`i`th 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.

`f`and

`v`for next level of subdivision>>

`k`th edge>>

`f`pointers for siblings>>

`f`pointers for neighbor children>> } }

The main loop of a subdivision step proceeds as follows: it creates
`vector`s 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
`i`th child face is next to the `i`th 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.

`i`th child face is adjacent to the

`i`th vertex of the original face and the fourth child face is in the center of the subdivided face. Three edge vertices need to be computed; they are numbered so that the

`i`th edge vertex is along the

`i`th edge of the original face.

`f`and

`v`for next level of subdivision>>=

`k`th edge>>

`f`pointers for siblings>>

`f`pointers for neighbor children>> } }

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.*

`k`th edge>>

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`.

`vert`one-ring in

`pRing`>> Point3f p = (1 - valence * beta) * vert->p; for (int i = 0; i < valence; ++i) p += beta * pRing[i]; return p; }

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

`vert`one-ring in

`pRing`>>=

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.

`vert`one-ring in

`pRing`>> Point3f p = (1 - 2 * beta) * vert->p; p += beta * pRing[0]; p += beta * pRing[valence - 1]; return p; }

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.

`k`th edge>>

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.

`k`th edge>>=

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.

`i`th child is adjacent to the

`i`th vertex of the original face, and such that the

`i`th child face’s

`i`th vertex is the child of the

`i`th vertex of the original face. The vertices of the center child are oriented such that the

`i`th vertex is the odd vertex along the

`i`th edge of the parent face.

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

- The odd vertices’
`SDVertex::startFace`pointers need to store a pointer to one of their adjacent faces. - Similarly, the even vertices’
`SDVertex::startFace`pointers must be set. - The new faces’ neighbor
`f[i]`pointers need to be set to point to the neighboring faces. - 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.

`f`pointers for siblings>>

`f`pointers for neighbor children>> } }

If a vertex is the `i`th vertex of its `startFace`, then it is
guaranteed that it will be adjacent to the `i`th 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.

`f`pointers for siblings>>

`f`pointers for neighbor children>> } }

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 `k`th edge of the interior face, and the interior
face is across the st edge of the `k`th face.

`f`pointers for siblings>>=

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 `k`th and
`PREV(k)`th edges of the `i`th child need to be set. To set the `k`th
edge of the `k`th child, we first find the `k`th 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.

`f`pointers for neighbor children>>=

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

For the `k`th child face (for ), the `k`th
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

Valence | ||
---|---|---|

2 | ||

3 | ||

4 (regular) |

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 `Triangle`s 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.