in KeealoopatdthaauaesLntncod.Erebod.Sab个 1 Introduction
Bounding the Distance between a Loop Subdivision Surface and Its Limit Mesh Zhangjin Huang and Guoping Wang School of Electronic Engineering and Computer Science, Peking University, Beijing 100871, China zhangjin.huang@gmail.com http://www.graphics.pku.edu.cn Abstract. Given a control mesh of a Loop subdivision surface, by pushing the control vertices to their limit positions, a limit mesh of the Loop surface is obtained. Compared with the control mesh, the limit mesh is a tighter linear approximation in general, which inscribes the limit surface. We derive an upper bound on the distance between a Loop subdivision surface patch and its limit triangle in terms of the maximum norm of the mixed second differences of the initial control vertices and a constant that depends only on the valence of the patch’s extraordinary vertex. A subdivision depth estimation formula for the limit mesh approximation is also proposed. Key words: Loop subdivision surfaces, Limit mesh, Error bound, Subdivision depth 1 Introduction A subdivision surface is defined as the limit of a sequence of successively refined control meshes. In practice such as surface rendering, surface trimming and surface intersection, a linear approximation (for example, a refined control mesh) is used to substitute the smooth surface. It is natural to ask the following questions: How does one estimate the error (distance) between a limit surface and its linear approximation (for instance, the control mesh)? How many steps of subdivision would be necessary to satify a user-specified error tolerance? Loop subdivision surfaces generalize the quartic three-directional box spline surfaces to triangular meshes of arbitrary topology [1]. Several works have been devoted to studying the distance between a Loop surface and its control mesh approximation [2–6]. These error estimation techniques can be classified into two classes. One is the vertex based method [2, 3], which measures the distance between the control vertices and their limit positions. In [2], Lanquetin et al. derived a wrong exponential bound which resulted in a improper subdivision depth formula. In [3], Wang et al. proposed a proper exponential bound with an awkward subdivision depth estimation technique. The vertex based bounds all suffer from
2 r be of the patch's d in 5e 2 Definition and Notation the initis 2.1 Distane (.)-F(e.l
2 Z. Huang and G. Wang one problem as pointed out in [3]: they may be smaller than the actual maximal distance between the limit surface and its control mesh in some cases. The other is the patch based method [4–6], which estimates the parametric distance between a Loop surface patch and its control mesh. Wu et al. presented an accurate error measure for adaptive subdivision and interference detection [4, 5]. But their estimate is dependent on recursive subdivision, thus can not be used to pre-compute the error bound after n steps of subdivision or predict the recursion depth within a user-specified tolerance. In [6], Huang derived a bound in terms of the maximum norm of the mixed second differences of the initial control vertices, and a subdivision depth estimation approach. The limit mesh of a Loop surface, obtained by pushing the control vertices of its control mesh to their limit positions, has been applied in surface reconstruction [8]. The limit mesh is usually considered a better linear approximation than the corresponding control mesh. But the approximation error for Loop surfaces has not been investigated yet. Recently, an effective approach has been proposed to bound the error of the limit mesh approximation to a Catmull-Clark subdivision surface [7]. With an analogous technique, in this paper we derive a bound on the maximum distance between a Loop surface patch and its limit triangle. The bound is expressed in terms of the maximum norm of the mixed second differences of the Loop patch’s initial control mesh and a constant that depends only on the valence of the patch’s extraordinary vertex. The paper is organized as follows. Section 2 introduces some definitions and notations. In Sections 3 and 4, we derive distance bounds for regular Loop patches and extraordinary Loop patches, respectively. And a subdivision depth estimation method for the limit mesh approximation is presented in Section 5. In Section 6, we compare the limit mesh approximation with the control mesh approximation. Finally we conclude the paper with some suggestions for future work. 2 Definition and Notation Without loss of generality, we assume the initial control mesh has been subdivided at least once, isolating the extraordinary vertices so that each face is a triangle and contains at most one extraordinary vertex. 2.1 Distance Given a control mesh of a Loop subdivision surface S˜, for each interior triangle F in the control mesh, there is a corresponding limit triangle F in the limit mesh, and a corresponding surface patch S in the limit surface S˜. Obviously, the limit triangle F is the triangle formed by connecting the three corner points of S. The distance between a Loop patch S and the corresponding limit triangle F is defined as the maximum distance between S(v, w) and F(v, w), that is, max (v,w)∈Ω kS(v, w) − F(v, w)k
2.2 Sccond Order Norm D +Pa-P P.IP.+P2-P-P✉ e+12nmwmkaP
Bounding the Distance between a Loop Surface and Its Limit Mesh 3 bc bc bc bc bc bc bc bc bc bc bc bc bc 5 6 n n+6 4 1 n+1 n+5 3 2 n+2 n+4 n+3 S = S 0 0 (a) bc bc bc bc bc bc bc bc bc bc bc bc bc b b b b b b b b b b b b b b 5 6 n b b b b b n+6 n+12 4 1 n+1 n+5 n+11 3 2 n+2 n+10 n+4 n+3 n+7 n+9 n+8 S1 0 S1 1 S1 2 S1 3 (b) Fig. 1. (a) Ordering of the control vertices of an extraordinary patch S of valence n. (b) Ordering of the new control vertices (solid dots) after one step of Loop subdivision. where Ω is the unit triangle, S(v, w) is the Stam’s parametrization [9] of S over Ω, and F(v, w) is the linear parametrization of F over Ω. 2.2 Second Order Norm The control mesh Π of a Loop patch S consists of n + 6 control vertices Π = {Pi : i = 1, 2, . . . , n+ 6}, where n is the valence of F’s only extraordinary vertex (if has, otherwise n = 6) and called the valence of the patch S (see Figure 1(a)). The second order norm of Π, denoted M = M0 = M0 0 , is defined as the maximum norm of the following n + 9 mixed second differences (MSDs) {αi : 1 ≤ i ≤ n + 9} of the n + 6 vertices of Π: M = max{kP1 + P2 − Pn+1 − P3k, {kP1 + Pi − Pi−1 − Pi+1k : 3 ≤ i ≤ n}, kP1 + Pn+1 − Pn − P2k, kP2 + Pn+1 − P1 − Pn+2k, kP2 + Pn+2 − Pn+1 − Pn+3k, kP2 + Pn+3 − Pn+2 − Pn+4k, kP2 + Pn+4 − Pn+3 − P3k, kP2 + P3 − Pn+4 − P1k, kPn+1 + Pn − P1 − Pn+6k, kPn+1 + Pn+6 − Pn − Pn+5k, kPn+1 + Pn+5 − Pn+6 − Pn+2k, kPn+1 + Pn+2 − Pn+5 − P2k} = max{kαik : i = 1, . . . , n + 9} . (1) M is also called the (level-0) second order norm of the patch S. For a regular patch (n = 6), there are 15 mixed second differences. Through subdivision we can generate n+12 new vertices P1 i , i = 1, . . . , n+12 (see Figure 1(b)), which are called the level-1 control vertices of S. All these level- 1 control vertices compose S’s level-1 control mesh Π1 = {P1 i : i = 1, 2, . . . , n + 12}. We use Pk i and Πk to represent the level-k control vertices and level-k control mesh of S, respectively, after k steps of subdivision on Π. The level-1 control vertices form four control vertex sets Π1 0 , Π1 1 , Π1 2 and Π1 3 , corresponding to the control meshes of the subpatches S 1 0 , S 1 1 , S 1 2 and S 1 3
4 3 Regular Patch Quartie Triangular Form pSa物sR罗ata 2 我k ((.)v0]and
4 Z. Huang and G. Wang respectively (see Figure 1b), where Π1 0 = {P 1 i : 1 ≤ i ≤ n + 6}. The subpatch S 1 0 is an extraordinary patch, but S 1 1 , S 1 2 and S 1 3 are regular triangular patches [9]. Following the definition in Equation (1), one can define the second order norms M1 i for S 1 i , i = 0, 1, 2, 3, respectively. M1 = max{M1 i : i = 0, 1, 2, 3} is defined as the second order norm of the level-1 control mesh Π1 . After k steps of subdivision on Π, one gets 4k control point sets Πk i : i = 0, 1, . . . , 4 k − 1 corresponding to the 4k subpatches S k i : i = 0, 1, . . . , 4 k − 1 of S, with S k 0 being the only level-k extraordinary patch (if n 6= 6). We denote Mk i and Mk as the second order norms of Πk i and Πk , respectively. 3 Regular Patch In this section, both the regular Loop patch S and its corresponding limit triangle F are first expressed in quartic triangular B´ezier form. Then we bound the distance between S and F by measuring the deviations between their corresponding B´ezier points. 3.1 Quartic Triangular B´ezier Form b b b b b b b b b 2 5 9 b b b 1 4 8 12 3 7 11 6 10 S b b b Ω (0,0) u=1 (0,1) (1,0) w=1 v=1 Fig. 2. Control vertices of a quartic box spline patch and their ordering. If S is a regular Loop patch, then S(u, v) can be expressed as a quartic box spline patch defined by 12 control vertices pi , 1 ≤ i ≤ 12 (see Figure 2): S(v, w) = X 12 i=1 piNi(v, w), (v, w) ∈ Ω , (2) where Ni(v, w), 1 ≤ i ≤ 12 are the quartic box spline basis functions. The expressions of Ni(v, w) refer to [9]. The surface is defined over the unit triangle: Ω = {(v, w)| v ∈ [0, 1] and w ∈ [0, 1 − v]}
in 12uw Caaa创 su=∑bae.een 3) 是 o程 F(c.w)=uba hearpa F,=∑5成c回 0 ,tL-F子，》，，k之0,0≤j+≤4 are the Rzir points. sic.)-Fc.)1-)s-1(c.).8
Bounding the Distance between a Loop Surface and Its Limit Mesh 5 b b b b b b b b b b b b b b b 1 3 6 10 15 5 9 14 8 13 12 11 7 4 2 (a) u 4 4u 3 v 4u 3w 6u 2 v 2 12u 2 vw 6u 2w 2 4uv3 12uv2w 12uvw2 4uw3 v 4 4v 3w 6v 2w 2 4vw3 w 4 (b) Fig. 3. (a) Ordering of the B´ezier points of a quartic triangular B´ezier patch. (b) Quartic Bernstein polynomials corresponding to the B´ezier points. We introduce the third parameter u = 1 − v − w such that (u, v, w) forms a barycentric system of coordinates for the unit triangle. S is a quartic triangular B´ezier patch, thus S(u, v) can be written in terms of Bernstein polynomials [11]: S(u, v) = X 15 i=1 biBi(v, w), (v, w) ∈ Ω , (3) where bi , 1 ≤ i ≤ 15 are the B´ezier points of S, and Bi(v, w), 1 ≤ i ≤ 15 are the quartic Bernstein polynomials (see Figure 3). The correspondence between the standard representation B4 ijk(u, v, w) = 4! i!j!k! u iv jw k , i, j, k ≥ 0, i+j +k = 4 and Bi(v, w), 1 ≤ i ≤ 15 is (i, j, k) 7→ (j+k)(j+k+1) 2 + k + 1. With the algorithm developed in [10], we can get the 15 × 12 matrix T which converts from the 12 control vertices (pi) to the 15 B´ezier points (bi) (see Appendix 1). Note that F = {b1, b11, b15} is the limit triangle corresponding to the center triangle F = {p4, p7, p8} (see Figures 2 and 3). The linear parametrization of F is F(v, w) = ub1 + vb11 + wb15 , where u = 1−v−w. By the linear precision property of the Bernstein polynomials [11], we can express F(v, w) as the following quartic B´ezier form: F(v, w) = X 15 i=1 biBi(v, w) , (4) where b(j+k)(j+k+1) 2 +k+1 = F( j 4 , k 4 ), j, k ≥ 0, 0 ≤ j + k ≤ 4 are the B´ezier points. 3.2 Bounding the Distance According to the previous analysis, for (v, w) ∈ Ω, it follows that kS(v, w) − F(v, w)k = X 15 i=1 (bi − bi)Bi(v, w) ≤ X 15 i=1 kbi − bikBi(v, w) . (5)
6 h:他置o llet poiblec1..5 Ib,-5≤6M hdon where-1,2.....15 are undetermined real coefficients.It follows that Ib-e2Ials24anls上hiw anonniwrra 香-mm工时 6) 4∑-h- -6--3 6-6--是 It follows that 5四-4De,四-+w-2-m- Thcorem e ance bet)and w 1S,四y-r,sB(M
6 Z. Huang and G. Wang In the following, we compute the smallest possible constants δi , i = 1, 2, . . . , 15 such that kbi − bik is bounded by kbi − bik ≤ δiM , where M is the second order norm of S. It is obvious that δ1 = δ11 = δ15 = 0. bi and bi , 1 ≤ i ≤ 15 are the convex combinations of the control vertices pi , i = 1, 2, . . . , 12, and bi−bi can be expressed as a linear combination of the 15 MSDs αl , l = 1, 2, . . . 15 defined in Equation (1) as follows: bi − bi = X 15 l=1 x i lαl , where x i l , l = 1, 2, . . . , 15 are undetermined real coefficients. It follows that kbi − bik ≤ X 15 l=1 kx i lαlk ≤ X 15 l=1 |x i l |kαlk ≤ X 15 l=1 |x i l |M . Therefore, to get a tight upper bound for kbi − bik, we solve the following constrained minimization problem: δi = minX 15 l=1 |x i l | s.t. X 15 l=1 x i lαl = bi − bi . (6) By symmetry, we only need to solve three constrained minimization problems. With the help of the symbolic computation of Mathematica, we have δ1 = δ11 = δ15 = 0 , δ2 = δ3 = δ7 = δ10 = δ12 = δ14 = 1 4 , δ4 = δ6 = δ13 = 1 3 , δ5 = δ8 = δ9 = 5 12 . It follows that B(v, w) = X 15 i=1 δiBi(v, w) = v + w − v 2 − vw − w 2 . We obtain a bound on the pointwise distance between S(v, w) and F(v, w): Theorem 1. For (v, w) ∈ Ω, we have kS(v, w) − F(v, w)k ≤ B(v, w)M , where B(v, w) = v + w − v 2 − vw − w 2 is called the distance bound function of S(v, w) with respect to F(v, w)
路的257E子u以ed 。-w2-3-63号 是pd (e.)-F(e. 4 Extraordinary Patch 4.1 Stam'Paramctrization Vig.Pteparameder domin }
Bounding the Distance between a Loop Surface and Its Limit Mesh 7 By symmetry, the maximum of B(u, v),(u, v) ∈ Ω must occur on the diagonal B(t, t) = 2t − 3t 2 , 0 ≤ t ≤ 1 2 . Since max 0≤t≤1/2 2t − 3t 2 = 1 3 = B( 1 3 , 1 3 ) , we have a bound on the maximal distance between S(u, v) and F(u, v) as stated in the following theorem: Theorem 2. The distance between a regular Loop patch S and the corresponding limit triangle F is bounded by max (v,w)∈Ω kS(v, w) − F(v, w)k ≤ 1 3 M . 4 Extraordinary Patch For an extraordinary patch, we first partition it into regular triangular subpatches with Stam’s parametrization [9], then derive a distance bound function for each regular subpatch with the technique developed in Section 3. 4.1 Stam’s Parametrization Ω 1 1 Ω 1 2 Ω 1 3 Ω 2 1 Ω 2 2 Ω 2 3 Ω3 1 Ω3 2 Ω3 3 (0,0) (1,0) (0,1) Fig. 4. Partition of the parameter domain Ω. Through subdivision an extraordinary Loop patch S of valence n can be partitioned into an infinite sequence of regular triangular patches {S k m}, k ≥ 1, m = 1, 2, 3. If we partition the unit triangle Ω into an infinite set of tiles {Ωk m}, k ≥ 1, m = 1, 2, 3 (see Figure 4), accordingly, with Ω k 1 = {(v, w)| v ∈ [2−k , 2 −k+1] and w ∈ [0, 2 −k+1 − v]} , Ω k 2 = {(v, w)| v ∈ [0, 2 −k ] and w ∈ [2−k − v, 2 −k ]} , Ω k 3 = {(v, w)| v ∈ [0, 2 −k ] and w ∈ [2−k , 2 −k+1 − v]}
8 .Hunng and G.Wang s,=s5=， P-r+s+mr,,ep。 rleFa be ptidtrdefnd F,时l选-t点c,》 oc晋他2T Bi-F(e.)-F(ea,wa).=F(rs.) 4.Douding the Distanc s6国)-F6,=s低，-偏训≤2-1以6，.9 Natir that is not the limit E( a-五引≤4M,1≤f≤15 m Equa (9)tha s,o-t训s5酷M,化，叫en
8 Z. Huang and G. Wang Each tile Ωk m corresponds to a box spline patch S k m. And S k m(v, w) is defined over the unit triangle with the form as Equation (2). Therefore, the parametrization for S(v, w) is constructed as follows [9]: S(v, w) |Ωkm = S k m(˜v, w˜) = S k m(t k m(v, w)) , where the transformation t k m maps the tile Ωk m onto the unit triangle Ω. The center triangle of S’s control mesh is F = {P1, P2, Pn+1} (see Figure 1), and the corresponding limit triangle is F = {P1, P2, Pn+1}, with Pi being the limit point of Pi , i = 1, 2, n + 1. Let F(u, v) be the linear parametrization of F: F(v, w) = uP1 + vP2 + wPn+1, (v, w) ∈ Ω . (7) The limit triangle F can be partitioned into sub-triangles defined over Ωk m as follows: F(v, w) |Ωkm = Fbk m(t k m(v, w))) . Here Fbk m is the linear patch defined as Fbk m(v, w) = ub1 + vb11 + wb15, (v, w) ∈ Ω . (8) If the three corners of Ωk m are (v1, w1),(v2, w2) and (v3, w3), which correspond to (0, 0),(1, 0) and (0, 1) in Ω via the transformation t k m, respectively. Then b1 = F(v1, w1) , b11 = F(v2, w2) , b15 = F(v3, w3) . 4.2 Bounding the Distance Similar to the analysis in Section 3, we can rewrite S k m(v, w) and Fbk m(v, w) into the quartic B´ezier forms as Equations (3) and (4), respectively. Thus for (v, w) ∈ Ωk m, we have kS(v, w) − F(v, w)k = kS k m(˜v, w˜) − Fbk m(˜v, w˜)k ≤ X 15 i=1 kbi − bikBi(˜v, w˜). (9) Notice that Fbk m is not the limit triangle of the triangular patch S k m but one portion of the extraordinary patch S’s limit triangle F. So we can not use the results for kbi − bik derived in Section 3 directly. By solving 15 constrained minimization problems with the form similar to Equation (6), we have kbi − bik ≤ δiM, 1 ≤ i ≤ 15 . Consequently, it follows from Equation (9) that kS k m(v, w) − Fbk m(v, w)k ≤ Bk m(v, w)M, (v, w) ∈ Ω
where the quartic Bezic function 院em-∑iae,e时en (an) 器款告 8el=.≥1.m=123 sEC m点置a而er5对aen Ise,o)-F化，训sM， in the domsin元-Pn,)l+w≥】 m照o4ea制 1 ble 2.F 5 Subdivision Depth Estimation 1+与≤，)A货，≥0 (12
Bounding the Distance between a Loop Surface and Its Limit Mesh 9 where the quartic B´ezier function B k m(v, w) = X 15 i=1 δiBi(v, w), (v, w) ∈ Ω (10) is the distance bound function of S k m(v, w) with respect to Fbk m(v, w). Then the distance bound function of S(v, w) with respect to F(v, w), B(v, w),(v, w) ∈ Ω, can be defined as follows: B(v, w) |Ωkm = B k m(t k m(v, w)), k ≥ 1, m = 1, 2, 3 . It is obvious that B(0, 0) = B(1, 0) = B(0, 1) = 0, and B(v, w) is a piecewise quartic triangular B´ezier function over Ω away from (0, 0). Let β(n) = max(v,w)∈Ω B(v, w), we have the following theorem on the maximal distance between S(v, w) and F(v, w): Theorem 3. The distance between an extraordinary Loop patch S of valence n and the corresponding limit face F is bounded by max (v,w)∈Ω kS(v, w) − F(v, w)k ≤ β(n)M , (11) where β(n) is a constant that depends only on n, the valence of S. For 3 ≤ n ≤ 50, we investigate the maximums of the quartic B´ezier functions B k m(v, w), k ≥ 1, m = 1, 2, 3 over Ω assisted by plotting the graph of B(v, w). The following facts are found: 1. B(v, w) attains its maximum in the domain Ω = Ω ∩ {(v, w)| v + w ≥ 1 4 } (shaded region in Figure 4). 2. B(v, w) attains its maximum either on the diagonal B(t, t), t ∈ [0, 1 2 ] or on the borders B(t, 0) and B(0, t), t ∈ [0, 1]. By symmetry, to compute the value of β(n), at most four distance bound functions corresponding to S 1 1 , S 1 2 , S 2 1 , S 2 2 are needed to be analyzed. Numerical results for β(n), 3 ≤ n ≤ 10 are given in Table 2. For regular Loop patches, the constant β(n) = 1 3 is optimum. But for extraordinary Loop patches, the constants can be improved through further subdividing the subpatches S k m. 5 Subdivision Depth Estimation Before estimating subdivision depth, we investigate the recurrence relation between the second order norms of the control meshes of S at different levels. 5.1 Convergence Rate of Second Order Norm If the second order norms Mk+j 0 and Mj 0 satisfy the following recurrence inequality Mk+j 0 ≤ rj (n)Mk 0 , j ≥ 0 , (12)
10 Wang n,the v "-宫城 aa s=血 (13) g=:1st≤+明atiamle be n n创- -P红 rn1(n=c1(n以.n>3 +≤rn.j1 thn th ve nhem
10 Z. Huang and G. Wang where rj (n) is a constant which depends on n, the valence of the extraordinary vertex, and r0(n) ≡ 1. We call rj (n) the j-step convergence rate of second order norm. The convergence rate reflects how fast the control mesh converges to the limit surface. The smaller the convergence rate is, the faster the control mesh converges. It is obvious that rj+k ≤ rj rk. Let α k i , i = 1, 2, . . . , n + 9 be the MSDs of Πk 0 , k ≥ 0 defined as in Equation (1). For each l = 1, 2, . . . , n + 9, we express α k+1 l as a linear combination of α k i : α k+1 l = nX +9 i=1 x l iα k i , where x l i , i = 1, 2, . . . , n + 9 are undetermined real coefficients. Then we can bound kα k+1 l k by cl(n)Mk 0 , where cl(n) is the solution of the following constrained minimization problem cl(n) = min nX +9 i=1 |x l i | , s.t. nX +9 i=1 x l iα k i = α k+1 l . (13) Since Mk+1 0 = max{kα k+1 l k : 1 ≤ l ≤ n + 9}, we get an estimate for r1(n) as follows r1(n) = max 1≤l≤n+9 cl(n) . By symmetry, we only need to solve at most four constrained minimization problems corresponding to α k+1 1 = P k+1 1 + P k+1 2 − P k+1 n+1 − P k+1 3 , α k+1 n+1 = P2 + Pn+1 − P1 − Pn+2, α k+1 n+2 = P2 + Pn+2 − Pn+1 − Pn+3, and α k+1 n+3 = P2 + Pn+3 − Pn+2 − Pn+4, respectively. Since P k+1 1 is the extraordinary vertex, it is not surprising to find out that c1(n) is the maximum for n > 3. The special case is r1(3) = c4(3) = 0.4375 > c1(3) = 0.3125. Then it follows that r1(n) = c1(n), n > 3 . Similarly, we can estimate rj (n), n ≥ 3, j > 1 by solving only one constrained minimization problem (13) with α k+1 1 replaced by α k+j 1 . Then we have Lemma 1. If Mk 0 represents the second order norm of the level-k extraordinary subpatch S k 0 , k ≥ 0 of valence n, then it follows that Mk+j 0 ≤ rj (n)Mk 0 , j ≥ 1 . The above lemma works in a more general sense, that is, if Mk 0 is replaced with Mk , the second order norm of the level-k control mesh Πk , the estimates for rj (n) still work. Though r2(n) can be roughly estimated as r1(n) 2 , in practice r2(n) may derive better results than r1(n) 2 as shown in the next subsection. Table 1 shows
