- Email: [email protected]
Non-stationary subdivision schemes for surface interpolation based on exponential polynomials
Applied Numerical Mathematics 60 (2010) 130–141
Contents lists available at ScienceDirect
Applied Numerical Mathematics www.elsevier.com/locate/apnum
Non-stationary subdivision schemes for surface interpolation based on exponential polynomials ✩ Yeon Ju Lee a , Jungho Yoon b,∗ a b
Department of Mathematical Sciences, KAIST, 373-1, Yuseong-gu, Daejeon, 305-701, Republic of Korea Department of Mathematics, Ewha Womans University, Seoul, 120-750, Republic of Korea
a r t i c l e
i n f o
Article history: Received 12 November 2008 Received in revised form 9 September 2009 Accepted 14 October 2009 Available online 23 October 2009 MSC: 41A05 41A25 41A30 65D10 65D17
a b s t r a c t This paper is concerned with non-stationary interpolatory subdivision schemes that can reproduce a large class of (complex) exponential polynomials. It enables our scheme to exactly reproduce the parametric surfaces such as torus and spheres. The subdivision rules are obtained by using the reproducing property of exponential polynomials which constitute a shift-invariant space S. In this study, we are particularly interested in the schemes based on the known butterfly-shaped stencils, proving that these schemes have the same smoothness and approximation order as the classical Butterfly interpolatory scheme. © 2009 IMACS. Published by Elsevier B.V. All rights reserved.
Keywords: Non-stationary subdivision Exponential polynomial Interpolation Asymptotical equivalence Smoothness Approximation order
1. Introduction Subdivision scheme is an efficient method for generating smooth curves and surfaces from a finite set of control points in Rd , d 1. A non-stationary subdivision scheme consists of recursive refinements of an initial sparse sequence with the use of rules that may vary from level to level but are the same everywhere on the same level. Starting with the control points P 0 = { pn0 : n ∈ Zd }, we define new sets of points P k = { pnk : n ∈ Zd } generated by the relation
P k = S a[k−1] · · · S a[0] P 0 . Non-stationary subdivision schemes are useful because they can provide design flexibility and the masks can be adapted to the geometrical configuration of the given data. Exponential (or trigonometric) polynomials and B-splines were used as a key ingredient for the construction of non-stationary subdivision schemes in [1,12,13,16,17,19–21]. Gaussian-based interpo✩ Yeon Ju Lee was supported by the Korea Research Foundation Grant funded by the Korean Government (MOEHRD) (KRF-2007-357-C00007). Jungho Yoon was supported by Priority Research Centers Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education, Science and Technology (2009-0093827). Corresponding author. E-mail addresses: [email protected] (Y.J. Lee), [email protected] (J. Yoon).
*
0168-9274/$30.00 © 2009 IMACS. Published by Elsevier B.V. All rights reserved. doi:10.1016/j.apnum.2009.10.005
Y.J. Lee, J. Yoon / Applied Numerical Mathematics 60 (2010) 130–141
131
Fig. 1. Exact reconstruction of a sphere and a torus.
latory schemes were studied in [14]. A general treatment of stationary schemes can be found in the selected references in [3,4,6,7,9,8]. In this study, instead of thinking about generating new vertices in R3 by taking local averages of vertices in the control polyhedron, we are thinking of scalar values which are assigned to the vertices of a triangulation in R2 . Therefore, starting with a given initial data f 0 = { f n0 ∈ R: n ∈ Z2 }, we consider f k = { f nk ∈ R: n ∈ Z2 } generated by the relation
f kj +1 =
[k]
a j −2n f nk ,
k ∈ Z+ , j ∈ Z2 .
(1)
n∈Z2
[k]
The set of coefficients a[k] := {an } is termed the mask of the rule at level k. We denote this rule by S a[k] and the corresponding non-stationary scheme by { S a[k] }. It is common to assume that for each level k, only a finite number of coefficients [k]
an are non-zero so that changes in a control point affect only its local neighborhood. This property clearly facilitates the practical implementation. A subdivision scheme is called stationary when the mask is independent of the levels; then we use the notation a := {an }. We denote the rule at each level by S a and the corresponding stationary scheme by { S a }. Also, a subdivision scheme is said to be uniformly convergent if for the given initial data δ = { f n0 = δn,0 : n ∈ Z2 }, there exists a function φ0 ∈ C 0 (R2 ), φ0 ≡ 0, satisfying
lim sup f nk − φ0 2−k n = 0,
k→∞ n∈Z2
(2)
where δn,0 denotes the Kronecker delta. In particular, φ0 is called the basic limit function of { S a[k] }, formally defined by
φ0 = lim S a[k] · · · S a[0] δ.
(3)
k→∞
Interpolatory subdivision schemes refine data by inserting values corresponding to intermediate points, using linear combinations of neighboring points. They are particularly attractive since they can control the shape of the resulting subdivision curves/surfaces in a more intuitive way than approximating schemes. The general form of their refinement rules (based on a planar parametric domain) is given as follows: for any j ∈ Z2 and k ∈ Z+ , 1 k f 2k+ j = fj, 1 f 2k+ = j +γ
[k]
aγ +2n f kj−n ,
γ ∈ E 2 \ 0,
n∈Z2
where E 2 is the extreme points of [0, 1]2 , i.e., E 2 = {(0, 0), (0, 1), (1, 0), (1, 1)}. A typical example of such stationary scheme is the four-point scheme by Dyn, Gregory, and Levin [10], and it is extended to Deslauriers–Dubuc schemes [8], where finer level points are determined by local polynomial interpolation of the coarse level points. The well-known generalization of the stationary 4-point scheme to the triangular mesh is the butterfly scheme [11]. It provides C 1 surfaces under regular topology (i.e., if every vertex in the triangular mesh has valence 6) and reproduces polynomials up to degree 3 by using a minimal support. Non-stationary 4-point interpolatory schemes were developed by Beccari et al. [1,2]. Also, in the paper [18], Jena et al. studied non-stationary 4-point interpolatory scheme reproducing elements in the linear space spanned by {1, x, cos(ax), sin(ax)}. All these schemes were developed for generating smooth curves but the corresponding non-stationary butterfly scheme and their analytical properties have not been investigated. Construction of non-stationary subdivision rules based on butterfly stencils are rather simple (see [5] or Section 2). However, fundamental questions on non-stationary subdivision schemes such as convergence, smoothness and approximation orders are yet to be answered. For this reason, in this study, we are concerned with a family of non-stationary interpolatory subdivision schemes based on the known butterfly-shaped stencils. Each scheme in this family is exact on a certain shift-invariant space (say, S) consisting of (complex) exponential polynomials. It enables our scheme to exactly reproduce the parametric surfaces such as torus and sphere (see Fig. 1). Our specific contribution is as follows: Under a certain condition of S, we will see that this new scheme provides the same smoothness C 1 as the original butterfly scheme. Moreover, we show that the non-stationary scheme provides the same approximation order four as the original butterfly scheme with the additional advantage that there are flexibilities in the choice of S and frequency factors.
132
Y.J. Lee, J. Yoon / Applied Numerical Mathematics 60 (2010) 130–141
Fig. 2. Stencils of the butterfly scheme. Larger dots indicate the insertion point.
In the following we use the two-dimensional index notation. First, let
Z2+ := (α1 , α2 ) ∈ Z2 : αi 0, i = 1, 2 .
The set of all complex numbers is denoted by C. For α = (α1 , α2 ) ∈ Z2+ , we set α ! = α1 !α2 ! and |α |1 := α1 + α2 . For α α x = (x1 , x2 ) ∈ R2 , xα = x1 1 x2 2 . The space of algebraic polynomials of degree n is denoted by Πn . For a matrix M, M∞ indicates its ∞-norm. Also, for a given smooth function f , denote by T f ,τ the Taylor polynomial f of degree 3 around τ and let R f ,τ be its remainder, i.e.,
T f ,τ (x) =
(x − τ )ν f (ν ) (τ )/ν !,
R f ,τ (x) =
(x − τ )ν f (ν ) (ξ )/ν !,
(4)
|ν |1 =4
|ν |1 3
for some ξ between x and
τ.
2. Construction of non-stationary schemes The non-stationary subdivision schemes proposed in this paper are constructed in a way of reproducing the exponential (or trigonometric) polynomials of the form
ϕ (x) = xα eβ·x , x ∈ R2 , α ∈ Z2+ , β ∈ C2 .
(5)
In practice, xα ’s are chosen to be low degree polynomials. Let
S = span{ϕ | = 1, . . . , 8}
(6)
be a shift-invariant space with linearly independent functions ϕ , = 1, . . . , 8, of the form as in (5). Refining the triangular mesh data based on the butterfly-shaped stencils involves three groups of vertices to evaluate three types of new vertices as shown in Fig. 2. Here and in the sequel, X j with j = 1, 2, 3 indicate the three types of stencils and p j /2 the insertion points corresponding to X j at level 0. Letting
p 1 = (1, 0),
p 2 = (0, 1),
p 3 = (1, 1),
(7)
for each level k = 0, 1, . . . , and stencil X j , the non-stationary subdivision rule is constructed by solving the linear system
ϕ p j 2−k−1 =
n∈X j
[k]
a p −2n ϕ n2−k , j
= 1, . . . , 8.
(8)
This linear system can be written in the matrix form
a[k] = B[k]
−1 [k]
b
(9)
,
where
[k]
a[k] = a p −2n : n ∈ X j , j B[k] ( , n) =
b[k] =
ϕ p j 2−k−1 : = 1, . . . , 8 ,
ϕ n2−k : n ∈ X j , = 1, . . . , 8 .
For the unique solution of the linear system (9), it is required that dim(S|2−k X j ) = dim S for each stencil X j . We will see later that the mask a[k] converges to the mask of the classical butterfly scheme as k → ∞. Example 2.1 (Sphere and torus). Let us define S as follows: For x = (u , v ),
S := span{sin u , cos u , sin v , cos v , sin u sin v , cos u cos v , sin u cos v , cos u sin v }.
Y.J. Lee, J. Yoon / Applied Numerical Mathematics 60 (2010) 130–141
133
Fig. 3. Reconstruction of a torus. The images from left to right show the subdivision results using our scheme and using the classical butterfly scheme.
Fig. 4. Surface metamorphosis. The images from left to right show the subdivision of a sphere while modulating the frequency of the trigonometric basis functions.
Then the three types of refinement rule are obtained by solving the linear system (9) using the butterfly stencils. Recursive application of this non-stationary subdivision rule exactly generates a sphere whose parametric equations are of the form
f (u , v ) = r sin u cos v ,
g (u , v ) = r sin u sin v ,
h(u , v ) = r cos u ,
(10)
where 0 v 2π and 0 u π ; see Fig. 1. It seems that there are two vertices (at poles) which look like extraordinary points with valences not equal to 6 but this is not the case. Because these functions in (10) are 2π -periodic, we use larger extended domain which are mapped into the same position. Eventually, no extraordinary stencil is necessary. The flowerlike shape obtained in Fig. 4 is due to this extended parameterization domain whose border is mapped onto R3 with this particular spinning covering (see Example 2.2). Similarly, the non-stationary scheme reconstructs a torus whose parametric equations are of the form
f (u , v ) = (1 + cos u ) cos v ,
g (u , v ) = (1 + cos u ) sin v ,
h(u , v ) = sin u ,
where 0 u , v 2π ; see Fig. 1. Also, a subdivision result by using the classical butterfly scheme is given in Fig. 3. Example 2.2. As shown in the images of Fig. 4, a flower-shaped model can be created from an initial, spherical mesh by modulating the frequencies of the trigonometric polynomials (see Example 2.1). In this example, our subdivision rules are obtained in a way of reproducing functions in the space
S := span{sin qu , cos qu , sin qv , cos qv , sin qu sin qv , cos qu cos qv , sin qu cos qv , cos qu sin qv } with q = 1.8. Trigonometric functions with higher frequencies have more oscillations and higher curvatures. Because our scheme reproduces these functions accurately, the limit surface becomes more oscillatory in order to capture these features. Since our scheme is implemented by local interpolation using exponential polynomials, we may perform the metamorphosis of a given initial coarse mesh into an interesting yet controllable shape by adopting different frequencies. The idea of surface modulation is to model a predictable shape using only a simple shape like a sphere. However, in order for such a metamorphosis to be useful, a deep study is necessary in the future. 3. Asymptotic equivalence Typical tools used for the analysis of the stationary schemes are Fourier transform, eigen analysis and Laurent polynomials [3,9]. However, these techniques are not applicable for the non-stationary case. Thus, for the analysis of non-stationary schemes, we adopt the notion of asymptotical equivalence [12]: a (non-stationary) subdivision scheme { S a[k] } is asymptotically equivalent to { S a } if
k∈Z+
where
S a[k] − S a ∞ < ∞,
(11)
134
Y.J. Lee, J. Yoon / Applied Numerical Mathematics 60 (2010) 130–141
S a[k] ∞ = max
[k] a
α −2n : α ∈ E 2
n∈Z2
with E 2 the extreme points of [0, 1]2 . 4. Analysis of convergence The proposed non-stationary butterfly scheme reproduces functions in the space S with dim S = 8. Hence, we choose an index set
T = ν ∈ Z2+ : |ν |1 3, ν = (1, 2), (2, 1) .
(12)
Note that #T = 8. Here, we assume that the set T is sorted by the polynomial ordering of a sequence of multi-indices, that is, an index α := (α1 , α2 ) comes before β := (β1 , β2 ) if |α |1 < |β|1 , or if |α |1 = |β|1 and α1 < β1 . Moreover, for notational convenience, we set
{ϕβ : β ∈ T} := {ϕ1 , . . . , ϕ8 }.
ϕβ satisfies the condition ϕβ(β) (0) = 0 with β ∈ T. Then, let us define the function
It is immediate that
W [ϕ1 , . . . , ϕ8 ](x) := det W(x) with the matrix W(x) given by
W(x) =
(ν )
ϕ j (x): j = 1, . . . , 8, ν ∈ T .
(13)
This function W [ϕ1 , . . . , ϕ8 ](x) is, in fact, the Wronskian of assume that
ϕ1 , . . . , ϕ8 on two-dimensional case. Throughout this paper, we
W [ϕ1 , . . . , ϕ8 ](0) = 0. We cite a basic result about the asymptotically equivalent subdivision schemes. Theorem 4.1. (See [12].) Let { S a } be a stationary subdivision scheme, and let { S a[k] } be a non-stationary subdivision which is asymp[k]
totically equivalent to { S a }. Assume that supp({an }), supp({an }) ⊂ [0, N ]2 , k ∈ Z+ , N < ∞. Then { S a[k] } is C 0 if { S a } is C 0 . Further, if
S a[k] − S a ∞ c2−k ,
k ∈ Z+ ,
then the basic limit function φ0 of { S a[k] } in (3) is Hölder continuous of exponent s > 0, i.e., φ0 ∈ H s . We are now ready to show that the non-stationary subdivision scheme { S a[k] } is asymptotically equivalent to the classical butterfly scheme { S a }, which implies that { S a[k] } is a convergent scheme. Recall that the mask {a p j −2n : n ∈ X j } of the butterfly scheme reproduce polynomials in Π3 in the sense that
a p j −2n p n2−k = p 2−k−1 p j ,
∀ p ∈ Π3 ,
(14)
n∈X j
with p j ( j = 1, 2, 3) in (7). [k]
Theorem 4.2. Let {an } be the mask of the butterfly scheme { S a }, and let {an } be the mask of the non-stationary scheme { S a[k] } reproducing functions in the space S = span{ϕβ : β ∈ T}. Then, there exists a constant c > 0 such that
[k]
maxan − an c2−k , n∈Z2
k K ∈ Z+ . [k]
[k]
Proof. Due to the fact that a2n = a2n = δn,0 with n, 0 ∈ Z2 , it suffices to estimate only the difference |a p j −2n − a p −2n | j with p j ( j = 1, 2, 3) in (7). For each β ∈ T, define a function Φβ by
Φβ (x) :=
c β (γ )ϕγ (x)
(15)
γ ∈T
so that the coefficients c β (γ ), (γ )
(γ )
γ ∈ T, are obtained by solving the linear system
Φβ (0) = ϕβ (0)(1 − δβ,γ ).
(16)
Y.J. Lee, J. Yoon / Applied Numerical Mathematics 60 (2010) 130–141
135
The existence of the unique solution of this linear system is clear from the condition det W(0) = 0. Further, since the mask {an[k] } reproduces ϕβ with β ∈ T (see (8)), it also reproduces Φβ such that we can get
ϕβ 2−k−1 p j − Φβ 2−k−1 p j =
n∈X j
[k]
a p −2n j
ϕβ n2−k − Φβ n2−k .
(17)
Now, for a smooth function g, put T g = T g ,0 , that is, the Taylor polynomial of g of degree 3 around 0 (see (4)). The condition (16) implies that
T ϕβ (x) − T Φβ (x) =
xβ
β!
ϕβ(β) (0) +
xν (ν )
ν! ν =(1,2),(2,1)
ϕβ (0) − Φβ(ν ) (0) .
Then, for any β ∈ T, define Pβ (x) by
Pβ (x) :=
(β)
xβ ϕβ (0)/β! (β) xβ ϕβ (0)/β! +
(2−k ·) = 2−|β|1 k P
Note that Pβ Eq. (17) becomes
if |β|1 2,
xν
(ν ) (ν ) ν =(1,2),(2,1) ν ! (ϕβ (0) − Φβ (0))
β ∈ T. Then, there exists a (smooth) function Rβ,k such that after some calculation,
β for any
Pβ ( p j /2) + Rβ,k ( p j /2) =
otherwise.
[k]
n∈X j
a p −2n Pβ (n) + Rβ,k (n) , j
= 2−k ,
(18)
with β ∈ T. Here, it is obvious that Rβ,k ( p j /2) and Rβ,k (n), n ∈ X j , are bounded by a constant c > 0 independent of k ∈ N and n ∈ X j . On the other hand, using the fact that the mask {a p j −2n : n ∈ X j } of the classical butterfly scheme reproduces Π3 in the sense of (14), we have
a p j −2n Pβ (n) = Pβ ( p j /2).
n∈X j
This linear system can be written in the matrix form
T·a=b
(19)
with
T = Pβ (n): n ∈ X j , β ∈ T , a = (a p j −2n : n ∈ X j ),
b = Pβ ( p j /2): β ∈ T .
(20) [k]
[k]
Then, in view of (18), (19) and (20), one can build a linear system which uniquely defines the mask a j := (a p −2n : n ∈ X j ), j in the form
(T + Fk )a[jk] = b + Gk (β)
with Fk ∞ , Gk ∞ c 1 < ∞. Since ϕβ (0) = 0, T is invertible and hence, T + Fk is also invertible for a sufficient small > 0. Note that the O ( ) perturbation of the non-singular matrix T causes the O ( ) perturbation of its inverse [15]. Thus, we can write [k]
a j = T−1 + F˜ k (b + Gk )
= T−1 b + T−1 Gk + F˜ k (b + Gk ) ,
with F˜ k ∞ c 2 < ∞, where T−1 b = a. Since T−1 Gk + F˜ k (b + Gk )∞ c 3 < ∞, we conclude that there is a constant c > 0 [k] such that a j − a∞ c2−k for any k K . 2 5. Smoothness 5.1. Sufficient condition for smoothness An analysis of the smoothness of univariate non-stationary subdivision schemes is discussed in [14]. In this paper, the result of [14] is extended to the case of non-stationary butterfly scheme. Specifically, our goal is to prove that the non-stationary scheme reproducing the space S in (6) has the same smoothness as the original butterfly scheme. To this end, we first show that a factor (1 + w k zλ ), λ ∈ Z2 , in the Laurent polynomial a[k] ( z) plays a role of smoothing factor if |1 − w k | c2−k . We then infer results on the smoothness of the suggested interpolatory non-stationary scheme from the smoothness of a stationary scheme, which is asymptotically equivalent to it.
136
Y.J. Lee, J. Yoon / Applied Numerical Mathematics 60 (2010) 130–141
Lemma 5.1. Let the Laurent polynomial of the non-stationary scheme { S a[k] } be of the form
a[k] ( z) =
1 2
1 + w k zλ b[k] ( z),
k > K ∈ Z+ ,
where λ ∈ Z2 and { S b[k] } is of compact support and C 0 . Suppose that
|1 − w k | c2−k ,
k K ∈ Z+ .
(21)
Then the basic limit function φ0 in (3) satisfies φ0 , ∂λ φ0 ∈ C where δ := {δn,0 : n ∈ Z }. 0
2
Remark. For the butterfly scheme in regular triangulation, we need to consider only the three directions: λ = (1, 0), (0, 1), and (1, 1). Proof. Due to Theorem 4.1, we can find that { S a[k] } is C 0 with the basic limit function φ0 of the form
φ0 =
φb (· − λt )h(t ) dt , R
where φb and h are the basic limit functions of { S b[k] } and { S 1+ w k z } respectively, and where h is bounded and supp{h} = [0, 1). Here, by virtue of Theorem 4.1, Theorem 4.2 and (21), φb is Hölder continuous of exponent s > 0. Then, it follows that φ0 ∈ C 0 . In order to prove ∂λ φ0 ∈ C 0 , it is sufficient to show that φ0 = g ∗ h ∈ C 1 (R) for any bounded and Hölder continuous function g of exponent s > 0, i.e., g ∈ H s . In fact, the property φ0 = g ∗ h ∈ C 1 (R) can be proved directly by using Lemma 2.5 in [14], 2 In order to guarantee the condition (21), we require a[k] to satisfy the following stronger condition: Condition A. A non-stationary subdivision scheme { S a[k] } satisfies Condition A if the corresponding Laurent polynomial a[k] ( z) is of the form (26) and for = 1, 2,
m ∂ [k] −(4−m)k a (θ ) , c2 ∂ zm
m = 0, . . . , 3 ,
with θ1 = (−1, 1) and θ2 = (1, −1). In what follows, we will prove that the Laurent polynomial a[k] ( z) associated with the scheme { S a[k] } satisfies the Condition A. It in fact implies that the corresponding scheme { S a[k] } has the same smoothness C 1 as the classical butterfly scheme. 5.2. Laurent polynomials of butterfly schemes To simplify the presentation of a subdivision scheme and its analysis, it is convenient to assign the Laurent polynomial to each subdivision mask. Denoting z = ( z1 , z2 ), the bivariate Laurent polynomial corresponding to the original butterfly scheme can be put in the following factored form [9]:
a( z ) =
an zn = 2−1 1 + z1−1 1 + z2−1 1 + z1−1 z2−1 z1 z2 1 − 16−1 q( z1 , z2 ) ,
(22)
n∈Z2
where
q( z) = q( z1 , z2 ) = 2z1−2 z2−1 + 2z1−1 z2−2 − 4z1−1 z2−1 − 4z1−1 − 4z2−1
+ 2z1−1 z2 + 2z1 z2−1 + 12 − 4z1 − 4z2 − 4z1 z2 + 2z12 z2 + 2z1 z22 . The Laurent polynomial a( z) with z = ( z1 , z2 ) has single zeros along the curves z1 = −1, z2 = −1 and z1 z2 = −1. In particular, we will see that a( z) has zeros of order 4 at the points (−1, 1) and (1, −1). In the following, for convenience, we use the notation
e 1 = (1, 0)
and e 2 = (0, 1).
(23)
Lemma 5.2. Let a( z) in (22) be the Laurent polynomial of the classical butterfly subdivision scheme { S a }. Then, for any 0 m 3,
∂m a(θ ) = 0, ∂ zm
= 1, 2,
with
θ1 = (−1, 1) and θ2 = (1, −1).
(24)
Y.J. Lee, J. Yoon / Applied Numerical Mathematics 60 (2010) 130–141
137
Proof. In this proof, we consider only the case = 1 because the case = 2 can be proved by the same way. For 0 m 3, we can write
∂m cm,r an nre1 θ1n , m a(θ1 ) = ∂ z1 2 0r m
(25)
n∈Z
pj
where e 1 = (1, 0), for some suitable constants cm,r . Seeing that θ1 = (−1) j for j = 1, 2, 3 and a2n = δ0,n with n ∈ Z2 , we obtain from the polynomial reproducing property in (14) that
2−r
an nre1 θ1n = δr ,0 +
3
p j −2n
a p j −2n ( p j /2 − n)re1 θ1
j =1 n∈X j
n∈Z2
= δr ,0 +
3
p
θ1 j
a p j −2n ( p j /2 − n)re1
n∈X j
j =1
= δr ,0 +
3 (−1) j δr ,0 = 0. j =1
Combining it with (25), the lemma is proved.
2
According to Theorem 4.2, the non-stationary scheme { S a[k] } is asymptotically equivalent to the original butterfly scheme [k]
{ S a } with the convergence property |an − an | = O (2−k ). Recalling that the Laurent polynomial a(z), z = (z1 , z2 ), in (22) has zeros along the lines z1 = −1, z2 = −1, and z1 z2 = −1, we see that a[k] (−1, z2 ) = a(−1, z2 ) + O (2−k ) = O (2−k ) for any fixed z2 . Similarly, we can show that a[k] ( z1 , −1) = a[k] ( z1 , − z2 ) = O (2−k ). Thus, it is not difficult to deduce that the Laurent polynomial a[k] ( z) associated with S a[k] is of the form a[k] ( z) :=
[k]
an zn
n∈Z
:= 2−1 (1 + w 1,k z1 )(1 + w 2,k z2 )(1 + w 3,k z1 z2 )c [k] ( z1 , z2 )
(26)
with w j ,k ∈ C and
|1 − w j ,k | = o(1),
j = 1, 2, 3, as k → ∞,
for some suitable Laurent polynomial c [k] ( z). 5.3. Analysis of smoothness Let a( z1 , z2 ) be the Laurent polynomial of the classical butterfly scheme as in (22). It is known that the subdivision scheme corresponding to the Laurent polynomial
b j ( z1 , z2 ) = 2(1 + z j )−1 a( z1 , z2 ),
j = 1, 2,
is C 0 (see [9]). Since the non-stationary scheme { S a[k] } is asymptotically equivalent to the butterfly scheme { S a }, the scheme corresponding to [k]
b j ( z1 , z2 ) = 2(1 + w j ,k z j )−1 a[k] ( z1 , z2 ),
j = 1, 2,
(27)
with a[k] ( z1 , z2 ) in (26) is also C 0 if |1 − w j ,k | c2−k for any j = 1, 2. Thus, according to Lemma 5.1, the C 1 -smoothness of the non-stationary scheme { S a[k] } can be verified by showing that |1 − w j ,k | c2−k . Therefore, for this proof, we show that Condition A on { S a[k] } implies the condition (21) for all the factors in the representation (26). Lemma 5.3. Suppose that the Condition A holds for the Laurent polynomial a[k] ( z) in (26) which is associated with the non-stationary subdivision scheme { S a[k] }. Then, for each j = 1, 2,
|1 − w j ,k | c2−k ,
k K ∈ Z+ .
Proof. For simplicity, let
ak,1 ( z1 ) := a[k] ( z1 , 1)
and ak,2 ( z2 ) := a[k] (1, z2 ),
(28)
138
Y.J. Lee, J. Yoon / Applied Numerical Mathematics 60 (2010) 130–141
i.e., ak, j , j = 1, 2, are univariate Laurent polynomials. Then, since { S a[k] } is asymptotically equivalent to { S a }, we deduce from Lemma 5.2 that ak, j ( z j ) can be written as follows:
ak,1 ( z1 ) = ck,1 ( z1 )
4
(1 + rk,n z1 ),
n =1
ak,2 ( z2 ) = ck,2 ( z2 )
4
(1 + sk,n z2 ),
(29)
n =1
where ck, j (−1) = c j + o(1) for some constants c j = 0 ( j = 1, 2). Now, we will show that
|1 − rk,n |, |1 − sk,n | c2−k ,
k K ∈ Z+ ,
for any n = 1, . . . , 4. This is indeed enough for the proof of |1 − w j ,k | c2−k for k K ∈ Z+ . For this purpose, observe that
m ∂ m [k] d a ) = a (θ j k, j (−1) ∂ zmj dzm j
with θ j in (24). Thus, the Condition A is reduced to the univariate case as follows:
m d c2−(4−m)k . (− a 1 ) dzm k, j j
It enables us to apply the known result in [14, Lemma 2.6] by which we can conclude that |1 − rk,n | c2−k . Also, the bound |1 − sk,n | c2−k can be proved by applying exactly the same technique. Therefore, the relation (28) is proved. 2 We now prove that the Laurent polynomial a[k] ( z) of the scheme { S a[k] } satisfies the Condition A. It in fact implies that the corresponding scheme { S a[k] } has the smoothness C 1 .
[k]
Theorem 5.4. Let a[k] ( z) = n∈Z2 an zn be the Laurent polynomial associated with the non-stationary interpolatory scheme { S a[k] }. Then, for any 0 m 3, we have
m ∂ [k] −k(4−m) a (θ ) , c2 ∂ zm
= 1, 2,
(30)
with θ1 = (−1, 1) and θ2 = (1, −1). Proof. In this proof, we consider only the case = 1 because the case = 2 can be done by exactly the same way. Note that
[k] ∂ m [k] μm,r an nre1 θ1n , m a (θ1 ) = ∂ z1 2 0r m
n∈Z
where e 1 = (1, 0), for some constants
M r ,k :=
μm,r . Hence, in order to prove (30), it is sufficient to show that for any 0 r m, [k] an nre1 θ1n = O 2−k(4−r ) , as k → ∞.
n∈Z2
[k]
[k]
Since {an } is the mask of an interpolatory scheme, a2n = δn,0 with n, 0 ∈ Z2 . Hence,
[k]
a2n (2n)re1 = δr ,0 .
n∈Z2
It induces that
2−r (k+1) M r ,k = δr ,0 +
3 [k] re1 (−1) j a p −2n ( p j /2 − n)2−k j
(31)
n∈X j
j =1
where X j are the butterfly stencils of the type j = 1, 2, 3 (see Fig. 1) and p j /2 is the corresponding insertion point (see (7)). In order to estimate the term in (31), for the given k ∈ N and j = 1, 2, 3, define [k, j ]
Φre1 (x) := Φre1 (x) :=
8
=1
[k, j ]
cre1 ( )ϕ (x)
(32)
Y.J. Lee, J. Yoon / Applied Numerical Mathematics 60 (2010) 130–141
[k, j ]
[k, j ]
(ν ) Φre p j 2−k−1 = δre1 ,ν (−1)ν ν !, 1
ν ∈ T.
139
so that the coefficient vector ck := cre1 := (cre1 ( ): = 1, . . . , 8) T is obtained by solving the linear system
(33)
These equations can be written in the matrix form
Wk · ck = r with
Wk := W p j 2−k−1 :=
(ν )
ϕ
r := rre1 := δre1 ,ν (−1)|ν |1 ν !:
p j 2−k−1 :
ν ∈ T, = 1, . . . , 8 ,
ν ∈T .
Since W [ϕ1 , . . . , ϕ8 ](x) = W(x) is continuous and W [ϕ1 , . . . , ϕ8 ](0) is non-zero, there exists K ∈ N so that for any k K , the matrix Wk is non-singular and Wk−1 ∞ c for some constant c > 0. It implies that for any compact set M in R2 and α ∈ Z2+ ,
(α ) Φre 1
L∞
6 (α ) −1 ϕ W r∞ k (M ) ∞ L
=1
∞ (M )
cα,M ,
k K,
(34)
[k]
for some constant c α , M > 0. On the other hand, since the mask {an } reproduces duces Φre1 . Hence, by (33), we get
ϕβ with β ∈ T (see (8)), it also repro-
[k] δr ,0 = Φre1 2−k−1 p j = a p −2n Φre1 n2−k . j n∈X j
This together with (31) leads to
2−r (k+1) M r ,k =
3 [k] re1 . (−1) j −1 a p −2n Φre1 n2−k − ( p j /2 − n)2−k j j =1
(35)
n∈X j
τ = 2−k−1 p j , write Φre1 n2−k = T Φre1 ,τ n2−k + R Φre1 ,τ n2−k
Now, letting
with T Φre1 ,τ and R Φre1 ,τ in (4). Then, invoking (33), we see that
T Φre1 ,τ n2−k − ( p j /2 − n)2−k
re1
=
ν (ν ) −k−1 (n − p j /2)2−k Φre 2 p j /ν !. 1
ν =(1,2),(2,1)
[k]
Due to Lemma 4.2, an = an + O (2−k ). Using the polynomial reproducing property of the mask {an } of the butterfly scheme (14), it is immediate that for each ν = (1, 2), (2, 1),
[k] −k ν −k ν −4k ( ( a n − p / 2 ) 2 a n − p / 2 ) 2 = c2−4k p − 2n j j j p j −2n + c2 n∈X j
(36)
n∈X j
for a constant c > 0. Next, by applying (34), it is obvious that
−k [k] c2−4k , a R n2 p j −2n Φre1 ,τ
k K.
n∈X j
This together with (35) and (36), we obtain the required result.
2
We are now ready to provide the main theorem of this section. Theorem 5.5. If { S a[k] } is the non-stationary interpolatory subdivision scheme reproducing the shift-invariant space S in (6), then { S a[k] } is C 1 , i.e., it has the same smoothness as the stationary butterfly subdivision scheme { S a }. Proof. From Lemmas 5.1, 5.3, and Theorem 5.4, the proof is immediate.
2
140
Y.J. Lee, J. Yoon / Applied Numerical Mathematics 60 (2010) 130–141
6. Approximation order An important issue in the implementation of subdivision algorithm is how to actually attain the original function as close as possible if the given initial data f 0 is sampled from an underlying function. A high quality reconstruction scheme should guarantee that the approximation error decreases when the quality of the sample increases. This approximation power is usually quantified by the approximation order. Suppose that f 0 is of the form f 0 := { f n0 = f (2−k0 n): n ∈ Z2 } for some positive integer k0 with a smooth function f . Then our concern is to find the largest exponent m > 0 such that
∞ f − f
L∞ (M )
C 2−mk0
with a constant C > 0 independent of k0 , where M is a compact set in R2 . The exponent m is called the approximation order of the subdivision scheme. It is well known that the original butterfly scheme reproduces cubic polynomials, and hence, it provides the approximation order of four. The goal of this section is to show that the non-stationary scheme { S a[k] } provides the same approximation order four as the original butterfly scheme. It is of interest to see that the butterfly-shaped stencils consist of just 8 vertices. Usually, it is expected to have 10 vertices to obtain the approximation order 4. Since the scheme { S a[k] } is uniformly convergent, its limit function can be written as
f ∞ = lim S a[k0 + ] · · · S a[k0 ] f 0 = →∞
φk0 2k0 · −n f 2−k0 n
(37)
n∈Z2
where φk0 is the basic limit function defined as in (3). Indeed, since { S a[k] } reproduces functions in S, it is clear that
ϕ =
ϕ 2−k0 n φk0 2k0 · −n , ∀ϕ ∈ S.
n∈Z2
This observation leads to the proof of the approximation order of the subdivision scheme { S a[k] }. In this paper, we are particularly interested in approximating functions f in the Sobolev space
s W∞ ( M ) = g : R2 → R:
g (α ) |α |1 s
L∞
< ∞ , (M )
s ∈ Z+ .
4 Theorem 6.1. Let M be a compact set in R2 . Assume that f ∈ W ∞ ( M ) and f 0 := { f n0 = f (2−k0 n): n ∈ Z2 }. Then the non-stationary 4 interpolatory scheme { S a[k] } has the approximation order four with respect to functions in W ∞ ( M ).
Proof. Let x ∈ M be a fixed point and define a function
Φ :=
8
μ ϕ (· − x)
=1
so that the coefficient vector {μ : = 1, . . . , 8} is obtained by solving the linear system
Φ (ν ) (x) = f (ν ) (x),
ν ∈T
with T in (12). The existence of the unique solution of this linear system is clear because det W(0) = 0. Since S is shiftinvariant and φ ∈ S, ϕ (· − x) belongs also to S. Also, since Φ is a linear combination of φ , = 1, . . . , 8, and the scheme { S a[k] } is exact for such functions, we easily get the identity
Φ=
φk0 2k0 · −n Φ n2−k0 .
n∈Z
Then, using the property f (x) = Φ(x), it follows that
f (x) − f ∞ (x) = Φ(x) −
=
φk0 2k0 x − n f n2−k0
n∈Z2
φk0 2k0 x − n Φ n2−k0 − f n2−k0 .
n∈Z2
Recall from (4) that T g ,x indicates the Taylor polynomial of a smooth function g of degree 3 around x and R g ,x is its remainder. Due to the fact Φ (ν ) (x) = f (ν ) (x) with ν ∈ T, we obtain
T Φ,x n2−k0 − T f ,x n2−k0 =
ν =(1,2),(2,1)
It follows that
(n2−k0 − x)ν (Φ − f )(ν ) (x). ν!
Y.J. Lee, J. Yoon / Applied Numerical Mathematics 60 (2010) 130–141
f (x) − f ∞ (x) =
141
(n2−k0 − x)ν (Φ − f )(ν ) (x) ν ! ν =(1,2),(2,1) n∈Z2 k 0 + φk0 2 x − n R Φ,x n2−k0 − R f ,x n2−k0 . φk0 2k0 x − n
(38)
n∈Z2
Letting φ be the basic limit function of the classical butterfly subdivision scheme S a , i.e., φ = S a∞ δ with δ = {δn,0 : n ∈ Z2 }, we find from [12, Lemma 15] and Theorem 4.2 that φk0 − φ∞ c2−k0 . Moreover, applying the polynomial reproducing property of φ in (14), we obtain
ν φ 2k0 x − n 2k0 x − n = 0,
ν ∈ (1, 2), (2, 1) .
n∈Z2
Note that since φk0 is compactly supported, #{n ∈ Z2 : φk0 (2k0 x − n) = 0} c for any x. It implies that
k −k ν k k ν −3k0 0 0 0 0 φ 2 x − n n2 − x 2 φ 2 x − n 2 x − n = k0 k0 n∈Z2
n∈Z2
ν 2−3k0 φ 2k0 x − n 2k0 x − n + c2−k0 n∈Z2
c2
−4k0
for some constant c > 0. Also, since Φ (ν ) is bounded on any compact set M and f term in (38) has the property O (2−4k0 ). It finishes the proof. 2
4 ∈ W∞ ( M ), it is obvious that the second
References [1] C. Beccari, G. Casciola, L. Romani, A non-stationary uniform tension controlled interpolating 4-point scheme reproducing conics, Comput. Aided Geom. Design 24 (2007) 1–9. [2] C. Beccari, G. Casciola, L. Romani, An interpolating 4-point C 2 ternary non-stationary subdivision scheme with tension control, Comput. Aided Geom. Design 24 (2007) 210–219. [3] A. Cavaretta, W. Dahmen, C.A. Micchelli, Stationary subdivision, Mem. Amer. Math. Soc. 93 (1991) 1–186. [4] G. Chaikin, An algorithm for high speed curve generation, Computer Graphics and Image Processing 3 (1974) 346–349. [5] Y.-J. Choi, Y.-J. Lee, J. Yoon, B.-G. Lee, Y.-J. Kim, A new class of non-stationary interpolatory subdivision schemes based on exponential polynomials, in: Geometric Modeling and Processing – GMP 2006, in: Lecture Notes in Computer Science, vol. 4077, 2006, pp. 563–570 (SpringerLink, PDF-file). [6] E. Cohen, T. Lyche, R. Riesenfeld, Discrete B-spline and subdivision techniques in computer-aided geometric design and computer graphics, Computer Graphics and Image Processing 14 (1980) 87–111. [7] D. Doo, M. Sabin, Behaviour of recursively division surfaces near extraordinary points, Comput.-Aided Design 10 (1978) 356–360. [8] G. Deslauriers, S. Dubuc, Symmetric iterative interpolation, Constr. Approx. 5 (1989) 49–68. [9] N. Dyn, Subdivision schemes in computer-aided geometric design, in: W.A. Light (Ed.), Advances in Numerical Analysis Vol. II: Wavelets, Subdivision Algorithms, Radial Basis Functions, Oxford University Press, 1992, pp. 36–104. [10] N. Dyn, J.A. Gregory, D. Levin, A four-point interpolatory subdivision scheme for curve design, Comput. Aided Geom. Design 4 (1987) 257–268. [11] N. Dyn, J.A. Gregory, D. Levin, A butterfly subdivision scheme for surface interpolation with tension control, ACM Trans. Graph. 9 (1990) 160–169. [12] N. Dyn, D. Levin, Analysis of asymptotically equivalent binary subdivision schemes, J. of Math. Anal. Appl. 193 (1995) 594–621. [13] N. Dyn, D. Levin, A. Luzzatto, Exponential reproducing subdivision scheme, Found. Comp. Math. 3 (2003) 187–206. [14] N. Dyn, D. Levin, J. Yoon, Analysis of univariate non-stationary subdivision schemes with application to Gaussian-based interpolatory schemes, SIAM J. Math. Anal. 39 (2007) 470–488. [15] G.H. Golub, C.F. Van Loan, Matrix Computations, John Hopkins University Press, Baltimore, 1996. [16] M.J. Jena, P. Shunmugaraj, P.J. Das, A subdivision algorithm for trigonometric spline curves, Comput. Aided Geom. Design 19 (2002) 71–88. [17] M.J. Jena, P. Shunmugaraj, P.J. Das, A non-stationary subdivision scheme for generalizing trigonometric surfaces to arbitrary meshes, Comput. Aided Geom. Design 20 (2003) 61–77. [18] M.K. Jena, P. Shunmugaraj, P.C. Das, A non-stationary subdivision scheme for curve interpolation, Anziam J. 44 (E) (2003) 216–235. [19] G. Morin, J. Warren, H. Weimer, A subdivision scheme for surfaces of revolution, Comput. Aided Geom. Design 18 (2001) 483–502. [20] L. Romani, From approximating subdivision schemes for exponential splines to high-performance interpolating algorithms, J. of Comp. and Appl. Math. 224 (2009) 383–396. [21] J. Warren, H. Weimer, Subdivision Methods for Geometric Design, Academic Press, New York, 2002.
Contents lists available at ScienceDirect
Applied Numerical Mathematics www.elsevier.com/locate/apnum
Non-stationary subdivision schemes for surface interpolation based on exponential polynomials ✩ Yeon Ju Lee a , Jungho Yoon b,∗ a b
Department of Mathematical Sciences, KAIST, 373-1, Yuseong-gu, Daejeon, 305-701, Republic of Korea Department of Mathematics, Ewha Womans University, Seoul, 120-750, Republic of Korea
a r t i c l e
i n f o
Article history: Received 12 November 2008 Received in revised form 9 September 2009 Accepted 14 October 2009 Available online 23 October 2009 MSC: 41A05 41A25 41A30 65D10 65D17
a b s t r a c t This paper is concerned with non-stationary interpolatory subdivision schemes that can reproduce a large class of (complex) exponential polynomials. It enables our scheme to exactly reproduce the parametric surfaces such as torus and spheres. The subdivision rules are obtained by using the reproducing property of exponential polynomials which constitute a shift-invariant space S. In this study, we are particularly interested in the schemes based on the known butterfly-shaped stencils, proving that these schemes have the same smoothness and approximation order as the classical Butterfly interpolatory scheme. © 2009 IMACS. Published by Elsevier B.V. All rights reserved.
Keywords: Non-stationary subdivision Exponential polynomial Interpolation Asymptotical equivalence Smoothness Approximation order
1. Introduction Subdivision scheme is an efficient method for generating smooth curves and surfaces from a finite set of control points in Rd , d 1. A non-stationary subdivision scheme consists of recursive refinements of an initial sparse sequence with the use of rules that may vary from level to level but are the same everywhere on the same level. Starting with the control points P 0 = { pn0 : n ∈ Zd }, we define new sets of points P k = { pnk : n ∈ Zd } generated by the relation
P k = S a[k−1] · · · S a[0] P 0 . Non-stationary subdivision schemes are useful because they can provide design flexibility and the masks can be adapted to the geometrical configuration of the given data. Exponential (or trigonometric) polynomials and B-splines were used as a key ingredient for the construction of non-stationary subdivision schemes in [1,12,13,16,17,19–21]. Gaussian-based interpo✩ Yeon Ju Lee was supported by the Korea Research Foundation Grant funded by the Korean Government (MOEHRD) (KRF-2007-357-C00007). Jungho Yoon was supported by Priority Research Centers Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education, Science and Technology (2009-0093827). Corresponding author. E-mail addresses: [email protected] (Y.J. Lee), [email protected] (J. Yoon).
*
0168-9274/$30.00 © 2009 IMACS. Published by Elsevier B.V. All rights reserved. doi:10.1016/j.apnum.2009.10.005
Y.J. Lee, J. Yoon / Applied Numerical Mathematics 60 (2010) 130–141
131
Fig. 1. Exact reconstruction of a sphere and a torus.
latory schemes were studied in [14]. A general treatment of stationary schemes can be found in the selected references in [3,4,6,7,9,8]. In this study, instead of thinking about generating new vertices in R3 by taking local averages of vertices in the control polyhedron, we are thinking of scalar values which are assigned to the vertices of a triangulation in R2 . Therefore, starting with a given initial data f 0 = { f n0 ∈ R: n ∈ Z2 }, we consider f k = { f nk ∈ R: n ∈ Z2 } generated by the relation
f kj +1 =
[k]
a j −2n f nk ,
k ∈ Z+ , j ∈ Z2 .
(1)
n∈Z2
[k]
The set of coefficients a[k] := {an } is termed the mask of the rule at level k. We denote this rule by S a[k] and the corresponding non-stationary scheme by { S a[k] }. It is common to assume that for each level k, only a finite number of coefficients [k]
an are non-zero so that changes in a control point affect only its local neighborhood. This property clearly facilitates the practical implementation. A subdivision scheme is called stationary when the mask is independent of the levels; then we use the notation a := {an }. We denote the rule at each level by S a and the corresponding stationary scheme by { S a }. Also, a subdivision scheme is said to be uniformly convergent if for the given initial data δ = { f n0 = δn,0 : n ∈ Z2 }, there exists a function φ0 ∈ C 0 (R2 ), φ0 ≡ 0, satisfying
lim sup f nk − φ0 2−k n = 0,
k→∞ n∈Z2
(2)
where δn,0 denotes the Kronecker delta. In particular, φ0 is called the basic limit function of { S a[k] }, formally defined by
φ0 = lim S a[k] · · · S a[0] δ.
(3)
k→∞
Interpolatory subdivision schemes refine data by inserting values corresponding to intermediate points, using linear combinations of neighboring points. They are particularly attractive since they can control the shape of the resulting subdivision curves/surfaces in a more intuitive way than approximating schemes. The general form of their refinement rules (based on a planar parametric domain) is given as follows: for any j ∈ Z2 and k ∈ Z+ , 1 k f 2k+ j = fj, 1 f 2k+ = j +γ
[k]
aγ +2n f kj−n ,
γ ∈ E 2 \ 0,
n∈Z2
where E 2 is the extreme points of [0, 1]2 , i.e., E 2 = {(0, 0), (0, 1), (1, 0), (1, 1)}. A typical example of such stationary scheme is the four-point scheme by Dyn, Gregory, and Levin [10], and it is extended to Deslauriers–Dubuc schemes [8], where finer level points are determined by local polynomial interpolation of the coarse level points. The well-known generalization of the stationary 4-point scheme to the triangular mesh is the butterfly scheme [11]. It provides C 1 surfaces under regular topology (i.e., if every vertex in the triangular mesh has valence 6) and reproduces polynomials up to degree 3 by using a minimal support. Non-stationary 4-point interpolatory schemes were developed by Beccari et al. [1,2]. Also, in the paper [18], Jena et al. studied non-stationary 4-point interpolatory scheme reproducing elements in the linear space spanned by {1, x, cos(ax), sin(ax)}. All these schemes were developed for generating smooth curves but the corresponding non-stationary butterfly scheme and their analytical properties have not been investigated. Construction of non-stationary subdivision rules based on butterfly stencils are rather simple (see [5] or Section 2). However, fundamental questions on non-stationary subdivision schemes such as convergence, smoothness and approximation orders are yet to be answered. For this reason, in this study, we are concerned with a family of non-stationary interpolatory subdivision schemes based on the known butterfly-shaped stencils. Each scheme in this family is exact on a certain shift-invariant space (say, S) consisting of (complex) exponential polynomials. It enables our scheme to exactly reproduce the parametric surfaces such as torus and sphere (see Fig. 1). Our specific contribution is as follows: Under a certain condition of S, we will see that this new scheme provides the same smoothness C 1 as the original butterfly scheme. Moreover, we show that the non-stationary scheme provides the same approximation order four as the original butterfly scheme with the additional advantage that there are flexibilities in the choice of S and frequency factors.
132
Y.J. Lee, J. Yoon / Applied Numerical Mathematics 60 (2010) 130–141
Fig. 2. Stencils of the butterfly scheme. Larger dots indicate the insertion point.
In the following we use the two-dimensional index notation. First, let
Z2+ := (α1 , α2 ) ∈ Z2 : αi 0, i = 1, 2 .
The set of all complex numbers is denoted by C. For α = (α1 , α2 ) ∈ Z2+ , we set α ! = α1 !α2 ! and |α |1 := α1 + α2 . For α α x = (x1 , x2 ) ∈ R2 , xα = x1 1 x2 2 . The space of algebraic polynomials of degree n is denoted by Πn . For a matrix M, M∞ indicates its ∞-norm. Also, for a given smooth function f , denote by T f ,τ the Taylor polynomial f of degree 3 around τ and let R f ,τ be its remainder, i.e.,
T f ,τ (x) =
(x − τ )ν f (ν ) (τ )/ν !,
R f ,τ (x) =
(x − τ )ν f (ν ) (ξ )/ν !,
(4)
|ν |1 =4
|ν |1 3
for some ξ between x and
τ.
2. Construction of non-stationary schemes The non-stationary subdivision schemes proposed in this paper are constructed in a way of reproducing the exponential (or trigonometric) polynomials of the form
ϕ (x) = xα eβ·x , x ∈ R2 , α ∈ Z2+ , β ∈ C2 .
(5)
In practice, xα ’s are chosen to be low degree polynomials. Let
S = span{ϕ | = 1, . . . , 8}
(6)
be a shift-invariant space with linearly independent functions ϕ , = 1, . . . , 8, of the form as in (5). Refining the triangular mesh data based on the butterfly-shaped stencils involves three groups of vertices to evaluate three types of new vertices as shown in Fig. 2. Here and in the sequel, X j with j = 1, 2, 3 indicate the three types of stencils and p j /2 the insertion points corresponding to X j at level 0. Letting
p 1 = (1, 0),
p 2 = (0, 1),
p 3 = (1, 1),
(7)
for each level k = 0, 1, . . . , and stencil X j , the non-stationary subdivision rule is constructed by solving the linear system
ϕ p j 2−k−1 =
n∈X j
[k]
a p −2n ϕ n2−k , j
= 1, . . . , 8.
(8)
This linear system can be written in the matrix form
a[k] = B[k]
−1 [k]
b
(9)
,
where
[k]
a[k] = a p −2n : n ∈ X j , j B[k] ( , n) =
b[k] =
ϕ p j 2−k−1 : = 1, . . . , 8 ,
ϕ n2−k : n ∈ X j , = 1, . . . , 8 .
For the unique solution of the linear system (9), it is required that dim(S|2−k X j ) = dim S for each stencil X j . We will see later that the mask a[k] converges to the mask of the classical butterfly scheme as k → ∞. Example 2.1 (Sphere and torus). Let us define S as follows: For x = (u , v ),
S := span{sin u , cos u , sin v , cos v , sin u sin v , cos u cos v , sin u cos v , cos u sin v }.
Y.J. Lee, J. Yoon / Applied Numerical Mathematics 60 (2010) 130–141
133
Fig. 3. Reconstruction of a torus. The images from left to right show the subdivision results using our scheme and using the classical butterfly scheme.
Fig. 4. Surface metamorphosis. The images from left to right show the subdivision of a sphere while modulating the frequency of the trigonometric basis functions.
Then the three types of refinement rule are obtained by solving the linear system (9) using the butterfly stencils. Recursive application of this non-stationary subdivision rule exactly generates a sphere whose parametric equations are of the form
f (u , v ) = r sin u cos v ,
g (u , v ) = r sin u sin v ,
h(u , v ) = r cos u ,
(10)
where 0 v 2π and 0 u π ; see Fig. 1. It seems that there are two vertices (at poles) which look like extraordinary points with valences not equal to 6 but this is not the case. Because these functions in (10) are 2π -periodic, we use larger extended domain which are mapped into the same position. Eventually, no extraordinary stencil is necessary. The flowerlike shape obtained in Fig. 4 is due to this extended parameterization domain whose border is mapped onto R3 with this particular spinning covering (see Example 2.2). Similarly, the non-stationary scheme reconstructs a torus whose parametric equations are of the form
f (u , v ) = (1 + cos u ) cos v ,
g (u , v ) = (1 + cos u ) sin v ,
h(u , v ) = sin u ,
where 0 u , v 2π ; see Fig. 1. Also, a subdivision result by using the classical butterfly scheme is given in Fig. 3. Example 2.2. As shown in the images of Fig. 4, a flower-shaped model can be created from an initial, spherical mesh by modulating the frequencies of the trigonometric polynomials (see Example 2.1). In this example, our subdivision rules are obtained in a way of reproducing functions in the space
S := span{sin qu , cos qu , sin qv , cos qv , sin qu sin qv , cos qu cos qv , sin qu cos qv , cos qu sin qv } with q = 1.8. Trigonometric functions with higher frequencies have more oscillations and higher curvatures. Because our scheme reproduces these functions accurately, the limit surface becomes more oscillatory in order to capture these features. Since our scheme is implemented by local interpolation using exponential polynomials, we may perform the metamorphosis of a given initial coarse mesh into an interesting yet controllable shape by adopting different frequencies. The idea of surface modulation is to model a predictable shape using only a simple shape like a sphere. However, in order for such a metamorphosis to be useful, a deep study is necessary in the future. 3. Asymptotic equivalence Typical tools used for the analysis of the stationary schemes are Fourier transform, eigen analysis and Laurent polynomials [3,9]. However, these techniques are not applicable for the non-stationary case. Thus, for the analysis of non-stationary schemes, we adopt the notion of asymptotical equivalence [12]: a (non-stationary) subdivision scheme { S a[k] } is asymptotically equivalent to { S a } if
k∈Z+
where
S a[k] − S a ∞ < ∞,
(11)
134
Y.J. Lee, J. Yoon / Applied Numerical Mathematics 60 (2010) 130–141
S a[k] ∞ = max
[k] a
α −2n : α ∈ E 2
n∈Z2
with E 2 the extreme points of [0, 1]2 . 4. Analysis of convergence The proposed non-stationary butterfly scheme reproduces functions in the space S with dim S = 8. Hence, we choose an index set
T = ν ∈ Z2+ : |ν |1 3, ν = (1, 2), (2, 1) .
(12)
Note that #T = 8. Here, we assume that the set T is sorted by the polynomial ordering of a sequence of multi-indices, that is, an index α := (α1 , α2 ) comes before β := (β1 , β2 ) if |α |1 < |β|1 , or if |α |1 = |β|1 and α1 < β1 . Moreover, for notational convenience, we set
{ϕβ : β ∈ T} := {ϕ1 , . . . , ϕ8 }.
ϕβ satisfies the condition ϕβ(β) (0) = 0 with β ∈ T. Then, let us define the function
It is immediate that
W [ϕ1 , . . . , ϕ8 ](x) := det W(x) with the matrix W(x) given by
W(x) =
(ν )
ϕ j (x): j = 1, . . . , 8, ν ∈ T .
(13)
This function W [ϕ1 , . . . , ϕ8 ](x) is, in fact, the Wronskian of assume that
ϕ1 , . . . , ϕ8 on two-dimensional case. Throughout this paper, we
W [ϕ1 , . . . , ϕ8 ](0) = 0. We cite a basic result about the asymptotically equivalent subdivision schemes. Theorem 4.1. (See [12].) Let { S a } be a stationary subdivision scheme, and let { S a[k] } be a non-stationary subdivision which is asymp[k]
totically equivalent to { S a }. Assume that supp({an }), supp({an }) ⊂ [0, N ]2 , k ∈ Z+ , N < ∞. Then { S a[k] } is C 0 if { S a } is C 0 . Further, if
S a[k] − S a ∞ c2−k ,
k ∈ Z+ ,
then the basic limit function φ0 of { S a[k] } in (3) is Hölder continuous of exponent s > 0, i.e., φ0 ∈ H s . We are now ready to show that the non-stationary subdivision scheme { S a[k] } is asymptotically equivalent to the classical butterfly scheme { S a }, which implies that { S a[k] } is a convergent scheme. Recall that the mask {a p j −2n : n ∈ X j } of the butterfly scheme reproduce polynomials in Π3 in the sense that
a p j −2n p n2−k = p 2−k−1 p j ,
∀ p ∈ Π3 ,
(14)
n∈X j
with p j ( j = 1, 2, 3) in (7). [k]
Theorem 4.2. Let {an } be the mask of the butterfly scheme { S a }, and let {an } be the mask of the non-stationary scheme { S a[k] } reproducing functions in the space S = span{ϕβ : β ∈ T}. Then, there exists a constant c > 0 such that
[k]
maxan − an c2−k , n∈Z2
k K ∈ Z+ . [k]
[k]
Proof. Due to the fact that a2n = a2n = δn,0 with n, 0 ∈ Z2 , it suffices to estimate only the difference |a p j −2n − a p −2n | j with p j ( j = 1, 2, 3) in (7). For each β ∈ T, define a function Φβ by
Φβ (x) :=
c β (γ )ϕγ (x)
(15)
γ ∈T
so that the coefficients c β (γ ), (γ )
(γ )
γ ∈ T, are obtained by solving the linear system
Φβ (0) = ϕβ (0)(1 − δβ,γ ).
(16)
Y.J. Lee, J. Yoon / Applied Numerical Mathematics 60 (2010) 130–141
135
The existence of the unique solution of this linear system is clear from the condition det W(0) = 0. Further, since the mask {an[k] } reproduces ϕβ with β ∈ T (see (8)), it also reproduces Φβ such that we can get
ϕβ 2−k−1 p j − Φβ 2−k−1 p j =
n∈X j
[k]
a p −2n j
ϕβ n2−k − Φβ n2−k .
(17)
Now, for a smooth function g, put T g = T g ,0 , that is, the Taylor polynomial of g of degree 3 around 0 (see (4)). The condition (16) implies that
T ϕβ (x) − T Φβ (x) =
xβ
β!
ϕβ(β) (0) +
xν (ν )
ν! ν =(1,2),(2,1)
ϕβ (0) − Φβ(ν ) (0) .
Then, for any β ∈ T, define Pβ (x) by
Pβ (x) :=
(β)
xβ ϕβ (0)/β! (β) xβ ϕβ (0)/β! +
(2−k ·) = 2−|β|1 k P
Note that Pβ Eq. (17) becomes
if |β|1 2,
xν
(ν ) (ν ) ν =(1,2),(2,1) ν ! (ϕβ (0) − Φβ (0))
β ∈ T. Then, there exists a (smooth) function Rβ,k such that after some calculation,
β for any
Pβ ( p j /2) + Rβ,k ( p j /2) =
otherwise.
[k]
n∈X j
a p −2n Pβ (n) + Rβ,k (n) , j
= 2−k ,
(18)
with β ∈ T. Here, it is obvious that Rβ,k ( p j /2) and Rβ,k (n), n ∈ X j , are bounded by a constant c > 0 independent of k ∈ N and n ∈ X j . On the other hand, using the fact that the mask {a p j −2n : n ∈ X j } of the classical butterfly scheme reproduces Π3 in the sense of (14), we have
a p j −2n Pβ (n) = Pβ ( p j /2).
n∈X j
This linear system can be written in the matrix form
T·a=b
(19)
with
T = Pβ (n): n ∈ X j , β ∈ T , a = (a p j −2n : n ∈ X j ),
b = Pβ ( p j /2): β ∈ T .
(20) [k]
[k]
Then, in view of (18), (19) and (20), one can build a linear system which uniquely defines the mask a j := (a p −2n : n ∈ X j ), j in the form
(T + Fk )a[jk] = b + Gk (β)
with Fk ∞ , Gk ∞ c 1 < ∞. Since ϕβ (0) = 0, T is invertible and hence, T + Fk is also invertible for a sufficient small > 0. Note that the O ( ) perturbation of the non-singular matrix T causes the O ( ) perturbation of its inverse [15]. Thus, we can write [k]
a j = T−1 + F˜ k (b + Gk )
= T−1 b + T−1 Gk + F˜ k (b + Gk ) ,
with F˜ k ∞ c 2 < ∞, where T−1 b = a. Since T−1 Gk + F˜ k (b + Gk )∞ c 3 < ∞, we conclude that there is a constant c > 0 [k] such that a j − a∞ c2−k for any k K . 2 5. Smoothness 5.1. Sufficient condition for smoothness An analysis of the smoothness of univariate non-stationary subdivision schemes is discussed in [14]. In this paper, the result of [14] is extended to the case of non-stationary butterfly scheme. Specifically, our goal is to prove that the non-stationary scheme reproducing the space S in (6) has the same smoothness as the original butterfly scheme. To this end, we first show that a factor (1 + w k zλ ), λ ∈ Z2 , in the Laurent polynomial a[k] ( z) plays a role of smoothing factor if |1 − w k | c2−k . We then infer results on the smoothness of the suggested interpolatory non-stationary scheme from the smoothness of a stationary scheme, which is asymptotically equivalent to it.
136
Y.J. Lee, J. Yoon / Applied Numerical Mathematics 60 (2010) 130–141
Lemma 5.1. Let the Laurent polynomial of the non-stationary scheme { S a[k] } be of the form
a[k] ( z) =
1 2
1 + w k zλ b[k] ( z),
k > K ∈ Z+ ,
where λ ∈ Z2 and { S b[k] } is of compact support and C 0 . Suppose that
|1 − w k | c2−k ,
k K ∈ Z+ .
(21)
Then the basic limit function φ0 in (3) satisfies φ0 , ∂λ φ0 ∈ C where δ := {δn,0 : n ∈ Z }. 0
2
Remark. For the butterfly scheme in regular triangulation, we need to consider only the three directions: λ = (1, 0), (0, 1), and (1, 1). Proof. Due to Theorem 4.1, we can find that { S a[k] } is C 0 with the basic limit function φ0 of the form
φ0 =
φb (· − λt )h(t ) dt , R
where φb and h are the basic limit functions of { S b[k] } and { S 1+ w k z } respectively, and where h is bounded and supp{h} = [0, 1). Here, by virtue of Theorem 4.1, Theorem 4.2 and (21), φb is Hölder continuous of exponent s > 0. Then, it follows that φ0 ∈ C 0 . In order to prove ∂λ φ0 ∈ C 0 , it is sufficient to show that φ0 = g ∗ h ∈ C 1 (R) for any bounded and Hölder continuous function g of exponent s > 0, i.e., g ∈ H s . In fact, the property φ0 = g ∗ h ∈ C 1 (R) can be proved directly by using Lemma 2.5 in [14], 2 In order to guarantee the condition (21), we require a[k] to satisfy the following stronger condition: Condition A. A non-stationary subdivision scheme { S a[k] } satisfies Condition A if the corresponding Laurent polynomial a[k] ( z) is of the form (26) and for = 1, 2,
m ∂ [k] −(4−m)k a (θ ) , c2 ∂ zm
m = 0, . . . , 3 ,
with θ1 = (−1, 1) and θ2 = (1, −1). In what follows, we will prove that the Laurent polynomial a[k] ( z) associated with the scheme { S a[k] } satisfies the Condition A. It in fact implies that the corresponding scheme { S a[k] } has the same smoothness C 1 as the classical butterfly scheme. 5.2. Laurent polynomials of butterfly schemes To simplify the presentation of a subdivision scheme and its analysis, it is convenient to assign the Laurent polynomial to each subdivision mask. Denoting z = ( z1 , z2 ), the bivariate Laurent polynomial corresponding to the original butterfly scheme can be put in the following factored form [9]:
a( z ) =
an zn = 2−1 1 + z1−1 1 + z2−1 1 + z1−1 z2−1 z1 z2 1 − 16−1 q( z1 , z2 ) ,
(22)
n∈Z2
where
q( z) = q( z1 , z2 ) = 2z1−2 z2−1 + 2z1−1 z2−2 − 4z1−1 z2−1 − 4z1−1 − 4z2−1
+ 2z1−1 z2 + 2z1 z2−1 + 12 − 4z1 − 4z2 − 4z1 z2 + 2z12 z2 + 2z1 z22 . The Laurent polynomial a( z) with z = ( z1 , z2 ) has single zeros along the curves z1 = −1, z2 = −1 and z1 z2 = −1. In particular, we will see that a( z) has zeros of order 4 at the points (−1, 1) and (1, −1). In the following, for convenience, we use the notation
e 1 = (1, 0)
and e 2 = (0, 1).
(23)
Lemma 5.2. Let a( z) in (22) be the Laurent polynomial of the classical butterfly subdivision scheme { S a }. Then, for any 0 m 3,
∂m a(θ ) = 0, ∂ zm
= 1, 2,
with
θ1 = (−1, 1) and θ2 = (1, −1).
(24)
Y.J. Lee, J. Yoon / Applied Numerical Mathematics 60 (2010) 130–141
137
Proof. In this proof, we consider only the case = 1 because the case = 2 can be proved by the same way. For 0 m 3, we can write
∂m cm,r an nre1 θ1n , m a(θ1 ) = ∂ z1 2 0r m
(25)
n∈Z
pj
where e 1 = (1, 0), for some suitable constants cm,r . Seeing that θ1 = (−1) j for j = 1, 2, 3 and a2n = δ0,n with n ∈ Z2 , we obtain from the polynomial reproducing property in (14) that
2−r
an nre1 θ1n = δr ,0 +
3
p j −2n
a p j −2n ( p j /2 − n)re1 θ1
j =1 n∈X j
n∈Z2
= δr ,0 +
3
p
θ1 j
a p j −2n ( p j /2 − n)re1
n∈X j
j =1
= δr ,0 +
3 (−1) j δr ,0 = 0. j =1
Combining it with (25), the lemma is proved.
2
According to Theorem 4.2, the non-stationary scheme { S a[k] } is asymptotically equivalent to the original butterfly scheme [k]
{ S a } with the convergence property |an − an | = O (2−k ). Recalling that the Laurent polynomial a(z), z = (z1 , z2 ), in (22) has zeros along the lines z1 = −1, z2 = −1, and z1 z2 = −1, we see that a[k] (−1, z2 ) = a(−1, z2 ) + O (2−k ) = O (2−k ) for any fixed z2 . Similarly, we can show that a[k] ( z1 , −1) = a[k] ( z1 , − z2 ) = O (2−k ). Thus, it is not difficult to deduce that the Laurent polynomial a[k] ( z) associated with S a[k] is of the form a[k] ( z) :=
[k]
an zn
n∈Z
:= 2−1 (1 + w 1,k z1 )(1 + w 2,k z2 )(1 + w 3,k z1 z2 )c [k] ( z1 , z2 )
(26)
with w j ,k ∈ C and
|1 − w j ,k | = o(1),
j = 1, 2, 3, as k → ∞,
for some suitable Laurent polynomial c [k] ( z). 5.3. Analysis of smoothness Let a( z1 , z2 ) be the Laurent polynomial of the classical butterfly scheme as in (22). It is known that the subdivision scheme corresponding to the Laurent polynomial
b j ( z1 , z2 ) = 2(1 + z j )−1 a( z1 , z2 ),
j = 1, 2,
is C 0 (see [9]). Since the non-stationary scheme { S a[k] } is asymptotically equivalent to the butterfly scheme { S a }, the scheme corresponding to [k]
b j ( z1 , z2 ) = 2(1 + w j ,k z j )−1 a[k] ( z1 , z2 ),
j = 1, 2,
(27)
with a[k] ( z1 , z2 ) in (26) is also C 0 if |1 − w j ,k | c2−k for any j = 1, 2. Thus, according to Lemma 5.1, the C 1 -smoothness of the non-stationary scheme { S a[k] } can be verified by showing that |1 − w j ,k | c2−k . Therefore, for this proof, we show that Condition A on { S a[k] } implies the condition (21) for all the factors in the representation (26). Lemma 5.3. Suppose that the Condition A holds for the Laurent polynomial a[k] ( z) in (26) which is associated with the non-stationary subdivision scheme { S a[k] }. Then, for each j = 1, 2,
|1 − w j ,k | c2−k ,
k K ∈ Z+ .
Proof. For simplicity, let
ak,1 ( z1 ) := a[k] ( z1 , 1)
and ak,2 ( z2 ) := a[k] (1, z2 ),
(28)
138
Y.J. Lee, J. Yoon / Applied Numerical Mathematics 60 (2010) 130–141
i.e., ak, j , j = 1, 2, are univariate Laurent polynomials. Then, since { S a[k] } is asymptotically equivalent to { S a }, we deduce from Lemma 5.2 that ak, j ( z j ) can be written as follows:
ak,1 ( z1 ) = ck,1 ( z1 )
4
(1 + rk,n z1 ),
n =1
ak,2 ( z2 ) = ck,2 ( z2 )
4
(1 + sk,n z2 ),
(29)
n =1
where ck, j (−1) = c j + o(1) for some constants c j = 0 ( j = 1, 2). Now, we will show that
|1 − rk,n |, |1 − sk,n | c2−k ,
k K ∈ Z+ ,
for any n = 1, . . . , 4. This is indeed enough for the proof of |1 − w j ,k | c2−k for k K ∈ Z+ . For this purpose, observe that
m ∂ m [k] d a ) = a (θ j k, j (−1) ∂ zmj dzm j
with θ j in (24). Thus, the Condition A is reduced to the univariate case as follows:
m d c2−(4−m)k . (− a 1 ) dzm k, j j
It enables us to apply the known result in [14, Lemma 2.6] by which we can conclude that |1 − rk,n | c2−k . Also, the bound |1 − sk,n | c2−k can be proved by applying exactly the same technique. Therefore, the relation (28) is proved. 2 We now prove that the Laurent polynomial a[k] ( z) of the scheme { S a[k] } satisfies the Condition A. It in fact implies that the corresponding scheme { S a[k] } has the smoothness C 1 .
[k]
Theorem 5.4. Let a[k] ( z) = n∈Z2 an zn be the Laurent polynomial associated with the non-stationary interpolatory scheme { S a[k] }. Then, for any 0 m 3, we have
m ∂ [k] −k(4−m) a (θ ) , c2 ∂ zm
= 1, 2,
(30)
with θ1 = (−1, 1) and θ2 = (1, −1). Proof. In this proof, we consider only the case = 1 because the case = 2 can be done by exactly the same way. Note that
[k] ∂ m [k] μm,r an nre1 θ1n , m a (θ1 ) = ∂ z1 2 0r m
n∈Z
where e 1 = (1, 0), for some constants
M r ,k :=
μm,r . Hence, in order to prove (30), it is sufficient to show that for any 0 r m, [k] an nre1 θ1n = O 2−k(4−r ) , as k → ∞.
n∈Z2
[k]
[k]
Since {an } is the mask of an interpolatory scheme, a2n = δn,0 with n, 0 ∈ Z2 . Hence,
[k]
a2n (2n)re1 = δr ,0 .
n∈Z2
It induces that
2−r (k+1) M r ,k = δr ,0 +
3 [k] re1 (−1) j a p −2n ( p j /2 − n)2−k j
(31)
n∈X j
j =1
where X j are the butterfly stencils of the type j = 1, 2, 3 (see Fig. 1) and p j /2 is the corresponding insertion point (see (7)). In order to estimate the term in (31), for the given k ∈ N and j = 1, 2, 3, define [k, j ]
Φre1 (x) := Φre1 (x) :=
8
=1
[k, j ]
cre1 ( )ϕ (x)
(32)
Y.J. Lee, J. Yoon / Applied Numerical Mathematics 60 (2010) 130–141
[k, j ]
[k, j ]
(ν ) Φre p j 2−k−1 = δre1 ,ν (−1)ν ν !, 1
ν ∈ T.
139
so that the coefficient vector ck := cre1 := (cre1 ( ): = 1, . . . , 8) T is obtained by solving the linear system
(33)
These equations can be written in the matrix form
Wk · ck = r with
Wk := W p j 2−k−1 :=
(ν )
ϕ
r := rre1 := δre1 ,ν (−1)|ν |1 ν !:
p j 2−k−1 :
ν ∈ T, = 1, . . . , 8 ,
ν ∈T .
Since W [ϕ1 , . . . , ϕ8 ](x) = W(x) is continuous and W [ϕ1 , . . . , ϕ8 ](0) is non-zero, there exists K ∈ N so that for any k K , the matrix Wk is non-singular and Wk−1 ∞ c for some constant c > 0. It implies that for any compact set M in R2 and α ∈ Z2+ ,
(α ) Φre 1
L∞
6 (α ) −1 ϕ W r∞ k (M ) ∞ L
=1
∞ (M )
cα,M ,
k K,
(34)
[k]
for some constant c α , M > 0. On the other hand, since the mask {an } reproduces duces Φre1 . Hence, by (33), we get
ϕβ with β ∈ T (see (8)), it also repro-
[k] δr ,0 = Φre1 2−k−1 p j = a p −2n Φre1 n2−k . j n∈X j
This together with (31) leads to
2−r (k+1) M r ,k =
3 [k] re1 . (−1) j −1 a p −2n Φre1 n2−k − ( p j /2 − n)2−k j j =1
(35)
n∈X j
τ = 2−k−1 p j , write Φre1 n2−k = T Φre1 ,τ n2−k + R Φre1 ,τ n2−k
Now, letting
with T Φre1 ,τ and R Φre1 ,τ in (4). Then, invoking (33), we see that
T Φre1 ,τ n2−k − ( p j /2 − n)2−k
re1
=
ν (ν ) −k−1 (n − p j /2)2−k Φre 2 p j /ν !. 1
ν =(1,2),(2,1)
[k]
Due to Lemma 4.2, an = an + O (2−k ). Using the polynomial reproducing property of the mask {an } of the butterfly scheme (14), it is immediate that for each ν = (1, 2), (2, 1),
[k] −k ν −k ν −4k ( ( a n − p / 2 ) 2 a n − p / 2 ) 2 = c2−4k p − 2n j j j p j −2n + c2 n∈X j
(36)
n∈X j
for a constant c > 0. Next, by applying (34), it is obvious that
−k [k] c2−4k , a R n2 p j −2n Φre1 ,τ
k K.
n∈X j
This together with (35) and (36), we obtain the required result.
2
We are now ready to provide the main theorem of this section. Theorem 5.5. If { S a[k] } is the non-stationary interpolatory subdivision scheme reproducing the shift-invariant space S in (6), then { S a[k] } is C 1 , i.e., it has the same smoothness as the stationary butterfly subdivision scheme { S a }. Proof. From Lemmas 5.1, 5.3, and Theorem 5.4, the proof is immediate.
2
140
Y.J. Lee, J. Yoon / Applied Numerical Mathematics 60 (2010) 130–141
6. Approximation order An important issue in the implementation of subdivision algorithm is how to actually attain the original function as close as possible if the given initial data f 0 is sampled from an underlying function. A high quality reconstruction scheme should guarantee that the approximation error decreases when the quality of the sample increases. This approximation power is usually quantified by the approximation order. Suppose that f 0 is of the form f 0 := { f n0 = f (2−k0 n): n ∈ Z2 } for some positive integer k0 with a smooth function f . Then our concern is to find the largest exponent m > 0 such that
∞ f − f
L∞ (M )
C 2−mk0
with a constant C > 0 independent of k0 , where M is a compact set in R2 . The exponent m is called the approximation order of the subdivision scheme. It is well known that the original butterfly scheme reproduces cubic polynomials, and hence, it provides the approximation order of four. The goal of this section is to show that the non-stationary scheme { S a[k] } provides the same approximation order four as the original butterfly scheme. It is of interest to see that the butterfly-shaped stencils consist of just 8 vertices. Usually, it is expected to have 10 vertices to obtain the approximation order 4. Since the scheme { S a[k] } is uniformly convergent, its limit function can be written as
f ∞ = lim S a[k0 + ] · · · S a[k0 ] f 0 = →∞
φk0 2k0 · −n f 2−k0 n
(37)
n∈Z2
where φk0 is the basic limit function defined as in (3). Indeed, since { S a[k] } reproduces functions in S, it is clear that
ϕ =
ϕ 2−k0 n φk0 2k0 · −n , ∀ϕ ∈ S.
n∈Z2
This observation leads to the proof of the approximation order of the subdivision scheme { S a[k] }. In this paper, we are particularly interested in approximating functions f in the Sobolev space
s W∞ ( M ) = g : R2 → R:
g (α ) |α |1 s
L∞
< ∞ , (M )
s ∈ Z+ .
4 Theorem 6.1. Let M be a compact set in R2 . Assume that f ∈ W ∞ ( M ) and f 0 := { f n0 = f (2−k0 n): n ∈ Z2 }. Then the non-stationary 4 interpolatory scheme { S a[k] } has the approximation order four with respect to functions in W ∞ ( M ).
Proof. Let x ∈ M be a fixed point and define a function
Φ :=
8
μ ϕ (· − x)
=1
so that the coefficient vector {μ : = 1, . . . , 8} is obtained by solving the linear system
Φ (ν ) (x) = f (ν ) (x),
ν ∈T
with T in (12). The existence of the unique solution of this linear system is clear because det W(0) = 0. Since S is shiftinvariant and φ ∈ S, ϕ (· − x) belongs also to S. Also, since Φ is a linear combination of φ , = 1, . . . , 8, and the scheme { S a[k] } is exact for such functions, we easily get the identity
Φ=
φk0 2k0 · −n Φ n2−k0 .
n∈Z
Then, using the property f (x) = Φ(x), it follows that
f (x) − f ∞ (x) = Φ(x) −
=
φk0 2k0 x − n f n2−k0
n∈Z2
φk0 2k0 x − n Φ n2−k0 − f n2−k0 .
n∈Z2
Recall from (4) that T g ,x indicates the Taylor polynomial of a smooth function g of degree 3 around x and R g ,x is its remainder. Due to the fact Φ (ν ) (x) = f (ν ) (x) with ν ∈ T, we obtain
T Φ,x n2−k0 − T f ,x n2−k0 =
ν =(1,2),(2,1)
It follows that
(n2−k0 − x)ν (Φ − f )(ν ) (x). ν!
Y.J. Lee, J. Yoon / Applied Numerical Mathematics 60 (2010) 130–141
f (x) − f ∞ (x) =
141
(n2−k0 − x)ν (Φ − f )(ν ) (x) ν ! ν =(1,2),(2,1) n∈Z2 k 0 + φk0 2 x − n R Φ,x n2−k0 − R f ,x n2−k0 . φk0 2k0 x − n
(38)
n∈Z2
Letting φ be the basic limit function of the classical butterfly subdivision scheme S a , i.e., φ = S a∞ δ with δ = {δn,0 : n ∈ Z2 }, we find from [12, Lemma 15] and Theorem 4.2 that φk0 − φ∞ c2−k0 . Moreover, applying the polynomial reproducing property of φ in (14), we obtain
ν φ 2k0 x − n 2k0 x − n = 0,
ν ∈ (1, 2), (2, 1) .
n∈Z2
Note that since φk0 is compactly supported, #{n ∈ Z2 : φk0 (2k0 x − n) = 0} c for any x. It implies that
k −k ν k k ν −3k0 0 0 0 0 φ 2 x − n n2 − x 2 φ 2 x − n 2 x − n = k0 k0 n∈Z2
n∈Z2
ν 2−3k0 φ 2k0 x − n 2k0 x − n + c2−k0 n∈Z2
c2
−4k0
for some constant c > 0. Also, since Φ (ν ) is bounded on any compact set M and f term in (38) has the property O (2−4k0 ). It finishes the proof. 2
4 ∈ W∞ ( M ), it is obvious that the second
References [1] C. Beccari, G. Casciola, L. Romani, A non-stationary uniform tension controlled interpolating 4-point scheme reproducing conics, Comput. Aided Geom. Design 24 (2007) 1–9. [2] C. Beccari, G. Casciola, L. Romani, An interpolating 4-point C 2 ternary non-stationary subdivision scheme with tension control, Comput. Aided Geom. Design 24 (2007) 210–219. [3] A. Cavaretta, W. Dahmen, C.A. Micchelli, Stationary subdivision, Mem. Amer. Math. Soc. 93 (1991) 1–186. [4] G. Chaikin, An algorithm for high speed curve generation, Computer Graphics and Image Processing 3 (1974) 346–349. [5] Y.-J. Choi, Y.-J. Lee, J. Yoon, B.-G. Lee, Y.-J. Kim, A new class of non-stationary interpolatory subdivision schemes based on exponential polynomials, in: Geometric Modeling and Processing – GMP 2006, in: Lecture Notes in Computer Science, vol. 4077, 2006, pp. 563–570 (SpringerLink, PDF-file). [6] E. Cohen, T. Lyche, R. Riesenfeld, Discrete B-spline and subdivision techniques in computer-aided geometric design and computer graphics, Computer Graphics and Image Processing 14 (1980) 87–111. [7] D. Doo, M. Sabin, Behaviour of recursively division surfaces near extraordinary points, Comput.-Aided Design 10 (1978) 356–360. [8] G. Deslauriers, S. Dubuc, Symmetric iterative interpolation, Constr. Approx. 5 (1989) 49–68. [9] N. Dyn, Subdivision schemes in computer-aided geometric design, in: W.A. Light (Ed.), Advances in Numerical Analysis Vol. II: Wavelets, Subdivision Algorithms, Radial Basis Functions, Oxford University Press, 1992, pp. 36–104. [10] N. Dyn, J.A. Gregory, D. Levin, A four-point interpolatory subdivision scheme for curve design, Comput. Aided Geom. Design 4 (1987) 257–268. [11] N. Dyn, J.A. Gregory, D. Levin, A butterfly subdivision scheme for surface interpolation with tension control, ACM Trans. Graph. 9 (1990) 160–169. [12] N. Dyn, D. Levin, Analysis of asymptotically equivalent binary subdivision schemes, J. of Math. Anal. Appl. 193 (1995) 594–621. [13] N. Dyn, D. Levin, A. Luzzatto, Exponential reproducing subdivision scheme, Found. Comp. Math. 3 (2003) 187–206. [14] N. Dyn, D. Levin, J. Yoon, Analysis of univariate non-stationary subdivision schemes with application to Gaussian-based interpolatory schemes, SIAM J. Math. Anal. 39 (2007) 470–488. [15] G.H. Golub, C.F. Van Loan, Matrix Computations, John Hopkins University Press, Baltimore, 1996. [16] M.J. Jena, P. Shunmugaraj, P.J. Das, A subdivision algorithm for trigonometric spline curves, Comput. Aided Geom. Design 19 (2002) 71–88. [17] M.J. Jena, P. Shunmugaraj, P.J. Das, A non-stationary subdivision scheme for generalizing trigonometric surfaces to arbitrary meshes, Comput. Aided Geom. Design 20 (2003) 61–77. [18] M.K. Jena, P. Shunmugaraj, P.C. Das, A non-stationary subdivision scheme for curve interpolation, Anziam J. 44 (E) (2003) 216–235. [19] G. Morin, J. Warren, H. Weimer, A subdivision scheme for surfaces of revolution, Comput. Aided Geom. Design 18 (2001) 483–502. [20] L. Romani, From approximating subdivision schemes for exponential splines to high-performance interpolating algorithms, J. of Comp. and Appl. Math. 224 (2009) 383–396. [21] J. Warren, H. Weimer, Subdivision Methods for Geometric Design, Academic Press, New York, 2002.