Fig. 1. (from left to right): input mesh; mesh with maximal error e below 0.5%, 0.2%, 0.1% and the resulting surface for e=0.1%.
Xiaobin Wu and Jörg Peters
University of Florida, USA
xwu@cise.ufl.edu
Abstracts:
A tight estimate on the maximum distance between a subdivision surface and its linear approximation is introduced to guide adaptive subdivision with guaranteed accuracy.
Subdivision provides a simple and powerful method for modeling free-form surfaces: given a polygonal input mesh, a sequence of refinements generates an ever denser mesh with a generically smooth limit surface. For modern graphics applications, input meshes can consist of thousands of faces. If, at each step of refinement, every mesh facet is split into a fixed number of new faces, the number of facets increases exponentially and the complexity of the mesh quickly exceeds the memory and processing limitations. The obvious answer, adaptive refinement, requires a good bound on how well planar triangles approximate the limit surface. Loose bounds waste resources and overly aggressive approximations can miss surface features. Due to the procedural definition of the subdivision surface, adaptive refinement is trickier than for standard spline surfaces.
In this paper, we leverage a new bound on the maximum distance between the limit surface and its linear approximation. The bound can be computed locally and efficiently, and yields a tight estimate with an error converging to zero under subdivision. Figure 1 shows adaptive meshes for a given threshold e and Figure 2 compares the new bound to conventional ones.
Fig. 1. (from left to right): input mesh; mesh with maximal error e below 0.5%, 0.2%, 0.1% and the resulting surface for e=0.1%.
Fig. 2. Comparison of bounds (in purple). (A) A subdivision patch (shaded) with its local control mesh (black)
(B) Axis-aligned bounding box (C) Kobbelt’s estimate (D) Our estimate
1.1. Previous work
A number of error estimates have been used for adaptive subdivision: sampling the distance between mesh node and its limit [8], oriented bounding boxes (directional distance between the limit surface and the interpolation of its three corners) [6, 5], axis aligned bounding boxes (to detect self intersection) [3]. Also various planarity tests have been used to guide adaptivity. M¨uller et al. used the angle between the limit normal of a vertex and the normals of its adjoining faces [9, 10], Zorin et al. measure planarity at corners and edge midpoints of the next subdivision level [15], Xu and Kondo [14] compute the angle difference between the normal of a face and that of all its neighboring faces.
A comparison of the surface max-norm bounds and our new bound is shown in Figure 2. Note that sampling and planarity tests do not yield explicit guaranteed maximal bounds between the mesh-polyhedron and the limit surface. Such an error guarantee, however, is crucial for manufacturing applications.
1.2. Overview
Our algorithm is as follows. For each facet of the input mesh, for each coordinate x, y, z, say the x-component, we compute two linear functions px and mx such that px (u, v) ≤ x(u, v) ≤ mx(u, v) . The x-error is ex =px−mx and the maximum error ||(ex, ey, ez)|| between the subdivision patch and the facet guides the refinement.
In Section 2, we give the details of each step above except for the important detail of parametrization discussed in Section 3. The implementation and the data structures are summarized in Section 4.
This section describes the central contribution of the paper: efficiently determining the maximum distance between for subdivision mesh, viewed as a polyhedron, and the limit surface.
We take advantage of locality. Just like NURBS patches, a mesh node (control point) only locally influences the surface shape. Conversely, any point on the surface also depends only on a small submesh. We will focus on Loop’s subdivision scheme [7] – it is not difficult to apply the ideas to other subdivision schemes. Loop subdivision is based on a triangle mesh. Each triangle, together with all its direct neighbor triangles, defines a curved triangular surface patch as depicted in 2 (A). Our goal is to approximate the maximum error between the patch and the triangle.
The difficulty is a two-fold. First, it is hard to compute the Hausdorff distance. Second, the limit surface is defined by recursive procedure and lacks a closed form representation. The idea is to parametrize the patch, bound each of the component functions: x(u, v), y(u, v), z(u, v) and combine the component bounds to an overall bound.
We can defer the detailed discussion of the (u, v) parametrization until the next section. For now, we assume that we have a description of x(u, v), y(u, v), z(u, v). In particular, after at most one (local) subdivision, we may assume that at most one of the three vertices of the center triangle (corresponding to the patch) has a valence n ≠ 6.
2.1. Bounding the Distance
To measure the maximum error between x(u, v) and its linear approximation, we decompose x(u, v) into a linear combination of n + 6 functions bi(u, v).
Each bi(u, v) is the limit, under Loop subdivision, of a local control net (vj , xj) 
R3 with vj = (uj , vj) abscissae laid out and indexed as in Figure 3 and all xj = 0 except for xj = 1 if j = i. The values of (uj , vj) we chose and two alternatives are discussed in Section 3.
Fig. 3. Indices of the vertices vj of the (u, v) domain mesh for n = 5 and n = 7; v0 is the vertex with valence n.
v1 and v2 are the other two vertices of the center triangle and have 6 neighbors.
By linear precision of Loop subdivision, the maximum distance between the lower bound mx and the upper bound px is zero for linear functions. Therefore we can remove a linear function l from the above equation without changing the maximum error:
where l interpolates (v0, x0), (v1, x1), (v2, x2) and di is the vertical difference between l and (vi, xi). By interval arithmetic, we obtain an upper and a lower bound of x(u, v) as the follows:
The error ex is bounded by the maximum difference:
In order to efficiently compute ex, we tabulate the error bound Ei := e(bi) for each basis function bi (the linear bound is defined by three scalar coefficients), as well as the barycentric coordinates (si, ti) of vi, i = 3..n + 5 with respect to the vertices v0, v1, v2 of the central triangle, i.e.
vi = (siv0 + tiv1 + (1 − si − ti)v2).
Then the distance di can be computed efficiently as
di = xi − (six0 + tix1 + (1 − si − ti)x2).
The following short algorithm sums up our process of computing the maximum error. Here input cf is a one dimensional array, containing one component of the n+6 control points. Variables s, t, E are pre-tabulated 1-dimensional floating point arrays of the values si, ti and Ei described above.
For each patch, we call this function three times to get ex, ey and ez. If ||< ex, ey, ez >|| < ε, the patch is rendered using a single triangle. Otherwise, we subdivide this patch into four subpatches and continue with the subpatches.
3. Parametrization
In this section, we discuss and compare three (u, v) parametrizations of the surface patches, i.e. the layout of the mesh points vi defining the parameters of each component. The extra-ordinary node with valence n ≠ 6 is v0. The comparison in Table 1 establishes that the carefully constructed, exact parametrization is slightly more efficient than the uniform parametrization and dramatically better than the binary parametrization.
3.1. Exact Parametrization
The parametrization suggested in [13] reproduces linear functions and is defined by the following construction (cf. Figure 3):
1. Set v0 = 0, the origin of the (u, v) plane.
2. The direct neighbors vi of v0 form a regular unit n−gon.
3. Extend the edge v0v1 and v0v2 by kn to get v4 and v6.
4. v5 is the average of v4 and v6.
5. v3, v7 are the reflection of v5, across v0v4 and v0v6, respectively.
kn is defined by the following formula where c := cos (π/n) :
We can now obtain any patch x(u, v) by subdividing the control net with vertices (vi, xi). The shaded areas in Figure 3, which depend on the valence of v0, is the limit of the vi-mesh under Loop subdivision. Note that the limit of the mesh with vertices vj are the shaded areas in Figure 3. Note that these areas fit exactly into the center triangle.
3.2. Uniform Parametrization
With the same construction as above, but kn=1 for all n, we obtain the uniform parametrization. The bottom boundary of the domain will either, for n < 6, pull back from the boundary of the center triangle or, for n > 6, push out of the triangle (see Figure 4) and therefore requires careful extrapolation to safeguard the bound.
Fig. 4. The domains of uniform parametrization for n=3 and n=10.
3.2. Binary Parametrization
The parametrization proposed in [11, 12] associates the vertices of the mesh under subdivision with a binarily refined grid. While this allows deducing the number of subdivision steps from the (u, v) position, it does not yield linear precision! This means, we cannot pull out the linear term l from the expansion of x(u, v) and cancel it when we compute the distance between px and mx. Therefore, asfor axis-aligned bounding boxes, the relative position of the object in space influences the local error and the error estimates can be dramatically worse as illustrated in Table 1. The vertex based method [8] by measuring the distance between the vertices and their limit positions is also listed in the Table 1 as ”v-error”.
Table1. Numbers of resulting triangles for different error measurements. Rows ”exact”, ”uniform”, ”binary” indicate three different parametrizations.
Row ”v-error” is for the method based on measuring the distance between the vertices and their limit points. ‘n/a’ indicates out of memory.
We support the adaptive subdivision process by a forest of balanced quad trees. Balancing according to [1], assures that neighboring faces differ only by one step of subdivision. The quad trees of this balanced adaptive subdivision, one tree per original facet, are linked so that we can directly access the neighbors for any given face at any level and without ascending the tree to the common parent (herein lies the difference to [15]).
The data structure has the following entries (cf. Figure5).
Fig. 5. Data structure. (left) The indices of the parent face and its 4 children. (right) Halfedge pointers for fast neighboring access.
Now, for any leaf face f, we can access its neighbors via the halfedge [n, i] in constant time. If the face n is a leaf, then it is the neighbor of f on side i. Otherwise, its children i and (i + 1 mod 3) are the two neighbors to f on side i. Due to tree balancing, they must be leaf nodes. In C code, our data structure is as follows:
The field entries can be filled by a depth first traversal of the quadtree forest. The time to build the data structure is therefore O(N) where N is number of leaf faces.
4.1. Gap Prevention
When two neighboring faces have different subdivision levels, a gap appears between them as illustrated in Figure 6 left. Rendering such a mesh results in irritating visual artifacts. Moreover, gaps spell serious trouble for the mesh processing and finite element computations.
Fig. 6. (left) A gap between adjacent faces at different levels of subdivision. (right) The gap is removed by splitting the coarse triangle.
We follow the standard recipe for removing the gaps by splitting the patch on the coarser side as illustrated in Figure 6 right. Such a process is done top to down recursively for each pair of original neighboring faces. The time complexity of this step is also O(N) where N is number of leaf faces.
5.1. View dependent adaptation
So far our focus was on true error control in 3 dimensions. For rendering, we can project our error vectors (ex, ey, ez) into the view plane and compute the projected size as adaptive refinement criterion.
5.1. Creases and Boundaries
Subdivision surfaces become a lot more exciting with Pixar’s semi-sharp creases ([2, 4]): a few steps of smoothing in only one direction followed by smooth subdivision. A simple way to avoid special bounds for the different crease configurations is to locally perform the uni-directional subdivision steps regardless of the error and then apply the regular bounds. Alternatively, since we only use upper and lower bounds, we need not generate new tables for every combination of two subdivision rules: if one rule results consistently in higher values than the other, say a unidirectional rule and a generic subdivision rule, we bound the upper function above and the lower function below to enclose the range of combinations. The same technique applies to boundaries that follow a spline curve.
In Figure 7, the 500-facet deer model is subdivided to meet error bounds of 0.5% and 0.1%, respectively. The total times to subdivide the model, including the error estimation, are 172ms and 906ms, respectively. The number of triangles generated are 4834 and 25980.
Fig. 7. (from left to right) input; e=0.5%; e=0.1%; and the surface with e=0.1%.
In Figure 8, we raytrace the subdivided mesh, on a P4 2.8G PC with 1G RAM. The total time to raytrace the input cat model of 671 triangles is 8 seconds. The time to adaptively subdivide the model to within a 0.2% bound and raytrace is 9 seconds. The total time to uniformly subdivide the model to the same bound and raytrace is 12 seconds.
Fig. 8. Ray-traced images at 800x600 resolution. (left) Input mesh (671 triangles). Ray-tracing time: 8s.
(middle) Adaptive subdivision with e=0.2%. Resulting number of triangles = 9662 and maximumsubdivision level: 4. Ray-tracing time: 9s.
(right) uniformly subdivided 4 times. Resulting number of triangles = 171776. Ray-tracing time: 12s.
We presented a tight estimate on the maximum distance between a subdivision surface and its linear approximation and applied it to adaptive subdivision. The computation is efficient and the algorithm is easy to implement. We are currently applying the approach to Catmull-Clark surfaces and develop a similar tight bound for patch normals to help silhouette and self-interference detection.
Appendix A. Convergence of the error estimate
We need to show that ex→0 under subdivision where 
. Since the terms e(bi) are constant, it suffices to prove di → 0 where di := si(xi − x0) + ti(xi − x1) + (1 − si − ti)(xi − x2). This follows from (xi − xj) → 0 under subdivision.
[1] M. Berg, M. Kreveld, O. Overmars, and O. Schwarzkopf. Computational Geometry: Algorithms and Applications. Springer-Verlag, Berlin, Heidelberg, New York, 1997.
[2] T. DeRose, M. Kass, and T. Truong. Subdivision surfaces in character animation. In Siggraph 98, pages 85–94, 1998.
[3] E. Grinspun and P. Schröder. Normal bounds for subdivisionsurface interference detection. In Proc Visualization, pages 333–340. IEEE, 2001.
[4] H. Hoppe, T. DeRose, T. Duchamp, M. Halstead, H. Jin, J. McDonald, J. Schweitzer, and W. Stuetzle. Piecewise smooth surface reconstruction. Computer Graphics, 28(Annual Conference Series):295–302, July 1994.
[5] L. Kobbelt. Tight bounding volumes for subdivision surfaces. In Pacific-Graphics’98, pages 17–26. IEEE, 1998.
[6] L. Kobbelt, K. Daubert, and H.-P. Seidel. Ray tracing of subdivision surfaces. In Rendering Techniques ’98 (Proceedings of the Eurographics Workshop), pages 69–80. Springer-Verlag, 1998.
[7] C. T. Loop. Smooth subdivision surfaces based on triangles, 1987. Master’s Thesis, Department of Mathematics, University of Utah.
[8] H. M¨uller and R. Jaeschke. Adaptive subdivision curves and surfaces. In Proceedings of the Conference on Computer Graphics International 1998 (CGI-98), pages 48–58, Los Alamitos, California, June 22–26 1998. IEEE Computer Society.
[9] K. M¨uller and S. Havemann. Subdivision surface tessellation on the fly using a versatile mesh data structure. In Computer Graphics Forum (Eurographics 2000), volume 19(3), 2000.
[10] V. Settgast, K. M¨uller, F. Fuenfzig, and D. Fellner. Adaptive tesselation of subdivision surfaces. Computers and Graphics, 28(1):73–78, Feb. 2004.
[11] J. Stam. Exact evaluation of catmull-clark subdivision surfaces at arbitrary parameter values. In SIGGRAPH 98 Proceedings, pages 395–404. Addison Wesley, 1998.
[12] J. Stam. Evaluation of loop subdivision surfaces, Aug. 27 1999.
[13] X. Wu and J. Peters. Interference detection for subdivision surfaces. Computer Graphics Forum, Eurographics 2004, 23(3):577–585, 2004.
[14] Z. Xu and K. Kondo. Local subdivision process with Doo-Sabin subdivision surfaces. In ShapeModeling International, Proceedings, 2002.
[15] D. Zorin, P. Schröder, and S. Sweldens. Interactive multiresolution mesh editing. In SIGGRAPH 97 Conference Proceedings, Annual Conference Series, pages 259–268. ACM SIGGRAPH, Addison Wesley, Aug. 1997. ISBN 0-89791-896-7.