Ternary Subdivision for Quadrilateral Meshes - CiteSeerX
Ternary Subdivision for Quadrilateral Meshes Tianyun Ni1 , Ahmad H. Nasri2 , and J¨org Peters1 1
2
Dept. CISE, University of Florida Dept. of Computer Science American University of Beirut
Abstract. A well-documented problem of Catmull and Clark subdivision surfaces is that, in the neighborhood of extraordinary points, the curvature is unbounded and fluctuates. In fact, since one of the eigenvalues that determines elliptic shape is too small, the limit surface can have a saddle point when the designer’s input mesh suggests a convex shape. Here, we replace, near the extraordinary point, Catmull-Clark subdivision by another set of rules based on refining each bi-cubic B-spline into nine. This provides many localized degrees of freedom for special rules that need not reach out to neighbor vertices in order to tune the behavior. In this paper, we provide a strategy for setting such degrees of freedom and exhibit tuned ternary quad subdivision that yields surfaces with bounded curvature, nonnegative weights and full contribution of elliptic and hyperbolic shape components. Keywords: Subdivision, ternary, bounded curvature, convex hull
1 Introduction Tuning a subdivision scheme means adjusting the subdivision rules or stencils to obtain a refined mesh and surface with prescribed properties [1, 2, 4, 6, 7, 11, 12]. To date, no stationary subdivision algorithm on quadrilateral meshes is known to yield bounded curvature while guaranteeing the convex hull property at the extraordinary nodes. This paper explores ternary refinement to obtain such a scheme, because ternary subdivision offers more parameters in a close vicinity of each extraordinary node to tune the subdivision than binary subdivision does and thereby allows localizing the tuning [6]. Ternary subdivision inserts two knots in each knot interval, splitting each quad into nine. If all nodes of a quad are of valence 4, the rules for bi-cubic B-splines shown in Figure 3 are applied and we have C 2 continuity. The challenge is with nodes of valence other than 4, called extraordinary nodes. In particular, we can devise rules at the extraordinary nodes and their newly created direct and diagonal neighbors to achieve – eigenvalues in order of magnitude as 1, λ, λ, λ2 , λ2 , λ2 , λi , . . . and 0 ≤ |λi | < λ2 for i = 7, .., 2n + 1; and – nonnegative weights. It is possible to set, in addition, λ = 1/3. Yet, Figure 2 shows that the macroscopic shape of the new scheme and Catmull-Clark are very similar (Figure 2), but, of course, the mesh is denser. In the following, we first recall the rules for regular binary and ternary subdivisions (Section 1.1). In section 2, we state the properties that the ternary subdivision should match. Section 3 analyzes the spectrum in terms of free parameters and Section 4 derives weights by interpolating the regular case.
2
. Fig. 1. Catmull-Clark subdivision (top) and ternary subdivision (bottom). (left top) Control net with vertices (x, y, 50(1 − x2 − 5y 2 )), where (x, y) are the positions of the nodes of the characteristic map of Catmull-Clark subdivision. (middle three) control net at subdivision steps 3,4,5 and (right) shaded limit surface.
Fig. 2. Three steps of Catmull-Clark subdivision (left) and ternary subdivision (right). While the density differes, the macroscopic shape is identical.
3
1.1 Background and Regular Case Subdivision is a widely adopted tool in computer graphics and is also making inroads into geometric modeling, if only for conceptual modeling. However, [5] pointed out that Catmull-Clark subdivision in particular, is lacking the full set of elliptic and hyperbolic subsubdominant eigenfunctions and therefore typically generates saddle shapes in the limit at vertices with valence greater than four (see Figure 1 top). This implies that any high-quality (standard, symmetric) scheme needs to have a spectrum 1, λ, λ, λ2 , λ2 , λ2 followed by smaller terms. With the strategy explained in the following, it is possible to achieve such a spectrum by making subdivision stencils near the extraordinary point depend on the valence (Figure 1 bottom). A major challenge is technical: there are either too few or far too many free parameters occurring in nonlinear inequality and equality constraints. This makes any selection by general optimization [2] impractical. In [6,
1
6
1
1
1
16 1
6
36
6
6
1
6 Vertex
1
1
6 1 Edge
1
76
16
4
64
40
10
16 0
100
1 76
361
76
19
304
16
76 Vertex
16
4
64 Edge
190 16
256
160
16
10
1 Face
40
1
Face
Fig. 3. The regular bicubic B-spline refinement vertex, edge, and facet refinement stencils of Catmull-Clark (left) and ternary quad subdivision (right).
7], Loop proposed improvements to his well-known triangular scheme to achieve both bounded curvature for binary and ternary schemes. Based on fixed-degree interpolation of known weights, this approach sets many unknown weights and reduces the number of free parameters. It turns out to be far more difficult on Catmull-Clark case since the leading eigenvalues now come from a 2 by 2 diagonal block matrix rather than from a 1 by 1 block as the Loop subdivision. The eigenvalues are therefore the roots of quadratic polynomials while, in the case of triangular meshes, the eigenvalue is simple. The problem for n-sided facets is even more complex due to 4 × 4 blocks but can be avoided by applying one initial Catmull-Clark subdivision step. Most recently, two other approaches to tuning subdivision schemes at extraordinary points have been proposed. Augsdoerfer et al [1] choose an optimal subdominant eigenvalue of the subdivision scheme to minimize the curvature variation of the central surface [10] and Ginkel and Umlauf [4] propose new techniques to avoid hybrid configurations by tuning eigenvalues and eigencoefficients. These two and the present approach have in common that they focus on the curvature around the extraordinary point. In all cases, the localized improvement may, but is not guaranteed to, yield overall curvature distribution is of high-quality.
4
2 Problem Statement
We want to derive a ternary, quad subdivision scheme that is stationary, (rotational and p-mirror) symmetric and affine invariant. We denote the leading eigenvalues, in order of magnitude as, 1, λ, λ, µ1 , µ2 , µ3 . To achieve bounded curvature and the convex hull property, we enforce the following constraints on the eigenvalues and the weights αi , βi , γi , δi and vi of the ternary quad scheme (see Figure 4): We focus on the leading
(i) C 1 scheme 1 > λ > µi , and the characteristic map is regular and injective. (ii) Bounded Curvature |µ1 | = |µ2 | = |µ3 | = λ2 . (iii) Convex Hull αi ≥ 0, βi ≥ 0, γi ≥ 0, δi ≥ 0, i = 0..n − 1, v1 , v2 ≥ 0, P Pn−1 e0 := 1 − n−1 i=0 (αi + βi ) ≥ 0, f0 := 1 − i=0 (γi + δi ) ≥ 0, P v0 := 1 − n 2i=1 (vi ) ≥ 0. (iv) Symmetry αl = αn−l , βl = βn−1−l , γl = γ(n+1−l) mod n , and δl = δn−l Table 1. Constraints on the spectrum and the stencil weights. In addition, we will be able to enforce µ1 = µ2 = µ3 > λ7 > 0 and the subsubdominant eigenvectors corresponding to µi come from Fourier blocks with index 0, 2, n − 2.
eigenvalues and hence the 1-ring of control mesh around extraordinary point and its eigenstructure. As shown in Figure 4, each refinement step generates three types of points corresponding to vertices, edges and faces respectively. Let n be the valence of a generic extraordinary point. We have, offhand, 4n + 2 weights to determine. Symmetry reduces this to at most 2n + 4 free parameters.
α2
v1
v1
δ1
β1
v2
γ2
α1
1 0 v0 0 1
e0 v1
1 0 0 1
v1
v2 Vertex Mask
δ0 1 0 0 f0 1
α0
γ0
γn−2
αn−2
v1
γ1
β0
v2
v2
αn−1
βn−2
βn−1
γn−1 δn−2
Edge Mask
Fig. 4. Refinement stencils used at an extraordinary point .
Face Mask
δn−1
5
3 Spectral Analysis The 1-ring subdivision matrix S with the control vertices shown in Figure 4 is: v0 v1 v2 v1 v2 · · · v1 v2 e0 α0 β0 α1 β1 · · · αn−1 βn−1 f0 γ0 δ0 γ1 δ1 · · · γn−1 δn−1 e0 αn−1 βn−1 α0 β0 · · · αn−2 βn−2 S := f0 γn−1 δn−1 γ0 δ0 · · · γn−2 δn−2 ∈ R(2n+1)×(2n+1) . .. .. .. .. .. . . .. .. . . . . . . . . e0 α1 β1 α2 β2 · · · α0 β0 δ1 γ2 δ2 · · · γ0 δ0 f0 γ1 Now denote the discrete Fourier transform of a cyclic sequence {φi } by φˆi :=
n−1 X
φj ω i,j ,
wi,j := e
√
−1 2π n ij
.
j=0
Lemma 1 The spectrum of S are the diagonal entries of − + − + − Λ = diag(1, λ+ 0 , λ0 , λi , λi , · · · , λn−1 , λn−1 ), where λ± i :=
(ˆ αi + δˆi ) ±
q (ˆ αi − δˆi )2 + 4βˆi γˆi 2
,
i = 1 . . . n − 1.
(1)
Proof. Since the 2n × 2n lower block matrix of S is cyclic, we can apply the discrete Fourier transform 1 N0 F := 0 I2 (ω i,j ) to obtain
v0 e0 f 0 e0 Q := F −1 SF = f0 .. . e0 f0
nv1 α ˆ0 γˆ0 0 0 .. . 0 0
nv2 0 0 · · · βˆ0 0 0 · · · δˆ0 0 0 · · · 0 α ˆ 1 βˆ1 · · · 0 γˆ1 δˆ1 · · · .. .. .. . . . . . . 0 0
0 0 0 0 0 .. .
0 0 ··· α ˆ n−1 0 0 · · · γˆn−1
0 0 0 0 0 .. .
. βˆn−1 δˆn−1
− The eigenvalues 1, λ+ 0 , and λ0 come from the first 3 × 3 block. Each 2 × 2 diagonal block α ˆ i βˆi , i = 1 . . . n − 1, γˆi δˆi
contributes the pair of eigenvalues (1).
⊓ ⊔
6 + By (iv) and since λ+ i and λn−i are complex conjugate, we anticipate the following lemma. + Lemma 2 λ+ i = λn−i ,
− λ− i = λn−i ,
and
− λ+ i , λi ∈ R.
Proof. By symmetry (iv), ( P n−1 2 αj cos 2πij 2 j=1 if n is odd, n α ˆ i = α0 + P n2 −1 2πij 2 j=1 αj cos n + α n2 cos πi if n is even.
(2)
2π(n−i)j Since cos 2πij and cos(πi) = 1 when i is even, α ˆi = α ˆ n−i ∈ R. The n = cos n π π i − i ˆ ˆ same reasoning applies to δi , βi e n and γˆi e n since, e.g.
π βˆi e n i =
( P n−1 −1 2 2 j=0 βj cos (2πj+π)i + β n−1 cos πi n 2 n P 2 −1 (2πj+π)i 2 j=0 βj cos n
if n is odd,
(3)
if n is even.
± Therefore, λ± i = λn−i ∈ R as claimed.
⊓ ⊔
For n = 4 the subdivision stencils are shown in Figure 3. Based on this subdivision matrix, we can compute the matrix Q in Section 3 and therefore compute explicitly the expansions collected in Lemma 3. Lemma 3 In the regular setting (n = 4), 337 19 1 19 , , , ], 729 81 √ 9 81 √ π π π π 88 4 4 2 2 [βˆ04 e n i , βˆ14 e n i , βˆ24 e n i , βˆ34 e n i ] = [ , , 0, − ], 729 81√ 81 √ π π π π 352 16 2 16 2 [ˆ γ04 e− n i , γˆ14 e− n i , γˆ24 e− n i , γˆ34 e− n i ] = [ , , 0, − ], 729 81 81 121 11 1 11 [δˆ04 , δˆ14 , δˆ24 , δˆ34 ] = [ , , , ]. 729 81 9 81 ˆ 41 , α ˆ 42 , α ˆ 43 ] = [ [ˆ α40 , α
(4)
4 Deriving Weights In the Irregular Case To derive the extraordinary rules , we use a variation of the following heuristic generalizing the one in [7]: (a) Associate the values of the regular stencil (Figure 3) with the real components of n = 4 points regularly spaced on the unit circle, starting at (cos(0), sin(0)) for α and δ and starting at (cos( nπ ), sin( πn )) for β and (cos(− πn ), sin(− πn )) for γ. For example, for the edge rule in Figure 3, we associate β0n and β3n with cos πn , and β1n and β2n with cos 3π n .
7
(b) Find a polynomial interpolant of low degree to the data. Due to symmetry, we need only consider the first ⌈n/2⌉ data. For example, for the edge rule, the interpolant 4 40 and 729 for n = 4. since the β(t) needs only interpolate the first two weights 729 remaining two yield identical interpolation conditions. (c) The stencil weights for the irregular case with n 6= 4 is obtained by evaluating the polynomial at the real components of n points regularly spaced on the unit circle, starting at (cos(0), sin(0)) for α and δ and starting at (cos( πn ), sin( πn )) for β and (cos(− πn ), sin(− nπ )) for γ. While the heuristic evidently creates rules consistent with the regular case, it needs to be amended to create useful stencils. Firstly, the linear interpolant β satisfying 40 4 β n (cos nπ ) = 729 and β n (cos 3π n ) = 729 is negative on some subinterval of [−1, 1]. Therefore it cannot give suitable coefficients βin ≥ 0 for n > 4. Adding the constraint β n (−1) = 0 yields the quadratic interpolant √ √ √ 1 (44 − 18 2)(1 − t)2 + (44 − 9 2)2t(1 − t) + (36 2)t2 β 4 (t) := 729
which is nonnegative on [−1, 1]. Secondly, for α and δ, nonnegativity is best achieved by only retaining the structure of the polynomial interpolant and perturbing its coefficients. That is, we try the Ansatz an,2 4 4an,3 4 an,1 (t + bn,1 )2 + cn,1 , β n := β , γ n := β , αn (t) := n n n an,4 (5) δ n (t) := (t + bn,4 )2 + cn,4 . n Evaluating the polynomials at the projections of regularly spaced points on the unit circle yields two of the three types of stencils at the extraordinary point and leaves, as degrees of freedom for these two types, the eight coefficients depending on n in (5). The following lemma proves the Ansatz to be good. Lemma 4 The subdivision rules, for valence n and j = 0 . . . n − 1 αnj := αn (u), βjn := β n (cos
δjn := δ n (u), 2πj + π ), n
2πj , n 2πj − π γjn := γ n (cos ) n u := cos
(6) (7)
are consistent with the regular rules and with a truncated Fourier expansion for n > 4 if δˆn α ˆn (8) an,1 := 4α ˆ n2 , bn,1 := 1 , an,4 := 4δˆ2n , bn,4 := 1 . an,1 an,4 Proof. By symmetry, the weights are a discrete cosine transformation of their Fourier coefficients: αni =
n−1 2πij 1X n α ˆ cos , n j=0 j n
δin =
n−1 2πij 1 X ˆn δ cos . n j=0 j n
(9)
8
For n = 4, we verify the consistency of the choice (6) by substituting (4) into (9). α4i =
19 151 1 (u + )2 + , 18 18 324
δi4 =
1 11 41 (u + )2 − . 18 18 324
Next, we decompose each polynomial into a vector of basis functions and a vector of coefficients: an,1 2 an,1 1 an,1 an,1 an,1 1 an,1 T (u+bn,1 )2 +cn,1 = [1, u, 2u2−1][ bn,1 + + cn,1 , 2 bn,1 , ] . n n 2 n n n 2 n
We mimic this decomposition for n > 4, by choosing α ˆ ni = 0 for 2 < i < n − 2. Then αni =
n−1 1X n 2πij 2πi 4πi 1 n 2 n 2 n T α ˆ cos = [1, cos , cos ][ α ˆ , α ˆ , α ˆ ] . n j=0 j n n n n 0 n 1 n 2
4πi 2 Since cos 2πi n = u, we have cos n = 2u − 1 and comparison of terms yields
an,1 bn,1 = α ˆ n1 ,
and
1 an,1 = 2α ˆ n2 . 2
Solving for an,1 , bn,1 in terms of α ˆ n1 , α ˆ n2 , n, and similarly for an,4 , bn,4 , yields (8).
⊓ ⊔
The rules of Lemma 4 preserve rotational and p-mirror symmetry. For n > 4, they reduce the free parameters from 2n+4 to eight, namely an,i , bn,i , cn,i , i = 1, 4.andan,2, an,3 (For n = 3 we use the original 2n+4 = 10 free parameters to enforce λ+ λ+ 1 = λ, 2 = + − 2 2 2 λ , λ0 = λ , λ0 < λ .) It remains to determine the eight free parameters so that for n > 4 the weights αj , δj , e0 , f0 are nonnegative (βj and γj are nonnegative since β 4 ≥ 0 on [−1, 1]) and + + + + Λ = diag(1, λ+ 1 , λn−1 , λ2 , λn−2 , λ0 , ν1 , ν2 , . . . , ν2n−5 ) 2
2
2
= diag(1, λ, λ, λ , λ , λ , . . . , λi , . . .),
(10) 2
where 0 < |λi | < λ ,
and then check injectivity and regularity of the characteristic map. By Equations (1) and (8), the constraints on the subdominant and subsubdominant eigenvalues are q (an,1 bn,1 − an,4 bn,4 )2 + 4βˆ1 γˆ1 a b + a b + n,1 n,1 n,4 n,4 (11) λ = λ+ = 1 2 q (an,1 −an,4 )2 an,1 +an,4 + + 4βˆ2 γˆ2 4 16 + 2 . (12) λ = λ2 = 2 Since βˆ1 γˆ1 and βˆ2 γˆ2 are constant, each is an equation in the parameters an,1 , bn,1 , an,4 , − 2 2 bn,4 . The additional necessary constraint, λ+ 0 = λ , λ0 < λ can be enforced by choice of v1 , v2 in the first 3 by 3 block of Q. The result is an underconstrained system with linear inequality constraints. It is easy to pick closed-form solutions that satisfy (i)–(iv). One such, with the additional constraint λ = 1/3, is stated in the following Lemma 5.
9
Lemma 5 For n = 5..10, the following choice satisfies (i)-(iv) and additionally λ = 1/3: an,1 = .413,
bn,1 = .624,
cn,1 = .229,
an,2 = 4,
an,4 = .260,
bn,4 = .286,
cn,4 = .058,
an,3 = 4.
(13)
Proof. We need only verify (i), (ii) and (iii) listed in Section 2 since the scheme is symmetric by construction. Since αn , βn , γn and δn are just a multiple of α4 , β4 , γ4 and δ4 and the latter are non-negative on the interval of interest, all weights are nonnegative summing to 1. An explicit check shows all eigenvalues other than 1, 13 , 13 , 91 , 19 , 19 for n = 3 . . . 10 to be less than 19 . In particular, λ1 = λn−1 = 13 is the subdominant eigenvalue and 19 correspond to the appropriate Fourier indices. The u-differences, scaled to unit size, fall strictly into the lower right quadrant and, by symmetry, the v-differences fall strictly into the upper right quadrant (Figure 6). The partial derivatives are a convex combination of the differences and hence all pairwise crossproducts are strictly positive. By [12], the characteristic map is regular. In addition, as shown in Figure 5, the first half-segment of the control net does not intersect the negative x-axis. By the convex hull property and Theorem 21 of [10], injectivity of the characteristic map follows. ⊓ ⊔
Fig. 5. Control polyhedron of the characteristic map for n = 3..10. Injectivity test: each red fat segment does not intersect the negative x-axis.
While the above solution passes all constraints, the resulting surfaces exhibit, upon closer scrutiny, a clear defect. As illustrated in Figure 7, the surfaces often fall into the range of ‘hybrid’ curvature, i.e. each surface ring added by refinement has both region of positive and of negative Gauss curvature [5]. Closer study shows that the oscillations in the surface depend on the ratio of α over β in the edge stencil, and the ratio of γ over δ in the face stencil. Therefore, we make these ratios equal to the ratio of the known
10
Fig. 6. The normalized differences in the u direction (green) and v direction (blue).
Fig. 7. A surface generated by the rules of Lemma 5 exhibits hybrid curvature: (left) shaded image shows slight oscillations; (right) superimposed needles in the direction of the normal, with length of Gauss curvature (Gauss curvature needles) show both negative Gauss curvature ( blue downward pointing needles in the top center part of the ring) and positive curvature.
satisfactory case n = 4. (As one referee pointed out, fixing this ratio is a restriction that automatically occurs when mask values rather than stencil values are adjusted.) Setting α40 αn 190 0 β n = β 4 = 40 yields the constraint 0
0
and enforcing
γ0n δ0n
=
γ04 δ04
19 an,1 (1 + bn,1 )2 + cn,1 = βn, n 4 0
=
160 100
(14)
is equivalent to
5 an,4 (1 + bn,4 )2 + cn,4 = γ0n . (15) n 8 To apply a variant of the technique in [4] and move an elliptic surface out of the hybrid region, we define the variable s that scales the relative contribution of the central control point to its next iterate. n−1 n−1 X X 337 an,1 α4i = s αni = s + an,1 b2n,1 + an,1 cn,1 = , 2 729 j=0 i=0
(16)
11 n−1 n−1 X X 121 an,4 δi4 = s δin = s + an,4 b2n,4 + an,4 cn,4 = . 2 729 j=0 i=0
(17)
Pn−1 Pn−1 Pn−1 Pn−1 Similarly, i=0 βin = s j=0 βi4 and i=0 γin = s j=0 γi4 yield an,2 = an,3 = 4s. Then (14) and (16), respectively (15) and (17) imply bn,1 =
19 n 4 nβ0
337 − 729 s− 2an,1
an,1 2
and bn,4 =
5 n 8 nγ0
− 121 729 s − 2an,4
an,4 2
.
(18)
We can now derive weights satisfying (i)–(iv) by the following procedure: (1) choose s ∈ [0..1] and set an,2 = an,3 = 4s; (2) substitute (18) into (11), (12) to derive an,1 , an,4 and bn,1 , bn,4 ; (3) compute cn,1 , cn,4 from (16) and (17). Figure 8 illustrates the choice s = 0.9 for n = 8. The hybrid curvature behavior has been eliminated but an oscillation in the curvature (that curvature bounded schemes tend to be prone to) is still noticeable.
Fig. 8. Surfaces generated by steps (1)–(3). There are no visible oscillations in the surface compared to Figure 7. The Gauss curvature(top) is strictly positive and the Gauss curvature(bottom) is strictly negative, but oscillations in the Gauss curvature remain.
5 Discussion The main contribution of the paper is to provide a manageable set of parameters that make it easy to satisfy the formal constraints (i)–(iv). This set can be further pruned
12
by minimizing regions of meshes where the scheme exhibits hybrid behavior [12]. Schemes in this reduced space still do not guarantee surface fairness in the sense of preventing curvature oscillations. Acknowledgment: This work was supported by NSF Grant CCF-0430891 and URB grant DDF-111135-788129.
References 1. U. Augsdorfer, N. Dodgson and M. Sabin, Tuning Subdivision by Minimising Gaussian Curvature Variation Near Extraordinary Vertices. Proceedings of Eurographics (2006) 2. L. Barthe, and L. Kobbelt, Subdivision Scheme tuning around extraordinary vertices, Computer Aided Geometric Design, 21: 561-583. 3. E. Catmull, and J. Clark, Recursively Generated B-spline Surfaces on Arbitrary Topological Meshes. Computer Aided Design, 10(6): 350-355.(1978) 4. I. Ginkel and G. Umlauf, Loop Subdivision with Curvature Control. Proceedings of Eurographics (2006) 5. K. Karˇciauskas, J. Peters and U. Reif, Shape Characterization of Subdivision Surfaces – Case Studies, Comp. Aided Geom. Design, 21(6), 2004, 601–614. 6. C. Loop, Smooth Ternary Subdivision of Triangle Meshes Curve and Surface Fitting (SaintMalo, 2002). 7. C. Loop, Bounded curvature triangle Mesh subdivision with the convex hull property. The Visual Computer, 18 (5-6), 316-325. 8. A. Nasri, I. Hasbini, J. Zheng, T. Sederberg, Quad-based Ternary Subdivision, presentation Dagstuhl Seminar on Geometric Modeling, 29 May - 03 June, 2005 for scientists and engineers, 9. H. Prautzsch and G. Umlauf , A G2-Subdivision Algorithm Geometric Modelling, 217-224. (1996) 10. U. Reif and J. Peters. Structural Analysis of Subdivision Surfaces – A Summary Topics in Multivariate Approximation and Interpolation, Elsevier Science Ltd, K. Jetter et al., 2005, 149-190. 11. M. Sabin , Cubic recursive division with bounded curvature Curves and Surfaces, pages 411-414. (1991) 12. G. Umlauf, Analysis and Tuning of Subdivision Algorithms. Proceedings of the 21st spring conference on Computer Graphics, 33-40. (2005)
2
Dept. CISE, University of Florida Dept. of Computer Science American University of Beirut
Abstract. A well-documented problem of Catmull and Clark subdivision surfaces is that, in the neighborhood of extraordinary points, the curvature is unbounded and fluctuates. In fact, since one of the eigenvalues that determines elliptic shape is too small, the limit surface can have a saddle point when the designer’s input mesh suggests a convex shape. Here, we replace, near the extraordinary point, Catmull-Clark subdivision by another set of rules based on refining each bi-cubic B-spline into nine. This provides many localized degrees of freedom for special rules that need not reach out to neighbor vertices in order to tune the behavior. In this paper, we provide a strategy for setting such degrees of freedom and exhibit tuned ternary quad subdivision that yields surfaces with bounded curvature, nonnegative weights and full contribution of elliptic and hyperbolic shape components. Keywords: Subdivision, ternary, bounded curvature, convex hull
1 Introduction Tuning a subdivision scheme means adjusting the subdivision rules or stencils to obtain a refined mesh and surface with prescribed properties [1, 2, 4, 6, 7, 11, 12]. To date, no stationary subdivision algorithm on quadrilateral meshes is known to yield bounded curvature while guaranteeing the convex hull property at the extraordinary nodes. This paper explores ternary refinement to obtain such a scheme, because ternary subdivision offers more parameters in a close vicinity of each extraordinary node to tune the subdivision than binary subdivision does and thereby allows localizing the tuning [6]. Ternary subdivision inserts two knots in each knot interval, splitting each quad into nine. If all nodes of a quad are of valence 4, the rules for bi-cubic B-splines shown in Figure 3 are applied and we have C 2 continuity. The challenge is with nodes of valence other than 4, called extraordinary nodes. In particular, we can devise rules at the extraordinary nodes and their newly created direct and diagonal neighbors to achieve – eigenvalues in order of magnitude as 1, λ, λ, λ2 , λ2 , λ2 , λi , . . . and 0 ≤ |λi | < λ2 for i = 7, .., 2n + 1; and – nonnegative weights. It is possible to set, in addition, λ = 1/3. Yet, Figure 2 shows that the macroscopic shape of the new scheme and Catmull-Clark are very similar (Figure 2), but, of course, the mesh is denser. In the following, we first recall the rules for regular binary and ternary subdivisions (Section 1.1). In section 2, we state the properties that the ternary subdivision should match. Section 3 analyzes the spectrum in terms of free parameters and Section 4 derives weights by interpolating the regular case.
2
. Fig. 1. Catmull-Clark subdivision (top) and ternary subdivision (bottom). (left top) Control net with vertices (x, y, 50(1 − x2 − 5y 2 )), where (x, y) are the positions of the nodes of the characteristic map of Catmull-Clark subdivision. (middle three) control net at subdivision steps 3,4,5 and (right) shaded limit surface.
Fig. 2. Three steps of Catmull-Clark subdivision (left) and ternary subdivision (right). While the density differes, the macroscopic shape is identical.
3
1.1 Background and Regular Case Subdivision is a widely adopted tool in computer graphics and is also making inroads into geometric modeling, if only for conceptual modeling. However, [5] pointed out that Catmull-Clark subdivision in particular, is lacking the full set of elliptic and hyperbolic subsubdominant eigenfunctions and therefore typically generates saddle shapes in the limit at vertices with valence greater than four (see Figure 1 top). This implies that any high-quality (standard, symmetric) scheme needs to have a spectrum 1, λ, λ, λ2 , λ2 , λ2 followed by smaller terms. With the strategy explained in the following, it is possible to achieve such a spectrum by making subdivision stencils near the extraordinary point depend on the valence (Figure 1 bottom). A major challenge is technical: there are either too few or far too many free parameters occurring in nonlinear inequality and equality constraints. This makes any selection by general optimization [2] impractical. In [6,
1
6
1
1
1
16 1
6
36
6
6
1
6 Vertex
1
1
6 1 Edge
1
76
16
4
64
40
10
16 0
100
1 76
361
76
19
304
16
76 Vertex
16
4
64 Edge
190 16
256
160
16
10
1 Face
40
1
Face
Fig. 3. The regular bicubic B-spline refinement vertex, edge, and facet refinement stencils of Catmull-Clark (left) and ternary quad subdivision (right).
7], Loop proposed improvements to his well-known triangular scheme to achieve both bounded curvature for binary and ternary schemes. Based on fixed-degree interpolation of known weights, this approach sets many unknown weights and reduces the number of free parameters. It turns out to be far more difficult on Catmull-Clark case since the leading eigenvalues now come from a 2 by 2 diagonal block matrix rather than from a 1 by 1 block as the Loop subdivision. The eigenvalues are therefore the roots of quadratic polynomials while, in the case of triangular meshes, the eigenvalue is simple. The problem for n-sided facets is even more complex due to 4 × 4 blocks but can be avoided by applying one initial Catmull-Clark subdivision step. Most recently, two other approaches to tuning subdivision schemes at extraordinary points have been proposed. Augsdoerfer et al [1] choose an optimal subdominant eigenvalue of the subdivision scheme to minimize the curvature variation of the central surface [10] and Ginkel and Umlauf [4] propose new techniques to avoid hybrid configurations by tuning eigenvalues and eigencoefficients. These two and the present approach have in common that they focus on the curvature around the extraordinary point. In all cases, the localized improvement may, but is not guaranteed to, yield overall curvature distribution is of high-quality.
4
2 Problem Statement
We want to derive a ternary, quad subdivision scheme that is stationary, (rotational and p-mirror) symmetric and affine invariant. We denote the leading eigenvalues, in order of magnitude as, 1, λ, λ, µ1 , µ2 , µ3 . To achieve bounded curvature and the convex hull property, we enforce the following constraints on the eigenvalues and the weights αi , βi , γi , δi and vi of the ternary quad scheme (see Figure 4): We focus on the leading
(i) C 1 scheme 1 > λ > µi , and the characteristic map is regular and injective. (ii) Bounded Curvature |µ1 | = |µ2 | = |µ3 | = λ2 . (iii) Convex Hull αi ≥ 0, βi ≥ 0, γi ≥ 0, δi ≥ 0, i = 0..n − 1, v1 , v2 ≥ 0, P Pn−1 e0 := 1 − n−1 i=0 (αi + βi ) ≥ 0, f0 := 1 − i=0 (γi + δi ) ≥ 0, P v0 := 1 − n 2i=1 (vi ) ≥ 0. (iv) Symmetry αl = αn−l , βl = βn−1−l , γl = γ(n+1−l) mod n , and δl = δn−l Table 1. Constraints on the spectrum and the stencil weights. In addition, we will be able to enforce µ1 = µ2 = µ3 > λ7 > 0 and the subsubdominant eigenvectors corresponding to µi come from Fourier blocks with index 0, 2, n − 2.
eigenvalues and hence the 1-ring of control mesh around extraordinary point and its eigenstructure. As shown in Figure 4, each refinement step generates three types of points corresponding to vertices, edges and faces respectively. Let n be the valence of a generic extraordinary point. We have, offhand, 4n + 2 weights to determine. Symmetry reduces this to at most 2n + 4 free parameters.
α2
v1
v1
δ1
β1
v2
γ2
α1
1 0 v0 0 1
e0 v1
1 0 0 1
v1
v2 Vertex Mask
δ0 1 0 0 f0 1
α0
γ0
γn−2
αn−2
v1
γ1
β0
v2
v2
αn−1
βn−2
βn−1
γn−1 δn−2
Edge Mask
Fig. 4. Refinement stencils used at an extraordinary point .
Face Mask
δn−1
5
3 Spectral Analysis The 1-ring subdivision matrix S with the control vertices shown in Figure 4 is: v0 v1 v2 v1 v2 · · · v1 v2 e0 α0 β0 α1 β1 · · · αn−1 βn−1 f0 γ0 δ0 γ1 δ1 · · · γn−1 δn−1 e0 αn−1 βn−1 α0 β0 · · · αn−2 βn−2 S := f0 γn−1 δn−1 γ0 δ0 · · · γn−2 δn−2 ∈ R(2n+1)×(2n+1) . .. .. .. .. .. . . .. .. . . . . . . . . e0 α1 β1 α2 β2 · · · α0 β0 δ1 γ2 δ2 · · · γ0 δ0 f0 γ1 Now denote the discrete Fourier transform of a cyclic sequence {φi } by φˆi :=
n−1 X
φj ω i,j ,
wi,j := e
√
−1 2π n ij
.
j=0
Lemma 1 The spectrum of S are the diagonal entries of − + − + − Λ = diag(1, λ+ 0 , λ0 , λi , λi , · · · , λn−1 , λn−1 ), where λ± i :=
(ˆ αi + δˆi ) ±
q (ˆ αi − δˆi )2 + 4βˆi γˆi 2
,
i = 1 . . . n − 1.
(1)
Proof. Since the 2n × 2n lower block matrix of S is cyclic, we can apply the discrete Fourier transform 1 N0 F := 0 I2 (ω i,j ) to obtain
v0 e0 f 0 e0 Q := F −1 SF = f0 .. . e0 f0
nv1 α ˆ0 γˆ0 0 0 .. . 0 0
nv2 0 0 · · · βˆ0 0 0 · · · δˆ0 0 0 · · · 0 α ˆ 1 βˆ1 · · · 0 γˆ1 δˆ1 · · · .. .. .. . . . . . . 0 0
0 0 0 0 0 .. .
0 0 ··· α ˆ n−1 0 0 · · · γˆn−1
0 0 0 0 0 .. .
. βˆn−1 δˆn−1
− The eigenvalues 1, λ+ 0 , and λ0 come from the first 3 × 3 block. Each 2 × 2 diagonal block α ˆ i βˆi , i = 1 . . . n − 1, γˆi δˆi
contributes the pair of eigenvalues (1).
⊓ ⊔
6 + By (iv) and since λ+ i and λn−i are complex conjugate, we anticipate the following lemma. + Lemma 2 λ+ i = λn−i ,
− λ− i = λn−i ,
and
− λ+ i , λi ∈ R.
Proof. By symmetry (iv), ( P n−1 2 αj cos 2πij 2 j=1 if n is odd, n α ˆ i = α0 + P n2 −1 2πij 2 j=1 αj cos n + α n2 cos πi if n is even.
(2)
2π(n−i)j Since cos 2πij and cos(πi) = 1 when i is even, α ˆi = α ˆ n−i ∈ R. The n = cos n π π i − i ˆ ˆ same reasoning applies to δi , βi e n and γˆi e n since, e.g.
π βˆi e n i =
( P n−1 −1 2 2 j=0 βj cos (2πj+π)i + β n−1 cos πi n 2 n P 2 −1 (2πj+π)i 2 j=0 βj cos n
if n is odd,
(3)
if n is even.
± Therefore, λ± i = λn−i ∈ R as claimed.
⊓ ⊔
For n = 4 the subdivision stencils are shown in Figure 3. Based on this subdivision matrix, we can compute the matrix Q in Section 3 and therefore compute explicitly the expansions collected in Lemma 3. Lemma 3 In the regular setting (n = 4), 337 19 1 19 , , , ], 729 81 √ 9 81 √ π π π π 88 4 4 2 2 [βˆ04 e n i , βˆ14 e n i , βˆ24 e n i , βˆ34 e n i ] = [ , , 0, − ], 729 81√ 81 √ π π π π 352 16 2 16 2 [ˆ γ04 e− n i , γˆ14 e− n i , γˆ24 e− n i , γˆ34 e− n i ] = [ , , 0, − ], 729 81 81 121 11 1 11 [δˆ04 , δˆ14 , δˆ24 , δˆ34 ] = [ , , , ]. 729 81 9 81 ˆ 41 , α ˆ 42 , α ˆ 43 ] = [ [ˆ α40 , α
(4)
4 Deriving Weights In the Irregular Case To derive the extraordinary rules , we use a variation of the following heuristic generalizing the one in [7]: (a) Associate the values of the regular stencil (Figure 3) with the real components of n = 4 points regularly spaced on the unit circle, starting at (cos(0), sin(0)) for α and δ and starting at (cos( nπ ), sin( πn )) for β and (cos(− πn ), sin(− πn )) for γ. For example, for the edge rule in Figure 3, we associate β0n and β3n with cos πn , and β1n and β2n with cos 3π n .
7
(b) Find a polynomial interpolant of low degree to the data. Due to symmetry, we need only consider the first ⌈n/2⌉ data. For example, for the edge rule, the interpolant 4 40 and 729 for n = 4. since the β(t) needs only interpolate the first two weights 729 remaining two yield identical interpolation conditions. (c) The stencil weights for the irregular case with n 6= 4 is obtained by evaluating the polynomial at the real components of n points regularly spaced on the unit circle, starting at (cos(0), sin(0)) for α and δ and starting at (cos( πn ), sin( πn )) for β and (cos(− πn ), sin(− nπ )) for γ. While the heuristic evidently creates rules consistent with the regular case, it needs to be amended to create useful stencils. Firstly, the linear interpolant β satisfying 40 4 β n (cos nπ ) = 729 and β n (cos 3π n ) = 729 is negative on some subinterval of [−1, 1]. Therefore it cannot give suitable coefficients βin ≥ 0 for n > 4. Adding the constraint β n (−1) = 0 yields the quadratic interpolant √ √ √ 1 (44 − 18 2)(1 − t)2 + (44 − 9 2)2t(1 − t) + (36 2)t2 β 4 (t) := 729
which is nonnegative on [−1, 1]. Secondly, for α and δ, nonnegativity is best achieved by only retaining the structure of the polynomial interpolant and perturbing its coefficients. That is, we try the Ansatz an,2 4 4an,3 4 an,1 (t + bn,1 )2 + cn,1 , β n := β , γ n := β , αn (t) := n n n an,4 (5) δ n (t) := (t + bn,4 )2 + cn,4 . n Evaluating the polynomials at the projections of regularly spaced points on the unit circle yields two of the three types of stencils at the extraordinary point and leaves, as degrees of freedom for these two types, the eight coefficients depending on n in (5). The following lemma proves the Ansatz to be good. Lemma 4 The subdivision rules, for valence n and j = 0 . . . n − 1 αnj := αn (u), βjn := β n (cos
δjn := δ n (u), 2πj + π ), n
2πj , n 2πj − π γjn := γ n (cos ) n u := cos
(6) (7)
are consistent with the regular rules and with a truncated Fourier expansion for n > 4 if δˆn α ˆn (8) an,1 := 4α ˆ n2 , bn,1 := 1 , an,4 := 4δˆ2n , bn,4 := 1 . an,1 an,4 Proof. By symmetry, the weights are a discrete cosine transformation of their Fourier coefficients: αni =
n−1 2πij 1X n α ˆ cos , n j=0 j n
δin =
n−1 2πij 1 X ˆn δ cos . n j=0 j n
(9)
8
For n = 4, we verify the consistency of the choice (6) by substituting (4) into (9). α4i =
19 151 1 (u + )2 + , 18 18 324
δi4 =
1 11 41 (u + )2 − . 18 18 324
Next, we decompose each polynomial into a vector of basis functions and a vector of coefficients: an,1 2 an,1 1 an,1 an,1 an,1 1 an,1 T (u+bn,1 )2 +cn,1 = [1, u, 2u2−1][ bn,1 + + cn,1 , 2 bn,1 , ] . n n 2 n n n 2 n
We mimic this decomposition for n > 4, by choosing α ˆ ni = 0 for 2 < i < n − 2. Then αni =
n−1 1X n 2πij 2πi 4πi 1 n 2 n 2 n T α ˆ cos = [1, cos , cos ][ α ˆ , α ˆ , α ˆ ] . n j=0 j n n n n 0 n 1 n 2
4πi 2 Since cos 2πi n = u, we have cos n = 2u − 1 and comparison of terms yields
an,1 bn,1 = α ˆ n1 ,
and
1 an,1 = 2α ˆ n2 . 2
Solving for an,1 , bn,1 in terms of α ˆ n1 , α ˆ n2 , n, and similarly for an,4 , bn,4 , yields (8).
⊓ ⊔
The rules of Lemma 4 preserve rotational and p-mirror symmetry. For n > 4, they reduce the free parameters from 2n+4 to eight, namely an,i , bn,i , cn,i , i = 1, 4.andan,2, an,3 (For n = 3 we use the original 2n+4 = 10 free parameters to enforce λ+ λ+ 1 = λ, 2 = + − 2 2 2 λ , λ0 = λ , λ0 < λ .) It remains to determine the eight free parameters so that for n > 4 the weights αj , δj , e0 , f0 are nonnegative (βj and γj are nonnegative since β 4 ≥ 0 on [−1, 1]) and + + + + Λ = diag(1, λ+ 1 , λn−1 , λ2 , λn−2 , λ0 , ν1 , ν2 , . . . , ν2n−5 ) 2
2
2
= diag(1, λ, λ, λ , λ , λ , . . . , λi , . . .),
(10) 2
where 0 < |λi | < λ ,
and then check injectivity and regularity of the characteristic map. By Equations (1) and (8), the constraints on the subdominant and subsubdominant eigenvalues are q (an,1 bn,1 − an,4 bn,4 )2 + 4βˆ1 γˆ1 a b + a b + n,1 n,1 n,4 n,4 (11) λ = λ+ = 1 2 q (an,1 −an,4 )2 an,1 +an,4 + + 4βˆ2 γˆ2 4 16 + 2 . (12) λ = λ2 = 2 Since βˆ1 γˆ1 and βˆ2 γˆ2 are constant, each is an equation in the parameters an,1 , bn,1 , an,4 , − 2 2 bn,4 . The additional necessary constraint, λ+ 0 = λ , λ0 < λ can be enforced by choice of v1 , v2 in the first 3 by 3 block of Q. The result is an underconstrained system with linear inequality constraints. It is easy to pick closed-form solutions that satisfy (i)–(iv). One such, with the additional constraint λ = 1/3, is stated in the following Lemma 5.
9
Lemma 5 For n = 5..10, the following choice satisfies (i)-(iv) and additionally λ = 1/3: an,1 = .413,
bn,1 = .624,
cn,1 = .229,
an,2 = 4,
an,4 = .260,
bn,4 = .286,
cn,4 = .058,
an,3 = 4.
(13)
Proof. We need only verify (i), (ii) and (iii) listed in Section 2 since the scheme is symmetric by construction. Since αn , βn , γn and δn are just a multiple of α4 , β4 , γ4 and δ4 and the latter are non-negative on the interval of interest, all weights are nonnegative summing to 1. An explicit check shows all eigenvalues other than 1, 13 , 13 , 91 , 19 , 19 for n = 3 . . . 10 to be less than 19 . In particular, λ1 = λn−1 = 13 is the subdominant eigenvalue and 19 correspond to the appropriate Fourier indices. The u-differences, scaled to unit size, fall strictly into the lower right quadrant and, by symmetry, the v-differences fall strictly into the upper right quadrant (Figure 6). The partial derivatives are a convex combination of the differences and hence all pairwise crossproducts are strictly positive. By [12], the characteristic map is regular. In addition, as shown in Figure 5, the first half-segment of the control net does not intersect the negative x-axis. By the convex hull property and Theorem 21 of [10], injectivity of the characteristic map follows. ⊓ ⊔
Fig. 5. Control polyhedron of the characteristic map for n = 3..10. Injectivity test: each red fat segment does not intersect the negative x-axis.
While the above solution passes all constraints, the resulting surfaces exhibit, upon closer scrutiny, a clear defect. As illustrated in Figure 7, the surfaces often fall into the range of ‘hybrid’ curvature, i.e. each surface ring added by refinement has both region of positive and of negative Gauss curvature [5]. Closer study shows that the oscillations in the surface depend on the ratio of α over β in the edge stencil, and the ratio of γ over δ in the face stencil. Therefore, we make these ratios equal to the ratio of the known
10
Fig. 6. The normalized differences in the u direction (green) and v direction (blue).
Fig. 7. A surface generated by the rules of Lemma 5 exhibits hybrid curvature: (left) shaded image shows slight oscillations; (right) superimposed needles in the direction of the normal, with length of Gauss curvature (Gauss curvature needles) show both negative Gauss curvature ( blue downward pointing needles in the top center part of the ring) and positive curvature.
satisfactory case n = 4. (As one referee pointed out, fixing this ratio is a restriction that automatically occurs when mask values rather than stencil values are adjusted.) Setting α40 αn 190 0 β n = β 4 = 40 yields the constraint 0
0
and enforcing
γ0n δ0n
=
γ04 δ04
19 an,1 (1 + bn,1 )2 + cn,1 = βn, n 4 0
=
160 100
(14)
is equivalent to
5 an,4 (1 + bn,4 )2 + cn,4 = γ0n . (15) n 8 To apply a variant of the technique in [4] and move an elliptic surface out of the hybrid region, we define the variable s that scales the relative contribution of the central control point to its next iterate. n−1 n−1 X X 337 an,1 α4i = s αni = s + an,1 b2n,1 + an,1 cn,1 = , 2 729 j=0 i=0
(16)
11 n−1 n−1 X X 121 an,4 δi4 = s δin = s + an,4 b2n,4 + an,4 cn,4 = . 2 729 j=0 i=0
(17)
Pn−1 Pn−1 Pn−1 Pn−1 Similarly, i=0 βin = s j=0 βi4 and i=0 γin = s j=0 γi4 yield an,2 = an,3 = 4s. Then (14) and (16), respectively (15) and (17) imply bn,1 =
19 n 4 nβ0
337 − 729 s− 2an,1
an,1 2
and bn,4 =
5 n 8 nγ0
− 121 729 s − 2an,4
an,4 2
.
(18)
We can now derive weights satisfying (i)–(iv) by the following procedure: (1) choose s ∈ [0..1] and set an,2 = an,3 = 4s; (2) substitute (18) into (11), (12) to derive an,1 , an,4 and bn,1 , bn,4 ; (3) compute cn,1 , cn,4 from (16) and (17). Figure 8 illustrates the choice s = 0.9 for n = 8. The hybrid curvature behavior has been eliminated but an oscillation in the curvature (that curvature bounded schemes tend to be prone to) is still noticeable.
Fig. 8. Surfaces generated by steps (1)–(3). There are no visible oscillations in the surface compared to Figure 7. The Gauss curvature(top) is strictly positive and the Gauss curvature(bottom) is strictly negative, but oscillations in the Gauss curvature remain.
5 Discussion The main contribution of the paper is to provide a manageable set of parameters that make it easy to satisfy the formal constraints (i)–(iv). This set can be further pruned
12
by minimizing regions of meshes where the scheme exhibits hybrid behavior [12]. Schemes in this reduced space still do not guarantee surface fairness in the sense of preventing curvature oscillations. Acknowledgment: This work was supported by NSF Grant CCF-0430891 and URB grant DDF-111135-788129.
References 1. U. Augsdorfer, N. Dodgson and M. Sabin, Tuning Subdivision by Minimising Gaussian Curvature Variation Near Extraordinary Vertices. Proceedings of Eurographics (2006) 2. L. Barthe, and L. Kobbelt, Subdivision Scheme tuning around extraordinary vertices, Computer Aided Geometric Design, 21: 561-583. 3. E. Catmull, and J. Clark, Recursively Generated B-spline Surfaces on Arbitrary Topological Meshes. Computer Aided Design, 10(6): 350-355.(1978) 4. I. Ginkel and G. Umlauf, Loop Subdivision with Curvature Control. Proceedings of Eurographics (2006) 5. K. Karˇciauskas, J. Peters and U. Reif, Shape Characterization of Subdivision Surfaces – Case Studies, Comp. Aided Geom. Design, 21(6), 2004, 601–614. 6. C. Loop, Smooth Ternary Subdivision of Triangle Meshes Curve and Surface Fitting (SaintMalo, 2002). 7. C. Loop, Bounded curvature triangle Mesh subdivision with the convex hull property. The Visual Computer, 18 (5-6), 316-325. 8. A. Nasri, I. Hasbini, J. Zheng, T. Sederberg, Quad-based Ternary Subdivision, presentation Dagstuhl Seminar on Geometric Modeling, 29 May - 03 June, 2005 for scientists and engineers, 9. H. Prautzsch and G. Umlauf , A G2-Subdivision Algorithm Geometric Modelling, 217-224. (1996) 10. U. Reif and J. Peters. Structural Analysis of Subdivision Surfaces – A Summary Topics in Multivariate Approximation and Interpolation, Elsevier Science Ltd, K. Jetter et al., 2005, 149-190. 11. M. Sabin , Cubic recursive division with bounded curvature Curves and Surfaces, pages 411-414. (1991) 12. G. Umlauf, Analysis and Tuning of Subdivision Algorithms. Proceedings of the 21st spring conference on Computer Graphics, 33-40. (2005)
