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.24: Subdivision of a Tetrahedron. From left to right, zero, one, two, and six subdivision steps have been used. (At zero levels, the vertices are just moved to lie on the limit surface.) As more subdivision is done, the mesh approaches the limit surface, the smooth surface described by the original mesh. Notice how the specular highlights become progressively more accurate and the silhouette edges appear smoother as more levels of subdivision are performed.

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.

Figure 3.25: Subdivision Applied to the Killeroo Model. (1) The control mesh describes (2) the resulting subdivision surface. Subdivision is well suited to modeling shapes like this one, since it’s easy to add detail locally by refining the control mesh, and there are no limitations on the topology of the final surface. (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.

Figure 3.26: Basic Refinement Process for Loop Subdivision. (left) The control mesh before subdivision. (right) The new mesh after one subdivision step. Each triangular face of the mesh has been subdivided into four new faces by splitting each of the edges and connecting the new vertices with new edges.

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.

<<LoopSubdiv Function Definitions>>= 
std::vector<std::shared_ptr<Shape>> LoopSubdivide( const Transform *ObjectToWorld, const Transform *WorldToObject, bool reverseOrientation, int nLevels, int nIndices, const int *vertexIndices, int nVertices, const Point3f *p) { std::vector<SDVertex *> vertices; std::vector<SDFace *> faces; <<Allocate LoopSubdiv vertices and faces>> 
std::unique_ptr<SDVertex[]> verts(new SDVertex[nVertices]); for (int i = 0; i < nVertices; ++i) { verts[i] = SDVertex(p[i]); vertices.push_back(&verts[i]); } int nFaces = nIndices / 3; std::unique_ptr<SDFace[]> fs(new SDFace[nFaces]); for (int i = 0; i < nFaces; ++i) faces.push_back(&fs[i]);
<<Set face to vertex pointers>> 
const int *vp = vertexIndices; for (int i = 0; i < nFaces; ++i, vp += 3) { SDFace *f = faces[i]; for (int j = 0; j < 3; ++j) { SDVertex *v = vertices[vp[j]]; f->v[j] = v; v->startFace = f; } }
<<Set neighbor pointers in faces>> 
std::set<SDEdge> edges; for (int i = 0; i < nFaces; ++i) { SDFace *f = faces[i]; for (int edgeNum = 0; edgeNum < 3; ++edgeNum) { <<Update neighbor pointer for edgeNum>> 
int v0 = edgeNum, v1 = NEXT(edgeNum); SDEdge e(f->v[v0], f->v[v1]); if (edges.find(e) == edges.end()) { <<Handle new edge>> 
e.f[0] = f; e.f0edgeNum = edgeNum; edges.insert(e);
} else { <<Handle previously seen edge>> 
e = *edges.find(e); e.f[0]->f[e.f0edgeNum] = f; f->f[edgeNum] = e.f[0]; edges.erase(e);
}
} }
<<Finish vertex initialization>> 
for (int i = 0; i < nVertices; ++i) { SDVertex *v = vertices[i]; SDFace *f = v->startFace; do { f = f->nextFace(v); } while (f && f != v->startFace); v->boundary = (f == nullptr); if (!v->boundary && v->valence() == 6) v->regular = true; else if (v->boundary && v->valence() == 4) v->regular = true; else v->regular = false; }
<<Refine subdivision mesh into triangles>> 
std::vector<SDFace *> f = faces; std::vector<SDVertex *> v = vertices; MemoryArena arena; for (int i = 0; i < nLevels; ++i) { <<Update f and v for next level of subdivision>> 
std::vector<SDFace *> newFaces; std::vector<SDVertex *> newVertices; <<Allocate next level of children in mesh tree>> 
for (SDVertex *vertex : v) { vertex->child = arena.Alloc<SDVertex>(); vertex->child->regular = vertex->regular; vertex->child->boundary = vertex->boundary; newVertices.push_back(vertex->child); } for (SDFace *face : f) { for (int k = 0; k < 4; ++k) { face->children[k] = arena.Alloc<SDFace>(); newFaces.push_back(face->children[k]); } }
<<Update vertex positions and create new edge vertices>> 
<<Update vertex positions for even vertices>> 
for (SDVertex *vertex : v) { if (!vertex->boundary) { <<Apply one-ring rule for even vertex>> 
if (vertex->regular) vertex->child->p = weightOneRing(vertex, 1.f / 16.f); else vertex->child->p = weightOneRing(vertex, beta(vertex->valence()));
} else { <<Apply boundary rule for even vertex>> 
vertex->child->p = weightBoundary(vertex, 1.f / 8.f);
} }
<<Compute new odd edge vertices>> 
std::map<SDEdge, SDVertex *> edgeVerts; for (SDFace *face : f) { for (int k = 0; k < 3; ++k) { <<Compute odd vertex on kth edge>> 
SDEdge edge(face->v[k], face->v[NEXT(k)]); SDVertex *vert = edgeVerts[edge]; if (!vert) { <<Create and initialize new odd vertex>> 
vert = arena.Alloc<SDVertex>(); newVertices.push_back(vert); vert->regular = true; vert->boundary = (face->f[k] == nullptr); vert->startFace = face->children[3];
<<Apply edge rules to compute new vertex position>> 
if (vert->boundary) { vert->p = 0.5f * edge.v[0]->p; vert->p += 0.5f * edge.v[1]->p; } else { vert->p = 3.f/8.f * edge.v[0]->p; vert->p += 3.f/8.f * edge.v[1]->p; vert->p += 1.f/8.f * face->otherVert(edge.v[0], edge.v[1])->p; vert->p += 1.f/8.f * face->f[k]->otherVert(edge.v[0], edge.v[1])->p; }
edgeVerts[edge] = vert; }
} }
<<Update new mesh topology>> 
<<Update even vertex face pointers>> 
for (SDVertex *vertex : v) { int vertNum = vertex->startFace->vnum(vertex); vertex->child->startFace = vertex->startFace->children[vertNum]; }
<<Update face neighbor pointers>> 
for (SDFace *face : f) { for (int j = 0; j < 3; ++j) { <<Update children f pointers for siblings>> 
face->children[3]->f[j] = face->children[NEXT(j)]; face->children[j]->f[NEXT(j)] = face->children[3];
<<Update children f pointers for neighbor children>> 
SDFace *f2 = face->f[j]; face->children[j]->f[j] = f2 ? f2->children[f2->vnum(face->v[j])] : nullptr; f2 = face->f[PREV(j)]; face->children[j]->f[PREV(j)] = f2 ? f2->children[f2->vnum(face->v[j])] : nullptr;
} }
<<Update face vertex pointers>> 
for (SDFace *face : f) { for (int j = 0; j < 3; ++j) { <<Update child vertex pointer to new even vertex>> 
face->children[j]->v[j] = face->v[j]->child;
<<Update child vertex pointer to new odd vertex>> 
SDVertex *vert = edgeVerts[SDEdge(face->v[j], face->v[NEXT(j)])]; face->children[j]->v[NEXT(j)] = vert; face->children[NEXT(j)]->v[j] = vert; face->children[3]->v[j] = vert;
} }
<<Prepare for next level of subdivision>> 
f = newFaces; v = newVertices;
} <<Push vertices to limit surface>> 
std::unique_ptr<Point3f[]> pLimit(new Point3f[v.size()]); for (size_t i = 0; i < v.size(); ++i) { if (v[i]->boundary) pLimit[i] = weightBoundary(v[i], 1.f / 5.f); else pLimit[i] = weightOneRing(v[i], loopGamma(v[i]->valence())); } for (size_t i = 0; i < v.size(); ++i) v[i]->p = pLimit[i];
<<Compute vertex tangents on limit surface>> 
std::vector<Normal3f> Ns; Ns.reserve(v.size()); std::vector<Point3f> pRing(16, Point3f()); for (SDVertex *vertex : v) { Vector3f S(0,0,0), T(0,0,0); int valence = vertex->valence(); if (valence > (int)pRing.size()) pRing.resize(valence); vertex->oneRing(&pRing[0]); if (!vertex->boundary) { <<Compute tangents of interior face>> 
for (int j = 0; j < valence; ++j) { S += std::cos(2 * Pi * j / valence) * Vector3f(pRing[j]); T += std::sin(2 * Pi * j / valence) * Vector3f(pRing[j]); }
} else { <<Compute tangents of boundary face>> 
S = pRing[valence - 1] - pRing[0]; if (valence == 2) T = Vector3f(pRing[0] + pRing[1] - 2 * vertex->p); else if (valence == 3) T = pRing[1] - vertex->p; else if (valence == 4) // regular T = Vector3f(-1 * pRing[0] + 2 * pRing[1] + 2 * pRing[2] + -1 * pRing[3] + -2 * vertex->p); else { Float theta = Pi / float(valence - 1); T = Vector3f(std::sin(theta) * (pRing[0] + pRing[valence - 1])); for (int k = 1; k < valence - 1; ++k) { Float wt = (2 * std::cos(theta) - 2) * std::sin((k) * theta); T += Vector3f(wt * pRing[k]); } T = -T; }
} Ns.push_back(Normal3f(Cross(S, T))); }
<<Create triangle mesh from subdivision mesh>> 
{ size_t ntris = f.size(); std::unique_ptr<int[]> verts(new int[3*ntris]); int *vp = verts.get(); size_t totVerts = v.size(); std::map<SDVertex *, int> usedVerts; for (size_t i = 0; i < totVerts; ++i) usedVerts[v[i]] = i; for (size_t i = 0; i < ntris; ++i) { for (int j = 0; j < 3; ++j) { *vp = usedVerts[f[i]->v[j]]; ++vp; } } return CreateTriangleMesh(ObjectToWorld, WorldToObject, reverseOrientation, ntris, verts.get(), totVerts, pLimit.get(), nullptr, &Ns[0], nullptr, nullptr); }
}

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.

<<Allocate LoopSubdiv vertices and faces>>= 
std::unique_ptr<SDVertex[]> verts(new SDVertex[nVertices]); for (int i = 0; i < nVertices; ++i) { verts[i] = SDVertex(p[i]); vertices.push_back(&verts[i]); } int nFaces = nIndices / 3; std::unique_ptr<SDFace[]> fs(new SDFace[nFaces]); for (int i = 0; i < nFaces; ++i) faces.push_back(&fs[i]);

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.

<<LoopSubdiv Local Structures>>= 
struct SDVertex { <<SDVertex Constructor>> 
SDVertex(const Point3f &p = Point3f(0, 0, 0)) : p(p) { }
<<SDVertex Methods>> 
int valence(); void oneRing(Point3f *p);
Point3f p; SDFace *startFace = nullptr; SDVertex *child = nullptr; bool regular = false, boundary = false; };

<<SDVertex Constructor>>= 
SDVertex(const Point3f &p = Point3f(0, 0, 0)) : p(p) { }

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.

Figure 3.27: Each triangular face stores three pointers to SDVertex objects v[i] and three pointers to neighboring faces f[i]. Neighboring faces are indexed using the convention that the ith edge is the edge from v[i] to v[(i+1)%3], and the neighbor across the ith edge is in f[i].

<<LoopSubdiv Local Structures>>+=  
struct SDFace { <<SDFace Constructor>> 
SDFace() { for (int i = 0; i < 3; ++i) { v[i] = nullptr; f[i] = nullptr; } for (int i = 0; i < 4; ++i) children[i] = nullptr; }
<<SDFace Methods>> 
int vnum(SDVertex *vert) const { for (int i = 0; i < 3; ++i) if (v[i] == vert) return i; Severe("Basic logic error in SDFace::vnum()"); return -1; } SDFace *nextFace(SDVertex *vert) { return f[vnum(vert)]; } SDFace *prevFace(SDVertex *vert) { return f[PREV(vnum(vert))]; } SDVertex *nextVert(SDVertex *vert) { return v[NEXT(vnum(vert))]; } SDVertex *prevVert(SDVertex *vert) { return v[PREV(vnum(vert))]; } SDVertex *otherVert(SDVertex *v0, SDVertex *v1) { for (int i = 0; i < 3; ++i) if (v[i] != v0 && v[i] != v1) return v[i]; Severe("Basic logic error in SDVertex::otherVert()"); return nullptr; }
SDVertex *v[3]; SDFace *f[3]; SDFace *children[4]; };

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.

<<LoopSubdiv Macros>>= 
#define NEXT(i) (((i) + 1) % 3) #define PREV(i) (((i) + 2) % 3)

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, normal v 0 and normal v 1 , with an edge between them. We expect that one of the triangular faces that has this edge will specify its three vertices so that normal v 0 is before normal v 1 , and that the other face will specify its vertices so that normal v 1 is before normal v 0 (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.

Figure 3.28: All of the faces in the input mesh must be specified so that each shared edge is given no more than once in each direction. Here, the edge from normal v 0 to normal v 1 is traversed from normal v 0 to normal v 1 by one face and from normal v 1 to normal v 0 by the other. Another way to think of this is in terms of face orientation: all faces’ vertices should be given consistently in either clockwise or counterclockwise order, as seen from outside the mesh.

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.

<<Set face to vertex pointers>>= 
const int *vp = vertexIndices; for (int i = 0; i < nFaces; ++i, vp += 3) { SDFace *f = faces[i]; for (int j = 0; j < 3; ++j) { SDVertex *v = vertices[vp[j]]; f->v[j] = v; v->startFace = f; } }

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.

<<LoopSubdiv Local Structures>>+= 
struct SDEdge { <<SDEdge Constructor>> 
SDEdge(SDVertex *v0 = nullptr, SDVertex *v1 = nullptr) { v[0] = std::min(v0, v1); v[1] = std::max(v0, v1); f[0] = f[1] = nullptr; f0edgeNum = -1; }
<<SDEdge Comparison Function>> 
bool operator<(const SDEdge &e2) const { if (v[0] == e2.v[0]) return v[1] < e2.v[1]; return v[0] < e2.v[0]; }
SDVertex *v[2]; SDFace *f[2]; int f0edgeNum; };

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 Subscript normal a , v Subscript normal b ) is correctly recognized as the same as the edge (v Subscript normal b , v Subscript normal a ), regardless of what order the vertices are provided in.

<<SDEdge Constructor>>= 
SDEdge(SDVertex *v0 = nullptr, SDVertex *v1 = nullptr) { v[0] = std::min(v0, v1); v[1] = std::max(v0, v1); f[0] = f[1] = nullptr; f0edgeNum = -1; }

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.

<<SDEdge Comparison Function>>= 
bool operator<(const SDEdge &e2) const { if (v[0] == e2.v[0]) return v[1] < e2.v[1]; return v[0] < e2.v[0]; }

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 upper O left-parenthesis log n right-parenthesis time.

<<Set neighbor pointers in faces>>= 
std::set<SDEdge> edges; for (int i = 0; i < nFaces; ++i) { SDFace *f = faces[i]; for (int edgeNum = 0; edgeNum < 3; ++edgeNum) { <<Update neighbor pointer for edgeNum>> 
int v0 = edgeNum, v1 = NEXT(edgeNum); SDEdge e(f->v[v0], f->v[v1]); if (edges.find(e) == edges.end()) { <<Handle new edge>> 
e.f[0] = f; e.f0edgeNum = edgeNum; edges.insert(e);
} else { <<Handle previously seen edge>> 
e = *edges.find(e); e.f[0]->f[e.f0edgeNum] = f; f->f[edgeNum] = e.f[0]; edges.erase(e);
}
} }

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.

<<Update neighbor pointer for edgeNum>>= 
int v0 = edgeNum, v1 = NEXT(edgeNum); SDEdge e(f->v[v0], f->v[v1]); if (edges.find(e) == edges.end()) { <<Handle new edge>> 
e.f[0] = f; e.f0edgeNum = edgeNum; edges.insert(e);
} else { <<Handle previously seen edge>> 
e = *edges.find(e); e.f[0]->f[e.f0edgeNum] = f; f->f[edgeNum] = e.f[0]; edges.erase(e);
}

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.

<<Handle new edge>>= 
e.f[0] = f; e.f0edgeNum = edgeNum; edges.insert(e);

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.

<<Handle previously seen edge>>= 
e = *edges.find(e); e.f[0]->f[e.f0edgeNum] = f; f->f[edgeNum] = e.f[0]; edges.erase(e);

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

Figure 3.29: Given a vertex v[i] and a face that it is incident to, f, we define the next face as the face adjacent to f across the edge from v[i] to v[NEXT(i)]. The previous face is defined analogously.

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.

Figure 3.30: We can determine if a vertex is a boundary vertex by starting from the adjacent face startFace and following next face pointers around the vertex. If we come to a face that has no next neighbor face, then the vertex is on a boundary. If we return to startFace, it’s an interior vertex.

<<Finish vertex initialization>>= 
for (int i = 0; i < nVertices; ++i) { SDVertex *v = vertices[i]; SDFace *f = v->startFace; do { f = f->nextFace(v); } while (f && f != v->startFace); v->boundary = (f == nullptr); if (!v->boundary && v->valence() == 6) v->regular = true; else if (v->boundary && v->valence() == 4) v->regular = true; else v->regular = false; }

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

<<LoopSubdiv Inline Functions>>= 
inline int SDVertex::valence() { SDFace *f = startFace; if (!boundary) { <<Compute valence of interior vertex>> 
int nf = 1; while ((f = f->nextFace(this)) != startFace) ++nf; return nf;
} else { <<Compute valence of boundary vertex>> 
int nf = 1; while ((f = f->nextFace(this)) != nullptr) ++nf; f = startFace; while ((f = f->prevFace(this)) != nullptr) ++nf; return nf + 1;
} }

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.

<<Compute valence of interior vertex>>= 
int nf = 1; while ((f = f->nextFace(this)) != startFace) ++nf; return nf;

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.

<<Compute valence of boundary vertex>>= 
int nf = 1; while ((f = f->nextFace(this)) != nullptr) ++nf; f = startFace; while ((f = f->prevFace(this)) != nullptr) ++nf; return nf + 1;

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.

<<SDFace Methods>>= 
int vnum(SDVertex *vert) const { for (int i = 0; i < 3; ++i) if (v[i] == vert) return i; Severe("Basic logic error in SDFace::vnum()"); return -1; }

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.

<<SDFace Methods>>+=  
SDFace *nextFace(SDVertex *vert) { return f[vnum(vert)]; }

<<SDFace Methods>>+=  
SDFace *prevFace(SDVertex *vert) { return f[PREV(vnum(vert))]; }

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

Figure 3.31: Given a vertex v on a face f, the method f->prevVert(v) returns the previous vertex around the face from v, and f->nextVert(v) returns the next vertex, where “next” and “previous” are defined by the original ordering of vertices when this face was defined.

<<SDFace Methods>>+=  
SDVertex *nextVert(SDVertex *vert) { return v[NEXT(vnum(vert))]; }

<<SDFace Methods>>+=  
SDVertex *prevVert(SDVertex *vert) { return v[PREV(vnum(vert))]; }

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.

<<Refine subdivision mesh into triangles>>= 
std::vector<SDFace *> f = faces; std::vector<SDVertex *> v = vertices; MemoryArena arena; for (int i = 0; i < nLevels; ++i) { <<Update f and v for next level of subdivision>> 
std::vector<SDFace *> newFaces; std::vector<SDVertex *> newVertices; <<Allocate next level of children in mesh tree>> 
for (SDVertex *vertex : v) { vertex->child = arena.Alloc<SDVertex>(); vertex->child->regular = vertex->regular; vertex->child->boundary = vertex->boundary; newVertices.push_back(vertex->child); } for (SDFace *face : f) { for (int k = 0; k < 4; ++k) { face->children[k] = arena.Alloc<SDFace>(); newFaces.push_back(face->children[k]); } }
<<Update vertex positions and create new edge vertices>> 
<<Update vertex positions for even vertices>> 
for (SDVertex *vertex : v) { if (!vertex->boundary) { <<Apply one-ring rule for even vertex>> 
if (vertex->regular) vertex->child->p = weightOneRing(vertex, 1.f / 16.f); else vertex->child->p = weightOneRing(vertex, beta(vertex->valence()));
} else { <<Apply boundary rule for even vertex>> 
vertex->child->p = weightBoundary(vertex, 1.f / 8.f);
} }
<<Compute new odd edge vertices>> 
std::map<SDEdge, SDVertex *> edgeVerts; for (SDFace *face : f) { for (int k = 0; k < 3; ++k) { <<Compute odd vertex on kth edge>> 
SDEdge edge(face->v[k], face->v[NEXT(k)]); SDVertex *vert = edgeVerts[edge]; if (!vert) { <<Create and initialize new odd vertex>> 
vert = arena.Alloc<SDVertex>(); newVertices.push_back(vert); vert->regular = true; vert->boundary = (face->f[k] == nullptr); vert->startFace = face->children[3];
<<Apply edge rules to compute new vertex position>> 
if (vert->boundary) { vert->p = 0.5f * edge.v[0]->p; vert->p += 0.5f * edge.v[1]->p; } else { vert->p = 3.f/8.f * edge.v[0]->p; vert->p += 3.f/8.f * edge.v[1]->p; vert->p += 1.f/8.f * face->otherVert(edge.v[0], edge.v[1])->p; vert->p += 1.f/8.f * face->f[k]->otherVert(edge.v[0], edge.v[1])->p; }
edgeVerts[edge] = vert; }
} }
<<Update new mesh topology>> 
<<Update even vertex face pointers>> 
for (SDVertex *vertex : v) { int vertNum = vertex->startFace->vnum(vertex); vertex->child->startFace = vertex->startFace->children[vertNum]; }
<<Update face neighbor pointers>> 
for (SDFace *face : f) { for (int j = 0; j < 3; ++j) { <<Update children f pointers for siblings>> 
face->children[3]->f[j] = face->children[NEXT(j)]; face->children[j]->f[NEXT(j)] = face->children[3];
<<Update children f pointers for neighbor children>> 
SDFace *f2 = face->f[j]; face->children[j]->f[j] = f2 ? f2->children[f2->vnum(face->v[j])] : nullptr; f2 = face->f[PREV(j)]; face->children[j]->f[PREV(j)] = f2 ? f2->children[f2->vnum(face->v[j])] : nullptr;
} }
<<Update face vertex pointers>> 
for (SDFace *face : f) { for (int j = 0; j < 3; ++j) { <<Update child vertex pointer to new even vertex>> 
face->children[j]->v[j] = face->v[j]->child;
<<Update child vertex pointer to new odd vertex>> 
SDVertex *vert = edgeVerts[SDEdge(face->v[j], face->v[NEXT(j)])]; face->children[j]->v[NEXT(j)] = vert; face->children[NEXT(j)]->v[j] = vert; face->children[3]->v[j] = vert;
} }
<<Prepare for next level of subdivision>> 
f = newFaces; v = newVertices;
} <<Push vertices to limit surface>> 
std::unique_ptr<Point3f[]> pLimit(new Point3f[v.size()]); for (size_t i = 0; i < v.size(); ++i) { if (v[i]->boundary) pLimit[i] = weightBoundary(v[i], 1.f / 5.f); else pLimit[i] = weightOneRing(v[i], loopGamma(v[i]->valence())); } for (size_t i = 0; i < v.size(); ++i) v[i]->p = pLimit[i];
<<Compute vertex tangents on limit surface>> 
std::vector<Normal3f> Ns; Ns.reserve(v.size()); std::vector<Point3f> pRing(16, Point3f()); for (SDVertex *vertex : v) { Vector3f S(0,0,0), T(0,0,0); int valence = vertex->valence(); if (valence > (int)pRing.size()) pRing.resize(valence); vertex->oneRing(&pRing[0]); if (!vertex->boundary) { <<Compute tangents of interior face>> 
for (int j = 0; j < valence; ++j) { S += std::cos(2 * Pi * j / valence) * Vector3f(pRing[j]); T += std::sin(2 * Pi * j / valence) * Vector3f(pRing[j]); }
} else { <<Compute tangents of boundary face>> 
S = pRing[valence - 1] - pRing[0]; if (valence == 2) T = Vector3f(pRing[0] + pRing[1] - 2 * vertex->p); else if (valence == 3) T = pRing[1] - vertex->p; else if (valence == 4) // regular T = Vector3f(-1 * pRing[0] + 2 * pRing[1] + 2 * pRing[2] + -1 * pRing[3] + -2 * vertex->p); else { Float theta = Pi / float(valence - 1); T = Vector3f(std::sin(theta) * (pRing[0] + pRing[valence - 1])); for (int k = 1; k < valence - 1; ++k) { Float wt = (2 * std::cos(theta) - 2) * std::sin((k) * theta); T += Vector3f(wt * pRing[k]); } T = -T; }
} Ns.push_back(Normal3f(Cross(S, T))); }
<<Create triangle mesh from subdivision mesh>> 
{ size_t ntris = f.size(); std::unique_ptr<int[]> verts(new int[3*ntris]); int *vp = verts.get(); size_t totVerts = v.size(); std::map<SDVertex *, int> usedVerts; for (size_t i = 0; i < totVerts; ++i) usedVerts[v[i]] = i; for (size_t i = 0; i < ntris; ++i) { for (int j = 0; j < 3; ++j) { *vp = usedVerts[f[i]->v[j]]; ++vp; } } return CreateTriangleMesh(ObjectToWorld, WorldToObject, reverseOrientation, ntris, verts.get(), totVerts, pLimit.get(), nullptr, &Ns[0], nullptr, nullptr); }

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.

Figure 3.32: Basic Subdivision of a Single Triangular Face. Four child faces are created, ordered such that the ith child face is adjacent to the ith 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 ith edge vertex is along the ith edge of the original face.

<<Update f and v for next level of subdivision>>= 
std::vector<SDFace *> newFaces; std::vector<SDVertex *> newVertices; <<Allocate next level of children in mesh tree>> 
for (SDVertex *vertex : v) { vertex->child = arena.Alloc<SDVertex>(); vertex->child->regular = vertex->regular; vertex->child->boundary = vertex->boundary; newVertices.push_back(vertex->child); } for (SDFace *face : f) { for (int k = 0; k < 4; ++k) { face->children[k] = arena.Alloc<SDFace>(); newFaces.push_back(face->children[k]); } }
<<Update vertex positions and create new edge vertices>> 
<<Update vertex positions for even vertices>> 
for (SDVertex *vertex : v) { if (!vertex->boundary) { <<Apply one-ring rule for even vertex>> 
if (vertex->regular) vertex->child->p = weightOneRing(vertex, 1.f / 16.f); else vertex->child->p = weightOneRing(vertex, beta(vertex->valence()));
} else { <<Apply boundary rule for even vertex>> 
vertex->child->p = weightBoundary(vertex, 1.f / 8.f);
} }
<<Compute new odd edge vertices>> 
std::map<SDEdge, SDVertex *> edgeVerts; for (SDFace *face : f) { for (int k = 0; k < 3; ++k) { <<Compute odd vertex on kth edge>> 
SDEdge edge(face->v[k], face->v[NEXT(k)]); SDVertex *vert = edgeVerts[edge]; if (!vert) { <<Create and initialize new odd vertex>> 
vert = arena.Alloc<SDVertex>(); newVertices.push_back(vert); vert->regular = true; vert->boundary = (face->f[k] == nullptr); vert->startFace = face->children[3];
<<Apply edge rules to compute new vertex position>> 
if (vert->boundary) { vert->p = 0.5f * edge.v[0]->p; vert->p += 0.5f * edge.v[1]->p; } else { vert->p = 3.f/8.f * edge.v[0]->p; vert->p += 3.f/8.f * edge.v[1]->p; vert->p += 1.f/8.f * face->otherVert(edge.v[0], edge.v[1])->p; vert->p += 1.f/8.f * face->f[k]->otherVert(edge.v[0], edge.v[1])->p; }
edgeVerts[edge] = vert; }
} }
<<Update new mesh topology>> 
<<Update even vertex face pointers>> 
for (SDVertex *vertex : v) { int vertNum = vertex->startFace->vnum(vertex); vertex->child->startFace = vertex->startFace->children[vertNum]; }
<<Update face neighbor pointers>> 
for (SDFace *face : f) { for (int j = 0; j < 3; ++j) { <<Update children f pointers for siblings>> 
face->children[3]->f[j] = face->children[NEXT(j)]; face->children[j]->f[NEXT(j)] = face->children[3];
<<Update children f pointers for neighbor children>> 
SDFace *f2 = face->f[j]; face->children[j]->f[j] = f2 ? f2->children[f2->vnum(face->v[j])] : nullptr; f2 = face->f[PREV(j)]; face->children[j]->f[PREV(j)] = f2 ? f2->children[f2->vnum(face->v[j])] : nullptr;
} }
<<Update face vertex pointers>> 
for (SDFace *face : f) { for (int j = 0; j < 3; ++j) { <<Update child vertex pointer to new even vertex>> 
face->children[j]->v[j] = face->v[j]->child;
<<Update child vertex pointer to new odd vertex>> 
SDVertex *vert = edgeVerts[SDEdge(face->v[j], face->v[NEXT(j)])]; face->children[j]->v[NEXT(j)] = vert; face->children[NEXT(j)]->v[j] = vert; face->children[3]->v[j] = vert;
} }
<<Prepare for next level of subdivision>> 
f = newFaces; v = newVertices;

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.

<<Allocate next level of children in mesh tree>>= 
for (SDVertex *vertex : v) { vertex->child = arena.Alloc<SDVertex>(); vertex->child->regular = vertex->regular; vertex->child->boundary = vertex->boundary; newVertices.push_back(vertex->child); } for (SDFace *face : f) { for (int k = 0; k < 4; ++k) { face->children[k] = arena.Alloc<SDFace>(); newFaces.push_back(face->children[k]); } }

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.

<<Update vertex positions and create new edge vertices>>= 
<<Update vertex positions for even vertices>> 
for (SDVertex *vertex : v) { if (!vertex->boundary) { <<Apply one-ring rule for even vertex>> 
if (vertex->regular) vertex->child->p = weightOneRing(vertex, 1.f / 16.f); else vertex->child->p = weightOneRing(vertex, beta(vertex->valence()));
} else { <<Apply boundary rule for even vertex>> 
vertex->child->p = weightBoundary(vertex, 1.f / 8.f);
} }
<<Compute new odd edge vertices>> 
std::map<SDEdge, SDVertex *> edgeVerts; for (SDFace *face : f) { for (int k = 0; k < 3; ++k) { <<Compute odd vertex on kth edge>> 
SDEdge edge(face->v[k], face->v[NEXT(k)]); SDVertex *vert = edgeVerts[edge]; if (!vert) { <<Create and initialize new odd vertex>> 
vert = arena.Alloc<SDVertex>(); newVertices.push_back(vert); vert->regular = true; vert->boundary = (face->f[k] == nullptr); vert->startFace = face->children[3];
<<Apply edge rules to compute new vertex position>> 
if (vert->boundary) { vert->p = 0.5f * edge.v[0]->p; vert->p += 0.5f * edge.v[1]->p; } else { vert->p = 3.f/8.f * edge.v[0]->p; vert->p += 3.f/8.f * edge.v[1]->p; vert->p += 1.f/8.f * face->otherVert(edge.v[0], edge.v[1])->p; vert->p += 1.f/8.f * face->f[k]->otherVert(edge.v[0], edge.v[1])->p; }
edgeVerts[edge] = vert; }
} }

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.

<<Update vertex positions for even vertices>>= 
for (SDVertex *vertex : v) { if (!vertex->boundary) { <<Apply one-ring rule for even vertex>> 
if (vertex->regular) vertex->child->p = weightOneRing(vertex, 1.f / 16.f); else vertex->child->p = weightOneRing(vertex, beta(vertex->valence()));
} else { <<Apply boundary rule for even vertex>> 
vertex->child->p = weightBoundary(vertex, 1.f / 8.f);
} }

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 beta (Figure 3.33). The vertex we are updating, in the center, is weighted by 1 minus n beta , where n is the valence of the vertex. Thus, the new position v prime for a vertex v is

normal v prime equals left-parenthesis 1 minus n beta right-parenthesis normal v plus sigma-summation Underscript i equals 1 Overscript upper N Endscripts beta normal v Subscript i Baseline period

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.

Figure 3.33: The new position v prime for a vertex v is computed by weighting the adjacent vertices v Subscript i by a weight beta and weighting v by left-parenthesis 1 minus n beta right-parenthesis , where n is the valence of v. The adjacent vertices v Subscript i are collectively referred to as the one-ring around v.

The specific weight beta 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  beta value based on the vertex’s valence that ensures smoothness. For regular interior vertices, beta() returns 1 slash 16 . Since this is a common case, the implementation uses 1 slash 16 directly instead of calling beta() every time.

<<Apply one-ring rule for even vertex>>= 
if (vertex->regular) vertex->child->p = weightOneRing(vertex, 1.f / 16.f); else vertex->child->p = weightOneRing(vertex, beta(vertex->valence()));

<<LoopSubdiv Inline Functions>>+=  
inline Float beta(int valence) { if (valence == 3) return 3.f / 16.f; else return 3.f / (8.f * valence); }

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.

<<LoopSubdiv Function Definitions>>+=  
static Point3f weightOneRing(SDVertex *vert, Float beta) { <<Put vert one-ring in pRing>> 
int valence = vert->valence(); Point3f *pRing = ALLOCA(Point3f, valence); vert->oneRing(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.

<<Put vert one-ring in pRing>>= 
int valence = vert->valence(); Point3f *pRing = ALLOCA(Point3f, valence); vert->oneRing(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.

<<LoopSubdiv Function Definitions>>+=  
void SDVertex::oneRing(Point3f *p) { if (!boundary) { <<Get one-ring vertices for interior vertex>> 
SDFace *face = startFace; do { *p++ = face->nextVert(this)->p; face = face->nextFace(this); } while (face != startFace);
} else { <<Get one-ring vertices for boundary vertex>> 
SDFace *face = startFace, *f2; while ((f2 = face->nextFace(this)) != nullptr) face = f2; *p++ = face->nextVert(this)->p; do { *p++ = face->prevVert(this)->p; face = face->prevFace(this); } while (face != nullptr);
} }

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

<<Get one-ring vertices for interior vertex>>= 
SDFace *face = startFace; do { *p++ = face->nextVert(this)->p; face = face->nextFace(this); } while (face != startFace);

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.

<<Get one-ring vertices for boundary vertex>>= 
SDFace *face = startFace, *f2; while ((f2 = face->nextFace(this)) != nullptr) face = f2; *p++ = face->nextVert(this)->p; do { *p++ = face->prevVert(this)->p; face = face->prevFace(this); } while (face != nullptr);

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 Subscript 1 and v Subscript 2 to compute the new position v prime  as

normal v prime equals left-parenthesis 1 minus 2 beta right-parenthesis normal v plus beta normal v 1 plus beta normal v 2 period

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

Figure 3.34: Subdivision on a Boundary Edge. The new position for the vertex in the center is computed by weighting it and its two neighbor vertices by the weights shown.

<<Apply boundary rule for even vertex>>= 
vertex->child->p = weightBoundary(vertex, 1.f / 8.f);

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.

<<LoopSubdiv Function Definitions>>+= 
static Point3f weightBoundary(SDVertex *vert, Float beta) { <<Put vert one-ring in pRing>> 
int valence = vert->valence(); Point3f *pRing = ALLOCA(Point3f, valence); vert->oneRing(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.

Figure 3.35: Subdivision Rule for Edge Split. The position of the new odd vertex, marked with an open circle, is found by weighting the two vertices at the ends of the edge and the two vertices opposite it on the adjacent triangles. (a) The weights for an interior vertex; (b) the weights for a boundary vertex.

<<Compute new odd edge vertices>>= 
std::map<SDEdge, SDVertex *> edgeVerts; for (SDFace *face : f) { for (int k = 0; k < 3; ++k) { <<Compute odd vertex on kth edge>> 
SDEdge edge(face->v[k], face->v[NEXT(k)]); SDVertex *vert = edgeVerts[edge]; if (!vert) { <<Create and initialize new odd vertex>> 
vert = arena.Alloc<SDVertex>(); newVertices.push_back(vert); vert->regular = true; vert->boundary = (face->f[k] == nullptr); vert->startFace = face->children[3];
<<Apply edge rules to compute new vertex position>> 
if (vert->boundary) { vert->p = 0.5f * edge.v[0]->p; vert->p += 0.5f * edge.v[1]->p; } else { vert->p = 3.f/8.f * edge.v[0]->p; vert->p += 3.f/8.f * edge.v[1]->p; vert->p += 1.f/8.f * face->otherVert(edge.v[0], edge.v[1])->p; vert->p += 1.f/8.f * face->f[k]->otherVert(edge.v[0], edge.v[1])->p; }
edgeVerts[edge] = vert; }
} }

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.

<<Compute odd vertex on kth edge>>= 
SDEdge edge(face->v[k], face->v[NEXT(k)]); SDVertex *vert = edgeVerts[edge]; if (!vert) { <<Create and initialize new odd vertex>> 
vert = arena.Alloc<SDVertex>(); newVertices.push_back(vert); vert->regular = true; vert->boundary = (face->f[k] == nullptr); vert->startFace = face->children[3];
<<Apply edge rules to compute new vertex position>> 
if (vert->boundary) { vert->p = 0.5f * edge.v[0]->p; vert->p += 0.5f * edge.v[1]->p; } else { vert->p = 3.f/8.f * edge.v[0]->p; vert->p += 3.f/8.f * edge.v[1]->p; vert->p += 1.f/8.f * face->otherVert(edge.v[0], edge.v[1])->p; vert->p += 1.f/8.f * face->f[k]->otherVert(edge.v[0], edge.v[1])->p; }
edgeVerts[edge] = vert; }

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.

<<Create and initialize new odd vertex>>= 
vert = arena.Alloc<SDVertex>(); newVertices.push_back(vert); vert->regular = true; vert->boundary = (face->f[k] == nullptr); vert->startFace = face->children[3];

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 3 slash 8 , and the two vertices opposite the edge are given weight 1 slash 8 (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.

<<Apply edge rules to compute new vertex position>>= 
if (vert->boundary) { vert->p = 0.5f * edge.v[0]->p; vert->p += 0.5f * edge.v[1]->p; } else { vert->p = 3.f/8.f * edge.v[0]->p; vert->p += 3.f/8.f * edge.v[1]->p; vert->p += 1.f/8.f * face->otherVert(edge.v[0], edge.v[1])->p; vert->p += 1.f/8.f * face->f[k]->otherVert(edge.v[0], edge.v[1])->p; }

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

<<SDFace Methods>>+= 
SDVertex *otherVert(SDVertex *v0, SDVertex *v1) { for (int i = 0; i < 3; ++i) if (v[i] != v0 && v[i] != v1) return v[i]; Severe("Basic logic error in SDVertex::otherVert()"); return nullptr; }

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.

Figure 3.36: Each face is split into four child faces, such that the ith child is adjacent to the ith vertex of the original face, and such that the ith child face’s ith vertex is the child of the ith vertex of the original face. The vertices of the center child are oriented such that the ith vertex is the odd vertex along the ith edge of the parent face.

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.

<<Update new mesh topology>>= 
<<Update even vertex face pointers>> 
for (SDVertex *vertex : v) { int vertNum = vertex->startFace->vnum(vertex); vertex->child->startFace = vertex->startFace->children[vertNum]; }
<<Update face neighbor pointers>> 
for (SDFace *face : f) { for (int j = 0; j < 3; ++j) { <<Update children f pointers for siblings>> 
face->children[3]->f[j] = face->children[NEXT(j)]; face->children[j]->f[NEXT(j)] = face->children[3];
<<Update children f pointers for neighbor children>> 
SDFace *f2 = face->f[j]; face->children[j]->f[j] = f2 ? f2->children[f2->vnum(face->v[j])] : nullptr; f2 = face->f[PREV(j)]; face->children[j]->f[PREV(j)] = f2 ? f2->children[f2->vnum(face->v[j])] : nullptr;
} }
<<Update face vertex pointers>> 
for (SDFace *face : f) { for (int j = 0; j < 3; ++j) { <<Update child vertex pointer to new even vertex>> 
face->children[j]->v[j] = face->v[j]->child;
<<Update child vertex pointer to new odd vertex>> 
SDVertex *vert = edgeVerts[SDEdge(face->v[j], face->v[NEXT(j)])]; face->children[j]->v[NEXT(j)] = vert; face->children[NEXT(j)]->v[j] = vert; face->children[3]->v[j] = vert;
} }

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.

<<Update even vertex face pointers>>= 
for (SDVertex *vertex : v) { int vertNum = vertex->startFace->vnum(vertex); vertex->child->startFace = vertex->startFace->children[vertNum]; }

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.

<<Update face neighbor pointers>>= 
for (SDFace *face : f) { for (int j = 0; j < 3; ++j) { <<Update children f pointers for siblings>> 
face->children[3]->f[j] = face->children[NEXT(j)]; face->children[j]->f[NEXT(j)] = face->children[3];
<<Update children f pointers for neighbor children>> 
SDFace *f2 = face->f[j]; face->children[j]->f[j] = f2 ? f2->children[f2->vnum(face->v[j])] : nullptr; f2 = face->f[PREV(j)]; face->children[j]->f[PREV(j)] = f2 ? f2->children[f2->vnum(face->v[j])] : nullptr;
} }

For the first step, recall that the interior child face is always stored in children[3]. Furthermore, the monospace k plus 1 st child face (for monospace k equals 0 comma 1 comma 2 ) is across the kth edge of the interior face, and the interior face is across the monospace k plus 1 st edge of the kth face.

<<Update children f pointers for siblings>>= 
face->children[3]->f[j] = face->children[NEXT(j)]; face->children[j]->f[NEXT(j)] = face->children[3];

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.

<<Update children f pointers for neighbor children>>= 
SDFace *f2 = face->f[j]; face->children[j]->f[j] = f2 ? f2->children[f2->vnum(face->v[j])] : nullptr; f2 = face->f[PREV(j)]; face->children[j]->f[PREV(j)] = f2 ? f2->children[f2->vnum(face->v[j])] : nullptr;

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

<<Update face vertex pointers>>= 
for (SDFace *face : f) { for (int j = 0; j < 3; ++j) { <<Update child vertex pointer to new even vertex>> 
face->children[j]->v[j] = face->v[j]->child;
<<Update child vertex pointer to new odd vertex>> 
SDVertex *vert = edgeVerts[SDEdge(face->v[j], face->v[NEXT(j)])]; face->children[j]->v[NEXT(j)] = vert; face->children[NEXT(j)]->v[j] = vert; face->children[3]->v[j] = vert;
} }

For the kth child face (for monospace k equals 0 comma 1 comma 2 ), 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.

<<Update child vertex pointer to new even vertex>>= 
face->children[j]->v[j] = face->v[j]->child;

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.

<<Update child vertex pointer to new odd vertex>>= 
SDVertex *vert = edgeVerts[SDEdge(face->v[j], face->v[NEXT(j)])]; face->children[j]->v[NEXT(j)] = vert; face->children[NEXT(j)]->v[j] = vert; face->children[3]->v[j] = vert;

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:

<<Prepare for next level of subdivision>>= 
f = newFaces; v = newVertices;

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 1 slash 5 and the center vertex by 3 slash 5 . The rule for interior vertices is based on a function loopGamma(), which computes appropriate vertex weights based on the valence of the vertex.

<<Push vertices to limit surface>>= 
std::unique_ptr<Point3f[]> pLimit(new Point3f[v.size()]); for (size_t i = 0; i < v.size(); ++i) { if (v[i]->boundary) pLimit[i] = weightBoundary(v[i], 1.f / 5.f); else pLimit[i] = weightOneRing(v[i], loopGamma(v[i]->valence())); } for (size_t i = 0; i < v.size(); ++i) v[i]->p = pLimit[i];

<<LoopSubdiv Inline Functions>>+= 
inline Float loopGamma(int valence) { return 1.f / (valence + 3.f / (8.f * beta(valence))); }

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.

<<Compute vertex tangents on limit surface>>= 
std::vector<Normal3f> Ns; Ns.reserve(v.size()); std::vector<Point3f> pRing(16, Point3f()); for (SDVertex *vertex : v) { Vector3f S(0,0,0), T(0,0,0); int valence = vertex->valence(); if (valence > (int)pRing.size()) pRing.resize(valence); vertex->oneRing(&pRing[0]); if (!vertex->boundary) { <<Compute tangents of interior face>> 
for (int j = 0; j < valence; ++j) { S += std::cos(2 * Pi * j / valence) * Vector3f(pRing[j]); T += std::sin(2 * Pi * j / valence) * Vector3f(pRing[j]); }
} else { <<Compute tangents of boundary face>> 
S = pRing[valence - 1] - pRing[0]; if (valence == 2) T = Vector3f(pRing[0] + pRing[1] - 2 * vertex->p); else if (valence == 3) T = pRing[1] - vertex->p; else if (valence == 4) // regular T = Vector3f(-1 * pRing[0] + 2 * pRing[1] + 2 * pRing[2] + -1 * pRing[3] + -2 * vertex->p); else { Float theta = Pi / float(valence - 1); T = Vector3f(std::sin(theta) * (pRing[0] + pRing[valence - 1])); for (int k = 1; k < valence - 1; ++k) { Float wt = (2 * std::cos(theta) - 2) * std::sin((k) * theta); T += Vector3f(wt * pRing[k]); } T = -T; }
} Ns.push_back(Normal3f(Cross(S, T))); }

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 w Subscript i .

Figure 3.37: To compute tangents for interior vertices, the one-ring vertices are weighted with weights w Subscript i . The center vertex, where the tangent is being computed, always has a weight of 0.

To compute the first tangent vector bold s , the weights are

w Subscript i Baseline equals cosine left-parenthesis StartFraction 2 pi i Over n EndFraction right-parenthesis comma

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

w Subscript i Baseline equals sine left-parenthesis StartFraction 2 pi i Over n EndFraction right-parenthesis period

<<Compute tangents of interior face>>= 
for (int j = 0; j < valence; ++j) { S += std::cos(2 * Pi * j / valence) * Vector3f(pRing[j]); T += std::sin(2 * Pi * j / valence) * Vector3f(pRing[j]); }

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.

Figure 3.38: Tangents at boundary vertices are also computed as weighted averages of the adjacent vertices. However, some of the boundary tangent rules incorporate the value of the center vertex.

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

bold s equals normal v Subscript n minus 1 Baseline minus normal v 0 period

The second tangent, known as the transverse tangent, is computed based on the vertex’s valence. The center vertex is given a weight w Subscript normal c and the one-ring vertices are given weights specified by a vector left-parenthesis w 0 comma w 1 comma ellipsis comma w Subscript n minus 1 Baseline right-parenthesis . The transverse tangent rules we will use are

Valence w Subscript normal c w Subscript i
2 negative 2 left-parenthesis 1 comma 1 right-parenthesis
3 negative 1 left-parenthesis 0 comma 1 comma 0 right-parenthesis
4 (regular) negative 2 left-parenthesis negative 1 comma 2 comma 2 comma negative 1 right-parenthesis

For valences of 5 and higher, w Subscript normal c Baseline equals 0 and

StartLayout 1st Row 1st Column w 0 2nd Column equals w Subscript n minus 1 Baseline equals sine theta 2nd Row 1st Column w Subscript i 2nd Column equals left-parenthesis 2 cosine theta minus 2 right-parenthesis sine left-parenthesis theta i right-parenthesis comma EndLayout

where

theta equals StartFraction pi Over n minus 1 EndFraction period

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

<<Compute tangents of boundary face>>= 
S = pRing[valence - 1] - pRing[0]; if (valence == 2) T = Vector3f(pRing[0] + pRing[1] - 2 * vertex->p); else if (valence == 3) T = pRing[1] - vertex->p; else if (valence == 4) // regular T = Vector3f(-1 * pRing[0] + 2 * pRing[1] + 2 * pRing[2] + -1 * pRing[3] + -2 * vertex->p); else { Float theta = Pi / float(valence - 1); T = Vector3f(std::sin(theta) * (pRing[0] + pRing[valence - 1])); for (int k = 1; k < valence - 1; ++k) { Float wt = (2 * std::cos(theta) - 2) * std::sin((k) * theta); T += Vector3f(wt * pRing[k]); } T = -T; }

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.