Efﬁcient Mesh Generation Using Subdivision Surfaces

Polygonal meshes and particularly triangular meshes are the most used structure for 3D modelling. The ‘di-rect edges’ data structure is the most efﬁcient way to represent them and subdivision surfaces is an appropri-112


INTRODUCTION
We needed an adequate method to represent the surface of a 3D virtual dog skull that will be part of a bigger project [1]. Our implementation of a parametric triangular mesh generation algorithm helped us to validate the technique and to confirm our choice.
A surface separates the interior from the exterior by a boundary. In order to generate a mathematical surface we need to establish continuity and neighbourhood relationships between samples. Complex shapes require piecewise representations and are split in sub-regions each one of which is defined by an individual function. These representations are just approximations of the surface and it can be demonstrated from the Taylor theorem that the approximation error in an interval h of the surface by a polynomial of degree p is O(h p+1 ). In order to decrease the error one can increase the degree of the polynomial or use more segments or sub-regions. Operations required by surfaces include evaluation, query and modification. Evaluation refers to the sampling of the surface and its attributes. Query aims to identify if a given point p R 3 is inside or outside the solid bounded by S. Modification refers to geometry changes such as surface deformations or topology changes. There are two major classes of representation. Parametric representations are defined by a vector valued parameterization function (1). Implicit or volumetric representations are defined by a zero set of a scalar valued function (2) [2]. Parametric representations are better for evaluation and modification, implicit representations are better for query [2].
f: Ω S, where S is a surface, Ω R2 and S = f (Ω) R3. (1) Polygonal meshes approximate a surface by a mesh of planar polygonal facets. They are low complexity representations more efficient for rendering. In particular, triangular meshes remain the most used structure for 3D modelling. In [3] they argue that quadrilaterals are better than triangles capturing the symmetries of objects and that they are compatible with bi cubic patches used in commercial software. However, in [2] they state that triangles have become increasingly popular because they allow more flexible and efficient processing and avoid conversion errors. We decided to follow the approach from [2] and use triangular meshes. Additionally, we found that subdivision surfaces are an appropriate way to generate and create parametric polygonal meshes. They reduce the representation of a complex surface to a simpler control mesh that is able to generate refined surfaces.

METHODOLOGY
In order to find an appropriate representation, we reviewed literature on 3D modelling. From this, we identified the advantages of using triangular meshes. Also, we studied and created an algorithm to convert from the 3DS file format to the direct edges data structure which we found to be the most efficient one. Then, from the review of subdivision surfaces we chose the √3 subdivision method for generating the mesh. At this stage, we had to find the right formulas for counting elements in the mesh (i.e. vertices and faces) in order to create an efficient implementation. We created and downloaded several input meshes in the 3DS format, converted them and tested our algorithm. The input meshes included closed regular surfaces, planes with boundaries, surfaces with non subdividable polygons and the arbitrary mesh of a spaceship. We calculated memory usage, performance and element growth. The implementation will be further used in the creation of the virtual dog head. • Access and enumeration of individual vertices, edges, faces.

TRIANGULAR MESHES
• Oriented traversal of edges of a face (next edge in a face).
• Access to at least one face attached to a given vertex.
The direct edges data structure is the most efficient to deal with triangular meshes [2]. It is based on indices as references to each element where indexing follows rules that implicitly encode connectivity information (Table 1). Here, each edge is represented into two opposing halfedges consistently oriented counter clockwise ( Figure 2). This data structure is only useful for triangular meshes and provides no explicit representation of edges (though it can be specialised to other polygons). It groups the three halfedges belonging to a common triangle. Additional information is explicitly stored in arrays. For instance, for each vertex, the index of an outgoing halfedge is stored and; for each halfedge, the index of its opposite half edge and the index of the vertex the halfedge points to are stored. Boundaries are managed with special negative indices that indicate that the edge vertex is invalid.

SUBDIVISION SURFACES
Subdivision surfaces help in the creation and refinement of proper 3D mesh models. They are capable of representing an arbitrary geometrically unrestricted topology. They produce an efficient hierarchical structure and an object can be modelled as a low resolution control mesh from which we can generate new meshes by refining the previous one. Some subdivision methods such as Catmull and Clark [4] work on quadrilaterals or extend subdivision to a n-sided problem and are not restricted to triangles such as Doo and Sabin [5, 6]. Others such as butterfly [4], √3 Subdivision [6] and Loop [8] are specialised to triangles. [8] Has been widely used [9,10,11].
One can define: M 0 (P 0 , K 0 ) as the original coarse mesh (which can be used as the control mesh) and M j (P j , K j ) as the j times subdivided mesh.

Indices and operation
Description / calculation Here, P j are the mesh points at level j of subdivision and K j contains the description of the topology at level j. M j , with j is the approximation of a B spline limit surface. A subdivision scheme S takes the vertices from level j to level j+1 so that S(K j ) = K j + 1 . A subdivision matrix or stencil SM maps a mesh M j to a topologically equivalent refined mesh M j+1 . Eigenvalues or characteristic roots are a special set of scalars associated with linear systems of equations such as matrix equations. In a subdivision scheme S, the eigenvalues of the subdivision matrix are important to determine if the method converges to a limit smooth surface. For instance, all subdivision schemes must guarantee adequate design of the SM stencil so that eigenvalues have a certain distribution and a continuous surface approximating the limit surface can be generated [7,9]. In SM, every row is a rule to compute the position of a new vertex and every column tells how one old vertex contributes to the vertex positions in the refined mesh [2]. In stationary methods the refinement functions are the same for every subdivision level. However, different rules are applied to define sharp features, creases or to deal with extraordinary vertices [9]. Also, stop conditions can be used for adaptive refinement [7]. In some methods old vertices are repositioned for smoothing [7].

Subdivision techniques for triangular meshes
Butterfly subdivision [4,9] converts each face into four using edge vertices. Artefacts and discontinuities are produced in vertices with valence different from 6. For example, it never produces smooth surfaces on extraordinary vertices and incorrect smooth regions can appear near high valence vertices due to eigenclustering (more than one eigenvalue per matrix). It is based on dyadic split (divide every edge in two) and does not produce C 2 continuity. The modified butterfly deals with artefacts of the original butterfly [4]. It approximates the behaviour of a C 2 surface based on dyadic split. Here, eigenclustering is avoided. Loop subdivision [9] is one of the simplest C 1 subdivision schemes. Here, stencils are smaller and convergence rate is better than in modified butterfly. It is based on dyadic split and is not C 2 at extraordinary vertices. The scheme does not look smooth for large valences, due to eigenclustering. The modified Loop [9,12] increases the stencil in a minimal way around a vertex and avoids eigenclustering. Here, vertices have better structure at extraordinary points. However, ripples appear. The ternary Loop subdivision [13] uses three stencils. It achieves bounded curvature, manifold support, and convex Hull. It is C 2 continuous.
However, No rules have been defined to deal with boundaries or with sharp features.
√3 Subdivision [7] is a stationary subdivision scheme with slower topological refinement with trisection of every original edge (every two steps). It inserts a new vertex at the centre of every face (Figure 1a and 1b). Then, it creates the new faces ( Figure 1c) and flips every original edge (Figure 1d). To do that, it uses simple stencils of minimum size and maximum symmetry. Here, a new vertex is calculated as the average of the three old ones (a new vertex only affects one face). Then, the old ones are relaxed using (6). The scheme uses a generation index to perform adaptive refinement and allows sharp feature lines. It is C 2 continuous except at extraordinary points. All new vertices have valence 6 and the valence of the old ones is not changed. We chose this method because it produces more levels of subdivision with lower number of triangles and simpler rules.

Converting from 3DS to Direct Edges
Using a 3DS format loader [14] we created an algorithm that converts a 3DS mesh into a Direct Edges data structure. In our algorithm, only the vertex list and faces descriptions were processed preventing unnecessary duplication. Our application opens the 3DS file. Next, it reads the vertices and adds them to a vertex list. Then, it reads the faces and adds them to a face list assigning the proper vertex indices. It creates the halfedges list from the faces list. For each one it assigns the vertex it points to and, for each vertex it assigns one outgoing halfedge. Next, for each halfedge it finds and assigns its opposite. Finally, it detects boundary vertices and saves the data.

√3 subdivison surfaces using the Direct Edges Data Structure
In our implementation a main subdivision cycle calculates the subdivision surface until it reaches the number of levels defined by the user. At each level, the algorithm performs three tasks. First, it calculates the new vertices in each face (Table 2). Here, if it is a non boundary subdividable face or a boundary face in even level (Figure 2c), the mid point is obtained and three faces are created. If the face is not subdividable then one face is recreated (Figure 2a). If the face is a boundary face and the level of subdivision is odd, two new vertices are calculated and three faces are created (Figure 2d). After the faces have been subdivided a mesh similar to the shown in Figure 1c is produced.
The next task (Table 2) is to flip the edges in each face of the new subdivided mesh. Here, the edges are flipped only if the face and its mate are sub dividable and non boundary or, when the face belongs to a boundary at an odd level of subdivision.
After the edges have been flipped a mesh similar to the shown in Figure  1d is produced. The final task is to re-position the three old vertices in each face. Calculations are done taking care of avoiding boundary vertices, because these can cause visible discontinuities. Here, the calculation of the vertex's neighbourhood is required (Table 2).

Memory allocation
Every time the calculations for a new level of subdivision start, the memory is allocated according to the new needs. The result of the last subdivision step is stored in a file that is used as a starting point for the next level.
The use of a file eliminated the need of using heavy intermediate memory objects. At level 0, memory for the control mesh is allocated according to the number of vertices, faces and halfedges which were defined as in Table 3. Memory for subsequent levels of subdivisions is calculated according to (7) (8) and (9).

Assigning halfedges and counting
Our main challenge was to keep track of the halfedges after each task. Assigning proper halfedges numbers is very important before the flipping operation, otherwise the mesh will be corrupted. The naive solution is to create costly loops to find and assign halfedge numbers. However, counting helps to use the direct edges data structure efficiently. The first and easiest thing is to calculate the interior halfedges using formula (3). For example, a face F has three interior halfedges. These are assigned counter clock wise starting from vertex a. For instance: • Halfedge 0 goes from a to b.
• Halfedge 1 goes from b to c.
• Halfedge 2 goes from c to a.   Figure 2c and 1, 5, 4, 8 in Figure 2d) is straightforward. Here, halfedge 1 of one new face is next to halfedge 2 of the other one (11). Opposite halfedges to the 0 halfedge of new boundary faces are assigned -1 (non existent).

Swapping faces
In Figure 3a, face F [a, b, c] is going to be swapped with face M [b, a, d].
Let us assume that F is a regular face. Swapping cannot be achieved if: • M has already been swapped in the current subdivision level.
• M is a boundary face.
• M is a non sub dividable face.
If swapping is possible the old opposite halfedge values HE1, HE2, HE3 and HE4 have to be stored and reassigned in order to keep consistency and connectivity. Then, the two faces are swapped as shown with the dotted line and the new 0 halfedge is assigned along this line for each face in this new arrangement. The new faces are described in Table 4. A boundary face in an odd level is swapped only if the same conditions for regular swapping are met (except that M could have been swapped); only the leftmost and rightmost faces (Figure 2d) are candidates for swapping, the middle face remains the same. Here (Figure 2b)

Detecting boundaries and subdividable faces
A boundary face is detected by checking if it contains at least one -1 valued opposite halfedge. Boundary vertices are detected as a last step in the conversion from the 3DS to the direct edges format. Here, a flag is assigned to each boundary vertex. In Figure 2a: • a is a boundary vertex if the opposite of halfdegde 0 or 2 is next to a boundary face.
• b is a boundary vertex if the opposite of halfdegde 0 or 1 is next to a boundary face.
• c is a boundary vertex if the opposite of halfdegde 1 or 2 is next to a boundary face.
Not keeping proper track of the boundary vertices may produce distortion. The detection of a sub dividable face is yet to be implemented. Under certain flatness criteria (based on differential geometry) a face could be tagged as non sub dividable during the conversion process.

Calculating a neighbourhood
This calculation returns the number of neighbours of a vertex and their average coordinate value. In Figure  3b vertex c has 4 neighbours (a, b, e, d) and their average is obtained with (16). Our algorithm deals with closed manifolds, boundaries and distorted polygons. Let us assume that the first detected outgoing halfedge in Figure   3b is HE 2. First, we add vertex a to the average and increment the neighbourhood. Then, we find the opposite halfedge (in this case 10). If it is the 2 halfedge of the next face we subtract 2 in order to find the halfedge which points to the next vertex. If it is the 0 or 1 halfedge of the next face we add 1 in order to find the halfedge which points to the next vertex. Then, we add this vertex to the average and increment the number of neighbours. The algorithm continues until the halfedge pointing to the next vertex is equal to the first halfedge or until an edge or a distortion is detected ( Table 2).

RESULTS
Our direct edges implementation uses fewer elements. While in the 3DS format a cube contains 26 * vertices, in the direct edges format it only contains 8. While in the 3DS format a spaceship contains 649 * vertices, in the direct edges format it only contains 260. The number of faces remained the same in both cases. The size of one vertex, one polygon and one halfedge (Table 3) can be easily calculated as a float and an integer occupy 32 bits each. According to (7), (8) and (9)

DISCUSSION
Parametric triangular meshes allow efficient validation (used for display) and modification (especially when defined with subdivision surfaces). We decided to use them, because due to their popularity, many methods for processing them have been developed. Coarse meshes may bring distorted triangles and edges (i.e. "nonmanifold meshes are problematic for most algorithms, since around non-manifold configurations there exist no well-defined local geodesic neighbourhood"), [2, p.19]. We manage them with the definition of proper stop conditions. Our algorithm can be applied to an arbitrary mesh ( Figure 5). The understanding of the direct edges data structure allowed efficient implementation and reduced the number of elements used in the original 3DS format. Counting was essential in memory allocation and in keeping track of the halfedges for each new face. It allowed the effective implementation of the basic operations required by the subdivision algorithm.
Subdivision surfaces have a major advantage, which is providing a control mesh (which will facilitate deformation) and a refined triangular mesh (which will ease display). Their use will bring C 2 continuity to most regions of the mesh, producing a more natural look [2]. The benefits of the √3 subdivision brought to our attention a better and newer subdivision scheme which is supported by less complex stencils and has a simpler way to deal with sharp features, boundaries and with adaptive refinement. In √3 subdivision growth is slower than in other subdivision techniques. The inclusion of non sub dividable faces lowers the number of new generated faces. Apart from boundary faces, the growth is still exponential (Figure 3). However, the √3 subdivision produces fewer triangles than other methods providing similar quality and the mesh produced is visually appealing with a few levels of subdivision.
One simple way to define edges or sharp features is to pre assign -1 valued opposite halfedges to the region of interest. We want to propose preserving flat regions with smaller number of triangles. For instance, differently from the original algorithm we want to avoid subdividing well defined flat regions during adaptive refinement. To achieve this, our algorithm still requires a flatness function able to predetermine what faces are non sub dividable. Londra [1] will be a virtual dog capable of displaying facial expressions. She will use a skull model implemented with a polygonal mesh generation / representation technique suitable for deformation, parameterization and animation, because the skull shape will change according to conformational anatomical parameters. From our implementation and the tests made, we have found that √3 subdivision surfaces are an appropriate modelling technique for the skull. Additionally, our skull will include non penetration features, requiring the use of an alternate implicit model (implicit models are better for query operations). Here, the idea will be to create new dog breed heads from new skull shapes obtained altering the model through its control mesh that will produce the new refined meshes. To account for head type transitions, Londra will not require topological changes, only geometric ones. In consequence, the complexities involved in the former will be avoided.