Splines and Multiresolution Analysis - Springer Link
Splines and Multiresolution Analysis Brigitte Forster .
Introduction... ............. ............. ............ ............. ............. ..
.
Historical Background.. .. .. .. .. .. ... .. .. .. .. .. .. .. .. .. .. ... .. .. .. .. .. .. .. .. .. .
. .. ... ... ... ... ... .. ..
Mathematical Modeling and Application.. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. . Mathematical Foundations. .. .. .. .. .. .. .. .. .. .. .. .. .. .. . .. .. .. .. .. .. .. .. .. .. .. .. .. .. . Regularity and Decay Under the Fourier Transform. . . . . . . . . . . .. . . . . . . . . . . . . . . . . Criteria for Riesz Sequences and Multiresolution Analyses. . . . . . . . . . . . . . . . . . . . . Regularity of Multiresolution Analysis. ... ... ... ... ... ... .... ... ... ... ... ... ... .... . Order of Approximation. ... ........ ........ ........ ....... ........ ........ ........ .... Wavelets. .. .. ... .. ... .. ... .. ... .. ... .. ... .. ... .. ... .. ... .. ... .. ... .. ... .. ... .. ... .. ... .. . B-Splines. .... ...... ....... ...... ...... ...... ....... ...... ...... ...... ....... ...... ...... . Polyharmonic B-Splines. ....... ....... ....... ....... ....... ....... ....... ....... ......
. ..
Survey on Mathematical Analysis of Methods. . . . .. . . . . . .. . . . . . .. . . . . . .. . . . . . Schoenberg’s B-Splines for Image Analysis – the Tensor Product Approach... ........ ....... ....... ........ ....... ........ ....... ....... ........ ....... .... .. Fractional and Complex B-Splines. . . .. .. . .. .. . .. .. . .. .. . .. .. . .. .. . .. .. .. . .. .. . .. .. . .. Polyharmonic B-Splines and Variants. . .. . .. . .. . .. . .. . .. . .. . .. . .. . .. . .. . .. . .. . .. . .. .. Splines on Other Lattices. . .. .. . .. .. .. .. . .. .. .. . .. .. .. .. . .. .. .. . .. .. .. .. . .. .. .. . .. .. .. ... Splines on the Quincunx Lattice. .. .... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... . ... Splines on the Hexagonal Lattice. .. .. .. .. . .. .. .. .. .. .. .. .. .. . .. .. .. .. .. .. .. .. .. . .. .. .
Numerical Methods. ........... ................ ................ ................ .
.
Open Questions. ..... ..... ..... .... ..... ..... ..... ..... .... ..... ..... ..... .... ..
.
Conclusion.. .... .... ..... .... .... ..... .... .... ..... .... .... ..... .... .... .... ....
.
Cross-References. .. ..... ..... ..... ..... ...... ..... ..... ..... ..... ..... ...... ....
Otmar Scherzer (ed.), Handbook of Mathematical Methods in Imaging, DOI ./---_, © Springer Science+Business Media LLC
Splines and Multiresolution Analysis
Abstract: Splines and multiresolution are two independent concepts, which – considered together – yield a vast variety of bases for image processing and image analysis. The idea of a multiresolution analysis is to construct a ladder of nested spaces that operate as some sort of mathematical looking glass. It allows to separate coarse parts in a signal or in an image from the details of various sizes. Spline functions are piecewise or domainwise polynomials in one dimension (D) resp. nD. There is a variety of spline functions that generate multiresolution analyses. The viewpoint in this chapter is the modeling of such spline functions in frequency domain via Fourier decay to generate functions with specified smoothness in time domain resp. space domain. The mathematical foundations are presented and illustrated at the example of cardinal B-splines as generators of multiresolution analyses. Other spline models such as complex B-splines, polyharmonic splines, hexagonal splines, and others are considered. For all these spline families exist fast and stable multiresolution algorithms which can be elegantly implemented in frequency domain. The chapter closes with a look on open problems in the field. AMS Subject Classification (): A Spline approximation D Numerical analysis – Splines U Computing methodologies and applications – Image processing T Numerical methods in Fourier analysis Keywords Cardinal B-splines ⋅ image processing ⋅ multiresolution analysis ⋅ order of approximation ⋅ polyharmonic B-splines ⋅ riesz basis ⋅ scaling function two-scale relation
.
Introduction
This chapter deals with two originally independent concepts, which were recently combined, and since together have a strong impact on signal and image analysis: the concept of splines, i.e., of piecewise polynomials, and the concept of multiresolution analysis, i.e., splitting functions – or more general data – into coarse approximations and details of various sizes. Eversince the combination of the two concepts, they have led to a load of new applications in, e.g., signal and image analysis, as well as in signal and image reconstruction, computer vision, numerics of partial differential equations, and other numerical fields. An impression of the vast area of applications can be gained in, e.g., [, , , ]. Already the spline functions alone proved to be very useful for mathematical analysis as well as for signal and image processing, analysis and representation, for computer graphics and many more, see, e.g., [, , , , , , ]. An example for a family of spline functions are I. J. Schoenberg’s polynomial splines with uniform knots [, ]: β m (t) =
m+ k m+ ) (t − k)+m , ∑ (−) ( m! k= k
m ∈ N.
(.)
Here, t+m denotes the one-sided power function, i.e., t+m = for t < and t+m = t m for t ≥ .
Splines and Multiresolution Analysis
1 0.8 0.6 0.4 0.2
–1
1
2
4
3
5
⊡ Fig. -
Cardinal B-splines of degree m = , . . . ,
The B-splines β m can be easily generated by an iterative process. Let β (t) = χ[,) (t) be the characteristic function of the interval [, ). Then the B-spline of degree m is derived by the convolution product β m = β ∗ β m− = β ∗ . . . ∗ β
m+-times where
for m ∈ N,
β ∗ β m− (x) = ∫ β (y)β m− (x − y) d y = ∫ R
(.)
β m− (x − y) d y.
For an illustration of the cardinal B-splines, see > Fig. -. Splines had their second breakthrough as Battle [] and Lemarié [] discovered that B-splines generate multiresolution analyses. The simple form of the B-splines and their compact support, in particular, were convenient for designing multiresolution algorithms and fast implementations. Definition Let A : Rn → Rn be a linear mapping that leaves Zn invariant, i.e., A(Zn ) ⊂ Zn and that has (real or complex) eigenvalues with absolute values greater than . A multiresolution analysis associated with the dilation matrix A is a sequence of nested subspaces (Vj ) j∈Z of L (Rn ) such that the following conditions hold: (i) (ii) (iii) (iv) (v) (vi)
. . . ⊂ V− ⊂ V ⊂ V ⊂ . . ., ⋂ j∈Z Vj = {}, Span ⋃ j∈Z Vj is dense in L (Rn ), f ∈ Vj ⇐⇒ f (A− j ●) ∈ V , f ∈ V ⇐⇒ f (● − k) ∈ V for all k ∈ Zn . There exists a so-called scaling function φ ∈ V such that the family {φ(● − k)} k∈Zn of translates of φ forms a Riesz basis of V .
Here, L (Rn ) denotes the vector space of square-integrable functions f : Rn → C with norm ∥ f ∥ = (∫ ∣ f (x)∣ dx) Rn
Splines and Multiresolution Analysis
and corresponding inner product ⟨ f , g⟩ = ∫
Rn
f (x) g(x) dx,
where g denotes the complex conjugate of g. The elements in L (Rn ) are also called functions of finite energy. Riesz bases are a slightly more general concept than orthonormal bases. In fact, Riesz bases are equivalent to orthonormal bases and therefore can be generated by applying a topological isomorphism on some orthonormal basis. Definition A sequence of functions { f n }n∈Z in a Hilbert space V is called a Riesz sequence if there exist positive constants A and B, the Riesz bounds, such that A∥c∥ l ≤ ∥ ∑ c k f k ∥ ≤ B∥c∥ l k∈Z n
V
for all scalar sequences C = (c k ) k∈Zn ⊂ l (Zn ). A Riesz sequence is called a Riesz basis if it additionally spans the space V . A good introduction to Riesz bases, their properties, and their relation to orthonormal bases is given in the monography by Young []. Multiresolution constructions with splines are treated in numerous sources. As a starting point, there are, e.g., the books by Christensen [, ] and Wojtaszczyk []. The mathematical properties in Definition have intuitive interpretations. A function f ∈ L (Rn ), which is projected orthogonally on Vj , is approximated with the so-called resolution A j . In fact, let P j : L (Rn ) → Vj denote the orthogonal projection operator. Then (ii) gives that by going to lower resolutions, all details are lost: lim ∥P j f ∥ = .
j→−∞
Whereas when the resolution is increased, j → ∞, more and more details are added. By (iii), the projection then converges to the original function f : lim ∥ f − P j f ∥ = .
j→∞
Hereby, the rate of convergence depends on the regularity of f . The approximation spaces Vj are nested, which allows for computing coarser approximations in Vk for k < j for functions f ∈ Vj . The scaling Ak enlarges details. Property (iv) shows that the approximation spaces have a similar structure over the scales and emanate from one another. The translation invariance (v) ensures that the analysis of a function in Vj is independent of the starting time or location. And property (vi) finally ensures the beautiful and mighty property that the whole sequence of nested approximation spaces can be generated by translates and scalings of one
Splines and Multiresolution Analysis
single function – the scaling function. In fact, (vi) together with (iv) yields that {φ(A j ● −k), k ∈ Zn } is a Riesz basis for Vj . While moving from the coarser approximation space Vj to the finer, larger space Vj+ , information has to be added. In fact, there is an orthogonal complement Wj , j ∈ Z, such that Vj+ = Vj ⊕ Wj . These spaces are called detail spaces or wavelet spaces. It is well known that these spaces also possess a Riesz basis spanned by shifts of ∣ det A∣ − generators, the wavelets ψ , . . . , ψ∣ det A∣− . Here, A is the dilation matrix in Definition . The wavelets can be constructed from the scaling function. As a consequence, the knowledge of just the single function φ allows for the construction of the approximation spaces Vj and for the wavelet spaces Wj . Detailed information on the generation of wavelets and their properties can be found in various books, e.g., [, , , , ]. Example A simple example for a multiresolution analysis on L (R) is given by piecewise constant functions. Consider the characteristic function φ = χ[,) of the interval semiopen interval [, ). Then φ generates a dyadic multiresolution analysis, i.e., for A = . The approximation spaces are Vj = span {χ[,) ( j ● −k)} k∈Z
L (R)
.
They consist of functions constant on intervals of the form [k− j , (k + )− j ). The spaces are obviously nested and separate L (R) in the sense of Definition (ii). Since piecewise constant functions are dense in L (R), (iii) holds. (iv) – (vi) hold by construction. In fact, this multiresolution analysis is generated by the B-spline β as scaling function. The B-spline basis operates as mean-value operator over the support interval. The corresponding wavelet extracts the details, i.e., the deviation from the mean value. To this end, it operates as a difference operator. > Figure - shows the scaling function β and the corresponding wavelet, the so-called Haar wavelet. In > Fig. -, an example of a multiresolution is given.
1
1
0.8 0.5 0.6 0.2
0.4
0.4
0.6
0.8
–0.5
0.2
–1 0.2
0.4
0.6
0.8
1
⊡ Fig. -
Left: The B-spline β is a scaling function and operates as mean value. Right: The corresponding wavelet, the Haar wavelet, operates as a local difference operator
1
Splines and Multiresolution Analysis
1 0.8 0.6 0.4 0.2 2 a 1
4
6
8
Original step function 1 0.75 0.5 0.25
0.8 0.6 0.4 0.2 2
4
b 1
–0.25 –0.5 –0.75 –1
2
4
6
8
4
6
8
4
6
8
6 8 First iteration: approximation and details 1 0.75 0.5 0.25
0.8 0.6 0.4 0.2 2
4
c 1
–0.25 –0.5 –0.75 –1
2
6 8 Second iteration: approximation and details 1 0.75 0.5 0.25
0.8 0.6 0.4 0.2 2 d
4
–0.25 –0.5 –0.75 –1
2
6 8 Third iteration: approximation and details
⊡ Fig. - Multiresolution decomposition of the step function (a) with β as scaling function and with the Haar wavelet. The approximations (left column) are further iterated and decomposed into a coarser approximation and details (right column), until the coarsest approximation step, here the mean value, is reached. The sum of the coarsest approximation in the third iteration and of all details yields the original function (a). No information is lost in a multiresolution decomposition.
Splines and Multiresolution Analysis
.
Historical Background
The idea of piecewise polynomial functions and splines goes back to Schoenberg [, ]. In the s, when computer power started to be used for numerical procedures such as data fitting, interpolation, solving differential equations or computer aided geometric design, splines experienced an extreme upsurge. Schoenberg invented and strongly contributed to the concept of cardinal splines, which have equidistant nodes on the integers, see, e.g., [, , ] and many more. As a parallel development, in the s, the adaption of signal resolution to only process relevant details for a particular task evolved. For example, for computer vision, a multiresolution pyramid was introduced by Burt and Adelson []. It allowed to process an image first on a low resolution level and then selectively increase the resolution and add more detailed information when needed. The definition of a dyadic multiresolution analysis, i.e., A = Id, was contributed by Mallat [] and Meyer []. An interesting and in some parts historical collection on the most important articles in multiresolution and wavelet theory was assembled by Heil and Walnut []. The concepts of splines and multiresolution where joined by Lemarié [] and Battle [], when they showed that cardinal B-splines are scaling functions for multiresolution analyses. This led to many more developments of piecewise polynomial scaling functions for various settings and also multidimensions [], as, e.g., polyharmonic B-splines [, ] and other functions inspired from radial basis functions []. In , S. Mallat published his famous algorithm for multiresolution and wavelet analysis []. He had developed an efficient numerical method such that multiresolution decompositions could be calculated in a fast way. For the splines, M. Unser et al. proposed a fast implementation [, , ] which strongly contributed to the breakthrough of splines for signal and image analysis. In the last years, periodic, fractional, and complex versions of splines for multiresolution were developed, e.g., [, , , , ]. Many of them use a Fourier domain filter algorithm which allows for infinite impulse response filters. The former important feature of compact support of the cardinal B-splines and other functions is no longer a limiting criterion. Therefore, it can be expected that many new contributions on splines will still be made in future by modeling signal and image features in Fourier domain.
.
Mathematical Modeling and Application
..
Mathematical Foundations
...
Regularity and Decay Under the Fourier Transform
An important idea behind splines and multiresolution is the relation between regularity in time domain and decay in frequency domain, respectively between decay in time domain and regularity in frequency domain. To illustrate this, the notion of the Schwartz space is very useful [, , , ].
Splines and Multiresolution Analysis
Definition The subspace of functions f ∈ C ∞ (Rn ) with sup sup ( + ∥x∥ ) N ∣D α f (x)∣ < ∞
∣α∣≤N x∈R n
for all N = , , , . . .
is called the space of rapidly decreasing functions or Schwartz space S(Rn ). The norms induce a Fréchet space topology, i.e., the space S is complete and metrizable. α
αn
Here, D α = ( ∂x∂ ) ⋯ ( ∂x∂ n ) for every multi-index α = (α , . . . , α n ) ∈ Non . The dual space S ′ (Rn ), endowed with the weak-∗-topology, is called space of tempered distributions. The following famous linear transform relates the view points of the space domain and of the frequency domain: Definition The Fourier transform, defined by F f (ω) := ̂ f (ω) := ∫
Rn
f (x)e −i⟨ω,x⟩ dx,
ω ∈ Rn ,
is a topological isomorphism on L (Rn ) and on S(Rn ). Its inverse is given by f (x) =
i⟨ω,x⟩ ̂ dω ∫ f (ω)e (π)n Rn
in L (Rn ) resp. in S(Rn ).
The Fourier transform can be extended to the space of tempered distributions. For T ∈ S ′ (Rn ) the Fourier transform is defined in a weak sense as ̂ ̂) for all φ ∈ S(Rn ). F T(φ) := T(φ) := T(φ Also on S′ (Rn ), the Fourier transform is a topological isomorphism. The Fourier transform has the nice property to relate polynomials and differential operators. Theorem (i)
Let f ∈ S(R). Then for all k ∈ N F ( f (k) ) (ω) = (iω) k ̂ f (ω), and
(ii)
(iii)
̂ f (k) (ω) = F ((−i ●) k f ) (ω).
Let P be an algebraic polynomial in Rn , say P(x) = ∑α c α x α . . . x nα n , and let f ∈ S(Rn ). Then ̂f = P(iD) ̂ f and P f, F (P ( D) f ) = P ̂ i where P(iD) = ∑α c α i ∣α∣ D α . Part (ii) also holds for f ∈ S ′ (Rn ).
Splines and Multiresolution Analysis
Example The Fourier transform of the polynomial x k is the tempered distribution k i k ddx k δ, k ∈ N . For the construction of a multiresolution analysis, the scaling function can be used as a starting point. The idea is to choose a scaling function of a certain regularity, such that the generated multiresolution analysis inherits the smoothness properties. In particular, for the splines the idea is to model the regularity via decay in Fourier domain. The following theorem gives a motivation for this. The result can be deduced from the considerations above, and the fact that S(Rn ) is dense in L (Rn ): Theorem
Let f ∈ L (Rn ) and its Fourier transform decay as ∣̂ f (ω)∣ ≤ C( + ∥ω∥)−N−ε
for some ε > . Then all partial derivatives of order ≤ N − n are continuous and in L (Rn ). These results allow to construct scaling function with explicit regularity and decay properties, in space and in frequency domain. However, some criteria are needed to verify that the constructed function generates a multiresolution analysis.
...
Criteria for Riesz Sequences and Multiresolution Analyses
The following is an explicit criterion to verify whether some function φ is a scaling function. Theorem Let A be a dilation matrix and let φ ∈ L (Rn ) be some function satisfying the following properties: (i) (ii)
{φ(● − k)} k∈Zn is a Riesz sequence in L (Rn ). φ satisfies a scaling relation. I.e., there is a sequence of coefficients (a k ) k∈Zn such that φ(A− x) = ∑ a k φ(x + k)
in L (Rn ).
(.)
k∈Z n
(iii)
̂∣ is continuous at and φ ̂() ≠ . ∣φ
Then the spaces Vj = span {φ(A j ● −k)} k∈Zn ,
j ∈ Z,
form a multiresolution analysis of L (Rn ) with respect to the dilation matrix A. Proof
See, e.g., [, Theorem .] for the D case.
∎
For particular applications, the Riesz basis property (i) of {φ(● − k)} k∈Zn in V is not enough, but an orthonormal basis is needed. An example for such an application is the
Splines and Multiresolution Analysis
denoising of signals contaminated with Gaussian white noise [, Chap. X, Sect. ..]. However, there is an elegant mathematical method to orthonormalize Riesz bases generated by shifts of a single function. Theorem (i)
Let φ ∈ L (Rn ). Then the following holds:
{φ(● − k)} k∈Zn is a Riesz sequence in L (Rn ) if and only if there are constants c and C, such that ̂(ω + πk)∣ ≤ C < ∞ almost everywhere. < c ≤ ∑ ∣φ k∈Z n
(ii)
̂(ω + πk)∣ is strictly positive and I.e., the autocorrelation filter M(ω) := ∑ k∈Zn ∣φ bounded from above. {φ(● − k)} k∈Zn is an orthonormal sequence if and only if ̂(ω + πk)∣ = almost everywhere. ∑ ∣φ k∈Z n
(iii)
If {φ(● − k)} k∈Zn is Riesz basis of a subspace X of L (Rn ), then there exists a function Φ ∈ L (Rn ), namely ̂(ω) φ ̂ Φ(ω) =√ ̂(ω + πk)∣ ∑ k∈Zn ∣φ
(.)
such that {Φ(● − k)} k∈Zn is an orthonormal basis of X. Proof
∎
See, e.g., [] and [, Chap. VII].
Due to this theorem, every scaling function can be orthonormalized. Let φ ∈ L (Rn ) be some scaling function that generates a multiresolution analysis {Vj } j∈Z of L (Rn ). Then the family {Φ j,k } k∈Zn with Φ j,k (x) = − j/ Φ(− j (x − k)), and Φ as defined in (> .) is an orthonormal basis of the space Vj , j ∈ Z. Example A simple possibility to construct a dyadic multiresolution analysis in L (Rn ) is the tensor product approach. Let (Vj ) j∈Z be a dyadic (i.e., A = ) multiresolution analysis of L (R) with scaling function φ. Then (V j ) j∈Z with V j = Vj ⊗ . . . ⊗ Vj together with
n -times
the scaling function φ(x , . . . , x n ) = φ(x ) ⋅ . . . ⋅ φ(x n ) forms a multiresolution analysis of L (Rn ) and dilation matrix Id. In the same way, the scaling function φ(x , . . . , x n ) = φ (x ) ⋅ . . . ⋅ φ n (x n ) generates a multiresolution analysis of L (Rn ), if every φ k , k = , . . . , n, is a scaling function of some D multiresolution analysis with dilation factor a ∈ N/{}.
Splines and Multiresolution Analysis
...
Regularity of Multiresolution Analysis
In signal and image analysis the choice of an appropriate analysis basis is crucial. Here, appropriate means that the features of the basis such as smoothness should be in accordance with the properties of the functions to analyze. For example, analyzing a smooth signal or image with a fractal basis in general yields results that are difficult to interpret and to work with in practice. In this case, the signal resp. the image model does not match the model of the basis. The next section will show that the family of spline bases helps to avoid such difficulties, because the splines allow for a good adjustment due to their regularity parameter m, cf. (> .) and (> .). The following definition specifies the term “regular.” (See [].) Definition Denote C r the class of r-times continuously differentiable functions in Rn , C the class of continuous functions, and C − the class of measurable functions. (i)
Let r = −, , , . . . . A function f : Rn → C is called r-regular, if f ∈ C r and ∣
(ii)
∂α Ak f (x)∣ ≤ ∂x α ( + ∥x∥) k
for every k ∈ N , every multi-index α with ∣α∣ ≤ max(r, ) and constants A k . A multiresolution analysis of L (Rn ) is called r-regular if it generated by an r-regular scaling function.
It is important to note that the orthonormalization procedure (> .) does not affect the regularity of the corresponding basis. For the orthonormalized scaling function Φ of a multiresolution analysis, the same regularity properties hold. Proposition Let φ ∈ L (Rn ) be an r-regular function, such that {φ(● − k)} k∈Zn forms a Riesz sequence. Then the via (> .) orthonormalized function Φ is also r-regular [].
...
Order of Approximation
Having found a scaling function that generates a multiresolution analysis, how good do the corresponding approximation spaces Vj approximate some function f ∈ L (Rn ) of a certain regularity? Let H k (Rn ) = { f ∈ L (Rn ) : ∥ f ∥H k :=
f ∥L < ∞}, ∥( + ∥●∥Rn ) k ̂ (π)n
k ∈ N ,
denote the Sobolev spaces. The following criterion for the order of approximation turns out to be easy to verify for splines.
Splines and Multiresolution Analysis
Theorem Let φ ∈ L (Rn ) satisfy the following properties [, Theorem .]: (i) (ii) (iii)
̂ is bounded on some neighborhood of the origin. /φ Let B ε be some open ball centered at the origin and let E := B ε + (πZn /{}). For ̂ of order ≤ α are in L (E). some α > k + n/, all derivatives of φ ̂(ω) = for all ∣γ∣ < k and all ω ∈ πZd /{}. Dγ φ
Then V = span {φ(● − k)} k∈Zn provides approximation order k: For f ∈ H k (Rn ), min{∥ f − s(●/h)∥L , s ∈ V } ≤ const. h k ∥ f ∥H k
...
for all h > .
Wavelets
For the step from a coarser approximation space Vj to a finder one Vj+ information has to be added. It is contained in the wavelet or detail space Wj , which is the orthonormal complement of Vj in Vj+ : Vj+ = Vj ⊕ Wj . It follows that Vj+m = Vj ⊕ ⊕ k− l= Wj+l , and hence L (Rn ) = ⊕ Wj
(.)
j∈Z
can be decomposed in a direct sum of mutually orthogonal detail spaces. Moreover, the detail spaces Wj inherit the scaling property from Definition (iv) for the approximation spaces Vj . For all j ∈ Z, f ∈ Wj ⇐⇒ f (A− j ●) ∈ W . The question now is whether there is also a simple basis generated by the shifts of one or few functions, the wavelets. The following definition is motivated from > Eq. (.). Definition Let A be a dilation matrix, and let {ψ l } l=, ... ,s be a set of functions in L (Rn ), such that the family {∣ det A∣ j/ ψ l (A j ● −k) ∣ l = , . . . , s, j ∈ Z, k ∈ Zn } forms an orthonormal basis of L (Rn ). Then {ψ l } l=, ... ,s is called wavelet set associated with A. What qualitative properties do the wavelets have? The approximation spaces Vj are generated by the scaling function, which operates as a low pass filter. This can be seen from ̂() ≠ , resp. from Theorem (i): /φ ̂ is bounded in some neighborhood Theorem (iii) φ of the origin. Therefore, the added details and thus the wavelets have to carry the highfrequency information. In addition, the wavelets ψ in W are elements of V and therefore have the form (.) ψ(A− x) = ∑ a k φ(x − k) k∈Z n
Splines and Multiresolution Analysis
in L -norm, where {a k } k∈Zn are the Fourier coefficients of a certain πZn -periodic function. Proposition Let (Vj ) j∈Z be a multiresolution analysis of L (Rn ) with respect to the dilation matrix A and with scaling function φ. Then for a function f ∈ V if and only if ̂ ̂(ω) almost everywhere. f (AT ω) = m f (ω)φ Here, m f ∈ L ([ , π]n ) and ∥m f ∥L ([ ,π]n ) =
∥ f ∥L (Rn ) . ∣ det A∣
For a proof, see, e.g., [, ]. Note that for a wavelet ψ as in mψ (ω) =
>
Eq. (.) there holds
∑ a k e i⟨ω,k⟩ . ∣ det A∣ k∈Zn
How many wavelets, i.e., generators of W , are needed to span the space? The parameter s in Definition is yet unspecified. In fact, s depends on the scaling matrix. A leaves the lattice Zn invariant, AZn ⊂ Zn . The number of cosets is ∣ det A∣ = ∣Zn /AZn ∣ (see [, Proposition .]). It turns out that q = ∣ det A∣ − wavelets are needed to generate the space W . To motivate this, for a start, let f ∈ V be an arbitrary function. Denote γ , . . . , γ q representatives of the q + cosets of AZn in Zn . Then each coset can be written as γ m + AZn , m = , . . . , q. The function f has the representation f (A− x) = ∑ c k ( f )φ(x − k), ∣ det A∣/ k∈Z n
(.)
or in Fourier domain ̂ f (AT ω) =
̂(ω), c f (ω)φ ∣ det A∣/
(.)
in L -sense and with an appropriate πZn -periodic function c f (ω) with Fourier coefficients (c k ( f )) k∈Zn . Then c f (ω) can be decomposed with respect to the cosets: q
c f (ω) = ∑ c k ( f )e i⟨ω,k⟩ = ∑
∑
m= k∈γ m +AZ n
k∈Z n
c k ( f )e i⟨ω,k⟩
q
q
= ∑ e i⟨ω,γ m ⟩ ∑ c k+γ m ( f )e i⟨ω,k⟩ = ∑ c m f (ω), m=
m=
k∈AZ n
where i⟨ω,γ m ⟩ i⟨ω,k⟩ = e i⟨ω,γ m ⟩ ∑ c Ak+γ m ( f )e i⟨ω, Ak⟩ cm ∑ c k+γ m ( f )e f (ω) = e k∈AZ n
=e
i⟨ω,γ m ⟩
∑ c Ak+γ m ( f )e k∈Z n
k∈Z n i⟨AT ω,k⟩
=
T e i⟨ω,γ m ⟩ κ m f (A ω).
Splines and Multiresolution Analysis
This representation exists for all functions V , in particular for φ and the wavelets. The following theorem indicates, how many wavelets are needed to generate the space W , such that W ⊕ V = V . Theorem Let φ ∈ V be a scaling function and let ψ , . . . , ψ q ∈ V . Then the family {φ(● − k)} k∈Zn is an orthonormal system if and only if q
m ∑ ∣κ φ (ω)∣ = almost everywhere.
(.)
m= q
The system {φ(● − k)} k∈Zn ∪ ⋃m= {ψ m (● − k)} k∈Zn is an orthonormal basis in V if and only if the so-called polyphase matrix ⎛κ φ (ω) ⎜ ⋮ ⎜ ⎝κ φq (ω)
κψ (ω) ⋮ q κψ (ω)
⋯ ⋱ ⋯
κψ q (ω)⎞ ⋮ ⎟ ⎟ q κψ q (ω)⎠
is unitary for almost all ω ∈ Rn . The proof for a more general version of this theorem is given in [, Sect. .]. A summary and a condition for r-regular wavelets yields in the following theorem. Theorem Consider a multiresolution analysis on Rn associated with a dilation matrix A. (i) (ii)
Then there exists an associated wavelet set consisting of q = ∣ det A∣ − functions. If the multiresolution analysis is r-regular, and in addition q + > n, then there exists an associated wavelet set consisting of q r-regular functions.
The idea of the proof is that for an r-regular function φ on Rn and a πZn -periodic C ∞ ̂(ω) = η(ω)φ ̂(ω) is an r-regular function. function η(ω) the convolution ψ defined by ψ For an explicit proof, see again []. Example As a continuation of Example , the wavelet function corresponding to φ = χ[,) is derived. To this end, consider the space L (R) and the dilation A = . Then q = det A − = ; thus γ = and γ = are representatives of the cosets of A. That is, there is only a single wavelet needed to generate W . > Equation (.) yields x √ φ ( ) = √ φ(x) + √ φ(x − ) for the normalized generator of V . Thus c (φ) = κ φ ( ω )
=
κ φ ( ω )
=
√
. Then by
>
√
and c (φ) =
√
e i ω . This implies
Eq. (.) of Theorem the family {φ(● − k)} k∈Z is
Splines and Multiresolution Analysis
orthogonal, since ∣κ φ (ω)∣ + ∣κ φ (ω)∣ = . The polyphase matrix (
κ φ (ω) κ φ (ω)
κψ (ω) ) κψ (ω)
can be completed to a unitary matrix by choosing κψ (ω) = wavelet then has the representation
√
= −κψ (ω). The corresponding
x √ ψ ( ) = √ φ(x) − √ φ(x − ), corresponding to (> .). This yields the Haar wavelet ψ as illustrated in
..
>
Fig. -.
B-Splines
Several of the criteria for scaling functions and multiresolution analyses given in the previous section are based on the Fourier representation of the scaling function, e.g., the Riesz sequence criterion and the orthonormalization trick in Theorem , as well as the criterion for the order of approximation in Theorem . For this reason, the modeling of a scaling function in Fourier domain to achieve certain specific properties is promising. Aiming at constructing a scaling function φ ∈ L (R) of regularity r = −, , , . . ., this property is considered in Fourier domain: It is a decay property of the Fourier transform ̂ (compare with > Sect. ...): φ ̂(ω) = O ( φ
) ∥ω∥r+
for ∥ω∥ → ∞.
Taking into account Theorem a first model for the scaling function in Fourier domain is ̂(ω) = φ where the function relation (> .)
(ω) , ω r+
ω ∈ R,
(.)
still has to be specified. Since scaling functions satisfy a scaling x φ ( ) = ∑ h k φ(x − k) in L (R), k∈Z
the Fourier transform of this equation yields ̂(ω) = H(ω)φ ̂(ω), φ where (h k ) k∈Z is the sequence of Fourier coefficients of the π-periodic function H. For the ansatz (> .), H(ω) =
̂ (ω) φ(ω) (ω) ω r+ = = . ̂(ω) φ (ω)r+ (ω) r+ (ω)
(.)
Splines and Multiresolution Analysis
This gives criteria for the choice of the function : (i) (ii) (iii)
vanishes at the origin and there has a zero of order r+. This ensures that φ̂ ∈ L (R) and that Theorem (iii) is satisfied. (ω) is a trigonometric function, to ensure that H(ω), the so-called scaling filter, is π-periodic. has no other zeros in [ −π, π ], except at the origin. Otherwise, the autocorrelation ̂(ω + πk)∣ would vanish somewhere, and the shifts of the filter A(ω) = ∑ k∈Z ∣φ function φ would fail to generate a Riesz sequence, see Theorem (i).
A simple function ensuring all three requirements (i), (ii), and (iii) is r+
(ω) = (sin(ω/)θ(ω/))
,
where θ is a π-periodic phase factor such that ∣θ∣ = , i.e., a shift in time domain. Choosing θ(ω) = e −i ω yields the cardinal B-splines as given in (> .) resp. (> .): β̂ (ω) = ∫
e −i ωt dt =
− e −i ω sin(ω/) −i ω/ = e . iω ω/
Since β has compact support, β̂ ∈ C ∞ []. Due to the convolution formula (> .), −i ω
−e β̂m (ω) = ( iω
m+
)
sin(ω/) −i ω/ e =( ) ω/
m+
.
(.)
The β m are scaling functions of regularity r = m − , as the verification of the criteria in Theorem shows. In fact, the following holds. Let A = . Integrability: Since by (> .) the functions β̂m are L -integrable, so are the β m , m ∈ N . Riesz sequence property: The shifted characteristic functions β (x − k) = χ[,) (x − k), k ∈ Z, are clearly orthonormal. Theorem (ii) thus yields ∑ ∣ β̂ (ω + πk)∣ = almost everywhere. k∈Z
To verify the Riesz sequence property for β m , the autocorrelation filter must be bounded with strictly positive constants from above and from below. It is m m+ . ∑ ∣ β̂ (ω + πk)∣ = ∑ ∣ β̂ (ω + πk)∣ k∈Z
k∈Z
In [ −π, π ], ∣β̂ ∣ is clearly positive (cf. > Fig. -), which gives ∣β̂ (π)∣ = /π as a positive bound from below. There is a constant c, such that < c = (/π)m+ < ∣β̂ (ω)∣m+ ≤ ∑ ∣ β̂m (ω + πk)∣m+ . k∈Z
Splines and Multiresolution Analysis
1 0.8 0.6 0.4 0.2 p
–4p –3p –2p –p
2p 3p 4p
⊡ Fig. -
∧ The function ∣β ∣ is strictly positive in the interval [ −π, π ]
Since the sequence (∣ β̂ (ω + πk)∣) k∈Z ∈ l (Z) for all ω ∈ R, the same is true for the sequence ∣ β̂m (ω + πk)∣ = ∣β̂ (ω + πk)∣m+ . This yields the existence of the requested upper bound c < ∞. Thus { β̂m (● − k)} k∈Z forms a Riesz sequence in L (R). Scaling relation: The scaling filter (> .) H(ω) = −m
m+ m + −i ωk ( − e −iω )m+ −m −i ω m+ −m )e = ( + e ) = ∑( k ( − e −i ω )m+ k=
))k∈Z ∈ l (Z). Hence, the is obviously π-periodic and has Fourier coefficients (−m (m+ k B-splines satisfy the scaling relation (> .) m+
β m (x/) = ∑ −m ( k=
m+ m )β (x + k). k
For β this equation reads β (x/) = β (x) + β (x + ), which is true since β (x/) = χ[,) (x/) = χ[,) (x). This equation and examples for scaling relations of other B-splines are illustrated in > Fig. -. ̂ at the origin: From > Eq. (.), Continuity and positivity of φ m (ω)∣ = ∣ ∣ β̂
sin(ω/) ∣ ω/
m+
,
m () = . Thus we have proved which has a continuous continuation at the origin, and β̂ the following conclusion:
Theorem The cardinal B-spline β m , m ∈ N , is a scaling function of an m − -regular multiresolution analysis with dilation . The order of approximation is m + .
Splines and Multiresolution Analysis
1
1
0.8
0.8
0.6
0.6
0.4
0.4
0.2
0.2
–1
1 m=0
2
3
0.7
–1
1
2 m=1
3
4
5
0.6
0.6
0.5
0.5
0.4
0.4
0.3
0.3 0.2
0.2
0.1
0.1 2
4
6
2
m=2
4 m=3
6
8
⊡ Fig. -
Scaling relation for B-splines βm , m = , . . . , . The B-spline versions βm (x/) ∈ V− are displayed with solid lines, the scaled translates in V are depicted dashed. The sum of the dashed functions gives the B-spline at the lower scale β m (x/)
Note that the cardinal B-splines β m with Fourier transform of the form (> .) are scaling functions, but they are not yet orthonormalized: The family {β m (● − k)} k∈Z spans V and is a Riesz basis, but it is not an orthonormal basis of V . Orthonormality can be achieved with Theorem and > Eq. (.): β̂m (ω) ̂m (ω) := √ B . ∑ k∈Z ∣ β̂m (ω + πk)∣ Figure - shows some orthonormalized B-spline scaling functions and the corresponding wavelets. >
..
Polyharmonic B-Splines
The same approach to model scaling functions in Fourier domain can be done in higher dimensions. We aim at constructing a scaling function for a multiresolution analysis of
Splines and Multiresolution Analysis
1.2 0.5
1 0.8 –4
0.6
–4
–2
–2
0.4
–0.5
0.2
–1 2
–0.2
4
2
4
2
4
2
4
–1.5 m=1
1
1
0.8
0.5
0.6 0.4
–4
–2
0.2 –4
–0.5
–2
2
4
–1
–0.2 m=2 1
0.5
0.8 0.6 –4
0.4
–2 –0.5
0.2 –4
–2
2
–1
4
–0.2 m=3
⊡ Fig. -
Orthonormalized B-splines Bm (left column) and corresponding wavelets (right column) for m = , , in time domain. Note that the orthonormalized B-splines and wavelets do not have compact support. Due to the orthonormalization procedure (> .) the orthonormalized B-spline is an infinite series of shifted compactly supported B-splines
L (Rn ) of the form ̂(ω) = φ
(ω) , ∥ω∥r
r ∈ N,
r > n/,
x ∈ Rn .
With an appropriate trigonometric polynomial n
r
(ω) = ( ∑ sin (ω k /)) , k=
ω = (ω , . . . , ω n ),
Splines and Multiresolution Analysis
̂ is a nonseparable scaling function for a multiresolution analysis of L (Rn ) with respect to φ dilation matrices A that are scaled rotations. The corresponding function in space domain φ is called elementary r-harmonic cardinal B-spline, or short polyharmonic B-spline P r . This terminology can be justified as follows. The Fourier transform in the sense of tempered distributions of the function /∥ω∥r is indeed a polynomial – up to a logarithmic factor for r − n even. In fact, in S ′ (Rn ), F − (/∥ ● ∥r )(x) = ∥x∥r−n (A(n, r) ln ∥x∥ + B(n, r)) =: ρ(x), with constants A(n, r), B(n, r) as given in [, Chap. , Sect. ], and A(n, r) = except for r − n even. (Note that for r > n/ on the right-hand side the final parts have to be considered.) The term polyharmonic comes from the fact that ρ is the Green function of the r-iterated Laplace operator Δr . However, with these considerations φ(x) = P r (x) = ∑
k ρ(x
+ k) almost everywhere
k∈Z
becomes an nD spline. Here ( k ) k∈Z is the sequence of Fourier coefficients of . Due to the decay in Fourier domain, the polyharmonic B-spline P r has continuous derivatives D β for multi-indices ∣β∣ < r − n. In the same way as for the B-splines, it can be shown with the theorems given in > Sect. .. that φ forms indeed a scaling function with approximation order r [, , ]. > Figure - shows the polyharmonic B-spline scaling function in space domain and in frequency domain.
2 1 1.5
0.8 0.6
1
0.4
0.5
0.2 0 –2
2 1
0 2
2
0 –2
2
0
0
–1 –2
–2
⊡ Fig. -
Polyharmonic B-spline for r = in space domain (left) and frequency domain (right)
Splines and Multiresolution Analysis
.
Survey on Mathematical Analysis of Methods
There are many function families that consist of piecewise polynomials and that are called splines, which in addition fulfil a multiresolution condition in the one or the other sense. These families can be classified by various aspects, e.g., by their dimensionality, by the lattice which is invariant under the corresponding dilation matrix, by the geometries they are defined on, or whether they provide phase information or not, and so on. The following sections list some of these spline approaches and illustrates their mathematical properties and features.
..
Schoenberg’s B-Splines for Image Analysis – the Tensor Product Approach
As mentioned in Example , multiresolution analyses for L (Rn ) and dilation matrix Id can be generated from tensor products of D dyadic multiresolution analyses. To analyze images with B-splines, the tensor product β m (x)β m (y) of B-splines is a scaling function for L (R ) and the dilation matrix A = Id. Since in D the determinant det A = , the corresponding details space W is spanned by three wavelet functions: ψ(x)β m (y),
β m (x)ψ(y),
ψ(x)ψ(y),
x, y ∈ R.
(.)
A drawback of this approach is the fact that these wavelets prefer horizontal, vertical, and diagonal directional features and are not sensitive to other directions, see > Fig. -. For the analysis of images with many isotropic features, the use of isotropic or steerable wavelets is recommended. However, the tensor approach is a simple and widely used wavelet approach. For an illustration of the respective image decomposition in coarse approximations and details of various sizes, see > Fig. -.
⊡ Fig. -
The three B-spline tensor wavelets (> .) show a preference for horizontal, vertical, and diagonal directions. Here, m = . Minimal resp. maximal function values are given in black resp. white
Splines and Multiresolution Analysis
a
b
c
d
⊡ Fig. - Decomposition of an image [, Part of IM.tif] into coarse approximations and details of various sizes. (a) Original image. (b) Matrix of the absolute values of the multiresolution coefficients. Large coefficients are white. The wavelet coefficients are depicted in the lower two and the upper right band of each scale; the approximation coefficients in the upper left band. (c) Two steps of themultiresolution decomposition. From left to right: Finest details and second finest details, remaining approximation of the original image. (d) Second finest details (c, center) split into the contribution of the three wavelets. From left to right: The decomposition into horizontal, vertical and diagonal details
..
Fractional and Complex B-Splines
The B-splines as described up to now have a discrete order of smoothness, i.e., they are C n -functions with n ∈ {−, , , , . . .}. For some applications, e.g., in medical imaging, where the order of smoothness of certain image classes is fractional, it would be favorable
Splines and Multiresolution Analysis
to have a spline and wavelet basis that is adaptable with respect to this regularity [, , ]. A first step in this direction was done by T. Blu and M. Unser, who proposed B-splines and wavelets of fractional orders []. They defined two variants of fractional B-splines, the causal ones and the symmetrical ones. The causal fractional B-spline is generated by applying the (α+)th fractional difference operator to the one-sided power function t+α : β+α (t) :=
α+ )(t − k)+α . Δ α+ t α = ∑ (−) k ( k Γ(α + ) + + Γ(α + ) k≥
The Fourier-transform representation is similar to the one of the classical B-splines (cf. > Eq. .): −i ω
−e β̂+α (ω) = ( iω
α+
)
.
Here again, the smoothness property β+α ∈ C m,γ (R) in time domain is gained by the fractional polynomial decay of order O(∣ω∣∣α+ ) in frequency domain. Note that C m,γ (R) denotes the Hölder space with exponent m = ⌊α + ⌋ and γ = α + − m, i.e., the space of m-times continuously differentiable functions f that are Hölder-regular with exponent < γ ≤ such that there is a constant C > with ∣D m f (t) − D m f (s)∣ ≤ C∣t − s∣ ∀s, t ∈ R. Although the fractional B-splines are not compactly supported, they decay in the order O(∣t∣−(α+) ) as t → ∞. They are elements of L (R) for α > , of L (R) for α > − and of the Sobolev spaces Wr (R) for r < α + . They share many properties with their classical B-spline relatives, such as the convolution property, their relation to difference operators, i.e., they are integral kernels for fractional difference operators [], and they are scaling functions for dyadic multiresolution analyses. This can be verified by the procedure given in > Sect. ... The causal fractional B-spline is not symmetric. Since for some signal and image analysis tasks symmetrical bases are preferred, in [] the symmetrical fractional B-splines β∗α are proposed. They are defined in Fourier domain as follows: −i ω
−e β̂∗α (ω) = ( iω
α+
)
α+
+ eiω sin(ω/) ) =∣ ∣ ( −iω ω/
α+
,
(.)
and therefore obviously are symmetrical in time domain. The same regularity and decay properties apply as for the causal fractional B-splines. The symmetrical fractional B-splines are also piecewise polynomials, as long as α ∉ N . For even integer degrees, the singularity introduced through the absolute value in > Eq. (.) causes that β∗m is a sum of integer shifts of the logarithmic term ∣t∣m ln(t) for m ∈ N . For the explicit time-domain representation and further details on these splines cf. []. In [], Blu and Unser defined another variant, the generalized fractional B-spline or (α, τ)-fractional spline β τα with a parameter τ ∈ R. Also these splines are defined via their
Splines and Multiresolution Analysis
Fourier domain representation: −i ω
−e β̂τα (ω) = ( iω
α+ +τ
)
α+
− eiω ) ( −iω
−τ
.
As above, the parameter α > controls the regularity of the splines. The parameter τ, in contrast, controls the position of the splines with respect to the grid Z. This can be justified by the following fact. All variants of the B-splines considered in this section converge to optimally time-frequency localized functions in the sense of Heisenberg, i.e., to Gaussians or Gabor functions, if the degree α becomes large. For a proof for the classical cardinal B-splines, see []. In the case of the (α, τ)-fractional splines [],
β τα (t) = O (e − α+ (t−τ) )
for α → ∞.
This explains the notion “shift parameter” for τ. Moreover, the parameter τ allows to interpolate the spline family between the two “knots,” the symmetrical ones (τ =) and ), see > Fig. -. Both parameters α and τ can be tuned the causal ones (τ = α+ independently and therefore allow for an individual adjustment of the analyzing basis. Another generalization are the complex B-splines []. There are two variants, both defined via their Fourier domain representation. Let z = α + iγ ∈ C, α > − , γ ∈ R and y ∈ C. Then −i ω
−e β̂z (ω) = ( iω
−i ω
−e β̂zy (ω) = ( iω
z+
) )
, z+ −y
− eiω ) ( −iω
z+ +y
are complex B-splines of complex degree z. The functions are well defined, because the −i ω function Ω(ω) = −ei ω never touches the negative real axis such that Ω(ω)z is uniquely 1
1
0.5
0.5
0
0
–0.5 –2
–1
0
1
2
3
4
–0.5 –2
–1
0
1
2
3
4
⊡ Fig. - The fractional (α, τ)-splines interpolate the families of the causal and the symmetric k for k = , , , , from the most right (causal) to the most left fractional splines. τ = α+ (symmetrical) function in each image. Right: α = .. Left: α =
Splines and Multiresolution Analysis
defined. β z and β zy are elements of the Sobolev spaces Wr (R) for r < α + . β z has the time-domain representation β z (t) =
k z+ )(t − k)+z , ∑ (−) ( k Γ(z + ) k≥
i.e., β z is a piecewise polynomial of complex degree. For more details on the properties of these families of complex splines cf. []. The idea behind the complex degree is as follows: The real part Re z = α operates as regularity and decay parameter in the same way as for the fractional B-splines. The imaginary part, however, causes an enhancement resp. damping of positive or negative frequencies. In fact, β̂z (ω) = β̂+α (ω)e −iγ ln ∣Ω(ω)∣ e γ arg Ω(ω) . The imaginary part γ of the complex degree introduces a phase and a scaling factor in frequency domain. The frequency components on the negative and positive real axis are enhanced with different sign, because arg Ω(ω) ≥ for negative ω and arg Ω(ω) ≤ for positive ω. > Figure - illustrates this effect. With real-valued functions, only symmetric spectra can be analyzed. The complex B-splines, however, allow for an approximate analysis of the positive or the negative frequencies, because the respective symmetric bands are damped due to the complex exponent. However, the complex B-splines inherit many properties of their classical and fractional relatives. All of the generalized B-spline families mentioned in this section have in common that they are scaling functions of dyadic multiresolution analyses. They are one-dimensional functions, but with the tensor approach mentioned in Example and > Sect. .. they are also suitable for image processing tasks. Although the fractional and the complex splines, in general, to not have compact support, they allow for fast analysis and synthesis algorithms. Due to their closed form in Fourier domain, they invite for an implementation of these algorithms in Fourier domain.
1 0.8
1.25
0.6
1
0.4 0.2
–15
–10
–5
⊡ Fig. -
6
1.5
5
10
15
–15
–10
–5
5 4
0.75
3
0.5
2
0.25
1 5
10
15
–15
–10
–5
5
10
∧ The frequency spectrum ∣β z ∣ for z = + iγ, γ = , , (from left to right). The spectrum of β = β+ is symmetric (right), whereas the spectra of β +i (center) and β +i (left) show an enhancement of the positive frequency axis
15
..
Splines and Multiresolution Analysis
Polyharmonic B-Splines and Variants
In > Sect. .. the so-called polyharmonic B-splines in Rn were introduced. They are defined in Fourier domain by the representation r
n ̂r (ω) = ( ∑ k= sin (ω k /) ) , P n ∑ k= ω k
r > n/, ω = (ω , . . . , ω n ).
These polyharmonic B-splines satisfy many properties of the classical Schoenberg splines; e.g. they are piecewise polynomial functions, they satisfy a convolution relation P r +r = P r ∗P r , they are positive functions, etc. However, they do not share the property that they converge to optimally space-frequency localized Gaussians as r increases []. This is due to the fact that the trigonometric polynomial in the numerator regularizes insufficiently at the origin: The second order derivative of ∑nk= sin (ω k /) n ∑ k= ω k is not continuous. Van De Ville et al. [] therefore proposed another localizing trigonometric polynomial: n
μ(ω) = ∑ sin (ω k /) − k=
n− n ∑ ∑ sin (ω k /) sin (ω l /). k= l=k+
(.)
A new function family then is defined in Fourier domain via r
̂r (ω) = ( μ(ω) ) , Q ∥ω∥
r>
n .
(.)
Qr is called isotropic polyharmonic B-spline. The function is piecewise polynomial (except for r − n even, where a logarithmic factor has to be added, see below), and shares with the polyharmonic splines their decay properties in Fourier domain and their regularity properties in space domain. Qr converges to a Gaussian as r increases, which makes the function family better suitable for image analysis than P r , because of the better spacefrequency localization. This effect is due to higher order rotation invariance or isotropy of the localizing trigonometric polynomial (> .): (ω) = + O(∥ω∥ ) vs.
μ(ω) = +
∥ω∥ + O(∥ω∥ ) as ∥ω∥ → .
̂r has a second order moment, and thus the central limit theorem can be This causes that Q ̂r has a higher applied to proof the convergence to the Gaussian function. In addition, Q ̂r ; therefore Qr decays faster. For a complete proof of the localization regularity than P property, see []. An example of the isotropic polyharmonic spline is given in > Fig. -. The polyharmonic B-splines and the isotropic polyharmonic B-splines both are realvalued functions. The isotropic B-spline is approximately rotation-invariant and therefore is suited for image analysis of isotropic features. A complex-valued variant of these B-splines in D was introduced in []. The idea is to design a spline scaling function that
Splines and Multiresolution Analysis
2 1.5 1 0.5 0 –2 0 2
2
0
–2
⊡ Fig. -
The isotropic polyharmonic spline Q . Compare with polyharmonic B-spline P
>
Fig. - of the classical
is approximately rotation covariant, instead of rotation invariant. Rotation covariant here means that the function intertwines with rotations up to a phase factor. Again, the design of the scaling function is done in Fourier domain, now using a perfectly rotation-covariant function ̂ ρ r,N (ω , ω ) =
(ω
+
r ω )
, (ω − iω ) N
where r ≥ and N ∈ N. In fact, for some rotation matrix R θ ∈ GL(, R), ̂ ρ r,N (R θ ω) = ρ r,N (ω). For localizing the function ̂ ρ r,N , the same trigonometric polynomials and e −i N θ ̂ μ as above can be used. The corresponding complex polyharmonic B-splines are then defined in frequency domain as N
̂ r,N (ω , ω ) = R
( (ω , ω ))r+ , r (ω + ω ) (ω − iω ) N
̂ r,N R μ (ω , ω ) =
(μ(ω , ω ))r+ . r (ω + ω ) (ω − iω ) N
or as
N
(.)
The case N = yields the real-valued polyharmonic splines. There are also other trigonometric polynomials that are suitable as localizing numerator for the real and the complex polyharmonic B-splines. With an appropriate choice the features of the polyharmonic splines can be tuned []. However, for both the real and the complex variant, the localizing multiplier has to fulfil moderate conditions to make the respective polyharmonic B-splines being a scaling function. In D, the following result holds (cf. []):
Theorem such that
Splines and Multiresolution Analysis
Let r > and N ∈ N . Let η(ω , ω ) be a bounded, πZ -periodic function, N
∣
(η(ω , ω ))r+ ∣ r (ω + ω ) (ω − iω ) N
is bounded in a neighborhood of the origin, and such that η(ω , ω ) ≠ for all (ω , ω ) ∈ [−π, π] /{(, )}. N ̂ = η r+ ⋅ ̂ ρ is the Fourier transform of a scaling function φ which generates a Then φ a b multiresolution analysis . . . V− ⊂ V ⊂ V . . . of L (R ) with dilation matrix A = ( ), −b a a, b ∈ Z: Vj = span {∣ det A∣ j/ φ(A j ● −k), k ∈ Z }. From the Fourier domain representation immediately follows that φ̂ ∈ L (R ) for r+ N > ̂ decays as ∣φ ̂(ω , ω )∣ = O(∥(ω , ω )∥−r−N ) when ∥(ω , ω )∥ → ∞. and that φ As a result, for all three variants of the polyharmonic B-splines, the classical ones P r , the isotropic ones Qr , and the complex ones Rr,N , the following properties hold: They are scaling functions for multiresolution analysis. Their smoothness parameter r can be chosen fractional and must fulfill r + N > n/ for integrability reasons. Then the scaling function in space domain is element of the Sobolov space φ ∈ Ws (R ) for all s < r + N − . The explicit space domain representation is φ(x) = ∑ η k ρ r,N (x + k) k∈Z N
for almost all x ∈ R . Here (η k ) k∈Z denotes the Fourier coefficients of η r+ . ρ r,N is the inverse Fourier transform of the Hadamard final part P f (̂ ρ r,N ) ∈ S ′ (R ). In fact, for r ∉ N , r− ρ r,N (x , x ) = c (x + x ) (x + ix ) N and for r ∈ N , ρ r,N (x , x ) = c (x + x )
r−
(x + ix ) N (ln π
√
x + x ) + c )
with appropriate constants c , c , c ∈ C. This justifies the notion spline for the function families. They all have a closed form in frequency domain. As in the case of the D cardinal B-splines there is a fast analysis and synthesis algorithm using the frequency domain representation and the fast Fourier transform, cf. > Sect. ..
..
Splines on Other Lattices
...
Splines on the Quincunx Lattice
The tensor product of two D dyadic multiresolution analyses yields a D multiresolution analysis with dilation matrix A = Id, cf. Example . As a consequence, the scaling factor while moving from one approximation space V to the next coarser space V is ∣ det A∣ = .
Splines and Multiresolution Analysis
6
6
6
5
5
5
4
4
4
3
3
3
2
2
2
1
1
1
a
1
2
3
4
5
6
b
1
2
3
4
5
6
c
1
2
3
4
5
6
⊡ Fig. -
Three iterations of the quincunx subsampling scheme. (a) Z , (b) Aq Z , (c) Aq Z . The thinning of the Z -lattice using the dilation matrix Aq is finer than dilation with the dyadic matrix A = Id, which in one step leads from (a) to (c)
For some image processing applications, especially in medical imaging, this scaling step size is too large. A step size of as in the D case would be preferred. Moreover, the decomposition of the wavelet space into three subspaces then would be avoided, and the eventual problematic of the directionality of the three involved wavelets would not arise. An example of a dilation matrix satisfying these requirement is the scaled rotation matrix Aq = ( −
)
with det A q = . It leads to the so-called quincunx lattice. This lattice is generated by applying A q to the cartesian lattice. It holds A q Z ⊂ Z , see > Fig. -. Since A q falls into the class of scaled rotations, the polyharmonic B-Spline construction including all variants are applicable for this case, cf. Theorem . Note that the tensor product approach in general is not suitable for the quincunx subsampling scheme.
...
Splines on the Hexagonal Lattice
Images as D objects are normally processed on the cartesian lattice, i.e., the image pixels are arranged on a rectangular grid. For image processing this arrangement has the drawback that not all neighbors of a pixel have the same relation: The centers of the diagonal neighbors have a larger distance to the center pixel than the adjacent ones. A higher degree of symmetry has the hexagonal lattice. It is therefore ideal for isotropic feature representation. The hexagonal lattice gives an optimal tesselation of R in the sense of the classical honeycomb conjecture, which says that any partition of the plane into regions of equal area has a perimeter at least that of regular hexagonal tiling []. The better isotropy of the hexagonal lattice is attractive for image analysis and has led to a series of articles on image processing methods (e.g., [, , ]) as well as on applications (e.g., [, , , ]).
Splines and Multiresolution Analysis
The hexagonal lattice is generated by applying the matrix √ Rh =
√ (
√/ ) /
on the cartesian lattice Λ h = R h Z . A scaling function of a multiresolution analysis of L (R ) in the hexagonal lattice fulfils all properties of Definition , but the last two. They change to (v) (vi)
f ∈ V ⇐⇒ f (● − R h k) ∈ V for all k ∈ Z . There exists a scaling function φ ∈ V , such that the family {φ(● − R h k)} k∈Z of translates of φ forms a Riesz basis of V .
Let A be a dilation matrix which leaves the hexagonal lattice invariant AΛ h ⊂ Λ h . Then A is of the form [] A = R h BR − h with B ∈ GL(, R) having only integer entries and with eigenvalues strictly larger than one. > Fig. - gives an example of two subsampling steps on the hexagonal lattice. There are several possible approaches to define spline functions on the hexagonal lattice. Sablonnière and Sbibih [] proposed to convolve piecewise linear pyramids to generate higher order B-splines. Van De Ville et al. [] started with the characteristic function of one hexagon and also used an iterative convolution procedure to construct B-splines of higher degree. However, both approaches lead to discrete order hexagonal Bsplines. If A is a scaled rotation, then fractional and complex B-splines on the hexagonal lattice can be defined in an analog way as in > Sect. .. for polyharmonic splines and
6
6
6
5
5
5
4
4
4
3
3
3
2
2
2
1
1
1
a
1
2
3
4
5
6
b
1
2
3
4
5
6
c
1
2
3
4
5
⊡ Fig. - Three iterations of subsampling of the hexagonal grid with A = Rh BR− h , where B = (− ). (a) Λh , (b) AΛh = Rh BR− h Λh = Rh BZ , (c) A Λh
6
Splines and Multiresolution Analysis
their (complex) variants []. Consider again the perfectly rotation-covariant (or for N = rotation-invariant) function ̂ ρ r,N (ω , ω ) =
, r (ω + ω ) (ω − iω ) N
where r > and N ∈ N . The idea is now to use a hexagonal-periodic trigonometric polynomial for localizing this function and to eliminate the singularity at the origin. Condat et al. [] proposed √ √ η h (ω , ω ) = √ ( − (cos (/ (−ω / + ω ) / ) √ √ √ + cos (/ (ω / + ω ) / ) + cos (−/ ω ))) , and defined the elementary polyharmonic hexagonal rotation-covariant B-spline via its frequency domain representation as ̂ r,N (ω , ω ) = R h
N
(η h (ω , ω ))r+ . r (ω + ω ) (ω − iω ) N
The B-spline in space domain then has the representation Rr,N h (x) = ∑ η h,k ρ r,N (x − R h k). k∈Z
r+ N
Here, (η h,k ) k∈Z denotes the sequence of Fourier coefficients of η h . > Figure - shows the localizing trigonometric polynomial η and the Fourier spectra of two hexagonal splines. For N = , the functions are the elementary polyharmonic hexagonal rotation-invariant B-splines. They are real-valued functions that converge to a Gaussian, as r → ∞, and therefore are well localized in the space domain as well as in the frequency domain. For N ∈ N, the splines are complex-valued functions and approximately rotation covariant in a neighborhood of the origin: ̂ r,N (ω) = e i N arg(ω) ( + C∥ω∥ + O(∥ω∥ )) R h
for ω → ,
where ω = (ω , ω ) and C ∈ R a constant. The translates of the complex B-spline Rr,N h form a Riesz basis of the approximation spaces j Vj = span {Rr,N h (A x − R h k) : k ∈ Z }
L (R )
,
j ∈ Z.
The ladder of spaces (Vj ) j∈Z generates a multiresolution analysis of L (R ) for the hexagonal grid and for scaled rotations A. Also in this case, the implementation of the analysis and synthesis algorithm can be elegantly performed in frequency domain [, ].
Splines and Multiresolution Analysis
4
10
2 0 –10
5 0 –5 –5
0 5
a
10
1 0.75 0.5 0.25 0 –10
10 5 0 –5 5
b
1 0.75 0.5 0.25 0 –10
10
10 5 0 –5
–5
0
–10
–5
0 5
–10
c
10
–10
⊡ Fig. -
∧ Localization of ρ with (a) the hexagonal-periodic trigonometric polynomial ηh yields the elementary polyharmonic hexagonal rotation-covariant B-spline. Frequency spectrum of ∧ r,N ∧ r+N/, ), (b) for r + N = , and (c) for r + N = . Rh (or equivalently of Rh
.
Numerical Methods
For the illustration of the numerical method, we focus on the quincunx dilation matrix A=( −
),
(.)
and consider the polyharmonic spline variants in D as defined in > Sect. ... Since det A = , the generators of the multiresolution space and the corresponding wavelet space are two functions: The scaling function (here the variant of the polyharmonic spline), and the associated wavelet. We consider the scaling function φ(x) = Qrμ (x) resp. φ(x) = Rr,N μ (x). It spans the ladder of nested approximation spaces {Vj } j∈Z via Vj = span {∣ det A∣ j/ φ(A j ● −k), k ∈ Z }
L (R )
,
j ∈ Z.
Splines and Multiresolution Analysis
Denote ̂(ω + πk)∣ M(ω) := ∑ ∣φ
(.)
k∈Z
the autocorrelation filter. It is bounded < M(ω) ≤ C for some positive constant C []. The scaling functions can be orthonormalized applying the procedure given in Theorem (iii): ̂(ω) φ ̂ . Φ(ω) =√ M(ω) The scaled shifts of Φ span the same spaces {Vj } j∈Z . The B-splines as scaling functions φ satisfy a refinement relation φ(A− x) = ∑ h k φ(x − k) almost everywhere and in L (R ). k∈Z
This relation in fact is a discrete convolution. The Fourier transform yields a πZ -periodic function H ∈ L (T ) of the form H(e i ω ) = ∣ det A∣ ⋅ =
̂(AT ω) ρ̂r,N (AT ω)ζ(AT ω) φ = ∣ det A∣ ⋅ ̂(ω) ̂ φ ρ r,N (ω)ζ(ω)
ζ(AT ω) ⋅ , r− ζ(ω) (a + b ) (a − ib) N
ω ∈ R .
N
Here, ζ(ω) = (μ(ω))r+ is the localizing multiplier for the isotropic polyharmonic Bspline in the case N = , and for the rotation-covariant polyharmonic B-splines in the case N ∈ N, cf. (> .), (> .), and (> .). The wavelet function ψ spanning a Riesz basis for the orthogonal complement Wj = span { j/ ψ(A j ● −k), k ∈ Z }
L (R )
in Vj+ = Wj ⊕Vj can also be gained in frequency domain. For the quincunx dilation matrix A as in (> .) a wavelet (or sometimes called a prewavelet, since the functions are not yet orthonormalized) is given by ̂(ω) = G(e i ω )φ ̂(ω) = e −i ω H(ω + (π, π)T )M(ω + (π, π)T )φ ̂(ω), ψ compare with family
>
Sect. .... The (pre-)wavelet Riesz basis for Wj is then given by the {ψ j,k = j/ ψ(A j ● −k), k ∈ Z }.
This basis in general is not orthonormal: ⟨ψ j,k , ψ j,l ⟩ ≠ δ k,l . A function f ∈ L (R ) can then be represented by the series f =
̃j,k ⟩ψ j,k = ∑ ⟨f,ψ j∈Z,k∈Z
̃j,k ⟩ψ j,k , ∑ ⟨f,ψ j∈Z,k∈Z
̃ j,k } k∈Z denotes the dual basis for each j ∈ Z: ⟨ψ ̃ j,k , ψ j,l ⟩ = δ k,l . Its generator in where {ψ frequency domain is
Splines and Multiresolution Analysis
̂(ω) M(ω + (π, π)T ) φ ̂ ̃(ω) = e −i ω H(ω + (π, π)T ) . ψ T M(A ω) M(ω) In contrast, the formula M N M(ω + (π, π)T ) T ̂ ̂(ω) Ψ(ω) =N ψ M(AT ω) generates an orthonormal wavelet basis. It corresponds to the orthonormal basis of V generated by the integer shifts of the orthonormalized scaling function Φ. These considerations show that there are three variants of a multiresolution implementation: An “orthonormal” one with respect to the orthonormalized scaling functions and corresponding orthonormal wavelets, one with the B-splines on the analysis side, ̃ j,k ̃ j,k + ∑ ⟨ f , ψ j,k ⟩ψ f = ∑ ⟨ f , φ j,k ⟩φ k∈Z
for f ∈ Vj+ ,
k∈Z
and finally one with the B-splines on the synthesis side: ̃j,k ⟩ψ j,k ̃ j,k ⟩φ j,k + ∑ ⟨ f , ψ f = ∑ ⟨f, φ k∈Z
for f ∈ Vj+ .
k∈Z
Both, the scaling filters H(e i ω ) and the wavelet filters G(e i ω ) = e −i ω H(ω + (π, π)T )M(ω + (π, π)T ) as well as their orthogonal and dual variants in our case are nonseparable and infinitely supported. Therefore, a spatial implementation of the decomposition would lead to truncation errors due to the necessary restriction to a finite number of samples. However, because of the closed form of H and therefore of G, the corresponding multiresolution decomposition or wavelet transform can be efficiently implemented in frequency domain. The respective image first undergoes an FFT, then is filtered in frequency domain by multiplication with the scaling filter H and the wavelet filter G. This method automatically imposes periodic boundary conditions. The coefficients resulting from the high pass filtering with G are the detail coefficients. They are stored, whereas the coefficients resulting from the low pass filtering H are reconsidered for the next iteration step. ̃ j,k . ̃ j+,k = ∑ ⟨ f , φ j,k ⟩φ ̃ j,k + ∑ ⟨ f , ψ j,k ⟩ψ ∑ ⟨ f , φ j+,k ⟩φ k∈Z
k∈Z
k∈Z
For details and tricks of the frequency domain implementation, cf. [, , , ]. > Figure - shows the multiresolution decomposition for the scaling function φ = R, μ . There it was assumed that the image is bandlimited and projected on the space V , which has the advantage that the coefficients do not depend on the chosen flavour of the scaling function, i.e., orthogonal, B-spline or dual. Qualitatively the transform is very similar to a multiscale gradient with the real part corresponding to the x -derivative and the imaginary part corresponding to the x -derivative [].
Splines and Multiresolution Analysis
a
b
c
d
⊡ Fig. -
Decomposition of an image [, Part of IM.tif] into approximation and wavelet coefficients. (a) Original image. (b) Matrix of the absolute values of the multiresolution coefficients. Large coefficients are white. The approximation coefficients in the upper left band, the other bands are wavelet coefficients on six scales. (c) Real part of the coefficient matrix, (d) imaginary part of the decomposition matrix for φ = R, μ and the corresponding wavelets. The coefficients had their intensity rescaled for better contrast
.
Open Questions
In this chapter a method for the construction of spline multiresolution bases was described. It yields a nice variety of new bases with several parameters for adaption and tuning. In the last decade, the notion of compressive sampling or compressed sensing arose, which is footing on the existence of well-adaptable bases. In fact, the idea behind compressed sensing is that certain functions have a sparse representation, if the underlying basis is smartly chosen. In this case, the function can be reconstructed from very few samples because of
Splines and Multiresolution Analysis
the prior knowledge of sparsity in this underlying basis. As a consequence, the knowledge on sparsity allows to sample such a signal at a rate significantly under the Nyquist rate. (The Shannon–Nyquist sampling theorem says that a signal must be sampled at least two times faster than the signal’s bandwidth to avoid loss of information.) In the last years, and virtually explosively in the last years, many important theoretical results were proven in this field, in particular by D. Donoho, E. Candès, J. Romberg, and T. Tao. For an introduction and references on compressed sensing, see e.g., [, ] and the website []. Compressed sensing is based upon two fundamental concepts: that of incoherence and that of sparsity. Let {x i } i=, ... ,N be an orthonormal basis of the vector space V. Let f = N ∑ i= s i x i with s i = ⟨ f , x i ⟩. The signal f is called k-sparse, if only k of the coefficients are nonzero, k ∈ N. A general linear measurement process for signals consists in computing M < N inner products y j = ⟨ f , y j ⟩ for a collection of vectors {y j } j . In matrix form, g = Y f = Y Xs, where Y and X are the matrixes with {y i } i and {x j } j as columns, and Y X is an M × N-matrix. If the function families Y and X are incoherent, i.e., if the incoherence measure √ √ μ(Y, X) = N max ∣⟨y i , x j ⟩∣ ∈ [, N] ≤i, j≤N
is close to one, then under mild additional conditions the k-sparse signal f can be reconstructed from M > Cμ (Y, X)k ln N samples with overwhelming probability. Wavelet bases have proven to be very suitable for compressed sensing. It is an open question to classify the signals from certain applications, and to estimate in which appropriate B-spline basis they have a k-sparse representation. Then adequate bases and function families incoherent with the spline bases have to be identified. In the last years, the concept of sparsity entered image processing. It has proven to help immensely to accelerate the solution of inverse problems and reconstruction algorithms, e.g., in medical imaging, such as in magnetic resonance imaging [], computed tomography [], photo-acoustic tomography [], tomosynthesis [], and others. In this area, as well as in other fields of imaging, it can be expected that the combination of splines – due to their easy modeling and the fast frequency domain algorithms – multiresolution and wavelets, and sparsity will lead to novel impressing fast algorithms for image reconstruction.
.
Conclusion
In the design procedure for scaling functions of multiresolution analyses, regularity and decay features, as well as symmetry properties can be tuned by an appropriate modeling in frequency domain. The idea is to start in frequency domain with a polynomial function P that fulfils the required symmetry features, and that has a degree, such that /P decays
Splines and Multiresolution Analysis
sufficiently fast. This assures that the resulting scaling function has the desired regularity. However, /P in general is not an L -function and has to be multiplied with a localizing trigonometric polynomial that eliminates the zeros in the denominator such that P becomes square integrable. The choice of this trigonometric polynomial has to be taken carefully to be compatible with the required features modeled in /P. Then under mild additional conditions the fraction ̂= φ P is the scaling function of a multiresolution analysis. This construction can be performed for D and higher dimensional spaces likewise. In time domain, the resulting scaling function is a piecewise polynomial, thus a spline. This design procedure for scaling functions unites the concepts of splines and of multiresolution. Interestingly, the polynomial in the denominator can be of a fractional or a complex degree and therefore allows a fine tuning of the scaling function’s properties. However, the scaling function then becomes an infinite series of shifted (truncated) polynomials. The numerical calculation with the approximating basis of the multiresolution analysis in time domain would cause truncation errors, which is unfavorable. But due to the construction of φ in frequency domain, and due to the closed form there, the implementation in frequency domain with periodic boundary conditions yields a fast and stable multiresolution algorithm suitable for image analysis tasks.
.
Cross-References
Astronomy Compressive Sensing > Gabor Analysis for Imaging > Neighborhood Filters and Local D Recovery > Sampling Methods > >
References and Further Reading . Aldroubi A, Unser MA (eds) () Wavelets in medicine and biology. CRC Press, Boca Raton . Baraniuk R () Compressive sensing. IEEE Signal Process Mag ():–, . Bartels RH, Bealty JC, Beatty JC () An introduction to splines for use in computer graphics and geometric modeling. Morgan Kaufman, Los Altos . Battle G () A block spin construction of ondelettes. Part : Lemarié functions. Commun Math Phys :–
. Blu T, Unser M () The fractional spline wavelet transform: definition and implementation. In: proceedings of the th IEEE International Conference on Acoustics, Speech, and Signal Processing (ICASSP’), vol , Istanbul, Turkey, – June, , pp – . Blu T, Unser M (). A complete family of scaling functions: the (α, τ)-fractional splines. In: Proceedings of the th International Conference on Acoustics, Speech, and Signal Processing (ICASSP’), vol , Hong Kong SAR,
.
.
.
.
.
.
.
.
. .
. . . .
.
.
Splines and Multiresolution Analysis
People’s Republic of China, April –, , pp – Buhmann MD () Radial basis functions: theory and implementations. Cambridge monographs on applied and computational mathematics. Cambridge University Press, Cambridge de Boor C, Höllig K, Riemenschneider S () Box splines, vol of Applied mathematical sciences. Springer, New York de Boor C, Devore RA, Ron A () Approximation from shiftinvariant subspaces of L(Rd). Trans Am Math Soc ():– Burt PJ, Adelson EH (Apr ) The Laplacian pyramid as a compact image code. IEEE Trans Commun ():– Candès EJ, Wakin MB () An introduction to compressive sampling. IEEE Signal Process Mag ():– Champeney DC () A handbook of Fourier theorems. Cambridge University Press, Cambridge Chen H- () Complex harmonic splines, periodic quasi-wavelets, theory and applications. Kluwer Academic, Dordrecht Choi JY, Kim MW, Seong W, Ye JC () Compressed sensing metal artifact removal in dental ct. In: Proceedings of IEEE International Symposium on Biomedical Imaging (ISBI), June– July , Boston, pp – Christensen O () An introduction to frames and riesz bases. Birkhäuser, Boston Christensen O () Frames and bases: an introductory course (applied and numerical harmonic analysis). Birkhäuser, Boston Chui CK () Multivariate splines. SIAM, Philadelphia Chui C () Wavelets—a tutorial in theory and practice. Academic, San Diego Chui CK (ed) () Wavelets: a tutorial in theory and applications. Academic, Boston Condat L, Forster-Heinlein B, Van De Ville D () A new family of rotation-covariant wavelets on the hexagonal lattice. In: SPIE Wavelets XII, Aug , San Diego Condat L () Image database. Online Ressource. http://www.greyc.ensicaen.fr/_ lcondat/imagebase.html (Version of Apr ) Dahmen W, Kurdila A, Oswald P (eds) () Multiscale wavelet methods for partial differential equations, vol of Wavelet
.
. .
.
.
.
.
.
.
. .
. .
.
.
analysis and its applications. Academic, San Diego Daubechies I () Ten lectures on wavelets. Society for Industrial and Applied Mathematics, Philadelphia Dierckx P () Curve and surface fitting with splines. McGraw-Hill, New York Feilner M, Van De Ville D, Unser M () An orthogonal family of quincunx wavelets with continuously adjustable order. IEEE Trans Image Process ():– Forster B, Blu T, Unser M () Complex B-splines. Appl Comput Harmon Anal (): – Forster B, Blu T, Van De Ville D, Unser M () Shiftinvariant spaces from rotation-covariant functions. Appl Comput Harmon Anal (): – Forster B, Massopust P () Statistical encounters with complex B-splines. Constr Approx ():– Frikel J () A new framework for sparse regularization in limited angle x-ray tomography. In IEEE international symposium on biomedical imaging, Rotterdam Giles RC, Kotiuga PR, Mansuripur M () Parallel micromagnetic simulations using Fourier methods on a regular hexagonal lattice. IEEE Trans Magn ():– Grigoryan AM () Efficient algorithms for computing the -D hexagonal Fourier transforms. IEEE Trans Signal Process (): – Hales TC () The honeycomb conjecture. Discr Comput Geom :– Heil C, Walnut DF () Fundamental papers in wavelet theory. Princeton University Press, Princeton. New edition Jones DS () Generalised functions. McGrawHill, London Lai M-J, Schumaker LL () Spline functions on triangulations. Cambridge University Press, Cambridge Laine AF, Schuler S, Fan J, Huda W () Mammographic feature enhancement by multiscale analysis. IEEE Trans Med Imaging ():– Legrand P () Local regularity and multifractal methods for image and signal analysis. In: Abry P, Gonçalves P, Véhel L (eds) Scaling, fractals and wavelets, chap . Wiley-ISTE, London
Splines and Multiresolution Analysis
. Lemarié P-G () Ondelettes a localisation exponentielle. J Math pures et Appl :– . Lesage F, Provost J () The application of compressed sensing for photo-acoustic tomography. IEEE Trans Med Imaging ():– . Lipow PR, Schoenberg IJ () Cardinal interpolation and spline functions. III: Cardinal hermite interpolation. Linear Algebra Appl : – . Louis AK, Maaß P, Rieder A () Wavelets: theory and applications. Wiley, New York . Lustig M, Donoho D, Pauly JM () Sparse MRI: the application of compressed sensing for rapid MR imaging. Magn Reson Med (): – . Mallat S () Multiresolution approximations and wavelet orthonormal bases of L(R). Trans Am Math Soc :– . Mallat SG () A wavelet tour of signal processing. Academic, San Diego . Mersereau RM () The processing of hexagonally sampled two-dimensioanl signals. Proc IEEE ():– . Meyer Y () Wavelets and operators. Cambridge University Press, Cambridge . Middleton L, Sivaswamy J () Hexagonal image processing: a practical approach. Advances in pattern recognition. Springer, Berlin . Nicolier F, Laligant O, Truchetet F () B-spline quincunx wavelet transforms and implementation in Fourier domain. Proc SPIE :– . Nicolier F, Laligant O, Truchetet F () Discrete wavelet transform implementation in Fourier domain for multidimensional signal. J Electron Imaging :– . Nürnberger G () Approximation by Spline functions. Springer, Berlin . Plonka G, Tasche M () On the computation of periodic splione wavelets. Appl Comput Harmon Anal :– . Püschel M, Rötteler M () Algebraic signal processing theory: D spatail hexagonal lattice. IEEE Trans Image Proc ():– . Rabut C (a) Elementary m-harmonic cardinal B-splines. Numer Algorithms :– . Rabut C (b) High level m-harmonic cardinal B-splines. Numer Algorithms :– . Rice University () Compressive sensing resources. Online Ressource. http://dsp. rice.edu/cs (Version of Apr )
. Rudin W () Functional analysis. International series in pure and applied mathematics. McGraw-Hill, New York . Sablonnière P, Sbibih D () B-splines with hexagonal support on a uniform three-direction mesh of the plane. C R Acad Sci Paris Série I :– . Schempp W () Complex contour integral representation of cardinal spline functions, vol of Contemporary mathematics. American Mathematical Society, Providence . Schoenberg IJ () Contributions to the problem of approximation of equidistant data by analytic functions. Part a.–on the problem of osculatory interpolation. A second class of analytic approximation formulae. Quart Appl Math :– . Schoenberg IJ () Contributions to the problem of approximation of equidistant data by analytic functions. Part a.–on the problem of smoothing or graduation. A first class of analytic approximation formulae. Quart Appl Math :– . Schoenberg IJ () Cardinal interpolation and spline functions. J Approx Theory :– . Schoenberg IJ () Cardinal interpolation and spline functions. II: Interpolation of data of power growth. J Approx Theory :– . Schwartz L () Théorie des distributions. Hermann, Paris . Unser M () Splines: A perfect fit for medical imaging. In: Sonka M., Fitzpatrick JM (eds) Progress in biomedical optics and imaging, vol , no. , vol , Part I of Proceedings of the SPIE international symposium on medical imaging: image processing (MI’), San Diego, – Feb, pp – . Unser M, Blu T (Mar ) Fractional splines and wavelets. SIAM Rev ():– . Unser M, Aldroubi A, Eden M (Mar ) Fast B-spline transforms for continuous image representation and interpolation. IEEE Trans Pattern Anal Mach Intell ():– . Unser M, Aldroubi A, Eden M (Mar ) On the asymptotic convergence of B-spline wavelets to gabor functions. IEEE Trans Info Theory : – . Unser M, Aldroubi A, Eden M (Feb a) Bspline signal processing: Part I—Theory. IEEE Trans Signal Process ():–
Splines and Multiresolution Analysis
. Unser M, Aldroubi A, Eden M (Feb b) Bspline signal processing: Part II—Efficient design and applications. IEEE Trans Signal Process ():– . Van De Ville D, Blu T, Unser M, Philips W, Lemahieu I, Van de Walle R () Hex-splines: a novel spline family for hexagonal lattices. IEEE Trans Image Process ():– . Van De Ville D, Blu T, Unser M (Nov ) Isotropic polyharmonic B-splines: scaling functions and wavelets. IEEE Trans Image Process ():– . Watson AB, Ahumuda AJ, Jr () Hexagonal orthogonal-oriented pyramid as a model of
image representation in visual cortex. IEEE Trans Biomed Eng ():– . Wendt H, Roux SG, Jaffard S, Abry P () Wavelet leaders and bootstrap for multifractal analysis of images. Signal Process :– . Wojtaszczyk P () A mathematical introduction to wavelets, vol of London mathematical society student texts. Cambridge University Press, Cambridge . Young R () An introduction to nonharmonic Fourier series. Academic, New York (revised first edition )
Introduction... ............. ............. ............ ............. ............. ..
.
Historical Background.. .. .. .. .. .. ... .. .. .. .. .. .. .. .. .. .. ... .. .. .. .. .. .. .. .. .. .
. .. ... ... ... ... ... .. ..
Mathematical Modeling and Application.. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. . Mathematical Foundations. .. .. .. .. .. .. .. .. .. .. .. .. .. .. . .. .. .. .. .. .. .. .. .. .. .. .. .. .. . Regularity and Decay Under the Fourier Transform. . . . . . . . . . . .. . . . . . . . . . . . . . . . . Criteria for Riesz Sequences and Multiresolution Analyses. . . . . . . . . . . . . . . . . . . . . Regularity of Multiresolution Analysis. ... ... ... ... ... ... .... ... ... ... ... ... ... .... . Order of Approximation. ... ........ ........ ........ ....... ........ ........ ........ .... Wavelets. .. .. ... .. ... .. ... .. ... .. ... .. ... .. ... .. ... .. ... .. ... .. ... .. ... .. ... .. ... .. ... .. . B-Splines. .... ...... ....... ...... ...... ...... ....... ...... ...... ...... ....... ...... ...... . Polyharmonic B-Splines. ....... ....... ....... ....... ....... ....... ....... ....... ......
. ..
Survey on Mathematical Analysis of Methods. . . . .. . . . . . .. . . . . . .. . . . . . .. . . . . . Schoenberg’s B-Splines for Image Analysis – the Tensor Product Approach... ........ ....... ....... ........ ....... ........ ....... ....... ........ ....... .... .. Fractional and Complex B-Splines. . . .. .. . .. .. . .. .. . .. .. . .. .. . .. .. . .. .. .. . .. .. . .. .. . .. Polyharmonic B-Splines and Variants. . .. . .. . .. . .. . .. . .. . .. . .. . .. . .. . .. . .. . .. . .. . .. .. Splines on Other Lattices. . .. .. . .. .. .. .. . .. .. .. . .. .. .. .. . .. .. .. . .. .. .. .. . .. .. .. . .. .. .. ... Splines on the Quincunx Lattice. .. .... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... . ... Splines on the Hexagonal Lattice. .. .. .. .. . .. .. .. .. .. .. .. .. .. . .. .. .. .. .. .. .. .. .. . .. .. .
Numerical Methods. ........... ................ ................ ................ .
.
Open Questions. ..... ..... ..... .... ..... ..... ..... ..... .... ..... ..... ..... .... ..
.
Conclusion.. .... .... ..... .... .... ..... .... .... ..... .... .... ..... .... .... .... ....
.
Cross-References. .. ..... ..... ..... ..... ...... ..... ..... ..... ..... ..... ...... ....
Otmar Scherzer (ed.), Handbook of Mathematical Methods in Imaging, DOI ./---_, © Springer Science+Business Media LLC
Splines and Multiresolution Analysis
Abstract: Splines and multiresolution are two independent concepts, which – considered together – yield a vast variety of bases for image processing and image analysis. The idea of a multiresolution analysis is to construct a ladder of nested spaces that operate as some sort of mathematical looking glass. It allows to separate coarse parts in a signal or in an image from the details of various sizes. Spline functions are piecewise or domainwise polynomials in one dimension (D) resp. nD. There is a variety of spline functions that generate multiresolution analyses. The viewpoint in this chapter is the modeling of such spline functions in frequency domain via Fourier decay to generate functions with specified smoothness in time domain resp. space domain. The mathematical foundations are presented and illustrated at the example of cardinal B-splines as generators of multiresolution analyses. Other spline models such as complex B-splines, polyharmonic splines, hexagonal splines, and others are considered. For all these spline families exist fast and stable multiresolution algorithms which can be elegantly implemented in frequency domain. The chapter closes with a look on open problems in the field. AMS Subject Classification (): A Spline approximation D Numerical analysis – Splines U Computing methodologies and applications – Image processing T Numerical methods in Fourier analysis Keywords Cardinal B-splines ⋅ image processing ⋅ multiresolution analysis ⋅ order of approximation ⋅ polyharmonic B-splines ⋅ riesz basis ⋅ scaling function two-scale relation
.
Introduction
This chapter deals with two originally independent concepts, which were recently combined, and since together have a strong impact on signal and image analysis: the concept of splines, i.e., of piecewise polynomials, and the concept of multiresolution analysis, i.e., splitting functions – or more general data – into coarse approximations and details of various sizes. Eversince the combination of the two concepts, they have led to a load of new applications in, e.g., signal and image analysis, as well as in signal and image reconstruction, computer vision, numerics of partial differential equations, and other numerical fields. An impression of the vast area of applications can be gained in, e.g., [, , , ]. Already the spline functions alone proved to be very useful for mathematical analysis as well as for signal and image processing, analysis and representation, for computer graphics and many more, see, e.g., [, , , , , , ]. An example for a family of spline functions are I. J. Schoenberg’s polynomial splines with uniform knots [, ]: β m (t) =
m+ k m+ ) (t − k)+m , ∑ (−) ( m! k= k
m ∈ N.
(.)
Here, t+m denotes the one-sided power function, i.e., t+m = for t < and t+m = t m for t ≥ .
Splines and Multiresolution Analysis
1 0.8 0.6 0.4 0.2
–1
1
2
4
3
5
⊡ Fig. -
Cardinal B-splines of degree m = , . . . ,
The B-splines β m can be easily generated by an iterative process. Let β (t) = χ[,) (t) be the characteristic function of the interval [, ). Then the B-spline of degree m is derived by the convolution product β m = β ∗ β m− = β ∗ . . . ∗ β
m+-times where
for m ∈ N,
β ∗ β m− (x) = ∫ β (y)β m− (x − y) d y = ∫ R
(.)
β m− (x − y) d y.
For an illustration of the cardinal B-splines, see > Fig. -. Splines had their second breakthrough as Battle [] and Lemarié [] discovered that B-splines generate multiresolution analyses. The simple form of the B-splines and their compact support, in particular, were convenient for designing multiresolution algorithms and fast implementations. Definition Let A : Rn → Rn be a linear mapping that leaves Zn invariant, i.e., A(Zn ) ⊂ Zn and that has (real or complex) eigenvalues with absolute values greater than . A multiresolution analysis associated with the dilation matrix A is a sequence of nested subspaces (Vj ) j∈Z of L (Rn ) such that the following conditions hold: (i) (ii) (iii) (iv) (v) (vi)
. . . ⊂ V− ⊂ V ⊂ V ⊂ . . ., ⋂ j∈Z Vj = {}, Span ⋃ j∈Z Vj is dense in L (Rn ), f ∈ Vj ⇐⇒ f (A− j ●) ∈ V , f ∈ V ⇐⇒ f (● − k) ∈ V for all k ∈ Zn . There exists a so-called scaling function φ ∈ V such that the family {φ(● − k)} k∈Zn of translates of φ forms a Riesz basis of V .
Here, L (Rn ) denotes the vector space of square-integrable functions f : Rn → C with norm ∥ f ∥ = (∫ ∣ f (x)∣ dx) Rn
Splines and Multiresolution Analysis
and corresponding inner product ⟨ f , g⟩ = ∫
Rn
f (x) g(x) dx,
where g denotes the complex conjugate of g. The elements in L (Rn ) are also called functions of finite energy. Riesz bases are a slightly more general concept than orthonormal bases. In fact, Riesz bases are equivalent to orthonormal bases and therefore can be generated by applying a topological isomorphism on some orthonormal basis. Definition A sequence of functions { f n }n∈Z in a Hilbert space V is called a Riesz sequence if there exist positive constants A and B, the Riesz bounds, such that A∥c∥ l ≤ ∥ ∑ c k f k ∥ ≤ B∥c∥ l k∈Z n
V
for all scalar sequences C = (c k ) k∈Zn ⊂ l (Zn ). A Riesz sequence is called a Riesz basis if it additionally spans the space V . A good introduction to Riesz bases, their properties, and their relation to orthonormal bases is given in the monography by Young []. Multiresolution constructions with splines are treated in numerous sources. As a starting point, there are, e.g., the books by Christensen [, ] and Wojtaszczyk []. The mathematical properties in Definition have intuitive interpretations. A function f ∈ L (Rn ), which is projected orthogonally on Vj , is approximated with the so-called resolution A j . In fact, let P j : L (Rn ) → Vj denote the orthogonal projection operator. Then (ii) gives that by going to lower resolutions, all details are lost: lim ∥P j f ∥ = .
j→−∞
Whereas when the resolution is increased, j → ∞, more and more details are added. By (iii), the projection then converges to the original function f : lim ∥ f − P j f ∥ = .
j→∞
Hereby, the rate of convergence depends on the regularity of f . The approximation spaces Vj are nested, which allows for computing coarser approximations in Vk for k < j for functions f ∈ Vj . The scaling Ak enlarges details. Property (iv) shows that the approximation spaces have a similar structure over the scales and emanate from one another. The translation invariance (v) ensures that the analysis of a function in Vj is independent of the starting time or location. And property (vi) finally ensures the beautiful and mighty property that the whole sequence of nested approximation spaces can be generated by translates and scalings of one
Splines and Multiresolution Analysis
single function – the scaling function. In fact, (vi) together with (iv) yields that {φ(A j ● −k), k ∈ Zn } is a Riesz basis for Vj . While moving from the coarser approximation space Vj to the finer, larger space Vj+ , information has to be added. In fact, there is an orthogonal complement Wj , j ∈ Z, such that Vj+ = Vj ⊕ Wj . These spaces are called detail spaces or wavelet spaces. It is well known that these spaces also possess a Riesz basis spanned by shifts of ∣ det A∣ − generators, the wavelets ψ , . . . , ψ∣ det A∣− . Here, A is the dilation matrix in Definition . The wavelets can be constructed from the scaling function. As a consequence, the knowledge of just the single function φ allows for the construction of the approximation spaces Vj and for the wavelet spaces Wj . Detailed information on the generation of wavelets and their properties can be found in various books, e.g., [, , , , ]. Example A simple example for a multiresolution analysis on L (R) is given by piecewise constant functions. Consider the characteristic function φ = χ[,) of the interval semiopen interval [, ). Then φ generates a dyadic multiresolution analysis, i.e., for A = . The approximation spaces are Vj = span {χ[,) ( j ● −k)} k∈Z
L (R)
.
They consist of functions constant on intervals of the form [k− j , (k + )− j ). The spaces are obviously nested and separate L (R) in the sense of Definition (ii). Since piecewise constant functions are dense in L (R), (iii) holds. (iv) – (vi) hold by construction. In fact, this multiresolution analysis is generated by the B-spline β as scaling function. The B-spline basis operates as mean-value operator over the support interval. The corresponding wavelet extracts the details, i.e., the deviation from the mean value. To this end, it operates as a difference operator. > Figure - shows the scaling function β and the corresponding wavelet, the so-called Haar wavelet. In > Fig. -, an example of a multiresolution is given.
1
1
0.8 0.5 0.6 0.2
0.4
0.4
0.6
0.8
–0.5
0.2
–1 0.2
0.4
0.6
0.8
1
⊡ Fig. -
Left: The B-spline β is a scaling function and operates as mean value. Right: The corresponding wavelet, the Haar wavelet, operates as a local difference operator
1
Splines and Multiresolution Analysis
1 0.8 0.6 0.4 0.2 2 a 1
4
6
8
Original step function 1 0.75 0.5 0.25
0.8 0.6 0.4 0.2 2
4
b 1
–0.25 –0.5 –0.75 –1
2
4
6
8
4
6
8
4
6
8
6 8 First iteration: approximation and details 1 0.75 0.5 0.25
0.8 0.6 0.4 0.2 2
4
c 1
–0.25 –0.5 –0.75 –1
2
6 8 Second iteration: approximation and details 1 0.75 0.5 0.25
0.8 0.6 0.4 0.2 2 d
4
–0.25 –0.5 –0.75 –1
2
6 8 Third iteration: approximation and details
⊡ Fig. - Multiresolution decomposition of the step function (a) with β as scaling function and with the Haar wavelet. The approximations (left column) are further iterated and decomposed into a coarser approximation and details (right column), until the coarsest approximation step, here the mean value, is reached. The sum of the coarsest approximation in the third iteration and of all details yields the original function (a). No information is lost in a multiresolution decomposition.
Splines and Multiresolution Analysis
.
Historical Background
The idea of piecewise polynomial functions and splines goes back to Schoenberg [, ]. In the s, when computer power started to be used for numerical procedures such as data fitting, interpolation, solving differential equations or computer aided geometric design, splines experienced an extreme upsurge. Schoenberg invented and strongly contributed to the concept of cardinal splines, which have equidistant nodes on the integers, see, e.g., [, , ] and many more. As a parallel development, in the s, the adaption of signal resolution to only process relevant details for a particular task evolved. For example, for computer vision, a multiresolution pyramid was introduced by Burt and Adelson []. It allowed to process an image first on a low resolution level and then selectively increase the resolution and add more detailed information when needed. The definition of a dyadic multiresolution analysis, i.e., A = Id, was contributed by Mallat [] and Meyer []. An interesting and in some parts historical collection on the most important articles in multiresolution and wavelet theory was assembled by Heil and Walnut []. The concepts of splines and multiresolution where joined by Lemarié [] and Battle [], when they showed that cardinal B-splines are scaling functions for multiresolution analyses. This led to many more developments of piecewise polynomial scaling functions for various settings and also multidimensions [], as, e.g., polyharmonic B-splines [, ] and other functions inspired from radial basis functions []. In , S. Mallat published his famous algorithm for multiresolution and wavelet analysis []. He had developed an efficient numerical method such that multiresolution decompositions could be calculated in a fast way. For the splines, M. Unser et al. proposed a fast implementation [, , ] which strongly contributed to the breakthrough of splines for signal and image analysis. In the last years, periodic, fractional, and complex versions of splines for multiresolution were developed, e.g., [, , , , ]. Many of them use a Fourier domain filter algorithm which allows for infinite impulse response filters. The former important feature of compact support of the cardinal B-splines and other functions is no longer a limiting criterion. Therefore, it can be expected that many new contributions on splines will still be made in future by modeling signal and image features in Fourier domain.
.
Mathematical Modeling and Application
..
Mathematical Foundations
...
Regularity and Decay Under the Fourier Transform
An important idea behind splines and multiresolution is the relation between regularity in time domain and decay in frequency domain, respectively between decay in time domain and regularity in frequency domain. To illustrate this, the notion of the Schwartz space is very useful [, , , ].
Splines and Multiresolution Analysis
Definition The subspace of functions f ∈ C ∞ (Rn ) with sup sup ( + ∥x∥ ) N ∣D α f (x)∣ < ∞
∣α∣≤N x∈R n
for all N = , , , . . .
is called the space of rapidly decreasing functions or Schwartz space S(Rn ). The norms induce a Fréchet space topology, i.e., the space S is complete and metrizable. α
αn
Here, D α = ( ∂x∂ ) ⋯ ( ∂x∂ n ) for every multi-index α = (α , . . . , α n ) ∈ Non . The dual space S ′ (Rn ), endowed with the weak-∗-topology, is called space of tempered distributions. The following famous linear transform relates the view points of the space domain and of the frequency domain: Definition The Fourier transform, defined by F f (ω) := ̂ f (ω) := ∫
Rn
f (x)e −i⟨ω,x⟩ dx,
ω ∈ Rn ,
is a topological isomorphism on L (Rn ) and on S(Rn ). Its inverse is given by f (x) =
i⟨ω,x⟩ ̂ dω ∫ f (ω)e (π)n Rn
in L (Rn ) resp. in S(Rn ).
The Fourier transform can be extended to the space of tempered distributions. For T ∈ S ′ (Rn ) the Fourier transform is defined in a weak sense as ̂ ̂) for all φ ∈ S(Rn ). F T(φ) := T(φ) := T(φ Also on S′ (Rn ), the Fourier transform is a topological isomorphism. The Fourier transform has the nice property to relate polynomials and differential operators. Theorem (i)
Let f ∈ S(R). Then for all k ∈ N F ( f (k) ) (ω) = (iω) k ̂ f (ω), and
(ii)
(iii)
̂ f (k) (ω) = F ((−i ●) k f ) (ω).
Let P be an algebraic polynomial in Rn , say P(x) = ∑α c α x α . . . x nα n , and let f ∈ S(Rn ). Then ̂f = P(iD) ̂ f and P f, F (P ( D) f ) = P ̂ i where P(iD) = ∑α c α i ∣α∣ D α . Part (ii) also holds for f ∈ S ′ (Rn ).
Splines and Multiresolution Analysis
Example The Fourier transform of the polynomial x k is the tempered distribution k i k ddx k δ, k ∈ N . For the construction of a multiresolution analysis, the scaling function can be used as a starting point. The idea is to choose a scaling function of a certain regularity, such that the generated multiresolution analysis inherits the smoothness properties. In particular, for the splines the idea is to model the regularity via decay in Fourier domain. The following theorem gives a motivation for this. The result can be deduced from the considerations above, and the fact that S(Rn ) is dense in L (Rn ): Theorem
Let f ∈ L (Rn ) and its Fourier transform decay as ∣̂ f (ω)∣ ≤ C( + ∥ω∥)−N−ε
for some ε > . Then all partial derivatives of order ≤ N − n are continuous and in L (Rn ). These results allow to construct scaling function with explicit regularity and decay properties, in space and in frequency domain. However, some criteria are needed to verify that the constructed function generates a multiresolution analysis.
...
Criteria for Riesz Sequences and Multiresolution Analyses
The following is an explicit criterion to verify whether some function φ is a scaling function. Theorem Let A be a dilation matrix and let φ ∈ L (Rn ) be some function satisfying the following properties: (i) (ii)
{φ(● − k)} k∈Zn is a Riesz sequence in L (Rn ). φ satisfies a scaling relation. I.e., there is a sequence of coefficients (a k ) k∈Zn such that φ(A− x) = ∑ a k φ(x + k)
in L (Rn ).
(.)
k∈Z n
(iii)
̂∣ is continuous at and φ ̂() ≠ . ∣φ
Then the spaces Vj = span {φ(A j ● −k)} k∈Zn ,
j ∈ Z,
form a multiresolution analysis of L (Rn ) with respect to the dilation matrix A. Proof
See, e.g., [, Theorem .] for the D case.
∎
For particular applications, the Riesz basis property (i) of {φ(● − k)} k∈Zn in V is not enough, but an orthonormal basis is needed. An example for such an application is the
Splines and Multiresolution Analysis
denoising of signals contaminated with Gaussian white noise [, Chap. X, Sect. ..]. However, there is an elegant mathematical method to orthonormalize Riesz bases generated by shifts of a single function. Theorem (i)
Let φ ∈ L (Rn ). Then the following holds:
{φ(● − k)} k∈Zn is a Riesz sequence in L (Rn ) if and only if there are constants c and C, such that ̂(ω + πk)∣ ≤ C < ∞ almost everywhere. < c ≤ ∑ ∣φ k∈Z n
(ii)
̂(ω + πk)∣ is strictly positive and I.e., the autocorrelation filter M(ω) := ∑ k∈Zn ∣φ bounded from above. {φ(● − k)} k∈Zn is an orthonormal sequence if and only if ̂(ω + πk)∣ = almost everywhere. ∑ ∣φ k∈Z n
(iii)
If {φ(● − k)} k∈Zn is Riesz basis of a subspace X of L (Rn ), then there exists a function Φ ∈ L (Rn ), namely ̂(ω) φ ̂ Φ(ω) =√ ̂(ω + πk)∣ ∑ k∈Zn ∣φ
(.)
such that {Φ(● − k)} k∈Zn is an orthonormal basis of X. Proof
∎
See, e.g., [] and [, Chap. VII].
Due to this theorem, every scaling function can be orthonormalized. Let φ ∈ L (Rn ) be some scaling function that generates a multiresolution analysis {Vj } j∈Z of L (Rn ). Then the family {Φ j,k } k∈Zn with Φ j,k (x) = − j/ Φ(− j (x − k)), and Φ as defined in (> .) is an orthonormal basis of the space Vj , j ∈ Z. Example A simple possibility to construct a dyadic multiresolution analysis in L (Rn ) is the tensor product approach. Let (Vj ) j∈Z be a dyadic (i.e., A = ) multiresolution analysis of L (R) with scaling function φ. Then (V j ) j∈Z with V j = Vj ⊗ . . . ⊗ Vj together with
n -times
the scaling function φ(x , . . . , x n ) = φ(x ) ⋅ . . . ⋅ φ(x n ) forms a multiresolution analysis of L (Rn ) and dilation matrix Id. In the same way, the scaling function φ(x , . . . , x n ) = φ (x ) ⋅ . . . ⋅ φ n (x n ) generates a multiresolution analysis of L (Rn ), if every φ k , k = , . . . , n, is a scaling function of some D multiresolution analysis with dilation factor a ∈ N/{}.
Splines and Multiresolution Analysis
...
Regularity of Multiresolution Analysis
In signal and image analysis the choice of an appropriate analysis basis is crucial. Here, appropriate means that the features of the basis such as smoothness should be in accordance with the properties of the functions to analyze. For example, analyzing a smooth signal or image with a fractal basis in general yields results that are difficult to interpret and to work with in practice. In this case, the signal resp. the image model does not match the model of the basis. The next section will show that the family of spline bases helps to avoid such difficulties, because the splines allow for a good adjustment due to their regularity parameter m, cf. (> .) and (> .). The following definition specifies the term “regular.” (See [].) Definition Denote C r the class of r-times continuously differentiable functions in Rn , C the class of continuous functions, and C − the class of measurable functions. (i)
Let r = −, , , . . . . A function f : Rn → C is called r-regular, if f ∈ C r and ∣
(ii)
∂α Ak f (x)∣ ≤ ∂x α ( + ∥x∥) k
for every k ∈ N , every multi-index α with ∣α∣ ≤ max(r, ) and constants A k . A multiresolution analysis of L (Rn ) is called r-regular if it generated by an r-regular scaling function.
It is important to note that the orthonormalization procedure (> .) does not affect the regularity of the corresponding basis. For the orthonormalized scaling function Φ of a multiresolution analysis, the same regularity properties hold. Proposition Let φ ∈ L (Rn ) be an r-regular function, such that {φ(● − k)} k∈Zn forms a Riesz sequence. Then the via (> .) orthonormalized function Φ is also r-regular [].
...
Order of Approximation
Having found a scaling function that generates a multiresolution analysis, how good do the corresponding approximation spaces Vj approximate some function f ∈ L (Rn ) of a certain regularity? Let H k (Rn ) = { f ∈ L (Rn ) : ∥ f ∥H k :=
f ∥L < ∞}, ∥( + ∥●∥Rn ) k ̂ (π)n
k ∈ N ,
denote the Sobolev spaces. The following criterion for the order of approximation turns out to be easy to verify for splines.
Splines and Multiresolution Analysis
Theorem Let φ ∈ L (Rn ) satisfy the following properties [, Theorem .]: (i) (ii) (iii)
̂ is bounded on some neighborhood of the origin. /φ Let B ε be some open ball centered at the origin and let E := B ε + (πZn /{}). For ̂ of order ≤ α are in L (E). some α > k + n/, all derivatives of φ ̂(ω) = for all ∣γ∣ < k and all ω ∈ πZd /{}. Dγ φ
Then V = span {φ(● − k)} k∈Zn provides approximation order k: For f ∈ H k (Rn ), min{∥ f − s(●/h)∥L , s ∈ V } ≤ const. h k ∥ f ∥H k
...
for all h > .
Wavelets
For the step from a coarser approximation space Vj to a finder one Vj+ information has to be added. It is contained in the wavelet or detail space Wj , which is the orthonormal complement of Vj in Vj+ : Vj+ = Vj ⊕ Wj . It follows that Vj+m = Vj ⊕ ⊕ k− l= Wj+l , and hence L (Rn ) = ⊕ Wj
(.)
j∈Z
can be decomposed in a direct sum of mutually orthogonal detail spaces. Moreover, the detail spaces Wj inherit the scaling property from Definition (iv) for the approximation spaces Vj . For all j ∈ Z, f ∈ Wj ⇐⇒ f (A− j ●) ∈ W . The question now is whether there is also a simple basis generated by the shifts of one or few functions, the wavelets. The following definition is motivated from > Eq. (.). Definition Let A be a dilation matrix, and let {ψ l } l=, ... ,s be a set of functions in L (Rn ), such that the family {∣ det A∣ j/ ψ l (A j ● −k) ∣ l = , . . . , s, j ∈ Z, k ∈ Zn } forms an orthonormal basis of L (Rn ). Then {ψ l } l=, ... ,s is called wavelet set associated with A. What qualitative properties do the wavelets have? The approximation spaces Vj are generated by the scaling function, which operates as a low pass filter. This can be seen from ̂() ≠ , resp. from Theorem (i): /φ ̂ is bounded in some neighborhood Theorem (iii) φ of the origin. Therefore, the added details and thus the wavelets have to carry the highfrequency information. In addition, the wavelets ψ in W are elements of V and therefore have the form (.) ψ(A− x) = ∑ a k φ(x − k) k∈Z n
Splines and Multiresolution Analysis
in L -norm, where {a k } k∈Zn are the Fourier coefficients of a certain πZn -periodic function. Proposition Let (Vj ) j∈Z be a multiresolution analysis of L (Rn ) with respect to the dilation matrix A and with scaling function φ. Then for a function f ∈ V if and only if ̂ ̂(ω) almost everywhere. f (AT ω) = m f (ω)φ Here, m f ∈ L ([ , π]n ) and ∥m f ∥L ([ ,π]n ) =
∥ f ∥L (Rn ) . ∣ det A∣
For a proof, see, e.g., [, ]. Note that for a wavelet ψ as in mψ (ω) =
>
Eq. (.) there holds
∑ a k e i⟨ω,k⟩ . ∣ det A∣ k∈Zn
How many wavelets, i.e., generators of W , are needed to span the space? The parameter s in Definition is yet unspecified. In fact, s depends on the scaling matrix. A leaves the lattice Zn invariant, AZn ⊂ Zn . The number of cosets is ∣ det A∣ = ∣Zn /AZn ∣ (see [, Proposition .]). It turns out that q = ∣ det A∣ − wavelets are needed to generate the space W . To motivate this, for a start, let f ∈ V be an arbitrary function. Denote γ , . . . , γ q representatives of the q + cosets of AZn in Zn . Then each coset can be written as γ m + AZn , m = , . . . , q. The function f has the representation f (A− x) = ∑ c k ( f )φ(x − k), ∣ det A∣/ k∈Z n
(.)
or in Fourier domain ̂ f (AT ω) =
̂(ω), c f (ω)φ ∣ det A∣/
(.)
in L -sense and with an appropriate πZn -periodic function c f (ω) with Fourier coefficients (c k ( f )) k∈Zn . Then c f (ω) can be decomposed with respect to the cosets: q
c f (ω) = ∑ c k ( f )e i⟨ω,k⟩ = ∑
∑
m= k∈γ m +AZ n
k∈Z n
c k ( f )e i⟨ω,k⟩
q
q
= ∑ e i⟨ω,γ m ⟩ ∑ c k+γ m ( f )e i⟨ω,k⟩ = ∑ c m f (ω), m=
m=
k∈AZ n
where i⟨ω,γ m ⟩ i⟨ω,k⟩ = e i⟨ω,γ m ⟩ ∑ c Ak+γ m ( f )e i⟨ω, Ak⟩ cm ∑ c k+γ m ( f )e f (ω) = e k∈AZ n
=e
i⟨ω,γ m ⟩
∑ c Ak+γ m ( f )e k∈Z n
k∈Z n i⟨AT ω,k⟩
=
T e i⟨ω,γ m ⟩ κ m f (A ω).
Splines and Multiresolution Analysis
This representation exists for all functions V , in particular for φ and the wavelets. The following theorem indicates, how many wavelets are needed to generate the space W , such that W ⊕ V = V . Theorem Let φ ∈ V be a scaling function and let ψ , . . . , ψ q ∈ V . Then the family {φ(● − k)} k∈Zn is an orthonormal system if and only if q
m ∑ ∣κ φ (ω)∣ = almost everywhere.
(.)
m= q
The system {φ(● − k)} k∈Zn ∪ ⋃m= {ψ m (● − k)} k∈Zn is an orthonormal basis in V if and only if the so-called polyphase matrix ⎛κ φ (ω) ⎜ ⋮ ⎜ ⎝κ φq (ω)
κψ (ω) ⋮ q κψ (ω)
⋯ ⋱ ⋯
κψ q (ω)⎞ ⋮ ⎟ ⎟ q κψ q (ω)⎠
is unitary for almost all ω ∈ Rn . The proof for a more general version of this theorem is given in [, Sect. .]. A summary and a condition for r-regular wavelets yields in the following theorem. Theorem Consider a multiresolution analysis on Rn associated with a dilation matrix A. (i) (ii)
Then there exists an associated wavelet set consisting of q = ∣ det A∣ − functions. If the multiresolution analysis is r-regular, and in addition q + > n, then there exists an associated wavelet set consisting of q r-regular functions.
The idea of the proof is that for an r-regular function φ on Rn and a πZn -periodic C ∞ ̂(ω) = η(ω)φ ̂(ω) is an r-regular function. function η(ω) the convolution ψ defined by ψ For an explicit proof, see again []. Example As a continuation of Example , the wavelet function corresponding to φ = χ[,) is derived. To this end, consider the space L (R) and the dilation A = . Then q = det A − = ; thus γ = and γ = are representatives of the cosets of A. That is, there is only a single wavelet needed to generate W . > Equation (.) yields x √ φ ( ) = √ φ(x) + √ φ(x − ) for the normalized generator of V . Thus c (φ) = κ φ ( ω )
=
κ φ ( ω )
=
√
. Then by
>
√
and c (φ) =
√
e i ω . This implies
Eq. (.) of Theorem the family {φ(● − k)} k∈Z is
Splines and Multiresolution Analysis
orthogonal, since ∣κ φ (ω)∣ + ∣κ φ (ω)∣ = . The polyphase matrix (
κ φ (ω) κ φ (ω)
κψ (ω) ) κψ (ω)
can be completed to a unitary matrix by choosing κψ (ω) = wavelet then has the representation
√
= −κψ (ω). The corresponding
x √ ψ ( ) = √ φ(x) − √ φ(x − ), corresponding to (> .). This yields the Haar wavelet ψ as illustrated in
..
>
Fig. -.
B-Splines
Several of the criteria for scaling functions and multiresolution analyses given in the previous section are based on the Fourier representation of the scaling function, e.g., the Riesz sequence criterion and the orthonormalization trick in Theorem , as well as the criterion for the order of approximation in Theorem . For this reason, the modeling of a scaling function in Fourier domain to achieve certain specific properties is promising. Aiming at constructing a scaling function φ ∈ L (R) of regularity r = −, , , . . ., this property is considered in Fourier domain: It is a decay property of the Fourier transform ̂ (compare with > Sect. ...): φ ̂(ω) = O ( φ
) ∥ω∥r+
for ∥ω∥ → ∞.
Taking into account Theorem a first model for the scaling function in Fourier domain is ̂(ω) = φ where the function relation (> .)
(ω) , ω r+
ω ∈ R,
(.)
still has to be specified. Since scaling functions satisfy a scaling x φ ( ) = ∑ h k φ(x − k) in L (R), k∈Z
the Fourier transform of this equation yields ̂(ω) = H(ω)φ ̂(ω), φ where (h k ) k∈Z is the sequence of Fourier coefficients of the π-periodic function H. For the ansatz (> .), H(ω) =
̂ (ω) φ(ω) (ω) ω r+ = = . ̂(ω) φ (ω)r+ (ω) r+ (ω)
(.)
Splines and Multiresolution Analysis
This gives criteria for the choice of the function : (i) (ii) (iii)
vanishes at the origin and there has a zero of order r+. This ensures that φ̂ ∈ L (R) and that Theorem (iii) is satisfied. (ω) is a trigonometric function, to ensure that H(ω), the so-called scaling filter, is π-periodic. has no other zeros in [ −π, π ], except at the origin. Otherwise, the autocorrelation ̂(ω + πk)∣ would vanish somewhere, and the shifts of the filter A(ω) = ∑ k∈Z ∣φ function φ would fail to generate a Riesz sequence, see Theorem (i).
A simple function ensuring all three requirements (i), (ii), and (iii) is r+
(ω) = (sin(ω/)θ(ω/))
,
where θ is a π-periodic phase factor such that ∣θ∣ = , i.e., a shift in time domain. Choosing θ(ω) = e −i ω yields the cardinal B-splines as given in (> .) resp. (> .): β̂ (ω) = ∫
e −i ωt dt =
− e −i ω sin(ω/) −i ω/ = e . iω ω/
Since β has compact support, β̂ ∈ C ∞ []. Due to the convolution formula (> .), −i ω
−e β̂m (ω) = ( iω
m+
)
sin(ω/) −i ω/ e =( ) ω/
m+
.
(.)
The β m are scaling functions of regularity r = m − , as the verification of the criteria in Theorem shows. In fact, the following holds. Let A = . Integrability: Since by (> .) the functions β̂m are L -integrable, so are the β m , m ∈ N . Riesz sequence property: The shifted characteristic functions β (x − k) = χ[,) (x − k), k ∈ Z, are clearly orthonormal. Theorem (ii) thus yields ∑ ∣ β̂ (ω + πk)∣ = almost everywhere. k∈Z
To verify the Riesz sequence property for β m , the autocorrelation filter must be bounded with strictly positive constants from above and from below. It is m m+ . ∑ ∣ β̂ (ω + πk)∣ = ∑ ∣ β̂ (ω + πk)∣ k∈Z
k∈Z
In [ −π, π ], ∣β̂ ∣ is clearly positive (cf. > Fig. -), which gives ∣β̂ (π)∣ = /π as a positive bound from below. There is a constant c, such that < c = (/π)m+ < ∣β̂ (ω)∣m+ ≤ ∑ ∣ β̂m (ω + πk)∣m+ . k∈Z
Splines and Multiresolution Analysis
1 0.8 0.6 0.4 0.2 p
–4p –3p –2p –p
2p 3p 4p
⊡ Fig. -
∧ The function ∣β ∣ is strictly positive in the interval [ −π, π ]
Since the sequence (∣ β̂ (ω + πk)∣) k∈Z ∈ l (Z) for all ω ∈ R, the same is true for the sequence ∣ β̂m (ω + πk)∣ = ∣β̂ (ω + πk)∣m+ . This yields the existence of the requested upper bound c < ∞. Thus { β̂m (● − k)} k∈Z forms a Riesz sequence in L (R). Scaling relation: The scaling filter (> .) H(ω) = −m
m+ m + −i ωk ( − e −iω )m+ −m −i ω m+ −m )e = ( + e ) = ∑( k ( − e −i ω )m+ k=
))k∈Z ∈ l (Z). Hence, the is obviously π-periodic and has Fourier coefficients (−m (m+ k B-splines satisfy the scaling relation (> .) m+
β m (x/) = ∑ −m ( k=
m+ m )β (x + k). k
For β this equation reads β (x/) = β (x) + β (x + ), which is true since β (x/) = χ[,) (x/) = χ[,) (x). This equation and examples for scaling relations of other B-splines are illustrated in > Fig. -. ̂ at the origin: From > Eq. (.), Continuity and positivity of φ m (ω)∣ = ∣ ∣ β̂
sin(ω/) ∣ ω/
m+
,
m () = . Thus we have proved which has a continuous continuation at the origin, and β̂ the following conclusion:
Theorem The cardinal B-spline β m , m ∈ N , is a scaling function of an m − -regular multiresolution analysis with dilation . The order of approximation is m + .
Splines and Multiresolution Analysis
1
1
0.8
0.8
0.6
0.6
0.4
0.4
0.2
0.2
–1
1 m=0
2
3
0.7
–1
1
2 m=1
3
4
5
0.6
0.6
0.5
0.5
0.4
0.4
0.3
0.3 0.2
0.2
0.1
0.1 2
4
6
2
m=2
4 m=3
6
8
⊡ Fig. -
Scaling relation for B-splines βm , m = , . . . , . The B-spline versions βm (x/) ∈ V− are displayed with solid lines, the scaled translates in V are depicted dashed. The sum of the dashed functions gives the B-spline at the lower scale β m (x/)
Note that the cardinal B-splines β m with Fourier transform of the form (> .) are scaling functions, but they are not yet orthonormalized: The family {β m (● − k)} k∈Z spans V and is a Riesz basis, but it is not an orthonormal basis of V . Orthonormality can be achieved with Theorem and > Eq. (.): β̂m (ω) ̂m (ω) := √ B . ∑ k∈Z ∣ β̂m (ω + πk)∣ Figure - shows some orthonormalized B-spline scaling functions and the corresponding wavelets. >
..
Polyharmonic B-Splines
The same approach to model scaling functions in Fourier domain can be done in higher dimensions. We aim at constructing a scaling function for a multiresolution analysis of
Splines and Multiresolution Analysis
1.2 0.5
1 0.8 –4
0.6
–4
–2
–2
0.4
–0.5
0.2
–1 2
–0.2
4
2
4
2
4
2
4
–1.5 m=1
1
1
0.8
0.5
0.6 0.4
–4
–2
0.2 –4
–0.5
–2
2
4
–1
–0.2 m=2 1
0.5
0.8 0.6 –4
0.4
–2 –0.5
0.2 –4
–2
2
–1
4
–0.2 m=3
⊡ Fig. -
Orthonormalized B-splines Bm (left column) and corresponding wavelets (right column) for m = , , in time domain. Note that the orthonormalized B-splines and wavelets do not have compact support. Due to the orthonormalization procedure (> .) the orthonormalized B-spline is an infinite series of shifted compactly supported B-splines
L (Rn ) of the form ̂(ω) = φ
(ω) , ∥ω∥r
r ∈ N,
r > n/,
x ∈ Rn .
With an appropriate trigonometric polynomial n
r
(ω) = ( ∑ sin (ω k /)) , k=
ω = (ω , . . . , ω n ),
Splines and Multiresolution Analysis
̂ is a nonseparable scaling function for a multiresolution analysis of L (Rn ) with respect to φ dilation matrices A that are scaled rotations. The corresponding function in space domain φ is called elementary r-harmonic cardinal B-spline, or short polyharmonic B-spline P r . This terminology can be justified as follows. The Fourier transform in the sense of tempered distributions of the function /∥ω∥r is indeed a polynomial – up to a logarithmic factor for r − n even. In fact, in S ′ (Rn ), F − (/∥ ● ∥r )(x) = ∥x∥r−n (A(n, r) ln ∥x∥ + B(n, r)) =: ρ(x), with constants A(n, r), B(n, r) as given in [, Chap. , Sect. ], and A(n, r) = except for r − n even. (Note that for r > n/ on the right-hand side the final parts have to be considered.) The term polyharmonic comes from the fact that ρ is the Green function of the r-iterated Laplace operator Δr . However, with these considerations φ(x) = P r (x) = ∑
k ρ(x
+ k) almost everywhere
k∈Z
becomes an nD spline. Here ( k ) k∈Z is the sequence of Fourier coefficients of . Due to the decay in Fourier domain, the polyharmonic B-spline P r has continuous derivatives D β for multi-indices ∣β∣ < r − n. In the same way as for the B-splines, it can be shown with the theorems given in > Sect. .. that φ forms indeed a scaling function with approximation order r [, , ]. > Figure - shows the polyharmonic B-spline scaling function in space domain and in frequency domain.
2 1 1.5
0.8 0.6
1
0.4
0.5
0.2 0 –2
2 1
0 2
2
0 –2
2
0
0
–1 –2
–2
⊡ Fig. -
Polyharmonic B-spline for r = in space domain (left) and frequency domain (right)
Splines and Multiresolution Analysis
.
Survey on Mathematical Analysis of Methods
There are many function families that consist of piecewise polynomials and that are called splines, which in addition fulfil a multiresolution condition in the one or the other sense. These families can be classified by various aspects, e.g., by their dimensionality, by the lattice which is invariant under the corresponding dilation matrix, by the geometries they are defined on, or whether they provide phase information or not, and so on. The following sections list some of these spline approaches and illustrates their mathematical properties and features.
..
Schoenberg’s B-Splines for Image Analysis – the Tensor Product Approach
As mentioned in Example , multiresolution analyses for L (Rn ) and dilation matrix Id can be generated from tensor products of D dyadic multiresolution analyses. To analyze images with B-splines, the tensor product β m (x)β m (y) of B-splines is a scaling function for L (R ) and the dilation matrix A = Id. Since in D the determinant det A = , the corresponding details space W is spanned by three wavelet functions: ψ(x)β m (y),
β m (x)ψ(y),
ψ(x)ψ(y),
x, y ∈ R.
(.)
A drawback of this approach is the fact that these wavelets prefer horizontal, vertical, and diagonal directional features and are not sensitive to other directions, see > Fig. -. For the analysis of images with many isotropic features, the use of isotropic or steerable wavelets is recommended. However, the tensor approach is a simple and widely used wavelet approach. For an illustration of the respective image decomposition in coarse approximations and details of various sizes, see > Fig. -.
⊡ Fig. -
The three B-spline tensor wavelets (> .) show a preference for horizontal, vertical, and diagonal directions. Here, m = . Minimal resp. maximal function values are given in black resp. white
Splines and Multiresolution Analysis
a
b
c
d
⊡ Fig. - Decomposition of an image [, Part of IM.tif] into coarse approximations and details of various sizes. (a) Original image. (b) Matrix of the absolute values of the multiresolution coefficients. Large coefficients are white. The wavelet coefficients are depicted in the lower two and the upper right band of each scale; the approximation coefficients in the upper left band. (c) Two steps of themultiresolution decomposition. From left to right: Finest details and second finest details, remaining approximation of the original image. (d) Second finest details (c, center) split into the contribution of the three wavelets. From left to right: The decomposition into horizontal, vertical and diagonal details
..
Fractional and Complex B-Splines
The B-splines as described up to now have a discrete order of smoothness, i.e., they are C n -functions with n ∈ {−, , , , . . .}. For some applications, e.g., in medical imaging, where the order of smoothness of certain image classes is fractional, it would be favorable
Splines and Multiresolution Analysis
to have a spline and wavelet basis that is adaptable with respect to this regularity [, , ]. A first step in this direction was done by T. Blu and M. Unser, who proposed B-splines and wavelets of fractional orders []. They defined two variants of fractional B-splines, the causal ones and the symmetrical ones. The causal fractional B-spline is generated by applying the (α+)th fractional difference operator to the one-sided power function t+α : β+α (t) :=
α+ )(t − k)+α . Δ α+ t α = ∑ (−) k ( k Γ(α + ) + + Γ(α + ) k≥
The Fourier-transform representation is similar to the one of the classical B-splines (cf. > Eq. .): −i ω
−e β̂+α (ω) = ( iω
α+
)
.
Here again, the smoothness property β+α ∈ C m,γ (R) in time domain is gained by the fractional polynomial decay of order O(∣ω∣∣α+ ) in frequency domain. Note that C m,γ (R) denotes the Hölder space with exponent m = ⌊α + ⌋ and γ = α + − m, i.e., the space of m-times continuously differentiable functions f that are Hölder-regular with exponent < γ ≤ such that there is a constant C > with ∣D m f (t) − D m f (s)∣ ≤ C∣t − s∣ ∀s, t ∈ R. Although the fractional B-splines are not compactly supported, they decay in the order O(∣t∣−(α+) ) as t → ∞. They are elements of L (R) for α > , of L (R) for α > − and of the Sobolev spaces Wr (R) for r < α + . They share many properties with their classical B-spline relatives, such as the convolution property, their relation to difference operators, i.e., they are integral kernels for fractional difference operators [], and they are scaling functions for dyadic multiresolution analyses. This can be verified by the procedure given in > Sect. ... The causal fractional B-spline is not symmetric. Since for some signal and image analysis tasks symmetrical bases are preferred, in [] the symmetrical fractional B-splines β∗α are proposed. They are defined in Fourier domain as follows: −i ω
−e β̂∗α (ω) = ( iω
α+
)
α+
+ eiω sin(ω/) ) =∣ ∣ ( −iω ω/
α+
,
(.)
and therefore obviously are symmetrical in time domain. The same regularity and decay properties apply as for the causal fractional B-splines. The symmetrical fractional B-splines are also piecewise polynomials, as long as α ∉ N . For even integer degrees, the singularity introduced through the absolute value in > Eq. (.) causes that β∗m is a sum of integer shifts of the logarithmic term ∣t∣m ln(t) for m ∈ N . For the explicit time-domain representation and further details on these splines cf. []. In [], Blu and Unser defined another variant, the generalized fractional B-spline or (α, τ)-fractional spline β τα with a parameter τ ∈ R. Also these splines are defined via their
Splines and Multiresolution Analysis
Fourier domain representation: −i ω
−e β̂τα (ω) = ( iω
α+ +τ
)
α+
− eiω ) ( −iω
−τ
.
As above, the parameter α > controls the regularity of the splines. The parameter τ, in contrast, controls the position of the splines with respect to the grid Z. This can be justified by the following fact. All variants of the B-splines considered in this section converge to optimally time-frequency localized functions in the sense of Heisenberg, i.e., to Gaussians or Gabor functions, if the degree α becomes large. For a proof for the classical cardinal B-splines, see []. In the case of the (α, τ)-fractional splines [],
β τα (t) = O (e − α+ (t−τ) )
for α → ∞.
This explains the notion “shift parameter” for τ. Moreover, the parameter τ allows to interpolate the spline family between the two “knots,” the symmetrical ones (τ =) and ), see > Fig. -. Both parameters α and τ can be tuned the causal ones (τ = α+ independently and therefore allow for an individual adjustment of the analyzing basis. Another generalization are the complex B-splines []. There are two variants, both defined via their Fourier domain representation. Let z = α + iγ ∈ C, α > − , γ ∈ R and y ∈ C. Then −i ω
−e β̂z (ω) = ( iω
−i ω
−e β̂zy (ω) = ( iω
z+
) )
, z+ −y
− eiω ) ( −iω
z+ +y
are complex B-splines of complex degree z. The functions are well defined, because the −i ω function Ω(ω) = −ei ω never touches the negative real axis such that Ω(ω)z is uniquely 1
1
0.5
0.5
0
0
–0.5 –2
–1
0
1
2
3
4
–0.5 –2
–1
0
1
2
3
4
⊡ Fig. - The fractional (α, τ)-splines interpolate the families of the causal and the symmetric k for k = , , , , from the most right (causal) to the most left fractional splines. τ = α+ (symmetrical) function in each image. Right: α = .. Left: α =
Splines and Multiresolution Analysis
defined. β z and β zy are elements of the Sobolev spaces Wr (R) for r < α + . β z has the time-domain representation β z (t) =
k z+ )(t − k)+z , ∑ (−) ( k Γ(z + ) k≥
i.e., β z is a piecewise polynomial of complex degree. For more details on the properties of these families of complex splines cf. []. The idea behind the complex degree is as follows: The real part Re z = α operates as regularity and decay parameter in the same way as for the fractional B-splines. The imaginary part, however, causes an enhancement resp. damping of positive or negative frequencies. In fact, β̂z (ω) = β̂+α (ω)e −iγ ln ∣Ω(ω)∣ e γ arg Ω(ω) . The imaginary part γ of the complex degree introduces a phase and a scaling factor in frequency domain. The frequency components on the negative and positive real axis are enhanced with different sign, because arg Ω(ω) ≥ for negative ω and arg Ω(ω) ≤ for positive ω. > Figure - illustrates this effect. With real-valued functions, only symmetric spectra can be analyzed. The complex B-splines, however, allow for an approximate analysis of the positive or the negative frequencies, because the respective symmetric bands are damped due to the complex exponent. However, the complex B-splines inherit many properties of their classical and fractional relatives. All of the generalized B-spline families mentioned in this section have in common that they are scaling functions of dyadic multiresolution analyses. They are one-dimensional functions, but with the tensor approach mentioned in Example and > Sect. .. they are also suitable for image processing tasks. Although the fractional and the complex splines, in general, to not have compact support, they allow for fast analysis and synthesis algorithms. Due to their closed form in Fourier domain, they invite for an implementation of these algorithms in Fourier domain.
1 0.8
1.25
0.6
1
0.4 0.2
–15
–10
–5
⊡ Fig. -
6
1.5
5
10
15
–15
–10
–5
5 4
0.75
3
0.5
2
0.25
1 5
10
15
–15
–10
–5
5
10
∧ The frequency spectrum ∣β z ∣ for z = + iγ, γ = , , (from left to right). The spectrum of β = β+ is symmetric (right), whereas the spectra of β +i (center) and β +i (left) show an enhancement of the positive frequency axis
15
..
Splines and Multiresolution Analysis
Polyharmonic B-Splines and Variants
In > Sect. .. the so-called polyharmonic B-splines in Rn were introduced. They are defined in Fourier domain by the representation r
n ̂r (ω) = ( ∑ k= sin (ω k /) ) , P n ∑ k= ω k
r > n/, ω = (ω , . . . , ω n ).
These polyharmonic B-splines satisfy many properties of the classical Schoenberg splines; e.g. they are piecewise polynomial functions, they satisfy a convolution relation P r +r = P r ∗P r , they are positive functions, etc. However, they do not share the property that they converge to optimally space-frequency localized Gaussians as r increases []. This is due to the fact that the trigonometric polynomial in the numerator regularizes insufficiently at the origin: The second order derivative of ∑nk= sin (ω k /) n ∑ k= ω k is not continuous. Van De Ville et al. [] therefore proposed another localizing trigonometric polynomial: n
μ(ω) = ∑ sin (ω k /) − k=
n− n ∑ ∑ sin (ω k /) sin (ω l /). k= l=k+
(.)
A new function family then is defined in Fourier domain via r
̂r (ω) = ( μ(ω) ) , Q ∥ω∥
r>
n .
(.)
Qr is called isotropic polyharmonic B-spline. The function is piecewise polynomial (except for r − n even, where a logarithmic factor has to be added, see below), and shares with the polyharmonic splines their decay properties in Fourier domain and their regularity properties in space domain. Qr converges to a Gaussian as r increases, which makes the function family better suitable for image analysis than P r , because of the better spacefrequency localization. This effect is due to higher order rotation invariance or isotropy of the localizing trigonometric polynomial (> .): (ω) = + O(∥ω∥ ) vs.
μ(ω) = +
∥ω∥ + O(∥ω∥ ) as ∥ω∥ → .
̂r has a second order moment, and thus the central limit theorem can be This causes that Q ̂r has a higher applied to proof the convergence to the Gaussian function. In addition, Q ̂r ; therefore Qr decays faster. For a complete proof of the localization regularity than P property, see []. An example of the isotropic polyharmonic spline is given in > Fig. -. The polyharmonic B-splines and the isotropic polyharmonic B-splines both are realvalued functions. The isotropic B-spline is approximately rotation-invariant and therefore is suited for image analysis of isotropic features. A complex-valued variant of these B-splines in D was introduced in []. The idea is to design a spline scaling function that
Splines and Multiresolution Analysis
2 1.5 1 0.5 0 –2 0 2
2
0
–2
⊡ Fig. -
The isotropic polyharmonic spline Q . Compare with polyharmonic B-spline P
>
Fig. - of the classical
is approximately rotation covariant, instead of rotation invariant. Rotation covariant here means that the function intertwines with rotations up to a phase factor. Again, the design of the scaling function is done in Fourier domain, now using a perfectly rotation-covariant function ̂ ρ r,N (ω , ω ) =
(ω
+
r ω )
, (ω − iω ) N
where r ≥ and N ∈ N. In fact, for some rotation matrix R θ ∈ GL(, R), ̂ ρ r,N (R θ ω) = ρ r,N (ω). For localizing the function ̂ ρ r,N , the same trigonometric polynomials and e −i N θ ̂ μ as above can be used. The corresponding complex polyharmonic B-splines are then defined in frequency domain as N
̂ r,N (ω , ω ) = R
( (ω , ω ))r+ , r (ω + ω ) (ω − iω ) N
̂ r,N R μ (ω , ω ) =
(μ(ω , ω ))r+ . r (ω + ω ) (ω − iω ) N
or as
N
(.)
The case N = yields the real-valued polyharmonic splines. There are also other trigonometric polynomials that are suitable as localizing numerator for the real and the complex polyharmonic B-splines. With an appropriate choice the features of the polyharmonic splines can be tuned []. However, for both the real and the complex variant, the localizing multiplier has to fulfil moderate conditions to make the respective polyharmonic B-splines being a scaling function. In D, the following result holds (cf. []):
Theorem such that
Splines and Multiresolution Analysis
Let r > and N ∈ N . Let η(ω , ω ) be a bounded, πZ -periodic function, N
∣
(η(ω , ω ))r+ ∣ r (ω + ω ) (ω − iω ) N
is bounded in a neighborhood of the origin, and such that η(ω , ω ) ≠ for all (ω , ω ) ∈ [−π, π] /{(, )}. N ̂ = η r+ ⋅ ̂ ρ is the Fourier transform of a scaling function φ which generates a Then φ a b multiresolution analysis . . . V− ⊂ V ⊂ V . . . of L (R ) with dilation matrix A = ( ), −b a a, b ∈ Z: Vj = span {∣ det A∣ j/ φ(A j ● −k), k ∈ Z }. From the Fourier domain representation immediately follows that φ̂ ∈ L (R ) for r+ N > ̂ decays as ∣φ ̂(ω , ω )∣ = O(∥(ω , ω )∥−r−N ) when ∥(ω , ω )∥ → ∞. and that φ As a result, for all three variants of the polyharmonic B-splines, the classical ones P r , the isotropic ones Qr , and the complex ones Rr,N , the following properties hold: They are scaling functions for multiresolution analysis. Their smoothness parameter r can be chosen fractional and must fulfill r + N > n/ for integrability reasons. Then the scaling function in space domain is element of the Sobolov space φ ∈ Ws (R ) for all s < r + N − . The explicit space domain representation is φ(x) = ∑ η k ρ r,N (x + k) k∈Z N
for almost all x ∈ R . Here (η k ) k∈Z denotes the Fourier coefficients of η r+ . ρ r,N is the inverse Fourier transform of the Hadamard final part P f (̂ ρ r,N ) ∈ S ′ (R ). In fact, for r ∉ N , r− ρ r,N (x , x ) = c (x + x ) (x + ix ) N and for r ∈ N , ρ r,N (x , x ) = c (x + x )
r−
(x + ix ) N (ln π
√
x + x ) + c )
with appropriate constants c , c , c ∈ C. This justifies the notion spline for the function families. They all have a closed form in frequency domain. As in the case of the D cardinal B-splines there is a fast analysis and synthesis algorithm using the frequency domain representation and the fast Fourier transform, cf. > Sect. ..
..
Splines on Other Lattices
...
Splines on the Quincunx Lattice
The tensor product of two D dyadic multiresolution analyses yields a D multiresolution analysis with dilation matrix A = Id, cf. Example . As a consequence, the scaling factor while moving from one approximation space V to the next coarser space V is ∣ det A∣ = .
Splines and Multiresolution Analysis
6
6
6
5
5
5
4
4
4
3
3
3
2
2
2
1
1
1
a
1
2
3
4
5
6
b
1
2
3
4
5
6
c
1
2
3
4
5
6
⊡ Fig. -
Three iterations of the quincunx subsampling scheme. (a) Z , (b) Aq Z , (c) Aq Z . The thinning of the Z -lattice using the dilation matrix Aq is finer than dilation with the dyadic matrix A = Id, which in one step leads from (a) to (c)
For some image processing applications, especially in medical imaging, this scaling step size is too large. A step size of as in the D case would be preferred. Moreover, the decomposition of the wavelet space into three subspaces then would be avoided, and the eventual problematic of the directionality of the three involved wavelets would not arise. An example of a dilation matrix satisfying these requirement is the scaled rotation matrix Aq = ( −
)
with det A q = . It leads to the so-called quincunx lattice. This lattice is generated by applying A q to the cartesian lattice. It holds A q Z ⊂ Z , see > Fig. -. Since A q falls into the class of scaled rotations, the polyharmonic B-Spline construction including all variants are applicable for this case, cf. Theorem . Note that the tensor product approach in general is not suitable for the quincunx subsampling scheme.
...
Splines on the Hexagonal Lattice
Images as D objects are normally processed on the cartesian lattice, i.e., the image pixels are arranged on a rectangular grid. For image processing this arrangement has the drawback that not all neighbors of a pixel have the same relation: The centers of the diagonal neighbors have a larger distance to the center pixel than the adjacent ones. A higher degree of symmetry has the hexagonal lattice. It is therefore ideal for isotropic feature representation. The hexagonal lattice gives an optimal tesselation of R in the sense of the classical honeycomb conjecture, which says that any partition of the plane into regions of equal area has a perimeter at least that of regular hexagonal tiling []. The better isotropy of the hexagonal lattice is attractive for image analysis and has led to a series of articles on image processing methods (e.g., [, , ]) as well as on applications (e.g., [, , , ]).
Splines and Multiresolution Analysis
The hexagonal lattice is generated by applying the matrix √ Rh =
√ (
√/ ) /
on the cartesian lattice Λ h = R h Z . A scaling function of a multiresolution analysis of L (R ) in the hexagonal lattice fulfils all properties of Definition , but the last two. They change to (v) (vi)
f ∈ V ⇐⇒ f (● − R h k) ∈ V for all k ∈ Z . There exists a scaling function φ ∈ V , such that the family {φ(● − R h k)} k∈Z of translates of φ forms a Riesz basis of V .
Let A be a dilation matrix which leaves the hexagonal lattice invariant AΛ h ⊂ Λ h . Then A is of the form [] A = R h BR − h with B ∈ GL(, R) having only integer entries and with eigenvalues strictly larger than one. > Fig. - gives an example of two subsampling steps on the hexagonal lattice. There are several possible approaches to define spline functions on the hexagonal lattice. Sablonnière and Sbibih [] proposed to convolve piecewise linear pyramids to generate higher order B-splines. Van De Ville et al. [] started with the characteristic function of one hexagon and also used an iterative convolution procedure to construct B-splines of higher degree. However, both approaches lead to discrete order hexagonal Bsplines. If A is a scaled rotation, then fractional and complex B-splines on the hexagonal lattice can be defined in an analog way as in > Sect. .. for polyharmonic splines and
6
6
6
5
5
5
4
4
4
3
3
3
2
2
2
1
1
1
a
1
2
3
4
5
6
b
1
2
3
4
5
6
c
1
2
3
4
5
⊡ Fig. - Three iterations of subsampling of the hexagonal grid with A = Rh BR− h , where B = (− ). (a) Λh , (b) AΛh = Rh BR− h Λh = Rh BZ , (c) A Λh
6
Splines and Multiresolution Analysis
their (complex) variants []. Consider again the perfectly rotation-covariant (or for N = rotation-invariant) function ̂ ρ r,N (ω , ω ) =
, r (ω + ω ) (ω − iω ) N
where r > and N ∈ N . The idea is now to use a hexagonal-periodic trigonometric polynomial for localizing this function and to eliminate the singularity at the origin. Condat et al. [] proposed √ √ η h (ω , ω ) = √ ( − (cos (/ (−ω / + ω ) / ) √ √ √ + cos (/ (ω / + ω ) / ) + cos (−/ ω ))) , and defined the elementary polyharmonic hexagonal rotation-covariant B-spline via its frequency domain representation as ̂ r,N (ω , ω ) = R h
N
(η h (ω , ω ))r+ . r (ω + ω ) (ω − iω ) N
The B-spline in space domain then has the representation Rr,N h (x) = ∑ η h,k ρ r,N (x − R h k). k∈Z
r+ N
Here, (η h,k ) k∈Z denotes the sequence of Fourier coefficients of η h . > Figure - shows the localizing trigonometric polynomial η and the Fourier spectra of two hexagonal splines. For N = , the functions are the elementary polyharmonic hexagonal rotation-invariant B-splines. They are real-valued functions that converge to a Gaussian, as r → ∞, and therefore are well localized in the space domain as well as in the frequency domain. For N ∈ N, the splines are complex-valued functions and approximately rotation covariant in a neighborhood of the origin: ̂ r,N (ω) = e i N arg(ω) ( + C∥ω∥ + O(∥ω∥ )) R h
for ω → ,
where ω = (ω , ω ) and C ∈ R a constant. The translates of the complex B-spline Rr,N h form a Riesz basis of the approximation spaces j Vj = span {Rr,N h (A x − R h k) : k ∈ Z }
L (R )
,
j ∈ Z.
The ladder of spaces (Vj ) j∈Z generates a multiresolution analysis of L (R ) for the hexagonal grid and for scaled rotations A. Also in this case, the implementation of the analysis and synthesis algorithm can be elegantly performed in frequency domain [, ].
Splines and Multiresolution Analysis
4
10
2 0 –10
5 0 –5 –5
0 5
a
10
1 0.75 0.5 0.25 0 –10
10 5 0 –5 5
b
1 0.75 0.5 0.25 0 –10
10
10 5 0 –5
–5
0
–10
–5
0 5
–10
c
10
–10
⊡ Fig. -
∧ Localization of ρ with (a) the hexagonal-periodic trigonometric polynomial ηh yields the elementary polyharmonic hexagonal rotation-covariant B-spline. Frequency spectrum of ∧ r,N ∧ r+N/, ), (b) for r + N = , and (c) for r + N = . Rh (or equivalently of Rh
.
Numerical Methods
For the illustration of the numerical method, we focus on the quincunx dilation matrix A=( −
),
(.)
and consider the polyharmonic spline variants in D as defined in > Sect. ... Since det A = , the generators of the multiresolution space and the corresponding wavelet space are two functions: The scaling function (here the variant of the polyharmonic spline), and the associated wavelet. We consider the scaling function φ(x) = Qrμ (x) resp. φ(x) = Rr,N μ (x). It spans the ladder of nested approximation spaces {Vj } j∈Z via Vj = span {∣ det A∣ j/ φ(A j ● −k), k ∈ Z }
L (R )
,
j ∈ Z.
Splines and Multiresolution Analysis
Denote ̂(ω + πk)∣ M(ω) := ∑ ∣φ
(.)
k∈Z
the autocorrelation filter. It is bounded < M(ω) ≤ C for some positive constant C []. The scaling functions can be orthonormalized applying the procedure given in Theorem (iii): ̂(ω) φ ̂ . Φ(ω) =√ M(ω) The scaled shifts of Φ span the same spaces {Vj } j∈Z . The B-splines as scaling functions φ satisfy a refinement relation φ(A− x) = ∑ h k φ(x − k) almost everywhere and in L (R ). k∈Z
This relation in fact is a discrete convolution. The Fourier transform yields a πZ -periodic function H ∈ L (T ) of the form H(e i ω ) = ∣ det A∣ ⋅ =
̂(AT ω) ρ̂r,N (AT ω)ζ(AT ω) φ = ∣ det A∣ ⋅ ̂(ω) ̂ φ ρ r,N (ω)ζ(ω)
ζ(AT ω) ⋅ , r− ζ(ω) (a + b ) (a − ib) N
ω ∈ R .
N
Here, ζ(ω) = (μ(ω))r+ is the localizing multiplier for the isotropic polyharmonic Bspline in the case N = , and for the rotation-covariant polyharmonic B-splines in the case N ∈ N, cf. (> .), (> .), and (> .). The wavelet function ψ spanning a Riesz basis for the orthogonal complement Wj = span { j/ ψ(A j ● −k), k ∈ Z }
L (R )
in Vj+ = Wj ⊕Vj can also be gained in frequency domain. For the quincunx dilation matrix A as in (> .) a wavelet (or sometimes called a prewavelet, since the functions are not yet orthonormalized) is given by ̂(ω) = G(e i ω )φ ̂(ω) = e −i ω H(ω + (π, π)T )M(ω + (π, π)T )φ ̂(ω), ψ compare with family
>
Sect. .... The (pre-)wavelet Riesz basis for Wj is then given by the {ψ j,k = j/ ψ(A j ● −k), k ∈ Z }.
This basis in general is not orthonormal: ⟨ψ j,k , ψ j,l ⟩ ≠ δ k,l . A function f ∈ L (R ) can then be represented by the series f =
̃j,k ⟩ψ j,k = ∑ ⟨f,ψ j∈Z,k∈Z
̃j,k ⟩ψ j,k , ∑ ⟨f,ψ j∈Z,k∈Z
̃ j,k } k∈Z denotes the dual basis for each j ∈ Z: ⟨ψ ̃ j,k , ψ j,l ⟩ = δ k,l . Its generator in where {ψ frequency domain is
Splines and Multiresolution Analysis
̂(ω) M(ω + (π, π)T ) φ ̂ ̃(ω) = e −i ω H(ω + (π, π)T ) . ψ T M(A ω) M(ω) In contrast, the formula M N M(ω + (π, π)T ) T ̂ ̂(ω) Ψ(ω) =N ψ M(AT ω) generates an orthonormal wavelet basis. It corresponds to the orthonormal basis of V generated by the integer shifts of the orthonormalized scaling function Φ. These considerations show that there are three variants of a multiresolution implementation: An “orthonormal” one with respect to the orthonormalized scaling functions and corresponding orthonormal wavelets, one with the B-splines on the analysis side, ̃ j,k ̃ j,k + ∑ ⟨ f , ψ j,k ⟩ψ f = ∑ ⟨ f , φ j,k ⟩φ k∈Z
for f ∈ Vj+ ,
k∈Z
and finally one with the B-splines on the synthesis side: ̃j,k ⟩ψ j,k ̃ j,k ⟩φ j,k + ∑ ⟨ f , ψ f = ∑ ⟨f, φ k∈Z
for f ∈ Vj+ .
k∈Z
Both, the scaling filters H(e i ω ) and the wavelet filters G(e i ω ) = e −i ω H(ω + (π, π)T )M(ω + (π, π)T ) as well as their orthogonal and dual variants in our case are nonseparable and infinitely supported. Therefore, a spatial implementation of the decomposition would lead to truncation errors due to the necessary restriction to a finite number of samples. However, because of the closed form of H and therefore of G, the corresponding multiresolution decomposition or wavelet transform can be efficiently implemented in frequency domain. The respective image first undergoes an FFT, then is filtered in frequency domain by multiplication with the scaling filter H and the wavelet filter G. This method automatically imposes periodic boundary conditions. The coefficients resulting from the high pass filtering with G are the detail coefficients. They are stored, whereas the coefficients resulting from the low pass filtering H are reconsidered for the next iteration step. ̃ j,k . ̃ j+,k = ∑ ⟨ f , φ j,k ⟩φ ̃ j,k + ∑ ⟨ f , ψ j,k ⟩ψ ∑ ⟨ f , φ j+,k ⟩φ k∈Z
k∈Z
k∈Z
For details and tricks of the frequency domain implementation, cf. [, , , ]. > Figure - shows the multiresolution decomposition for the scaling function φ = R, μ . There it was assumed that the image is bandlimited and projected on the space V , which has the advantage that the coefficients do not depend on the chosen flavour of the scaling function, i.e., orthogonal, B-spline or dual. Qualitatively the transform is very similar to a multiscale gradient with the real part corresponding to the x -derivative and the imaginary part corresponding to the x -derivative [].
Splines and Multiresolution Analysis
a
b
c
d
⊡ Fig. -
Decomposition of an image [, Part of IM.tif] into approximation and wavelet coefficients. (a) Original image. (b) Matrix of the absolute values of the multiresolution coefficients. Large coefficients are white. The approximation coefficients in the upper left band, the other bands are wavelet coefficients on six scales. (c) Real part of the coefficient matrix, (d) imaginary part of the decomposition matrix for φ = R, μ and the corresponding wavelets. The coefficients had their intensity rescaled for better contrast
.
Open Questions
In this chapter a method for the construction of spline multiresolution bases was described. It yields a nice variety of new bases with several parameters for adaption and tuning. In the last decade, the notion of compressive sampling or compressed sensing arose, which is footing on the existence of well-adaptable bases. In fact, the idea behind compressed sensing is that certain functions have a sparse representation, if the underlying basis is smartly chosen. In this case, the function can be reconstructed from very few samples because of
Splines and Multiresolution Analysis
the prior knowledge of sparsity in this underlying basis. As a consequence, the knowledge on sparsity allows to sample such a signal at a rate significantly under the Nyquist rate. (The Shannon–Nyquist sampling theorem says that a signal must be sampled at least two times faster than the signal’s bandwidth to avoid loss of information.) In the last years, and virtually explosively in the last years, many important theoretical results were proven in this field, in particular by D. Donoho, E. Candès, J. Romberg, and T. Tao. For an introduction and references on compressed sensing, see e.g., [, ] and the website []. Compressed sensing is based upon two fundamental concepts: that of incoherence and that of sparsity. Let {x i } i=, ... ,N be an orthonormal basis of the vector space V. Let f = N ∑ i= s i x i with s i = ⟨ f , x i ⟩. The signal f is called k-sparse, if only k of the coefficients are nonzero, k ∈ N. A general linear measurement process for signals consists in computing M < N inner products y j = ⟨ f , y j ⟩ for a collection of vectors {y j } j . In matrix form, g = Y f = Y Xs, where Y and X are the matrixes with {y i } i and {x j } j as columns, and Y X is an M × N-matrix. If the function families Y and X are incoherent, i.e., if the incoherence measure √ √ μ(Y, X) = N max ∣⟨y i , x j ⟩∣ ∈ [, N] ≤i, j≤N
is close to one, then under mild additional conditions the k-sparse signal f can be reconstructed from M > Cμ (Y, X)k ln N samples with overwhelming probability. Wavelet bases have proven to be very suitable for compressed sensing. It is an open question to classify the signals from certain applications, and to estimate in which appropriate B-spline basis they have a k-sparse representation. Then adequate bases and function families incoherent with the spline bases have to be identified. In the last years, the concept of sparsity entered image processing. It has proven to help immensely to accelerate the solution of inverse problems and reconstruction algorithms, e.g., in medical imaging, such as in magnetic resonance imaging [], computed tomography [], photo-acoustic tomography [], tomosynthesis [], and others. In this area, as well as in other fields of imaging, it can be expected that the combination of splines – due to their easy modeling and the fast frequency domain algorithms – multiresolution and wavelets, and sparsity will lead to novel impressing fast algorithms for image reconstruction.
.
Conclusion
In the design procedure for scaling functions of multiresolution analyses, regularity and decay features, as well as symmetry properties can be tuned by an appropriate modeling in frequency domain. The idea is to start in frequency domain with a polynomial function P that fulfils the required symmetry features, and that has a degree, such that /P decays
Splines and Multiresolution Analysis
sufficiently fast. This assures that the resulting scaling function has the desired regularity. However, /P in general is not an L -function and has to be multiplied with a localizing trigonometric polynomial that eliminates the zeros in the denominator such that P becomes square integrable. The choice of this trigonometric polynomial has to be taken carefully to be compatible with the required features modeled in /P. Then under mild additional conditions the fraction ̂= φ P is the scaling function of a multiresolution analysis. This construction can be performed for D and higher dimensional spaces likewise. In time domain, the resulting scaling function is a piecewise polynomial, thus a spline. This design procedure for scaling functions unites the concepts of splines and of multiresolution. Interestingly, the polynomial in the denominator can be of a fractional or a complex degree and therefore allows a fine tuning of the scaling function’s properties. However, the scaling function then becomes an infinite series of shifted (truncated) polynomials. The numerical calculation with the approximating basis of the multiresolution analysis in time domain would cause truncation errors, which is unfavorable. But due to the construction of φ in frequency domain, and due to the closed form there, the implementation in frequency domain with periodic boundary conditions yields a fast and stable multiresolution algorithm suitable for image analysis tasks.
.
Cross-References
Astronomy Compressive Sensing > Gabor Analysis for Imaging > Neighborhood Filters and Local D Recovery > Sampling Methods > >
References and Further Reading . Aldroubi A, Unser MA (eds) () Wavelets in medicine and biology. CRC Press, Boca Raton . Baraniuk R () Compressive sensing. IEEE Signal Process Mag ():–, . Bartels RH, Bealty JC, Beatty JC () An introduction to splines for use in computer graphics and geometric modeling. Morgan Kaufman, Los Altos . Battle G () A block spin construction of ondelettes. Part : Lemarié functions. Commun Math Phys :–
. Blu T, Unser M () The fractional spline wavelet transform: definition and implementation. In: proceedings of the th IEEE International Conference on Acoustics, Speech, and Signal Processing (ICASSP’), vol , Istanbul, Turkey, – June, , pp – . Blu T, Unser M (). A complete family of scaling functions: the (α, τ)-fractional splines. In: Proceedings of the th International Conference on Acoustics, Speech, and Signal Processing (ICASSP’), vol , Hong Kong SAR,
.
.
.
.
.
.
.
.
. .
. . . .
.
.
Splines and Multiresolution Analysis
People’s Republic of China, April –, , pp – Buhmann MD () Radial basis functions: theory and implementations. Cambridge monographs on applied and computational mathematics. Cambridge University Press, Cambridge de Boor C, Höllig K, Riemenschneider S () Box splines, vol of Applied mathematical sciences. Springer, New York de Boor C, Devore RA, Ron A () Approximation from shiftinvariant subspaces of L(Rd). Trans Am Math Soc ():– Burt PJ, Adelson EH (Apr ) The Laplacian pyramid as a compact image code. IEEE Trans Commun ():– Candès EJ, Wakin MB () An introduction to compressive sampling. IEEE Signal Process Mag ():– Champeney DC () A handbook of Fourier theorems. Cambridge University Press, Cambridge Chen H- () Complex harmonic splines, periodic quasi-wavelets, theory and applications. Kluwer Academic, Dordrecht Choi JY, Kim MW, Seong W, Ye JC () Compressed sensing metal artifact removal in dental ct. In: Proceedings of IEEE International Symposium on Biomedical Imaging (ISBI), June– July , Boston, pp – Christensen O () An introduction to frames and riesz bases. Birkhäuser, Boston Christensen O () Frames and bases: an introductory course (applied and numerical harmonic analysis). Birkhäuser, Boston Chui CK () Multivariate splines. SIAM, Philadelphia Chui C () Wavelets—a tutorial in theory and practice. Academic, San Diego Chui CK (ed) () Wavelets: a tutorial in theory and applications. Academic, Boston Condat L, Forster-Heinlein B, Van De Ville D () A new family of rotation-covariant wavelets on the hexagonal lattice. In: SPIE Wavelets XII, Aug , San Diego Condat L () Image database. Online Ressource. http://www.greyc.ensicaen.fr/_ lcondat/imagebase.html (Version of Apr ) Dahmen W, Kurdila A, Oswald P (eds) () Multiscale wavelet methods for partial differential equations, vol of Wavelet
.
. .
.
.
.
.
.
.
. .
. .
.
.
analysis and its applications. Academic, San Diego Daubechies I () Ten lectures on wavelets. Society for Industrial and Applied Mathematics, Philadelphia Dierckx P () Curve and surface fitting with splines. McGraw-Hill, New York Feilner M, Van De Ville D, Unser M () An orthogonal family of quincunx wavelets with continuously adjustable order. IEEE Trans Image Process ():– Forster B, Blu T, Unser M () Complex B-splines. Appl Comput Harmon Anal (): – Forster B, Blu T, Van De Ville D, Unser M () Shiftinvariant spaces from rotation-covariant functions. Appl Comput Harmon Anal (): – Forster B, Massopust P () Statistical encounters with complex B-splines. Constr Approx ():– Frikel J () A new framework for sparse regularization in limited angle x-ray tomography. In IEEE international symposium on biomedical imaging, Rotterdam Giles RC, Kotiuga PR, Mansuripur M () Parallel micromagnetic simulations using Fourier methods on a regular hexagonal lattice. IEEE Trans Magn ():– Grigoryan AM () Efficient algorithms for computing the -D hexagonal Fourier transforms. IEEE Trans Signal Process (): – Hales TC () The honeycomb conjecture. Discr Comput Geom :– Heil C, Walnut DF () Fundamental papers in wavelet theory. Princeton University Press, Princeton. New edition Jones DS () Generalised functions. McGrawHill, London Lai M-J, Schumaker LL () Spline functions on triangulations. Cambridge University Press, Cambridge Laine AF, Schuler S, Fan J, Huda W () Mammographic feature enhancement by multiscale analysis. IEEE Trans Med Imaging ():– Legrand P () Local regularity and multifractal methods for image and signal analysis. In: Abry P, Gonçalves P, Véhel L (eds) Scaling, fractals and wavelets, chap . Wiley-ISTE, London
Splines and Multiresolution Analysis
. Lemarié P-G () Ondelettes a localisation exponentielle. J Math pures et Appl :– . Lesage F, Provost J () The application of compressed sensing for photo-acoustic tomography. IEEE Trans Med Imaging ():– . Lipow PR, Schoenberg IJ () Cardinal interpolation and spline functions. III: Cardinal hermite interpolation. Linear Algebra Appl : – . Louis AK, Maaß P, Rieder A () Wavelets: theory and applications. Wiley, New York . Lustig M, Donoho D, Pauly JM () Sparse MRI: the application of compressed sensing for rapid MR imaging. Magn Reson Med (): – . Mallat S () Multiresolution approximations and wavelet orthonormal bases of L(R). Trans Am Math Soc :– . Mallat SG () A wavelet tour of signal processing. Academic, San Diego . Mersereau RM () The processing of hexagonally sampled two-dimensioanl signals. Proc IEEE ():– . Meyer Y () Wavelets and operators. Cambridge University Press, Cambridge . Middleton L, Sivaswamy J () Hexagonal image processing: a practical approach. Advances in pattern recognition. Springer, Berlin . Nicolier F, Laligant O, Truchetet F () B-spline quincunx wavelet transforms and implementation in Fourier domain. Proc SPIE :– . Nicolier F, Laligant O, Truchetet F () Discrete wavelet transform implementation in Fourier domain for multidimensional signal. J Electron Imaging :– . Nürnberger G () Approximation by Spline functions. Springer, Berlin . Plonka G, Tasche M () On the computation of periodic splione wavelets. Appl Comput Harmon Anal :– . Püschel M, Rötteler M () Algebraic signal processing theory: D spatail hexagonal lattice. IEEE Trans Image Proc ():– . Rabut C (a) Elementary m-harmonic cardinal B-splines. Numer Algorithms :– . Rabut C (b) High level m-harmonic cardinal B-splines. Numer Algorithms :– . Rice University () Compressive sensing resources. Online Ressource. http://dsp. rice.edu/cs (Version of Apr )
. Rudin W () Functional analysis. International series in pure and applied mathematics. McGraw-Hill, New York . Sablonnière P, Sbibih D () B-splines with hexagonal support on a uniform three-direction mesh of the plane. C R Acad Sci Paris Série I :– . Schempp W () Complex contour integral representation of cardinal spline functions, vol of Contemporary mathematics. American Mathematical Society, Providence . Schoenberg IJ () Contributions to the problem of approximation of equidistant data by analytic functions. Part a.–on the problem of osculatory interpolation. A second class of analytic approximation formulae. Quart Appl Math :– . Schoenberg IJ () Contributions to the problem of approximation of equidistant data by analytic functions. Part a.–on the problem of smoothing or graduation. A first class of analytic approximation formulae. Quart Appl Math :– . Schoenberg IJ () Cardinal interpolation and spline functions. J Approx Theory :– . Schoenberg IJ () Cardinal interpolation and spline functions. II: Interpolation of data of power growth. J Approx Theory :– . Schwartz L () Théorie des distributions. Hermann, Paris . Unser M () Splines: A perfect fit for medical imaging. In: Sonka M., Fitzpatrick JM (eds) Progress in biomedical optics and imaging, vol , no. , vol , Part I of Proceedings of the SPIE international symposium on medical imaging: image processing (MI’), San Diego, – Feb, pp – . Unser M, Blu T (Mar ) Fractional splines and wavelets. SIAM Rev ():– . Unser M, Aldroubi A, Eden M (Mar ) Fast B-spline transforms for continuous image representation and interpolation. IEEE Trans Pattern Anal Mach Intell ():– . Unser M, Aldroubi A, Eden M (Mar ) On the asymptotic convergence of B-spline wavelets to gabor functions. IEEE Trans Info Theory : – . Unser M, Aldroubi A, Eden M (Feb a) Bspline signal processing: Part I—Theory. IEEE Trans Signal Process ():–
Splines and Multiresolution Analysis
. Unser M, Aldroubi A, Eden M (Feb b) Bspline signal processing: Part II—Efficient design and applications. IEEE Trans Signal Process ():– . Van De Ville D, Blu T, Unser M, Philips W, Lemahieu I, Van de Walle R () Hex-splines: a novel spline family for hexagonal lattices. IEEE Trans Image Process ():– . Van De Ville D, Blu T, Unser M (Nov ) Isotropic polyharmonic B-splines: scaling functions and wavelets. IEEE Trans Image Process ():– . Watson AB, Ahumuda AJ, Jr () Hexagonal orthogonal-oriented pyramid as a model of
image representation in visual cortex. IEEE Trans Biomed Eng ():– . Wendt H, Roux SG, Jaffard S, Abry P () Wavelet leaders and bootstrap for multifractal analysis of images. Signal Process :– . Wojtaszczyk P () A mathematical introduction to wavelets, vol of London mathematical society student texts. Cambridge University Press, Cambridge . Young R () An introduction to nonharmonic Fourier series. Academic, New York (revised first edition )
