- Email: [email protected]
Properties of generators of quasi-interpolation operators of high approximation orders in spaces of polyharmonic splines
Journal of Computational and Applied Mathematics 267 (2014) 96–106
Contents lists available at ScienceDirect
Journal of Computational and Applied Mathematics journal homepage: www.elsevier.com/locate/cam
Properties of generators of quasi-interpolation operators of high approximation orders in spaces of polyharmonic splines Mira Bozzini, Milvia Rossini ∗,1 University of Milano-Bicocca, Italy
article
abstract
info
Article history: Received 8 March 2013 Received in revised form 14 November 2013 Keywords: Polyharmonic splines Quasi-interpolation operators High degree polynomial reproduction Multiresolution analysis Scaling functions Subdivision
We have presented in Bozzini et al. (2011) a procedure in spaces of m-harmonic splines in Rd that starts from a simple generator φ0 and recursively defines generators φ1 , φ2 , . . . , φm−1 with corresponding quasi-interpolation operators reproducing polynomials of degrees 3, 5, . . . , 2m − 1 respectively. In this paper we study the properties of generators φj , and we prove that these new generators are positive definite functions, and are scaling functions whenever φ0 has those properties. Moreover φ0 and φj generate the same multiresolution analysis. We show that it is possible to define a convergent subdivision scheme, and to provide in this way a fast computation of the quasi-interpolant. © 2014 Elsevier B.V. All rights reserved.
1. Introduction Polyharmonic splines have been introduced in [1] as a variational extension to many variables of polynomial splines in one variable. Because of their good properties (see e.g. [2–4]) they are studied by several authors, and widely used in scattered data approximation and interpolation. An important amount of research work has been done also on regular grids which occur in many applications such as image processing and multiresolution analysis. In general it is important to construct quasi-interpolation operators capable of reproducing some important functions of the space, and to provide a fast evaluation of the quasi-interpolant function. The space of m-harmonic splines V2m , 2m > d, is defined as the subspace of S ′ (Rd ) (the space of d-dimensional tempered distributions)
V2m = g ∈ S ′ (Rd )
C 2m−d−1 (Rd ) : ∆m g = 0, on Rd \ Zd ,
(1)
2 2 m m m−1 where ∆ := g ). Recently (see [5]) they i=1 ∂ /∂ xi is the Laplace operator and ∆ is defined iteratively by ∆ g = ∆(∆ m have been generalized by considering the fundamental solutions of more general differential operators j=1 (−∆ + κj2 I ).
d
It is well-known (see e.g. [6]) that V2m contains Π2m−1 (the space of polynomials defined on Rd of degree not exceeding 2m − 1), and that for any n ≤ 2m − 1 it is possible to construct quasi-interpolation (QI) operators reproducing Πn (see e.g. [7,8]). The generator φ ∈ V2m of the QI operator Qφ (x, h, f ) =
f (hl)φ(x/h − l)
l∈Zd
∗
Corresponding author. Tel.: +39 0264485715. E-mail addresses: [email protected] (M. Bozzini), [email protected] (M. Rossini).
1 Sponsored by the PRIN Project Real and Complex Varieties: Geometry, Topology, Harmonic Analysis. http://dx.doi.org/10.1016/j.cam.2014.01.029 0377-0427/© 2014 Elsevier B.V. All rights reserved.
(2)
M. Bozzini, M. Rossini / Journal of Computational and Applied Mathematics 267 (2014) 96–106
97
is required to decay fast enough at infinity, so that Qφ (x, h, f ) is well defined for f growing at infinity not faster than a polynomial of degree n. The QI operator (2) approximates smooth enough functions with L∞ -error of order hn+1 [7]. Important bases of V2m are m-harmonic B-splines, [9,4]. In the following we denote them by φ0 . They decay sufficiently fast, at least as ∥x∥−d−2 when ∥x∥ → ∞, they generate quasi-interpolation operators reproducing linear polynomials and then, for f with bounded derivatives of order 1, 2, 3, the sum Q0 (x, h, f ) =
f (hl)φ0 (x/h − l)
l∈Zd
approximates f in Rd with L∞ -error of order h2 . In addition they are valid scaling function for generating a multiresolution analysis of L2 (Rd ) [4]. To obtain higher approximation orders, the function φ0 should be replaced by a function, with a stronger decay rate as ∥x∥ → ∞. For this reason we have presented in [10] a procedure which starts from a simple generator φ0 , and recursively defines generators φ1 , φ2 , . . . , φm−1 with corresponding QI operators reproducing Π3 , Π5 , . . . , Π2m−1 respectively. The aim of this paper, is to study other important properties of the generators φj defined in [10]. We prove that these new generators are positive definite functions, are scaling functions whenever φ0 has those properties. Moreover φ0 and φj generate the same multiresolution analysis, and we show that it is possible to define a convergent subdivision scheme and to provide a computational refinement which allows us a fast computation of Qj (x, h, f ) =
f (hl)φj (x/h − l).
l∈Zd
The paper is organized as follows. In Section 2 we give some definitions, preliminary results, and we briefly recall the construction of the new generators φj . In Section 3 we give the properties of φj , and in Section 4 we show some numerical experiments. 2. Preliminaries 2.1. Definitions, notations, useful results
d
We use in this paper the multi-index notation. In particular for k ∈ Zd , |k| = i=1 |ki |. The symbol ∗ denotes the convolution operator between functions or sequences and the semi-discrete convolution between a sequence and a function. T := [−π , π]d is the d-dimensional torus, L2 (Rd /2π Zd ) is the set of functions locally square integrable and 2π Zd -periodic. In the sequel we assume 2m > d. We say that a linear transformation A on Rd is an acceptable dilation for Zd if it leaves Zd invariant, i.e. AZd ⊂ Zd , and all the eigenvalues of A satisfy |λi | > 1. These properties imply that |det A| is an integer ≥ 2. Important acceptable dilations for Z d are similarities, i.e matrices of the form A = ρ A0 , where A0 is an orthogonal matrix, and ρ is a real number such that |det A| = |ρ d | is an integer ≥ 2. In the case d = 2 these matrices are all of the form
q −p
p q
or
q p
p , −q
with q, p ∈ Z. They include the classic dyadic matrix 2I and the quincunx case q = p = 1 which gives |det A| = 2. A multiresolution analysis (MRA) of dilation A ∈ Rd×d is a sequence V = {Vl , l ∈ Z} of closed subspaces of L2 (Rd ) satisfying the following properties (i) (ii) (iii) (iv)
V l ⊂ Vl+1 , 2 d l∈Z Vl is dense in L (R ) and l∈Z Vl = {0}, −1 f (·) ∈ Vl ⇔ f (A ·) ∈ Vl−1 , there is a function of V0 , called the scaling function, whose integer translates form a Riesz basis of V0 .
Condition (iv) is equivalent to saying that V0 is the closed linear span of Zd translates of one function φ , which satisfies the condition c1 ≤
ˆ − 2π k)|2 ≤ c2 |φ(ω
k∈Zd
for some positive constants c1 , c2 . The whole multiresolution analysis may be regarded as being generated by φ . Note that the scaling function is not unique (see e.g. [11]). It is useful to remember the following results stated in [11, Section 2.3].
ˆ 0) = 1 and φ( ˆ 2kπ ) = 0 for all k ∈ Zd \ {0}, then φ is a valid scaling Lemma 1. Let φ ∈ L1 (Rd ) L2 (Rd ) be such that φ( function for generating a MRA V associated with (Zd , A), if and only if φ enjoys the following properties: (1) there is a function S (ω) ∈ L2 (Rd /2π Zd ) such that
ˆ ˆ B−1 ω) where B = AT ; φ(ω) = S (B−1 ω)φ( (2) for almost all ω ˆ − 2π k)|2 > 0. |φ(ω k∈Zd
(3)
(4)
98
M. Bozzini, M. Rossini / Journal of Computational and Applied Mathematics 267 (2014) 96–106
ˆ 0) = 1 and φ( ˆ 2kπ ) = 0 for all k ∈ Zd \ {0}, and suppose that φ generates Lemma 2. Let φ ∈ L1 (Rd ) L2 (Rd ) be such that φ( a MRA V associated with (Zd , A). Let γ = {γk }k∈Zd ∈ l1 (Zd ) such that γˆ (ω) ̸= 0, ω ∈ T . Then also the function
ϕ =γ ∗φ generates the MRA V associated with (Zd , A). The Green’s function of ∆m is
vm (x) =
Cm ∥x∥2m−d Cm ∥x∥2m−d log ∥x∥
2m − d ̸∈ 2Z 2m − d ∈ 2Z,
(5)
where the constant Cm is such that the generalized Fourier transform is
vˆ m (ω) = (−1)m ∥ω∥−2m .
(6)
A function φ0 : R → R is called m-harmonic B-spline if its Fourier transform has the form d
φˆ 0 (ω) = e2m (ω)ˆvm (ω),
(7)
where e2m is a real, even, and 2π Z -periodic trigonometric polynomial satisfying d
e2m (ω) = (−1)m ∥ω∥2m + O(∥ω∥2m+2 ),
(8)
for ω in a neighborhood of the origin U (0). For d = 1, φ0 is the symmetric B-spline of odd degree 2m − 1 with integer knots. In this case
e2m (ω) = −4 sin2
ω m 2
.
(9)
For d = 2, important examples of φ0 are the elementary m-harmonic B-spline, having
e2m (ω1 , ω2 ) = −4 sin2
ω1 2
− 4 sin2
ω2 m 2
,
(10)
and the isotropic m-harmonic B-spline having (see e.g. [12,13])
e2m (ω1 , ω2 ) = −
m ω1 ω2 2 1 + 4 sin2 + 4 sin2 + cos ω1 cos ω2 . 3 2 2
(11)
Generally, the trigonometric polynomial e2m , as happens for the elementary and isotropic B-spline, is requested to satisfy e2m (ω) = 0,
ω∈T
ω = 0.
⇔
(12)
This guarantees that φ0 is a valid scaling function for generating a MRA V associated with (Z , A) whenever A is a similarity. d
2.2. Construction of φj , j = 1, . . . , m Let us consider a m-harmonic spline φ0 (7). Starting from φ0 we have constructed in [10] new generators
φ1 , φ2 , . . . , φm−1 of QI operators reproducing polynomials up to degree 3, 5, . . . , 2m − 1 respectively. The procedure defines φj recursively as a linear combination of φj−1 (·) and φj−1 (·/2) with explicitly known coefficients which are independent of m. For j = 1, 2, . . . , m − 1, we define
φj (x) = aj φj−1 (x) + bj φj−1 (x/2),
x ∈ Rd ,
(13)
and choose the coefficients aj , bj so that
φˆ j (ω) = 1 +
∞
[ j]
h2i (ω),
ω ∈ U (0),
(14)
i=j+1
[j]
with h2i a homogeneous function of order 2i. This is possible since
φˆ j (ω) = aj φˆ j−1 (ω) + 2d bj φˆ j−1 (2ω),
(15)
and by our inductive hypothesis
φˆ j−1 (ω) = 1 +
∞ i=j
[j−1]
h2i
(ω),
ω ∈ U (0).
(16)
M. Bozzini, M. Rossini / Journal of Computational and Applied Mathematics 267 (2014) 96–106
99
The coefficients aj , bj in (15) are the solution of the system
aj + 2 d b j = 1 aj + 2d+2j bj = 0,
(17)
that is aj =
22j 22j
,
−1
bj = −
1 2d (22j − 1)
,
j = 1, . . . , m − 1.
(18)
Note that the coefficients aj , bj depend on neither φ0 nor m and that bj = −
aj 2d+2j
.
(19)
In [10], it has been proved that φj has the correct decay at infinity needed for the quasi-interpolation based on it to be well defined for polynomials of degree up to 2j + 1, and that these polynomials are reproduced by this quasi-interpolation. Namely for j = 1, 2, . . . , m − 1
φj (x) ≤
C
∥x ∥
d+2j+2
,
∥x∥ → ∞,
(20)
and φj satisfies the Strang Fix conditions of order 2j + 1, i.e.
φˆ j (0) = 1,
Dα φˆ j (0) = 0,
α
D φˆ j (2kπ ) = 0,
1 ≤ |α| ≤ 2j + 1,
(21)
k ∈ Z \ {0}, |α| ≤ 2j + 1. d
(22)
Consequently for j = 1, . . . , m − 1, the sum
f (l)φj (x − l)
(23)
l∈Zd
is well defined for any f growing at infinity not faster than a polynomial of degree 2j + 1, and if f ∈ Π2j+1 , the sum (23) equals f . Moreover Qj (x, h, f ) =
f (hl)φj (x/h − l)
(24)
l∈Zd
approximates f having bounded derivatives of order 2j + 1, 2j + 2, 2j + 3, with L∞ -error of order h2j+2 . 3. Properties of φj In this section we study some additional properties of the generators φj . It is trivial to verify that φj ∈ C 2m−d−1 (Rd ), and that (see (20)) φj ∈ L1 (Rd ) L2 (Rd ), j = 0, . . . , m − 1. Now we want to investigate if φj , j = 1, . . . , m − 1, preserves the properties of φ0 . In particular we are interested in proving that a strictly positive definite function φ0 , under some assumptions on e2m , generates strictly positive definite functions φj , and that if φ0 is a scaling function, it generates functions which are scaling functions too. In addition, we will prove that φ0 and φj generate the same multiresolution analysis of L2 (Rd ). We start by recalling the definition of positive definite function and its characterization (see e.g. [14]). [0]
Definition 3. A complex valued function Φ is called positive definite on Rd if N N
cl c¯k Φ (xl − xk ) ≥ 0
(25)
l=0 k=0
for any pairwise different points x1 , . . . , xN ∈ Rd , and c = [c1 , . . . , cN ] ∈ CN . The function Φ is called strictly positive definite on Rd if (25) is zero only for c ≡ 0. It is known that (see e.g. [14]). Proposition 4. Let Φ be a continuous function in L1 (Rd ). Φ is strictly positive definite if and only if Φ is bounded and its Fourier transform is non negative and not identically equal to zero. In virtue of this fundamental result, φ0 strictly positive definite implies that
φˆ 0 (ω) = (−1)m e2m (ω)∥ω∥−2m ≥ 0
100
M. Bozzini, M. Rossini / Journal of Computational and Applied Mathematics 267 (2014) 96–106
and not identically equal to zero, i.e. e2m (ω) := (−1)m e2m (ω) ≥ 0. [0]
(26)
By construction φˆ j (ω) has the form j] φˆ j (ω) = e[2m (ω)∥ω∥−2m ,
(27)
where [j]
[j−1]
[j−1]
e2m (ω) = aj e2m (ω) + 2d−2m bj e2m (2ω),
j = 1, . . . , m − 1.
(28)
By using (19), we get
[j]
[j−1]
[j−1]
e2m (ω) = aj e2m (ω) − 2−2m−2j e2m (2ω) ,
j = 1, . . . , m − 1.
(29)
Proposition 5 gives a sufficient condition that ensures that φj is strictly positive definite and, as we shall see in Proposition 6, this condition is always satisfied when we start from the most known m-harmonic B-splines: the symmetric odd degree B-spline with integer knots (d = 1), the elementary and the isotropic B-spline (d = 2). Proposition 5. Let φ0 be strictly positive definite and let
αj := 2−2m−2
j −1
2−2r = 2−2m
1 − 2−2j
.
3
r =0
(30)
If for j = 1, . . . , m − 1 e2m (ω) − αj e2m (2ω) ≥ 0, [0]
[0]
(31)
then φj is strictly positive definite. [j]
Proof. By recursion, we express each e2m , j = 1, . . . , m − 1, in terms of e2m . Let us begin with e2m . By definition
[0]
e2m (ω) = a1 e2m (ω) − 2−2m−2 e2m (2ω) , [1]
[0]
[0]
[1]
a1 > 0.
Trivially, if e2m (ω) ≥ 2−2m−2 e2m (2ω) and e2m (ω) is not identically equal to 2−2m−2 e2m (2ω), also e2m is non negative and not identically equal to zero. Since φ1 ∈ L1 (Rd ), φ1 is strictly positive definite. For j = 2, we get [0]
[0]
[0]
[0]
e2m (ω) = a2 e2m (ω) − 2−2m−4 e2m (2ω) , [2]
[1]
[1]
= a1 a2 e2m (ω) − 2 [0]
−2m−2
1+
[1]
a2 > 0
1 22
e2m (2ω) + 2 [0]
e2m (4ω) .
−2m−2 −2m−4 [0]
2
[0] [0] [0] [2] Since e2m (·) ≥ 0 and not identically zero, it follows that if e2m (ω) ≥ 2−2m−2 (1 + 212 )e2m (2ω) also e2m is non negative and not [j]
[0]
identically equal to zero. Going on with this simple stuff, it is not difficult to see that e2m can be written in terms of e2m as
[j]
e2m (ω) =
j
e2m (ω) − 2 [0]
al
−2m−2
l =1
+
j −2
j −1
e2m (2ω)
−2r [0]
2
(32)
r =0
−4m−4r −6
2
r =0
j
j−(r +2)
ak
k=r +1
e2m (4ω) .
−2k [r ]
2
(33)
k=0
Since e2m (ω) ≥ 0, r = 0, . . . , j − 1, and not identically equal to zero, from assumption (31) we get the proof. [r ]
Proposition 6. When φ0 is the symmetric odd degree B-spline with integer knots (d = 1), the elementary or the isotropic B-spline j] is non negative and vanishes only at ω ∈ 2π Zd . (d = 2), the functions φj , j = 1, . . . , m are strictly positive definite. Namely e[2m Proof. We give the proof for elementary B-spline of order m > 1. Analogous arguments can be carried out for the odd degree B-spline with integer knots and the isotropic m-harmonic B-spline for which the trigonometric polynomials are (9) and (11) respectively. [0] The trigonometric polynomial e2m is (see (10) and (26))
e2m (ω1 , ω2 ) = 4m sin2 [0]
ω1 2
+ sin2
ω2 m 2
.
(34)
M. Bozzini, M. Rossini / Journal of Computational and Applied Mathematics 267 (2014) 96–106
101
To verify condition (31) means to verify if, for j = 1, . . . , m − 1
sin2
ω1 2
+ sin2
ω2 m 2
m − αj sin2 ω1 + sin2 ω2 ≥ 0.
(35)
We can factorize the left hand side as
sin2
ω1 2
+ sin2
ω2
−1 m m−1−k ω1 ω2 k √ √ m α (sin2 ω + sin2 ω ) m α (sin2 ω + sin2 ω ) . sin2 + sin2 j 1 2 j 1 2
−
2
k=0
2
(36)
2
The sum m−1
sin2
k=0
ω1 2
ω2 k √ m
+ sin2
2
m−1−k αj (sin2 ω1 + sin2 ω2 )
is always non negative and it is zero only for (ω1 , ω2 ) = (2lπ , 2r π ) l, r ∈ Z. So we have to study the sign of p(ω1 , ω2 ) := sin2
ω1 2
+ sin2
ω2 2
−
2 √ 2 m α j sin ω1 + sin ω2 ,
(37)
that can be rewritten as p(ω1 , ω2 ) = sin2
= sin2
ω1
ω2 √ αj sin2 ω1 + sin2 − m αj sin2 ω2 2 ω1 ω2 ω2 √ √ 1 − 4 m αj cos2 + sin2 1 − 4 m αj cos2 .
− 2 ω1 2
√ m
2
√
2
2
From definition (30), 4 m αj < 1 for any j = 1, . . . , m − 1, then p(ω1 , ω2 ) is non negative, and it is equal to zero only when (ω1 , ω2 ) = (2lπ , 2r π ) l, r ∈ Z. By recursion we get the proof. It is known that the functions φ0 of Proposition 6 are scaling functions of a MRA V of L2 (Rd ) associated with (Zd , A), A a similarity. In Proposition 7, we prove that also φj , j = 1, . . . , m − 1 is a scaling function, and that φ0 and φj generate the same multiresolution analysis. Proposition 7. When φ0 is the symmetric odd degree B-spline with integer knots (d = 1), the elementary or the isotropic B-splines (d = 2), φj , j = 1, . . . , m − 1, is a scaling function of a MRA V of L2 (Rd ) associated with (Zd , A), A a similarity. In addition φj , j = 1, . . . , m − 1, and φ0 generate the same MRA. [j]
Proof. We observe that whenever A is a similarity, there is a function S2m (ω) ∈ L2 (Rd /2π Zd ) such that [j] φˆ j (Bω) = S2m (ω)φˆ j (ω),
B = AT ,
(38)
where [j]
e (Bω) [j] . S2m (ω) = ρ −2m 2m [ j] e2m (ω) Then, we consider
Φj (ω) =
2 φˆ j (ω − 2π k) ,
j = 1, . . . , m − 1.
(39)
k∈Z2
Φj is well defined and bounded for any ω because of the decay of vˆ m . Moreover φˆ j (0) = 1 and φˆ j (2π k) = 0, k ∈ Zd k ̸= 0, then for Lemma 1 we can conclude that φj , j = 1, . . . , m − 1, generates a MRA. [j]
[j]
In addition, e2m (ω) and e2m (ω) are equal to zero only at ω = 2π k k ∈ Zd , and the ratio between e2m (ω) and e2m (ω) is [0]
[0]
[j]
well-defined for any ω and equal to one at ω = 2π k k ∈ Zd . Then, since φˆ j (ω) = e2m (ω)ˆvm (ω), we can trivially rewrite φˆ j as [j]
φˆ j (ω) =
[j]
e2m (ω) [0] e (ω) e2m (ω)ˆvm (ω) = 2m φˆ 0 (ω), [0] [0] e2m (ω) e2m (ω) [j]
where the function γˆj (ω) := e2m (ω)/e2m (ω) is equal to one for ω = 2kπ , k ∈ Zd , and positive everywhere. Moreover γˆj ∈ L2 (Rd /2π Zd ) and then, for Lemma 2, φj j = 1, . . . , m − 1, and φ0 generate the same MRA V associated with (Zd , A). [j]
[0]
S2m (ω) is the 2π Zd periodic function [j]
S2m (ω) = |ρ|−d
k∈Zd
[j] sk e−i⟨k,ω⟩ .
102
M. Bozzini, M. Rossini / Journal of Computational and Applied Mathematics 267 (2014) 96–106
[j]
The coefficients s[j] = {sk }k∈Zd are the scaling coefficients of φj , i.e
φj (·) =
[j]
sk φj (A · −k),
k∈Zd
and as we show in the following, they can be computed explicitly via convolutions. To this end, we first show that the m-harmonic Lagrangian function Λ2m ∈ V2m (Λ2m (l) = δ0l , l ∈ Zd ) can be expressed in terms of shifts of φj . [j]
Let us consider the vector β [j] = {φj (k)}k∈Zd . From (20), β [j] ∈ l1 (Zd ) and, at least |βk | ≤ C |k|−d−2−2j . The absolutely convergent trigonometric series
βˆ [j] (ω) =
φj (k)e−iω·k ̸= 0,
∀ω ∈ Rd .
k∈Zd
In fact, from an application of the Poisson summation formula
φj (k)e−iω·k =
k∈Zd
φˆ j (ω − 2π k),
k∈Zd
[j]
and since φˆ j is non negative, also βˆ [j] (ω) is non negative. Being e2m (ω) = 0 only for ω ∈ 2π Zd , there is no ω ∈ Rd such
that φˆ j (ω − 2π k) = 0 for all k ∈ Zd . Then by the Wiener’s lemma, there are unique absolutely summable coefficients a[j] =
{ak[j] }k∈Zd , j = 1, . . . , m − 1, such that [j] Λ2m (x) = ak φj (x − k).
(40)
k∈Zd
Those coefficients satisfy a[j] ∗ β [j] = δ, [ j]
and by Lemma 2 in [15], a[j] is such that |ak | ≤ C |k|−d−2−2j . The vector a[j] , j = 1, . . . , m − 1, can be explicitly computed via the iterative algorithm in [16]. By using this result, also the scaling coefficients s[j] , j = 1, . . . , m − 1 can be computed explicitly. Let us consider φj (A−1 x). We know that V−1 ⊂ V0 , and that the integer translates of Λ2m (· − k) are a basis for V0 , then φj (A−1 x) = [j]
where, bl = φj (A can write
−1
l∈Zd
[j]
bl Λ(x − l)
l). In addition (40) holds and since all the involved vectors are in l (Z ) and decay sufficiently fast, we 1
φj (A−1 x) =
=
d
[j]
bl Λ(x − l)
l∈Zd
l∈Zd
=
[j]
bl
[j]
ak φj (x − k − l)
k∈Zd
[j] [j]
bk ar −k φj (x − r )
r ∈Zd k∈Zd
=
s[rj] φj (x − r ),
r ∈Zd
where s = a ∗ b[j] . By [15, Lemma 2], the scaling coefficients have the same decay of a[j] and b[j] , then s[j] ∈ l1 (Zd ), j = 1, . . . , m − 1. For instance, if φ0 is the m-harmonic isotropic B-spline, d = 2, m = 2, φ1 , decays at infinity as ∥x∥−8 (see [10]), [ j] and then sk ≤ C |k|−8 . Let us restrict ourself to the case A = 2I. The MRA V associated with (Zd , 2I ) induces a dyadic subdivision scheme having mask s[j] . Namely, let λ ∈ l1 (Zd ), the stationary subdivision scheme S is defined by the rules [j]
[ j]
λ0 := λ
λk+1 := S λk ,
k≥0
where
(S λk )r =
[j]
sr −2l λkl .
l∈Zd
∈ l (Z ), the scheme converges to Qj fλ (x) = λl φj (x − l).
Since s
[j]
1
d
l∈Zd
Having the refinement coefficients is important in order to compute via dyadic subdivision the quasi-interpolant (24). In this way we significantly reduce the computational cost. The numerical examples show that convergence to the limit surface is fast, after 2 or three levels we are at the limit.
M. Bozzini, M. Rossini / Journal of Computational and Applied Mathematics 267 (2014) 96–106
103
Table 1 d = 1, maximum absolute approximation error for Q1 . f
e∞
1 x/2 x2/4 x3/8
2.3e−14 2.1e−13 2.2e−12 2.1e−11
Table 2 d = 2, m = 2, maximum absolute approximation error for Q1 . f
e∞
1 (x + y)/4 (x2 + y2 )/8 (x3 + y3 )/16
5.4e−15 4.3e−13 4.4e−6 1.0e−5
Table 3 d = 2, m = 2, maximum absolute approximation errors for Q1 . f
size of λ0
e∞
F1 F5
17 × 17 9×9
1.6e−2 2.3e−2
4. Numerical examples In this section we evaluate the performances of the subdivision scheme presented in Section 3 when applied to the computation of the quasi-interpolant Q1 (corresponding to the case m = 2) of some polynomial functions and smooth functions in d = 1 and d = 2. 4.1. The case d = 1 Here φ0 is the symmetric B-spline of degree 2m − 1 with integer knots. φj j = 1, . . . , m − 1 is a spline of degree 2m − 1 with support [−2j m, 2j m]. The mask s[j] has not compact support, but it decays fast (see Fig. 1). In the numerical examples its size is 13. We show in Table 1 the maximum errors e∞ we get when computing via subdivision the quasi-interpolant (24) with m = 2, j = 1, h = 1, and f (x) a polynomial of degree 0, 1, 2, 3. We consider λ0 = {λ0l , l ∈ [−9, 9] Z}. The error is evaluated on the interval [−2, 2] after three levels of refinement. 4.2. The case d = 2 In the two-dimensional examples φ0 is the m-harmonic isotropic B-spline with m = 2. φ1 decays at infinity as ∥x∥−8 (see [10]) and reproduces polynomials of degree up to 3. As in Section 4.1, we show the performance of the subdivision scheme by considering in Table 2 the maximum errors e∞ obtained when computing via subdivision the quasi-interpolant 2 (24) with m = 2, j = 1, h = 1, and f (x) a polynomial of degree 0, 1, 2, 3. We consider λ0 = {λ0l , l ∈ [−15, 15]2 Z }. The 2 maximum error e∞ is evaluated on [−2, 2] after three levels of refinement and it is shown in Table 2. In order to have a correct polynomial reproduction, as the mask decays algebraically, it is necessary to use a truncated mask of size 65 × 65. In addition, we consider the well known Franke’s function F 1 and the functions F 5 [17], F 5(x, y) =
1 3
exp −20.25((x − 0.5)2 + (y − 0.5)2 ) ,
x, y ∈ [0, 1]2 .
Both test functions are depicted in Fig. 2. In Table 3 we show, the size of the starting vector λ0 , the maximum absolute approximation errors for F 1 and F 5 after two and three levels of refinement respectively. In Figs. 3 and 4, the starting vectors and the final approximation are depicted. In these two experiments the size of the mask is 17 × 17. We conclude the section by considering an application which shows the potential of the method when applied to the resampling of an image onto a grid with higher density, for instance for zooming purposes. We apply the subdivision scheme for the magnification with a zooming factor 4 of a 32 × 32 detail of Lena’s image (see Fig. 5). In Fig. 6 we show the detail of the image (left) and its resampling after 4 levels of refinements (right). Note how well the oblique contours are rendered.
104
M. Bozzini, M. Rossini / Journal of Computational and Applied Mathematics 267 (2014) 96–106
Fig. 1. Mask s1 , m = 2.
Fig. 2. Left: F 1. Right: F 5.
Fig. 3. F 1. Left: λ0 17 × 17. Right: λ2 .
Acknowledgment The authors are grateful to an anonymous referee for her/his useful and constructive comments and remarks.
M. Bozzini, M. Rossini / Journal of Computational and Applied Mathematics 267 (2014) 96–106
105
Fig. 4. F 5. Left: λ0 9 × 9. Right: λ3 .
Fig. 5. The test image Lena. In the square, the 32 × 32 detail is depicted.
Fig. 6. Left: Lena’s detail: λ0 32 × 32. Right: reconstruction after 4 levels of refinement.
References [1] J. Duchon, Interpolation des fonctions de deux variables suivant le principe de la flexion des plaques minces, Rev. Française Automat. Informat. Rech. Opér. Anal. Numer. 10 (1976) 5–12. [2] M.D. Buhmann, Radial Basis Functions, in: Cambridge Monographs on Applied and Computational Mathematics, Cambridge University Press, 2004. [3] A. Iske, Multiresolution methods in scattered data modelling, in: Lecture Notes in Computational Science and Engineering, Springer-Verlag, Berlin, 2004. [4] C. Rabut, M. Rossini, Polyharmonic multiresolution analysis: an overview and some new results, 48 (1–3) (2008) 135–160. [5] M. Bozzini, M. Rossini, R. Schaback, Generalized Whittle-Matérn and polyharmonic kernels, Adv. Comput. Math. 39 (1) (2013) 129–141.
106
M. Bozzini, M. Rossini / Journal of Computational and Applied Mathematics 267 (2014) 96–106
[6] W.R. Madych, S.A. Nelson, Polyharmonic cardinal splines, J. Approx. Theory 60 (1990) 141–156. [7] N. Dyn, I.R.H. Jackson, D. Levin, A. Ron, On multivariate approximation by integer translates of basis function, Israel J. Math. 74 (1992) 95–130. [8] N. Dyn, Approximation by translates of a radial function, in: J.R. Higgins, R.L. Sten (Eds.), Sampling Theory in Fourier and Signal Analysis; Advanced Topics, Oxford University Press, 1999, pp. 187–208. [9] C. Rabut, Elementary polyharmonic cardinal B-splines, Numer. Algorithms 2 (1992) 39–62. [10] M. Bozzini, N. Dyn, M. Rossini, Construction of generators of quasi-interpolation operators of high approximation orders in spaces of polyharmonic splines, J. Comput. Appl. Math. 236 (4) (2011) 557–565. [11] W.R. Madych, Some elementary properties of multiresolution analyses of L2 (Rn ), in: C.K. Chui (Ed.), Wavelets: A Tutorial in Theory and Applications, Academic Press, 1992, pp. 259–294. [12] M. Bozzini, C. Rabut, M. Rossini, A multiresolution analysis with a new family of polyharmonic B-splines, in: A. Cohen, J.L. Merrien, L. Schumaker (Eds.), Curve and Surface Fitting: Avignon 2006, Nashboro Press, 2007, pp. 51–60. [13] M. Rossini, On the construction of polyharmonic B-splines, J. Comput. Appl. Math. 221 (2008) 437–446. [14] G.F. Fasshauer, Meshfree Approximation Methods with MATLAB, in: Interdisciplinary Mathematical Sciences, vol. 6, World Scientific Publishers, Singapore, 2007. [15] B. Bacchelli, M. Bozzini, C. Rabut, M.L. Varas, Decomposition and reconstruction of multidimensional signals using polyharmonic pre-wavelets, Appl. Comput. Harmon. Anal. 18 (3) (2005) 282–299. [16] B. Bacchelli, M. Bozzini, C. Rabut, A fast wavelet algorithm for multidimensional signals using polyharmonic splines, in: A. Choen, J.L. Merrien, L.L. Schumaker (Eds.), Curves and Surface Fitting: Saint Malo 2002, Nashboro Press, 2003, pp. 21–30. [17] R.J. Renka, R. Brown, Algorithm 792: accuracy tests of ACM algoritm for interpolation of scattered data in the plane, ACM TOMS 25 (1999) 78–94.
Contents lists available at ScienceDirect
Journal of Computational and Applied Mathematics journal homepage: www.elsevier.com/locate/cam
Properties of generators of quasi-interpolation operators of high approximation orders in spaces of polyharmonic splines Mira Bozzini, Milvia Rossini ∗,1 University of Milano-Bicocca, Italy
article
abstract
info
Article history: Received 8 March 2013 Received in revised form 14 November 2013 Keywords: Polyharmonic splines Quasi-interpolation operators High degree polynomial reproduction Multiresolution analysis Scaling functions Subdivision
We have presented in Bozzini et al. (2011) a procedure in spaces of m-harmonic splines in Rd that starts from a simple generator φ0 and recursively defines generators φ1 , φ2 , . . . , φm−1 with corresponding quasi-interpolation operators reproducing polynomials of degrees 3, 5, . . . , 2m − 1 respectively. In this paper we study the properties of generators φj , and we prove that these new generators are positive definite functions, and are scaling functions whenever φ0 has those properties. Moreover φ0 and φj generate the same multiresolution analysis. We show that it is possible to define a convergent subdivision scheme, and to provide in this way a fast computation of the quasi-interpolant. © 2014 Elsevier B.V. All rights reserved.
1. Introduction Polyharmonic splines have been introduced in [1] as a variational extension to many variables of polynomial splines in one variable. Because of their good properties (see e.g. [2–4]) they are studied by several authors, and widely used in scattered data approximation and interpolation. An important amount of research work has been done also on regular grids which occur in many applications such as image processing and multiresolution analysis. In general it is important to construct quasi-interpolation operators capable of reproducing some important functions of the space, and to provide a fast evaluation of the quasi-interpolant function. The space of m-harmonic splines V2m , 2m > d, is defined as the subspace of S ′ (Rd ) (the space of d-dimensional tempered distributions)
V2m = g ∈ S ′ (Rd )
C 2m−d−1 (Rd ) : ∆m g = 0, on Rd \ Zd ,
(1)
2 2 m m m−1 where ∆ := g ). Recently (see [5]) they i=1 ∂ /∂ xi is the Laplace operator and ∆ is defined iteratively by ∆ g = ∆(∆ m have been generalized by considering the fundamental solutions of more general differential operators j=1 (−∆ + κj2 I ).
d
It is well-known (see e.g. [6]) that V2m contains Π2m−1 (the space of polynomials defined on Rd of degree not exceeding 2m − 1), and that for any n ≤ 2m − 1 it is possible to construct quasi-interpolation (QI) operators reproducing Πn (see e.g. [7,8]). The generator φ ∈ V2m of the QI operator Qφ (x, h, f ) =
f (hl)φ(x/h − l)
l∈Zd
∗
Corresponding author. Tel.: +39 0264485715. E-mail addresses: [email protected] (M. Bozzini), [email protected] (M. Rossini).
1 Sponsored by the PRIN Project Real and Complex Varieties: Geometry, Topology, Harmonic Analysis. http://dx.doi.org/10.1016/j.cam.2014.01.029 0377-0427/© 2014 Elsevier B.V. All rights reserved.
(2)
M. Bozzini, M. Rossini / Journal of Computational and Applied Mathematics 267 (2014) 96–106
97
is required to decay fast enough at infinity, so that Qφ (x, h, f ) is well defined for f growing at infinity not faster than a polynomial of degree n. The QI operator (2) approximates smooth enough functions with L∞ -error of order hn+1 [7]. Important bases of V2m are m-harmonic B-splines, [9,4]. In the following we denote them by φ0 . They decay sufficiently fast, at least as ∥x∥−d−2 when ∥x∥ → ∞, they generate quasi-interpolation operators reproducing linear polynomials and then, for f with bounded derivatives of order 1, 2, 3, the sum Q0 (x, h, f ) =
f (hl)φ0 (x/h − l)
l∈Zd
approximates f in Rd with L∞ -error of order h2 . In addition they are valid scaling function for generating a multiresolution analysis of L2 (Rd ) [4]. To obtain higher approximation orders, the function φ0 should be replaced by a function, with a stronger decay rate as ∥x∥ → ∞. For this reason we have presented in [10] a procedure which starts from a simple generator φ0 , and recursively defines generators φ1 , φ2 , . . . , φm−1 with corresponding QI operators reproducing Π3 , Π5 , . . . , Π2m−1 respectively. The aim of this paper, is to study other important properties of the generators φj defined in [10]. We prove that these new generators are positive definite functions, are scaling functions whenever φ0 has those properties. Moreover φ0 and φj generate the same multiresolution analysis, and we show that it is possible to define a convergent subdivision scheme and to provide a computational refinement which allows us a fast computation of Qj (x, h, f ) =
f (hl)φj (x/h − l).
l∈Zd
The paper is organized as follows. In Section 2 we give some definitions, preliminary results, and we briefly recall the construction of the new generators φj . In Section 3 we give the properties of φj , and in Section 4 we show some numerical experiments. 2. Preliminaries 2.1. Definitions, notations, useful results
d
We use in this paper the multi-index notation. In particular for k ∈ Zd , |k| = i=1 |ki |. The symbol ∗ denotes the convolution operator between functions or sequences and the semi-discrete convolution between a sequence and a function. T := [−π , π]d is the d-dimensional torus, L2 (Rd /2π Zd ) is the set of functions locally square integrable and 2π Zd -periodic. In the sequel we assume 2m > d. We say that a linear transformation A on Rd is an acceptable dilation for Zd if it leaves Zd invariant, i.e. AZd ⊂ Zd , and all the eigenvalues of A satisfy |λi | > 1. These properties imply that |det A| is an integer ≥ 2. Important acceptable dilations for Z d are similarities, i.e matrices of the form A = ρ A0 , where A0 is an orthogonal matrix, and ρ is a real number such that |det A| = |ρ d | is an integer ≥ 2. In the case d = 2 these matrices are all of the form
q −p
p q
or
q p
p , −q
with q, p ∈ Z. They include the classic dyadic matrix 2I and the quincunx case q = p = 1 which gives |det A| = 2. A multiresolution analysis (MRA) of dilation A ∈ Rd×d is a sequence V = {Vl , l ∈ Z} of closed subspaces of L2 (Rd ) satisfying the following properties (i) (ii) (iii) (iv)
V l ⊂ Vl+1 , 2 d l∈Z Vl is dense in L (R ) and l∈Z Vl = {0}, −1 f (·) ∈ Vl ⇔ f (A ·) ∈ Vl−1 , there is a function of V0 , called the scaling function, whose integer translates form a Riesz basis of V0 .
Condition (iv) is equivalent to saying that V0 is the closed linear span of Zd translates of one function φ , which satisfies the condition c1 ≤
ˆ − 2π k)|2 ≤ c2 |φ(ω
k∈Zd
for some positive constants c1 , c2 . The whole multiresolution analysis may be regarded as being generated by φ . Note that the scaling function is not unique (see e.g. [11]). It is useful to remember the following results stated in [11, Section 2.3].
ˆ 0) = 1 and φ( ˆ 2kπ ) = 0 for all k ∈ Zd \ {0}, then φ is a valid scaling Lemma 1. Let φ ∈ L1 (Rd ) L2 (Rd ) be such that φ( function for generating a MRA V associated with (Zd , A), if and only if φ enjoys the following properties: (1) there is a function S (ω) ∈ L2 (Rd /2π Zd ) such that
ˆ ˆ B−1 ω) where B = AT ; φ(ω) = S (B−1 ω)φ( (2) for almost all ω ˆ − 2π k)|2 > 0. |φ(ω k∈Zd
(3)
(4)
98
M. Bozzini, M. Rossini / Journal of Computational and Applied Mathematics 267 (2014) 96–106
ˆ 0) = 1 and φ( ˆ 2kπ ) = 0 for all k ∈ Zd \ {0}, and suppose that φ generates Lemma 2. Let φ ∈ L1 (Rd ) L2 (Rd ) be such that φ( a MRA V associated with (Zd , A). Let γ = {γk }k∈Zd ∈ l1 (Zd ) such that γˆ (ω) ̸= 0, ω ∈ T . Then also the function
ϕ =γ ∗φ generates the MRA V associated with (Zd , A). The Green’s function of ∆m is
vm (x) =
Cm ∥x∥2m−d Cm ∥x∥2m−d log ∥x∥
2m − d ̸∈ 2Z 2m − d ∈ 2Z,
(5)
where the constant Cm is such that the generalized Fourier transform is
vˆ m (ω) = (−1)m ∥ω∥−2m .
(6)
A function φ0 : R → R is called m-harmonic B-spline if its Fourier transform has the form d
φˆ 0 (ω) = e2m (ω)ˆvm (ω),
(7)
where e2m is a real, even, and 2π Z -periodic trigonometric polynomial satisfying d
e2m (ω) = (−1)m ∥ω∥2m + O(∥ω∥2m+2 ),
(8)
for ω in a neighborhood of the origin U (0). For d = 1, φ0 is the symmetric B-spline of odd degree 2m − 1 with integer knots. In this case
e2m (ω) = −4 sin2
ω m 2
.
(9)
For d = 2, important examples of φ0 are the elementary m-harmonic B-spline, having
e2m (ω1 , ω2 ) = −4 sin2
ω1 2
− 4 sin2
ω2 m 2
,
(10)
and the isotropic m-harmonic B-spline having (see e.g. [12,13])
e2m (ω1 , ω2 ) = −
m ω1 ω2 2 1 + 4 sin2 + 4 sin2 + cos ω1 cos ω2 . 3 2 2
(11)
Generally, the trigonometric polynomial e2m , as happens for the elementary and isotropic B-spline, is requested to satisfy e2m (ω) = 0,
ω∈T
ω = 0.
⇔
(12)
This guarantees that φ0 is a valid scaling function for generating a MRA V associated with (Z , A) whenever A is a similarity. d
2.2. Construction of φj , j = 1, . . . , m Let us consider a m-harmonic spline φ0 (7). Starting from φ0 we have constructed in [10] new generators
φ1 , φ2 , . . . , φm−1 of QI operators reproducing polynomials up to degree 3, 5, . . . , 2m − 1 respectively. The procedure defines φj recursively as a linear combination of φj−1 (·) and φj−1 (·/2) with explicitly known coefficients which are independent of m. For j = 1, 2, . . . , m − 1, we define
φj (x) = aj φj−1 (x) + bj φj−1 (x/2),
x ∈ Rd ,
(13)
and choose the coefficients aj , bj so that
φˆ j (ω) = 1 +
∞
[ j]
h2i (ω),
ω ∈ U (0),
(14)
i=j+1
[j]
with h2i a homogeneous function of order 2i. This is possible since
φˆ j (ω) = aj φˆ j−1 (ω) + 2d bj φˆ j−1 (2ω),
(15)
and by our inductive hypothesis
φˆ j−1 (ω) = 1 +
∞ i=j
[j−1]
h2i
(ω),
ω ∈ U (0).
(16)
M. Bozzini, M. Rossini / Journal of Computational and Applied Mathematics 267 (2014) 96–106
99
The coefficients aj , bj in (15) are the solution of the system
aj + 2 d b j = 1 aj + 2d+2j bj = 0,
(17)
that is aj =
22j 22j
,
−1
bj = −
1 2d (22j − 1)
,
j = 1, . . . , m − 1.
(18)
Note that the coefficients aj , bj depend on neither φ0 nor m and that bj = −
aj 2d+2j
.
(19)
In [10], it has been proved that φj has the correct decay at infinity needed for the quasi-interpolation based on it to be well defined for polynomials of degree up to 2j + 1, and that these polynomials are reproduced by this quasi-interpolation. Namely for j = 1, 2, . . . , m − 1
φj (x) ≤
C
∥x ∥
d+2j+2
,
∥x∥ → ∞,
(20)
and φj satisfies the Strang Fix conditions of order 2j + 1, i.e.
φˆ j (0) = 1,
Dα φˆ j (0) = 0,
α
D φˆ j (2kπ ) = 0,
1 ≤ |α| ≤ 2j + 1,
(21)
k ∈ Z \ {0}, |α| ≤ 2j + 1. d
(22)
Consequently for j = 1, . . . , m − 1, the sum
f (l)φj (x − l)
(23)
l∈Zd
is well defined for any f growing at infinity not faster than a polynomial of degree 2j + 1, and if f ∈ Π2j+1 , the sum (23) equals f . Moreover Qj (x, h, f ) =
f (hl)φj (x/h − l)
(24)
l∈Zd
approximates f having bounded derivatives of order 2j + 1, 2j + 2, 2j + 3, with L∞ -error of order h2j+2 . 3. Properties of φj In this section we study some additional properties of the generators φj . It is trivial to verify that φj ∈ C 2m−d−1 (Rd ), and that (see (20)) φj ∈ L1 (Rd ) L2 (Rd ), j = 0, . . . , m − 1. Now we want to investigate if φj , j = 1, . . . , m − 1, preserves the properties of φ0 . In particular we are interested in proving that a strictly positive definite function φ0 , under some assumptions on e2m , generates strictly positive definite functions φj , and that if φ0 is a scaling function, it generates functions which are scaling functions too. In addition, we will prove that φ0 and φj generate the same multiresolution analysis of L2 (Rd ). We start by recalling the definition of positive definite function and its characterization (see e.g. [14]). [0]
Definition 3. A complex valued function Φ is called positive definite on Rd if N N
cl c¯k Φ (xl − xk ) ≥ 0
(25)
l=0 k=0
for any pairwise different points x1 , . . . , xN ∈ Rd , and c = [c1 , . . . , cN ] ∈ CN . The function Φ is called strictly positive definite on Rd if (25) is zero only for c ≡ 0. It is known that (see e.g. [14]). Proposition 4. Let Φ be a continuous function in L1 (Rd ). Φ is strictly positive definite if and only if Φ is bounded and its Fourier transform is non negative and not identically equal to zero. In virtue of this fundamental result, φ0 strictly positive definite implies that
φˆ 0 (ω) = (−1)m e2m (ω)∥ω∥−2m ≥ 0
100
M. Bozzini, M. Rossini / Journal of Computational and Applied Mathematics 267 (2014) 96–106
and not identically equal to zero, i.e. e2m (ω) := (−1)m e2m (ω) ≥ 0. [0]
(26)
By construction φˆ j (ω) has the form j] φˆ j (ω) = e[2m (ω)∥ω∥−2m ,
(27)
where [j]
[j−1]
[j−1]
e2m (ω) = aj e2m (ω) + 2d−2m bj e2m (2ω),
j = 1, . . . , m − 1.
(28)
By using (19), we get
[j]
[j−1]
[j−1]
e2m (ω) = aj e2m (ω) − 2−2m−2j e2m (2ω) ,
j = 1, . . . , m − 1.
(29)
Proposition 5 gives a sufficient condition that ensures that φj is strictly positive definite and, as we shall see in Proposition 6, this condition is always satisfied when we start from the most known m-harmonic B-splines: the symmetric odd degree B-spline with integer knots (d = 1), the elementary and the isotropic B-spline (d = 2). Proposition 5. Let φ0 be strictly positive definite and let
αj := 2−2m−2
j −1
2−2r = 2−2m
1 − 2−2j
.
3
r =0
(30)
If for j = 1, . . . , m − 1 e2m (ω) − αj e2m (2ω) ≥ 0, [0]
[0]
(31)
then φj is strictly positive definite. [j]
Proof. By recursion, we express each e2m , j = 1, . . . , m − 1, in terms of e2m . Let us begin with e2m . By definition
[0]
e2m (ω) = a1 e2m (ω) − 2−2m−2 e2m (2ω) , [1]
[0]
[0]
[1]
a1 > 0.
Trivially, if e2m (ω) ≥ 2−2m−2 e2m (2ω) and e2m (ω) is not identically equal to 2−2m−2 e2m (2ω), also e2m is non negative and not identically equal to zero. Since φ1 ∈ L1 (Rd ), φ1 is strictly positive definite. For j = 2, we get [0]
[0]
[0]
[0]
e2m (ω) = a2 e2m (ω) − 2−2m−4 e2m (2ω) , [2]
[1]
[1]
= a1 a2 e2m (ω) − 2 [0]
−2m−2
1+
[1]
a2 > 0
1 22
e2m (2ω) + 2 [0]
e2m (4ω) .
−2m−2 −2m−4 [0]
2
[0] [0] [0] [2] Since e2m (·) ≥ 0 and not identically zero, it follows that if e2m (ω) ≥ 2−2m−2 (1 + 212 )e2m (2ω) also e2m is non negative and not [j]
[0]
identically equal to zero. Going on with this simple stuff, it is not difficult to see that e2m can be written in terms of e2m as
[j]
e2m (ω) =
j
e2m (ω) − 2 [0]
al
−2m−2
l =1
+
j −2
j −1
e2m (2ω)
−2r [0]
2
(32)
r =0
−4m−4r −6
2
r =0
j
j−(r +2)
ak
k=r +1
e2m (4ω) .
−2k [r ]
2
(33)
k=0
Since e2m (ω) ≥ 0, r = 0, . . . , j − 1, and not identically equal to zero, from assumption (31) we get the proof. [r ]
Proposition 6. When φ0 is the symmetric odd degree B-spline with integer knots (d = 1), the elementary or the isotropic B-spline j] is non negative and vanishes only at ω ∈ 2π Zd . (d = 2), the functions φj , j = 1, . . . , m are strictly positive definite. Namely e[2m Proof. We give the proof for elementary B-spline of order m > 1. Analogous arguments can be carried out for the odd degree B-spline with integer knots and the isotropic m-harmonic B-spline for which the trigonometric polynomials are (9) and (11) respectively. [0] The trigonometric polynomial e2m is (see (10) and (26))
e2m (ω1 , ω2 ) = 4m sin2 [0]
ω1 2
+ sin2
ω2 m 2
.
(34)
M. Bozzini, M. Rossini / Journal of Computational and Applied Mathematics 267 (2014) 96–106
101
To verify condition (31) means to verify if, for j = 1, . . . , m − 1
sin2
ω1 2
+ sin2
ω2 m 2
m − αj sin2 ω1 + sin2 ω2 ≥ 0.
(35)
We can factorize the left hand side as
sin2
ω1 2
+ sin2
ω2
−1 m m−1−k ω1 ω2 k √ √ m α (sin2 ω + sin2 ω ) m α (sin2 ω + sin2 ω ) . sin2 + sin2 j 1 2 j 1 2
−
2
k=0
2
(36)
2
The sum m−1
sin2
k=0
ω1 2
ω2 k √ m
+ sin2
2
m−1−k αj (sin2 ω1 + sin2 ω2 )
is always non negative and it is zero only for (ω1 , ω2 ) = (2lπ , 2r π ) l, r ∈ Z. So we have to study the sign of p(ω1 , ω2 ) := sin2
ω1 2
+ sin2
ω2 2
−
2 √ 2 m α j sin ω1 + sin ω2 ,
(37)
that can be rewritten as p(ω1 , ω2 ) = sin2
= sin2
ω1
ω2 √ αj sin2 ω1 + sin2 − m αj sin2 ω2 2 ω1 ω2 ω2 √ √ 1 − 4 m αj cos2 + sin2 1 − 4 m αj cos2 .
− 2 ω1 2
√ m
2
√
2
2
From definition (30), 4 m αj < 1 for any j = 1, . . . , m − 1, then p(ω1 , ω2 ) is non negative, and it is equal to zero only when (ω1 , ω2 ) = (2lπ , 2r π ) l, r ∈ Z. By recursion we get the proof. It is known that the functions φ0 of Proposition 6 are scaling functions of a MRA V of L2 (Rd ) associated with (Zd , A), A a similarity. In Proposition 7, we prove that also φj , j = 1, . . . , m − 1 is a scaling function, and that φ0 and φj generate the same multiresolution analysis. Proposition 7. When φ0 is the symmetric odd degree B-spline with integer knots (d = 1), the elementary or the isotropic B-splines (d = 2), φj , j = 1, . . . , m − 1, is a scaling function of a MRA V of L2 (Rd ) associated with (Zd , A), A a similarity. In addition φj , j = 1, . . . , m − 1, and φ0 generate the same MRA. [j]
Proof. We observe that whenever A is a similarity, there is a function S2m (ω) ∈ L2 (Rd /2π Zd ) such that [j] φˆ j (Bω) = S2m (ω)φˆ j (ω),
B = AT ,
(38)
where [j]
e (Bω) [j] . S2m (ω) = ρ −2m 2m [ j] e2m (ω) Then, we consider
Φj (ω) =
2 φˆ j (ω − 2π k) ,
j = 1, . . . , m − 1.
(39)
k∈Z2
Φj is well defined and bounded for any ω because of the decay of vˆ m . Moreover φˆ j (0) = 1 and φˆ j (2π k) = 0, k ∈ Zd k ̸= 0, then for Lemma 1 we can conclude that φj , j = 1, . . . , m − 1, generates a MRA. [j]
[j]
In addition, e2m (ω) and e2m (ω) are equal to zero only at ω = 2π k k ∈ Zd , and the ratio between e2m (ω) and e2m (ω) is [0]
[0]
[j]
well-defined for any ω and equal to one at ω = 2π k k ∈ Zd . Then, since φˆ j (ω) = e2m (ω)ˆvm (ω), we can trivially rewrite φˆ j as [j]
φˆ j (ω) =
[j]
e2m (ω) [0] e (ω) e2m (ω)ˆvm (ω) = 2m φˆ 0 (ω), [0] [0] e2m (ω) e2m (ω) [j]
where the function γˆj (ω) := e2m (ω)/e2m (ω) is equal to one for ω = 2kπ , k ∈ Zd , and positive everywhere. Moreover γˆj ∈ L2 (Rd /2π Zd ) and then, for Lemma 2, φj j = 1, . . . , m − 1, and φ0 generate the same MRA V associated with (Zd , A). [j]
[0]
S2m (ω) is the 2π Zd periodic function [j]
S2m (ω) = |ρ|−d
k∈Zd
[j] sk e−i⟨k,ω⟩ .
102
M. Bozzini, M. Rossini / Journal of Computational and Applied Mathematics 267 (2014) 96–106
[j]
The coefficients s[j] = {sk }k∈Zd are the scaling coefficients of φj , i.e
φj (·) =
[j]
sk φj (A · −k),
k∈Zd
and as we show in the following, they can be computed explicitly via convolutions. To this end, we first show that the m-harmonic Lagrangian function Λ2m ∈ V2m (Λ2m (l) = δ0l , l ∈ Zd ) can be expressed in terms of shifts of φj . [j]
Let us consider the vector β [j] = {φj (k)}k∈Zd . From (20), β [j] ∈ l1 (Zd ) and, at least |βk | ≤ C |k|−d−2−2j . The absolutely convergent trigonometric series
βˆ [j] (ω) =
φj (k)e−iω·k ̸= 0,
∀ω ∈ Rd .
k∈Zd
In fact, from an application of the Poisson summation formula
φj (k)e−iω·k =
k∈Zd
φˆ j (ω − 2π k),
k∈Zd
[j]
and since φˆ j is non negative, also βˆ [j] (ω) is non negative. Being e2m (ω) = 0 only for ω ∈ 2π Zd , there is no ω ∈ Rd such
that φˆ j (ω − 2π k) = 0 for all k ∈ Zd . Then by the Wiener’s lemma, there are unique absolutely summable coefficients a[j] =
{ak[j] }k∈Zd , j = 1, . . . , m − 1, such that [j] Λ2m (x) = ak φj (x − k).
(40)
k∈Zd
Those coefficients satisfy a[j] ∗ β [j] = δ, [ j]
and by Lemma 2 in [15], a[j] is such that |ak | ≤ C |k|−d−2−2j . The vector a[j] , j = 1, . . . , m − 1, can be explicitly computed via the iterative algorithm in [16]. By using this result, also the scaling coefficients s[j] , j = 1, . . . , m − 1 can be computed explicitly. Let us consider φj (A−1 x). We know that V−1 ⊂ V0 , and that the integer translates of Λ2m (· − k) are a basis for V0 , then φj (A−1 x) = [j]
where, bl = φj (A can write
−1
l∈Zd
[j]
bl Λ(x − l)
l). In addition (40) holds and since all the involved vectors are in l (Z ) and decay sufficiently fast, we 1
φj (A−1 x) =
=
d
[j]
bl Λ(x − l)
l∈Zd
l∈Zd
=
[j]
bl
[j]
ak φj (x − k − l)
k∈Zd
[j] [j]
bk ar −k φj (x − r )
r ∈Zd k∈Zd
=
s[rj] φj (x − r ),
r ∈Zd
where s = a ∗ b[j] . By [15, Lemma 2], the scaling coefficients have the same decay of a[j] and b[j] , then s[j] ∈ l1 (Zd ), j = 1, . . . , m − 1. For instance, if φ0 is the m-harmonic isotropic B-spline, d = 2, m = 2, φ1 , decays at infinity as ∥x∥−8 (see [10]), [ j] and then sk ≤ C |k|−8 . Let us restrict ourself to the case A = 2I. The MRA V associated with (Zd , 2I ) induces a dyadic subdivision scheme having mask s[j] . Namely, let λ ∈ l1 (Zd ), the stationary subdivision scheme S is defined by the rules [j]
[ j]
λ0 := λ
λk+1 := S λk ,
k≥0
where
(S λk )r =
[j]
sr −2l λkl .
l∈Zd
∈ l (Z ), the scheme converges to Qj fλ (x) = λl φj (x − l).
Since s
[j]
1
d
l∈Zd
Having the refinement coefficients is important in order to compute via dyadic subdivision the quasi-interpolant (24). In this way we significantly reduce the computational cost. The numerical examples show that convergence to the limit surface is fast, after 2 or three levels we are at the limit.
M. Bozzini, M. Rossini / Journal of Computational and Applied Mathematics 267 (2014) 96–106
103
Table 1 d = 1, maximum absolute approximation error for Q1 . f
e∞
1 x/2 x2/4 x3/8
2.3e−14 2.1e−13 2.2e−12 2.1e−11
Table 2 d = 2, m = 2, maximum absolute approximation error for Q1 . f
e∞
1 (x + y)/4 (x2 + y2 )/8 (x3 + y3 )/16
5.4e−15 4.3e−13 4.4e−6 1.0e−5
Table 3 d = 2, m = 2, maximum absolute approximation errors for Q1 . f
size of λ0
e∞
F1 F5
17 × 17 9×9
1.6e−2 2.3e−2
4. Numerical examples In this section we evaluate the performances of the subdivision scheme presented in Section 3 when applied to the computation of the quasi-interpolant Q1 (corresponding to the case m = 2) of some polynomial functions and smooth functions in d = 1 and d = 2. 4.1. The case d = 1 Here φ0 is the symmetric B-spline of degree 2m − 1 with integer knots. φj j = 1, . . . , m − 1 is a spline of degree 2m − 1 with support [−2j m, 2j m]. The mask s[j] has not compact support, but it decays fast (see Fig. 1). In the numerical examples its size is 13. We show in Table 1 the maximum errors e∞ we get when computing via subdivision the quasi-interpolant (24) with m = 2, j = 1, h = 1, and f (x) a polynomial of degree 0, 1, 2, 3. We consider λ0 = {λ0l , l ∈ [−9, 9] Z}. The error is evaluated on the interval [−2, 2] after three levels of refinement. 4.2. The case d = 2 In the two-dimensional examples φ0 is the m-harmonic isotropic B-spline with m = 2. φ1 decays at infinity as ∥x∥−8 (see [10]) and reproduces polynomials of degree up to 3. As in Section 4.1, we show the performance of the subdivision scheme by considering in Table 2 the maximum errors e∞ obtained when computing via subdivision the quasi-interpolant 2 (24) with m = 2, j = 1, h = 1, and f (x) a polynomial of degree 0, 1, 2, 3. We consider λ0 = {λ0l , l ∈ [−15, 15]2 Z }. The 2 maximum error e∞ is evaluated on [−2, 2] after three levels of refinement and it is shown in Table 2. In order to have a correct polynomial reproduction, as the mask decays algebraically, it is necessary to use a truncated mask of size 65 × 65. In addition, we consider the well known Franke’s function F 1 and the functions F 5 [17], F 5(x, y) =
1 3
exp −20.25((x − 0.5)2 + (y − 0.5)2 ) ,
x, y ∈ [0, 1]2 .
Both test functions are depicted in Fig. 2. In Table 3 we show, the size of the starting vector λ0 , the maximum absolute approximation errors for F 1 and F 5 after two and three levels of refinement respectively. In Figs. 3 and 4, the starting vectors and the final approximation are depicted. In these two experiments the size of the mask is 17 × 17. We conclude the section by considering an application which shows the potential of the method when applied to the resampling of an image onto a grid with higher density, for instance for zooming purposes. We apply the subdivision scheme for the magnification with a zooming factor 4 of a 32 × 32 detail of Lena’s image (see Fig. 5). In Fig. 6 we show the detail of the image (left) and its resampling after 4 levels of refinements (right). Note how well the oblique contours are rendered.
104
M. Bozzini, M. Rossini / Journal of Computational and Applied Mathematics 267 (2014) 96–106
Fig. 1. Mask s1 , m = 2.
Fig. 2. Left: F 1. Right: F 5.
Fig. 3. F 1. Left: λ0 17 × 17. Right: λ2 .
Acknowledgment The authors are grateful to an anonymous referee for her/his useful and constructive comments and remarks.
M. Bozzini, M. Rossini / Journal of Computational and Applied Mathematics 267 (2014) 96–106
105
Fig. 4. F 5. Left: λ0 9 × 9. Right: λ3 .
Fig. 5. The test image Lena. In the square, the 32 × 32 detail is depicted.
Fig. 6. Left: Lena’s detail: λ0 32 × 32. Right: reconstruction after 4 levels of refinement.
References [1] J. Duchon, Interpolation des fonctions de deux variables suivant le principe de la flexion des plaques minces, Rev. Française Automat. Informat. Rech. Opér. Anal. Numer. 10 (1976) 5–12. [2] M.D. Buhmann, Radial Basis Functions, in: Cambridge Monographs on Applied and Computational Mathematics, Cambridge University Press, 2004. [3] A. Iske, Multiresolution methods in scattered data modelling, in: Lecture Notes in Computational Science and Engineering, Springer-Verlag, Berlin, 2004. [4] C. Rabut, M. Rossini, Polyharmonic multiresolution analysis: an overview and some new results, 48 (1–3) (2008) 135–160. [5] M. Bozzini, M. Rossini, R. Schaback, Generalized Whittle-Matérn and polyharmonic kernels, Adv. Comput. Math. 39 (1) (2013) 129–141.
106
M. Bozzini, M. Rossini / Journal of Computational and Applied Mathematics 267 (2014) 96–106
[6] W.R. Madych, S.A. Nelson, Polyharmonic cardinal splines, J. Approx. Theory 60 (1990) 141–156. [7] N. Dyn, I.R.H. Jackson, D. Levin, A. Ron, On multivariate approximation by integer translates of basis function, Israel J. Math. 74 (1992) 95–130. [8] N. Dyn, Approximation by translates of a radial function, in: J.R. Higgins, R.L. Sten (Eds.), Sampling Theory in Fourier and Signal Analysis; Advanced Topics, Oxford University Press, 1999, pp. 187–208. [9] C. Rabut, Elementary polyharmonic cardinal B-splines, Numer. Algorithms 2 (1992) 39–62. [10] M. Bozzini, N. Dyn, M. Rossini, Construction of generators of quasi-interpolation operators of high approximation orders in spaces of polyharmonic splines, J. Comput. Appl. Math. 236 (4) (2011) 557–565. [11] W.R. Madych, Some elementary properties of multiresolution analyses of L2 (Rn ), in: C.K. Chui (Ed.), Wavelets: A Tutorial in Theory and Applications, Academic Press, 1992, pp. 259–294. [12] M. Bozzini, C. Rabut, M. Rossini, A multiresolution analysis with a new family of polyharmonic B-splines, in: A. Cohen, J.L. Merrien, L. Schumaker (Eds.), Curve and Surface Fitting: Avignon 2006, Nashboro Press, 2007, pp. 51–60. [13] M. Rossini, On the construction of polyharmonic B-splines, J. Comput. Appl. Math. 221 (2008) 437–446. [14] G.F. Fasshauer, Meshfree Approximation Methods with MATLAB, in: Interdisciplinary Mathematical Sciences, vol. 6, World Scientific Publishers, Singapore, 2007. [15] B. Bacchelli, M. Bozzini, C. Rabut, M.L. Varas, Decomposition and reconstruction of multidimensional signals using polyharmonic pre-wavelets, Appl. Comput. Harmon. Anal. 18 (3) (2005) 282–299. [16] B. Bacchelli, M. Bozzini, C. Rabut, A fast wavelet algorithm for multidimensional signals using polyharmonic splines, in: A. Choen, J.L. Merrien, L.L. Schumaker (Eds.), Curves and Surface Fitting: Saint Malo 2002, Nashboro Press, 2003, pp. 21–30. [17] R.J. Renka, R. Brown, Algorithm 792: accuracy tests of ACM algoritm for interpolation of scattered data in the plane, ACM TOMS 25 (1999) 78–94.