中国科学技术大学：Estimating Error Bounds and Subdivision Depths for Loop Subdivision Surfaces

Estimating Error Bounds and Subdivision Depths for Loop Subdivision Surfaces Zhangjin Huang Scbool of Electronic Engineering and Computer Science, Peking University,Beijing 100871,China 0ct,2007 Abstract We derive a bound on the maximal distanoe betwreen a Loop sub- division surfnce patch and its control mesh in terms of the maximum norm of the tnixed second dilferenoes of the control points and a con- stant that pemb only on the lenx of the pateh.A sulxlivisim dopth formmla is also proposed Keywords:Loop subdivision surface,control mesh.error bound,sub- division depth 1 Introduction A subdivision surfaoc is defined as the limit of a finer and finer control mesh by subdividing the mesh recursively.The Loop subdivision surface geteral izes the quartic three-directional box spline surface to triangular meshes of arbitrary topology [1]. Previons erroe estimation techniques for Loop suhdivision surfnces can be clasified into two classes. One is the vertex hnsed method 2,3,which measures the distance between the vertloes and their limit positions.Lanquetin et al.derived a wrong exponential bound and consequently a wrong subdivision depth formula [2].Wang et al.propoeed a proper exponential bound with an inefficient subdivision depth estimation technique 3 As pointed out in 3, the vertex basei bounds all sulfer Irom one proble:they iay be sanaller than the nctnnl distance in some enses. The other is the patch based method [4.5].which estimates the para- metric distanoe between a Loop surfnce patch and its control mesh.Wu ct aL.presented an accurate error measure [4,5].But their estimate is de- pendent on recursive subdivision,thus can not be used to pre-compute the

Estimating Error Bounds and Subdivision Depths for Loop Subdivision Surfaces Zhangjin Huang School of Electronic Engineering and Computer Science, Peking University, Beijing 100871, China Oct., 2007 Abstract We derive a bound on the maximal distance between a Loop sub￾division surface patch and its control mesh in terms of the maximum norm of the mixed second differences of the control points and a con￾stant that depends only on the valence of the patch. A subdivision depth formula is also proposed. Keywords: Loop subdivision surface, control mesh, error bound, sub￾division depth 1 Introduction A subdivision surface is defined as the limit of a finer and finer control mesh by subdividing the mesh recursively. The Loop subdivision surface general￾izes the quartic three-directional box spline surface to triangular meshes of arbitrary topology [1]. Previous error estimation techniques for Loop subdivision surfaces can be classified into two classes. One is the vertex based method [2, 3], which measures the distance between the vertices and their limit positions. Lanquetin et al. derived a wrong exponential bound and consequently a wrong subdivision depth formula [2]. Wang et al. proposed a proper exponential bound with an inefficient subdivision depth estimation technique [3]. As pointed out in [3], the vertex based bounds all suffer from one problem: they may be smaller than the actual distance in some cases. The other is the patch based method [4, 5], which estimates the para￾metric distance between a Loop surface patch and its control mesh. Wu et al. presented an accurate error measure [4, 5]. But their estimate is de￾pendent on recursive subdivision, thus can not be used to pre-compute the 1

error bound after several steps of subdivision or predict the recursion depth within a user-specified tolerance. Our tochniqc helongs to the latter.Based on Stam's parametrizntion 6):we derive a conservative but safe bound on the maximum distance be- tween a Loop surface patch and its control mesh.An cificient subdivision depth formula is abo given 2 Second Order Norm and Convergence Rate Without loss of grnerality,we nssume the initial trianpular control mesh has been subdivided at least once,isolating the extraordinary vertioes so that cach face is n trinngle and contains at most one ctraordinary vertex. 2.1 Distance and Second Order Norm Given a ooutrol mesh of a Loop subdivision surface S.for each interior trinngle face F in the control mesh,there is a corresponding surface patch snt山e limit surface S. 4中 4中 + (a (b) Figure 1:(a)Ordering of the coutrol vertices of an extraordinary patch S of valence n.(b)Ordering of the new control vertices (solid dots)after one step of Loop subdivision. The distence between a Loop patch S and the corresponding triangle F is definexl as the maximnun distance between S(c,w)and F(v,w),that is, max lIS(t,w)-F(v,w), le.u)cn where is the unit triangle !{(e,w)e [0,1]and w e 0.1), S(e,w)is the Stam's parametrization of S over 52,and F(e,w)is the linear pxarametrization of E over 2

error bound after several steps of subdivision or predict the recursion depth within a user-specified tolerance. Our technique belongs to the latter. Based on Stam’s parametrization [6], we derive a conservative but safe bound on the maximum distance be￾tween a Loop surface patch and its control mesh. An efficient subdivision depth formula is also given. 2 Second Order Norm and Convergence Rate Without loss of generality, we assume the initial triangular 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 and Second Order Norm Given a control mesh of a Loop subdivision surface S˜, for each interior triangle face F in the control mesh, there is a corresponding surface patch S in the limit surface S˜. 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 b b b b b 5 6 n 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 S 1 0 S 1 1 S 1 2 S 1 3 (b) Figure 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. The distance between a Loop patch S and the corresponding 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 , where Ω is the unit triangle Ω = {(v, w)|v ∈ [0, 1] and w ∈ [0, 1 − v]} , S(v, w) is the Stam’s parametrization of S over Ω, and F(v, w) is the linear parametrization of F over Ω. 2

+12 rtires P!=1. he level-l coutro or of th f subdivis as the 2.2 Convergence Rate edoderoeashMgaadhfsihgyteolowitgeureeete 哈≤r,j≥0 3

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 extraordi￾nary 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 : i ≤ i ≤ n + 9} of the n + 6 vertices of Π: M = max{kP1 + P2 − Pn+1 − P3k, {kP1 + Pi − Pi−1 − P(i−1)%n+2k : 3 ≤ i ≤ n + 1}, 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 + 6} . (1) M is also called as 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 ver￾tices 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 , respectively (see Figure 1b), where Π1 0 = {P1 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 reg￾ular triangular patches [6]. Following the definition in Euatioin (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 con￾trol 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. 2.2 Convergence Rate If second order norms Mk+j 0 and Mj 0 satisfy the following recurrence in￾equality Mk+j 0 ≤ rj (n)Mk 0 , j ≥ 0 , (2) 3

ned ini e=h2.Rt9 tee折ata花6nm 4m=ia艺hl d芝a=t =a创.n>3 rder no 4

where rj (n) is a constant which depends on n, the valence of the extraor￾dinary vertex, and r0(n) ≡ 1. We call rj (n) the j-step convergence rate of second order norm. In the following, we estimate rj (n), j ≥ 1 by solving constrained min￾imization problems. 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 can 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. It follows that kα k+1 l k ≤ nX +9 i=1 kx l iα k i k ≤ nX +9 i=1 |x l i |kα k i k ≤ nX +9 i=1 |x l i |Mk 0 . 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 . (3) 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 mini￾mization 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 extraor￾dinary 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 (3) with α k+1 1 replaced by α k+j 1 . Table 1 shows the convergence rates of the second order norm for 3 ≤ n ≤ 10. The convergence rate rj (6) = ( 1 4 ) j is sharp for regular patches. 4

Table 1:Comvergeace rates r)12.3 8 0 3 Distance Bound 3.1 Regular Patch (a} b e quartic box spline basis functions.The 5,-∑h,m (5

Table 1: Convergence rates ri(n), i = 1, 2, 3 n 3 4 5 6 7 8 9 10 r1(n) 0.437500 0.382813 0.540907 0.250000 0.674025 0.705136 0.726355 0.741711 r2(n) 0.082031 0.142700 0.258367 0.062500 0.372582 0.402608 0.424000 0.439960 r3(n) 0.020752 0.053148 0.118899 0.015625 0.197695 0.219995 0.236377 0.248872 3 Distance Bound 3.1 Regular Patch 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 (a) b b b Ω (0,0) u=1 (0,1) (1,0) w=1 v=1 (b) Figure 2: (a) Control vertices of a quartic box spline patch and their order￾ing. (b) Parameter domain Ω = {(v, w)|v ∈ [0, 1] and w ∈ [0, 1 − v]} with the third parameter u = 1 − v − w. If S is a regular Loop patch, then S(v, w) can be expressed as a quartic box spline patch defined over the unit triangle Ω with control vertices pi , 1 ≤ i ≤ 12 (see Figure 2) as follows: S(v, w) = X 12 i=1 piNi(v, w) , (4) where Ni(v, w), 1 ≤ i ≤ 12 are the quartic box spline basis functions. The expressions of Ni(v, w) refer to [6]. S is also a quartic triangular B´ezier patch, thus S(v, w) can be written in terms of Bernstein polynomials [8]: S(v, w) = X 15 i=1 biBi(v, w) , (5) 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 5

between the standard representation B(u,w)= 0.i+j+k=4adB,(u,o).1≤i≤15is6,i.k)一 k主 2+k+1. u do ao thua dwu 6) ecan 、 rization of the center triangle F={pa,pr,ps}(sx F,=p4+p:+p 6

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. 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) Figure 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. The 15 × 12 matrix T which converts from the 12 control vertices to the 15 B´ezier points can be computed with the algorithm developed in [7]: T = 1 24                            2 2 2 12 2 0 2 2 0 0 0 0 1 0 3 12 1 0 4 3 0 0 0 0 0 1 1 12 3 0 3 4 0 0 0 0 0 0 4 8 0 0 8 4 0 0 0 0 0 0 1 10 1 0 6 6 0 0 0 0 0 0 0 8 4 0 4 8 0 0 0 0 0 0 3 4 0 1 12 3 0 0 1 0 0 0 1 6 0 0 10 6 0 0 1 0 0 0 0 6 1 0 6 10 0 0 1 0 0 0 0 4 3 0 3 12 1 0 1 0 0 0 2 2 0 2 12 2 0 2 2 0 0 0 1 3 0 0 12 4 0 1 3 0 0 0 0 4 0 0 8 8 0 0 4 0 0 0 0 3 1 0 4 12 0 0 3 1 0 0 0 2 2 0 2 12 2 0 2 2                            The linear parametrization of the center triangle F = {p4, p7, p8} (see Figure 4(a)) is F(v, w) = up4 + vp7 + wp8 , where u = 1 − v − w. By the linear precision property of the Bernstein polynomials [8], we can express F(v, w) as the following quartic B´ezier form: F(v, w) = X 15 i=1 biBi(v, w) , (6) 6

=Ff,)i.k≥0.0≤i+k≤4 are the Be2e oepu s三-EBe时 Ve n botd,-万in terms of,-万l≤ii,irc 6=bu=ius= 2=函==加=2=4= d==4= 瓜=函=M=司 natz出&dt te th anniedthmaaaeotetmmseado 3.2 Extraordinary Patch 时={0e,rlv∈2k,2 and we0,2kH-} =(e,w).2-*]and w e [2- 7

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. It is obvious that p4 = b1, p7 = b11, p8 = b15. Hence, 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) . (7) We can bound kbi − bik in terms of M as kbi − bik ≤ δiM. Here δ1 = δ11 = δ15 = 1 2 , δ2 = δ3 = δ7 = δ10 = δ12 = δ14 = 1 4 , δ4 = δ6 = δ13 = 1 6 , δ5 = δ8 = δ9 = 1 12 . It is easy to know that P15 i=1 δiBi(v, w) reaches its maximum at the three corners of Ω, that is, max (v,w)∈Ω X 15 i=1 δiBi(v, w) = 1 2 . Then we have a bound on the maximal distance between S(v, w) and F(v, w) as stated in the following theorem: Theorem 3.1 The distance between a regular Loop patch S and the corre￾sponding center triangle F is bounded by max (v,w)∈Ω kS(v, w) − F(v, w)k ≤ 1 2 M . 3.2 Extraordinary Patch An extraordinary Loop patch S 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 accord￾ingly, (see Figure 4), 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]} . 7

g =(w)Ic0 and wc0,2 he distan Is,四)-re,lss(c,j-上5（，o +E-时e1-时-四侧 Lemma3.2f(e,）∈疏，m=l,23,hea 时s{m 【nmmn.3f(,)e哈：then for any0≤i≤k-2 ur haxe g卡2学end时 8

Ω 1 1 Ω 1 2 Ω 1 3 Ω 2 1 Ω 2 2 Ω 2 3 Figure 4: Partition of the parameter domain Ω. Each tile Ωk m corresponds to a box spline patch S k m. And we denote the parameter space corresponding to the extraordinary subpatch S k 0 by Ω k 0 = {(v, w)| v ∈ [0, 2 −k ] and w ∈ [0, 2 −k − v]} . Let F(v, w), F k m(v, w) and F k 0 (v, w) be the linear parametrization of the center faces of the control meshes of S, S k m and S k 0 , respectively. Using the triangle inequality, for (v, w) ∈ Ω k m, m = 1, 2, 3, the distance between an extraordinary Loop patch S(v, w) and the corresponding triangle F(v, w) can be bounded as kS(v, w) − F(v, w)k ≤ kS k m(v, w) − F k m(v, w)k + kF k m(v, w) − F k−1 0 (v, w)k + X k−2 i=0 kF i+1 0 (v, w) − F i 0 (v, w)k . (8) With a proof analogous to the one of Lemma 3 in [9], we have the following two lemmas. Lemma 3.2 If (v, w) ∈ Ω k m, m = 1, 2, 3, then F k m(v, w) − F k−1 0 (v, w) ≤  3 8Mk−1 0 , m=1,3; 1 8Mk−1 0 , m=2. , where Mk−1 0 is the second order norm of S k−1 0 . Lemma 3.3 If (v, w) ∈ Ω k m, then for any 0 ≤ i ≤ k − 2 we have F i+1 0 (v, w) − F i 0 (v, w) ≤ ω(n)Mi 0 , where ω(n) = 5 8 − (3+2 cos(2π/n))2 64 , and Mi 0 is the second order norm of S i 0 and F 0 0 = F. 8

It follows that if (v.w)m=1.3. 5u-I≤+gg+fm2好 Ise,回)-Ft,ol≤e62nmla 时。 Se,)-F,l≤x(njA ahere ca=a2m。 yr Cm5rgonayepadsad IS(,叫)-F,w训≤C4,A≥1 o=M is the secored order zortn of s

It follows that if (v, w) ∈ Ω k m, m = 1, 3, kS(v, w) − F(v, w)k ≤ 1 2 Mk m + 3 8 Mk−1 0 + ω(n) X k−2 i=0 Mi 0 ≤ 1 2 Mk m + 3 8 Mk−1 0 + ω(n) X k−2 i=0 ri(n)M0 (9) Here, ri(n) is the i-step convergence rate of second order norm. Let m → ∞ in Equation (9), we get kS(v, w) − F(v, w)k ≤ ω(n) X∞ i=0 ri(n)M0 . Because {Ω k m}, k ≥ 1, m = 1, 2, 3, form a partition of Ω, we have the fol￾lowing theorem on the maximal distance between S(v, w) and F(v, w), (v, w) ∈ Ω: Theorem 3.4 The distance between an extraordinary Loop patch S and the corresponding triangle F is bounded by max (v,w)∈Ω kS(v, w) − F(v, w)k ≤ C∞(n)M0 , where C∞(n) = ω(n) X∞ i=0 ri(n) , and M0 = M is the second order norm of S. For a regular patch with n = 6, we have C∞(6) = ω(6)/(1 − r1(6)) = (3/8)/(1 − 1/4) = 1/2. The result in Theorem 3.1 is obtained again. How￾ever, there are no explicit expressions of ri(n) for general n, we have the following practical corollary for error estimation. Corollary 3.5 The distance between an extraordinary Loop patch S and the corresponding triangle F is bounded as max (v,w)∈Ω kS(v, w) − F(v, w)k ≤ Cλ(n)M0 , λ ≥ 1 , where Cλ(n) = ω(n) Pλ−1 i=0 ri(n) 1 − rλ(n) , and M0 = M is the second order norm of S. Table 2 gives the numerical results of the bound constants Cλ(n) (λ = 1, 2, 3). 9

Table 2:Comparison of C(m)(12.3) 88的照e 4 Subdivision Depth Estimation Ise,回-Pe,sCma A光a0空m.e ISe,回)-e,训sCrx ≤((,<t a1Cmenmaya8almmrnm k=,M+方 10

Table 2: Comparison of Cλ(n) (λ = 1, 2, 3) n 3 4 5 6 7 8 9 10 C1(n) 1.000000 0.784810 0.915862 0.5 1.052763 1.087084 1.111167 1.129655 C2(n) 0.880851 0.781290 0.873611 0.5 0.915630 0.914924 0.911327 0.907421 C3(n) 0.872850 0.780397 0.858623 0.5 0.875407 0.866176 0.856245 0.847476 4 Subdivision Depth Estimation Because the distance between a level-k control mesh and the surface patch S is dominated by the distance between the level-k extraordinary subpatch and its corresponding control mesh, which, according to Corollary 3.5, is kS(v, w) − F k (v, w)k ≤ Cλ(n)Mk , where Mk is the second order norm of S’s level-k control mesh. Assume ǫ > 0 and k = λlj + j, 0 ≤ j ≤ λ − 1, let kS(v, w) − F k (v, w)k ≤ Cλ(n)rk(n)M0 ≤ Cλ(n)(rλ(n))lj rj (n)M0 0, after k = min 0≤j≤λ−1 λlj + j steps of subdivision on the control mesh of S, the distance between S and its level-k control mesh is smaller than ǫ. Here, lj =  log 1 rλ(n)  rj (n)Cλ(n)M0 ǫ  , 0 ≤ j ≤ λ − 1, λ ≥ 1 , where rj (n) and Cλ(n) are the same as in Corollary 3.5, M0 = M is the second order norm of S. In particular, for a regular patch, k = l log4  M0 2ǫ m. Assume the second order norm M0 = 1, and the error tolerance ǫ = 0.01. Table 3 shows the results of subdivision depths computed with different λ = 1, 2, 3. 10