This page describes the subdivision algorithm that will be implemented in JPatch. It is “open for discussion”.
Catmull Clark Subdivision Surfaces
Here’s the algorithm for Catmull Clark subdivision:
When a polyhedron is subdivided with the Catmull-Clark method, new vertices (called “face points”) are placed at the center of each original face, new “edge points” are similarly placed at the center of each original edge, and then new edges are added to connect the new edge points to the new adjacent face points. The positions of the vertices are calculated as follows:
- The face points are positioned as the average of the positions of the face’s original vertices;
- The edge point locations are calculated as the average of the center point of the original edge and the average of the locations of the two new adjacent face points;
- The old vertices are repositioned according to the equation
, where
- Q is the average of the new face points surrounding the old vertex,
- R is the average of the midpoints of the edges that share the old vertex,
- S is the old vertex point, and
- n is the number of edges that share the old vertex.
A naive approach would be to use a plain old “wing” or “half-wing” data structure to represent the surface and recursively apply the subdivision algorithm on this structure. This approach has several problems:
- The entire mesh needs to be in memory. For higher levels of subdivision this is a severe problem.
- Since the algorithm works on the entire mesh, it is not possible to perform adaptive tessellation.
- “Wing” or “Half-wing” data structures are very flexible, but tend to be slow due to the many pointer operations involved.
The only way to solve the first problem is to apply the subdivision algorithm to individual faces. SDS have “local control”, that means the limit surface of each face is fully determined by the control points of the original face plus its 1-neighbors (the points surrounding the face).
But the algorithm would still suffer from the slow operations on the “half-wing” data structure. Additionally, at each step new vertices and faces are created (and thus have to be allocated in heap memory), which further decreases performance because of the memory allocation and eventual garbage collection.
So I was looking for an algorithm that uses pre-allocated arrays (up to a given maximum subdivision depth) that would be used during subdivision. This way there is no need for memory allocation or garbage when performing the subdivision. The primary idea for the following algorithm is taken from the paper Adaptive Tessellation of Subdivision Surfaces.
Observations
After one level of Catmull Clark subdivision, the entire mesh consists only of 4-sided faces. Each vertex, though, can be of different valence (the number faces adjacent to the vertex). The following image shows a typical situation:
The left image shows one face to be subdivided (orange) plus its surrounding faces (gray). The top left vertex (of the orange face) is of valence 5, the top right vertex is of valence 3 and the two bottom vertices are of valence 4. Vertices of valence 4 are called regular, all other vertices are called irregular.
The middle image shows the situation after one level of subdivision, the right image after 2 levels of subdivision.
These images show that the vertices can be arranged in a 2 dimensional grid, with extra arrays for the corners.
The size of the array for storing the grid is 
, where n is the subdivision level. This array would be sufficient if all vertices were regular. For irregular vertices, separate arrays are needed. Their size is 
, where v is the valence. To keep things simple, the arrays that are needed to store the coordinates for each subdivision level are pre-allocated. There is only one such geometry array per subdivision level. It starts with elements for the corner vertices (space is allocated for a given maximum valence), followed by elements for the quadratic grid.
With this setup, we now need an algorithm that takes the array containing the vertices of level n and fills the array containing the vertices of level n + 1. It turns out that much of the information needed to perform this subdivision can be precomputed as well:
Stencils and weights
Each face-, edge- and (regular) vertex-point of level n is defined by 4, 6 or 9 vertices of level n - 1, respectively. This information can be stored for each vertex as stencil.
This image shows the stencils for face-, edge and regular vertex-points, together with the weights needed to compute the position of the new vertex. Different stencils and weights are used to compute crease or corner points. There are also stencils that can be used to project a point onto the limit surface and stencils to compute the tangents.
So, in addition to the arrays needed to store the vertex position, some (precomputed) arrays containing the stencils are needed. A stencil is an integer array, containing the indices of the vertices in the geometry array (of the previous level). Thus, a stencil array is an array of integer arrays (and contains multiple stencils). For each vertex of level n it contains a stencil with pointers to the vertices of level n - 1 that are needed to compute the position of this vertex. For each subdivision level there is one large stencil array for the quadratic grid (called patch stencil array), four stencil arrays (one for each corner, called fan stencil arrays) and one for each patch-corner vertex (called corner stencil. This is necessary because the patch-corners might be irregular).
This image shows the arrangement of the patch stencil array and the four fan stencil arrays. The colors of the dots indicate the type of stencil (face, edge or vertex), a white dot is an unused stencil. Note that the four stencils in each corner of the patch stencil array are unused (these could be omitted, but keeping them in the patch stencil array ensures a simple mapping of stencil array index to geometry array index, which could make it easier to run this algorithm as a GPU kernel). The stencils and weights for the remaining four corner vertices are currently computed ad-hoc.
To perform subdivision, the algorithm first loops over the patch stencil array. For each stencil it looks up the geometry array of the previous level (using the indices from the stencils), multiplies the positions with the pre-defined weights and stores the calculated position in the geometry array of this level. After that it loops over each fan stencil array, and finally the positions of the 4 remaining (possibly irregular) vertices at the patch corners are computed.
Algorithm
- One level of subdivision is performed using the “half wing” representation of the mesh. This results in a number of quadrilateral facettes. The vertices of each facette store information of how they have been computed, thus if an original vertex was moved, the level-1 vertices can be re-evaluated quickly. Each facette also stores arrays with pointers to the surrounding vertices (needed for subdivision). This step is only performed when the topology has been changed.
- Given the current modelview and projection matrices, the desired subdivision depth for each facette is computed (it can depend on the estimated screen size of the facette, its curvature and on a global render-quality setting).
- Each facette is subdivided until its subdivision depth (computed in step 2) has been reached. After that, the resulting vertices are projected to the limit surface using the limit stencils and normals are computed using the tangent stencils and computing the vector cross product of the two tangents. The limit points and normals are copied into arrays that can be rendered using glDrawArrays or glDrawElements.
- Taking account of the subdivision levels of the neighbor facettes, the facettes are stitched together to produce a watertight mesh.
Future improvements
I’ve implemented the algorithm described above and so far it works well and is pretty fast (implemented in Java and running on the CPU).
The next steps are:
- Add stencils for corners and creases.
- Allow hierarchical subdivision meshes.
One way to subdivide hierarchical meshes would be to store level 2+ vertices in the facette they belong to. When subdividing, these points are copied into the geometry array after each regular subdivision step (and thus overrides the results of the regular subdivision). Similarly there could be crease and corner information stored for higher levels, these had to be copied into the stencil arrays. One problem is that information from neighbor facettes is needed as well, but I think this can be solved.
Eventually it would be nice to implement this algorithm as a GPU kernel. The geometry and stencil arrays would be represented as textures, and a fragment shader would create the geometry arrays for the next subdivision level (by rendering into an auxiliary buffer that would be used as “texture” in the next step). At the end the resulting limit-positions and normals would be “rendered” into a buffer that’s used as VertexBuffer (to actually draw the surface) in the final step.