DDigital Shearlet Transforms
Gitta Kutyniok, Wang-Q Lim, and Xiaosheng Zhuang
Abstract
Over the past years, various representation systems which sparsely ap-proximate functions governed by anisotropic features such as edges in images havebeen proposed. We exemplarily mention the systems of contourlets, curvelets, andshearlets. Alongside the theoretical development of these systems, algorithmic re-alizations of the associated transforms were provided. However, one of the mostcommon shortcomings of these frameworks is the lack of providing a unified treat-ment of the continuum and digital world, i.e., allowing a digital theory to be a naturaldigitization of the continuum theory. In fact, shearlet systems are the only systemsso far which satisfy this property, yet still deliver optimally sparse approximationsof cartoon-like images. In this chapter, we provide an introduction to digital shearlettheory with a particular focus on a unified treatment of the continuum and digitalrealm. In our survey we will present the implementations of two shearlet transforms,one based on band-limited shearlets and the other based on compactly supportedshearlets. We will moreover discuss various quantitative measures, which allow anobjective comparison with other directional transforms and an objective tuning ofparameters. The codes for both presented transforms as well as the framework forquantifying performance are provided in the Matlab toolbox
ShearLab . One key property of wavelets, which enabled their success as a universal method-ology for signal processing, is the unified treatment of the continuum and digitalworld. In fact, the wavelet transform can be implemented by a natural digitization ofthe continuum theory, thus providing a theoretical foundation for the digital trans-form. Lately, it was observed that wavelets are however suboptimal when sparseapproximations of 2D functions are seeked. The reason is that these functions aretypically governed by anisotropic features such as edges in images or evolving shockfronts in solutions of transport equations. However, Besov models – which wavelets a r X i v : . [ m a t h . NA ] J un Gitta Kutyniok, Wang-Q Lim, and Xiaosheng Zhuang optimally encode – are clearly deficient to capture these features. Within the modelof cartoon-like images, introduced by Donoho in [13] in 1999, the suboptimal be-havior of wavelets for such 2D functions was made mathematically precise; see alsoChapter [4].Among the various directional representation systems which have since thenbeen proposed such as contourlets [12], curvelets [9], and shearlets, the shearletsystem is in fact the only one which delivers optimally sparse approximations ofcartoon-like images and still also allows for a unified treatment of the continuumand digital world. One main reason in comparison to the other two mentioned sys-tems is the fact that shearlets are affine systems, thereby enabling an extensive theo-retical framework, but parameterize directions by slope (in contrast to angles) whichgreatly supports treating the digital setting. As a thought experiment just note thata shear matrix leaves the digital grid Z invariant, which is in general not true forrotation.This raises the following questions, which we will answer in this chapter:(P1) What are the main desiderata for a digital shearlet theory?(P2) Which approaches do exist to derive a natural digitization of the continuumshearlet theory?(P3) How can we measure the accuracy to which the desiderata from (P1) arematched?(P4) Can we even introduce a framework within which different directional trans-forms can be objectively compared?Before delving into a detailed discussion, let us first contemplate about thesequestions on a more intuitive level. Several desiderata come to one’s mind, which guarantee a unified framework forboth the continuum and digital world, and provide an answer to (P1). The followingare the choices of desiderata which were considered in [20, 14]: • Parseval Frame Property.
The transform shall ideally have the tight frame prop-erty, which enables taking the adjoint as inverse transform. This property can bebroken into the following two parts, which most, but not all, transforms admit: (cid:5)
Algebraic Exactness.
The transform should be based on a theory for digitaldata in the sense that the analyzing functions should be an exact digitizationof the continuum domain analyzing elements. (cid:5)
Isometry of Pseudo-Polar Fourier Transform.
If the image is first mappedinto a different domain – here the pseudo-polar domain –, then this mapshould be an isometry. • Space-Frequency-Localization.
The analyzing elements of the associated trans-form should ideally be highly localized in space and frequency – to the extent towhich uncertainty principles allow this. igital Shearlet Transforms 3 • True Shear Invariance.
Shearing naturally occurs in digital imaging, and it can– in contrast to rotation – be precisely realized in the digital domain. Thus thetransform should be shear invariant, i.e., a shearing of the input image should bemirrored in a simple shift of the transform coefficients. • Speed.
The transform should admit an algorithm of order O ( N log N ) flops,where N is the number of digital points of the input image. • Geometric Exactness.
The transform should preserve geometric properties par-allel to those of the continuum theory, for example, edges should be mapped toedges in transform domain. • Robustness.
The transform should be resilient against impacts such as (hard)thresholding.
In general, two different types of shearlet systems are utilized today: Band-limitedshearlet systems and compactly supported shearlet systems (see also Chapters [1]and [4]). Regarding those from an algorithmic viewpoint, both have their particularadvantages and disadvantages:Algorithmic realizations of the band-limited shearlet transform have on the onehand typically a higher computational complexity due to the fact that the window-ing takes place in frequency domain. However, on the other hand, they do allow ahigh localization in frequency domain which is important, for instance, for handlingseismic data. Even more, band-limited shearlets do admit a precise digitization ofthe continuum theory.In contrast to this, algorithmic realizations of the compactly supported shearlettransform are much faster and have the advantage of achieving a high accuracy inspatial domain. But for a precise digitization one has to lower one’s sights slightly.A more comprehensive answer to (P2) will be provided in the sequel of this chapter,where we will present the digital transform based on band-limited shearlets intro-duced in [20] and the digital transform based on compactly supported shearlets from[22].
Since the introduction of directional representation systems by many pioneer re-searchers ([8, 9, 10, 11, 12]), various numerical implementations of their directionalrepresentation systems have been proposed. Let us next briefly survey the main fea-tures of the two closest to shearlets, which are the contourlet and curvelet algo-rithms.
Gitta Kutyniok, Wang-Q Lim, and Xiaosheng Zhuang • Curvelets [7]. The discrete curvelet transform is implemented in the softwarepackage
CurveLab , which comprises two different approaches. One is based onunequispaced FFTs, which are used to interpolate the function in the frequencydomain on different tiles with respect to different orientations of curvelets. Theother is based on frequency wrapping, which wraps each subband indexed byscale and angle into a fixed rectangle around the origin. Both approaches canbe realized efficiently in O ( N log N ) flops with N being the image size. Thedisadvantage of this approach is the lack of an associated continuum domaintheory. • Contourlets [12]. The implementation of contourlets is based on a directional fil-ter bank, which produces a directional frequency partitioning similar to the onegenerated by curvelets. The main advantage of this approach is that it allows atree-structured filter bank implementation, in which aliasing due to subsamplingis allowed to exist. Consequently, one can achieve great efficiency in terms ofredundancy and good spatial localization. A drawback of this approach is thatvarious artifacts are introduced and that an associated continuum domain theoryis missing.Summarizing, all the above implementations of directional representation systemshave their own advantages and disadvantages; one of the most common shortcom-ings is the lack of providing a unified treatment of the continuum and digital world.Besides the shearlet implementations we will present in this chapter, we wouldlike to refer to Chapter [2] for a discussion of the algorithm in [16] based on theLaplacian pyramid scheme and directional filtering. It should be though noted thatthis implementation is not focussed on a natural digitization of the continuum theoryand that the code was not made publicly available, both of which are crucial aspectsof the work presented in the sequel. We further would like to draw the reader’sattention to Chapter [3] which is based on [21] aiming at introducing a shearletMRA from a subdivision perspective. Finally, we should mention that a differentapproach to a shearlet MRA was recently undertaken in [17].
A major problem with many computation-based results in applied mathematics isthe non-availability of an accompanying code, and the lack of a fair and objectivecomparison with other approaches. The first problem can be overcome by followingthe philosophy of ‘reproducible research’ [15] and making the code publicly avail-able with sufficient documentation. In this spirit, the shearlet transforms presentedin this chapter are all downloadable from . Oneapproach to overcome the second obstacle is the provision of a carefully selectedset of prescribed performance measures aiming to prohibit a biased comparison onisolated tasks such as denoising and compression of specific standard images like‘Lena’, ‘Barbara’, etc. It seems far better from an intellectual viewpoint to care-fully decompose performance according to a more insightful array of tests, each igital Shearlet Transforms 5 one motivated by a particular well-understood property we are trying to obtain. Inthis chapter we will present such a framework for quantifying performance specifi-cally of implementations of directional transforms, which was originally introducedin [20, 14]. We would like to emphasize that such a framework does not only pro-vide the possibility of a fair and thorough comparison, but also enables the tuningof the parameters of an algorithm in a rational way, thereby providing an answer toboth (P3) and (P4).
Following the philosophy of the previously detailed thoughts,
ShearLab was in-troduced by Donoho, Shahram, and the authors. This software package contains • An algorithm based on band-limited shearlets introduced in [20]. • An algorithm based on compactly supported shearlets introduced in [22]. • A comprehensive framework for quantifying performance of directional repre-sentations in general.This chapter is also devoted to provide an introduction to and discuss the mathemat-ical foundation of these components.
In Section 2, we introduce and analyze the fast digital shearlet transform FDST,which is based on band-limited shearlets. Section 3 is then devoted to the presenta-tion and discussion of the digital separable shearlet transform DSST and the digitalnon-separable shearlet transform DNST. The framework of performance measuresfor parabolic scaling based transforms is provided in Section 4. In the same sec-tion, we further discuss these measures for the special cases of the three previouslyintroduced transforms.
The first algorithmic realization of a digital shearlet transform we will present,coined
Fast Digital Shearlet Transform (FDST) , is based on band-limited shear-lets. Let us start by defining the class of shearlet systems we are interested in. Re-ferring to Chapter [1], we will consider the cone-adapted discrete shearlet system SH ( φ , ψ , ˜ ψ ; ∆ , Λ , ˜ Λ ) = Φ ( φ ; ∆ ) ∪ Ψ ( ψ ; Λ ) ∪ ˜ Ψ ( ˜ ψ ; ˜ Λ ) with ∆ = Z and ShearLab (Version 1.1) is available from . Gitta Kutyniok, Wang-Q Lim, and Xiaosheng Zhuang Λ = ˜ Λ = { ( j , k , m ) : j ≥ , | k | ≤ j , m ∈ Z } . We wish to emphasize that this choice relates to a scaling by 4 j yielding an integervalued parabolic scaling matrix, which is better adapted to the digital setting than ascaling by 2 j . We further let ψ be a classical shearlet ( ˜ ψ likewise with ˜ ψ ( ξ , ξ ) = ψ ( ξ , ξ ) ), i.e., ˆ ψ ( ξ ) = ˆ ψ ( ξ , ξ ) = ˆ ψ ( ξ ) ˆ ψ ( ξ ξ ) , (1)where ψ ∈ L ( R ) is a wavelet with ˆ ψ ∈ C ∞ ( R ) and supp ˆ ψ ⊆ [ − , − ] ∪ [ , ] ,and ψ ∈ L ( R ) a ‘bump’ function satisfying ˆ ψ ∈ C ∞ ( R ) and supp ˆ ψ ⊆ [ − , ] . Weremark that the chosen support deviates slightly from the choice in the introduction,which is however just a minor adaption again to prepare for the digitization. Further,recall the definition of the cones C – C from Chapter [1].The digitization of the associated discrete shearlet transform will be performedin the frequency domain. Focussing, on the cone C , say, the discrete shearlet trans-form is of the form f (cid:55)→ (cid:104) f , ψ η (cid:105) = (cid:104) ˆ f , ˆ ψ η (cid:105) = (cid:68) ˆ f , − j ˆ ψ ( S Tk A − j · ) e π i (cid:104) A − j S k m , · (cid:105) (cid:69) , (2)where η = ( j , k , m , ι ) indexes scale j , orientation k , position m , and cone ι . Con-sidering this shearlet transform for continuum domain data (taking all cones intoaccount) implicitly induces a trapezoidal tiling of frequency space which is evi-dently not cartesian. A digital grid perfectly adapted to this situation is the so-called‘pseudo-polar grid’, which we will introduce and discuss subsequently in detail. Letus for now mention that this viewpoint enables representation of the discrete shearlettransform as a cascade of three steps:1) Classical Fourier transformation and change of variables to pseudo-polar coor-dinates.2) Weighting by a radial ‘density compensation’ factor.3) Decomposition into rectangular tiles and inverse Fourier transform of each tiles.Before discussing these steps in detail, let us give an overview of how these stepswill be faithfully digitized. First, it will be shown in Subsection 2.1, that the twooperations in Step 1) can be combined to the so-called pseudo-polar Fourier trans-form. An oversampling in radial direction of the pseudo-polar grid, on which thepseudo-polar Fourier transform is computed, will then enable the design of ‘density-compensation-style’ weights on those grid points leading to Steps 1) & 2) being anisometry. This will be discussed in Subsection 2.2. Subsection 2.3 is then concernedwith the digitization of the discrete shearlets to subband windows. Notice that adigital analog of (2) moreover requires an additional 2D-iFFT. Thus, concludingthe digitization of the discrete shearlet transform will cascade the following steps,which is the exact analogy of the continuum domain shearlet transform (2):(S1) PPFT: Pseudo-polar Fourier transform with oversampling factor in the radialdirection. igital Shearlet Transforms 7 (S2) Weighting: Multiplication by ‘density-compensation-style’ weights.(S3) Windowing: Decomposing the pseudo-polar grid into rectangular subbandwindows with additional 2D-iFFT.With a careful choice of the weights and subband windows, this transform is anisometry. Then the inverse transform can be computed by merely taking the adjointin each step. A final discussion on the FDST will be presented in Subsection 2.4. We start by discussing Step (S1).
In [5], a fast pseudo-polar Fourier transform (PPFT) which evaluates the discreteFourier transform at points on a trapezoidal grid in frequency space, the so-calledpseudo-polar grid, was already developed. However, the direct use of the PPFT isproblematic, since it is – as defined in [5] – not an isometry. The main obstacleis the highly nonuniform arrangement of the points on the pseudo-polar grid. Thisintuitively suggests to downweight points in regions of very high density by usingweights which correspond roughly to the density compensation weights underly-ing the continuous change of variables. This will be enabled by a sufficient radialoversampling of the pseudo-polar grid.This new pseudo-polar grid, which we will denote in the sequel by Ω R to indicatethe oversampling rate R , is defined by Ω R = Ω R ∪ Ω R , (3)where Ω R = { ( − nR · (cid:96) N , nR ) : − N ≤ (cid:96) ≤ N , − RN ≤ n ≤ RN } , (4) Ω R = { ( nR , − nR · (cid:96) N ) : − N ≤ (cid:96) ≤ N , − RN ≤ n ≤ RN } . (5)This grid is illustrated in Fig. 1. We remark that the pseudo-polar grid introducedin [5] coincides with Ω R for the particular choice R =
2. It should be emphasizedthat Ω R = Ω R ∪ Ω R is not a disjoint partitioning, nor is the mapping ( n , (cid:96) ) (cid:55)→ ( − nR · (cid:96) N , nR ) or ( nR , − nR · (cid:96) N ) injective. In fact, the center C = { ( , ) } (6)appears N + Ω R as well as Ω R , and the points on the seam lines Gitta Kutyniok, Wang-Q Lim, and Xiaosheng Zhuang Ω R Ω R Ω R Fig. 1
The pseudo-polar grid Ω R = Ω R ∪ Ω R for N = R = S R = { ( − nR , nR ) : − RN ≤ n ≤ RN , n (cid:54) = } , S R = { ( nR , − nR ) : − RN ≤ n ≤ RN , n (cid:54) = } . appear in both Ω R and Ω R . Definition 1.
Let N , R be positive integer, and let Ω R be the pseudo-polar grid givenby (3). For an N × N image I : = { I ( u , v ) : − N ≤ u , v ≤ N − } , the pseudo-polarFourier transform (PFFT) ˆ I of I evaluated on Ω R is then defined to beˆ I ( ω , ω ) = N / − ∑ u , v = − N / I ( u , v ) e − π im ( u ω + v ω ) , ( ω , ω ) ∈ Ω R , where m ≥ N is an integer.We wish to mention that m ≥ N is typically set to be m = R ( RN + ) for com-putational reasons (see also [5]), but we for now allow more generality. It was shown in [5], that the PPFT can be realized in O ( N log N ) flops with N × N being the size of the input image. We will now discuss how the extended pseudo-polar Fourier transform as defined in Definition 1 can be computed with similarcomplexity.For this, let I be an image of size N × N . Also, m is set – but not restricted –to be m = R ( RN + ) ; we will elaborate on this choice at the end of this subsec-tion. We now focus on Ω R , and mention that the PPFT on the other cone can becomputed similarly. Rewriting the pseudo-polar Fourier transform from Definition1, for ( ω , ω ) = ( − nR · (cid:96) N , nR ) ∈ Ω R , we obtain igital Shearlet Transforms 9 ˆ I ( ω , ω ) = N / − ∑ u , v = − N / I ( u , v ) e − π im ( u ω + v ω ) = N / − ∑ u = − N / N / − ∑ v = − N / I ( u , v ) e − π im ( u − n (cid:96) RN + v nR ) = N / − ∑ u = − N / (cid:32) N / − ∑ v = − N / I ( u , v ) e − π ivnRN + (cid:33) e − π iu (cid:96) · − n ( RN + ) · N . (7)This rewritten form, i.e., (7), suggests that the pseudo-polar Fourier transform ˆ I of I on Ω R can be obtained by performing the 1D FFT on the extension of I alongdirection v and then applying a fractional Fourier transform (frFT) along direction u . To be more specific, we require the following operations: Fractional Fourier Transform.
For c ∈ C N + , the (unaliased) discrete fractionalFourier transform by α ∈ C is defined to be ( F α N + c )( k ) : = N / ∑ j = − N / c ( j ) e − π i · j · k · α , k = − N , . . . , N . It was shown in [6], that the fractional Fourier transform F α N + c can be computedusing O ( N log N ) operations. For the special case of α = / ( N + ) , the fractionalFourier transform becomes the (unaliased) 1D discrete Fourier Transform (1D FFT),which in the sequel will be denoted by F . Similarly, the 2D discrete Fourier Trans-form (2D FFT) will be denoted by F , and the inverse of the F by F − (2D iFFT). Padding Operator . For N even, m > N an odd integer, and c ∈ C N , the paddingoperator E m , n gives a symmetrically zero padding version of c in the sense that ( E m , N c )( k ) = (cid:40) c ( k ) k = − N , . . . , N − , k ∈ {− m , . . . , m } \ {− N , . . . , N − } . Using these operators, (7) can be computed byˆ I ( ω , ω ) = N / − ∑ u = − N / F ◦ E RN + , N ◦ I ( u , n ) e − π iu (cid:96) · − n ( RN + ) · N / = N / ∑ u = − N / E N + , N ◦ F ◦ E RN + , N ◦ I ( u , n ) e − π iu (cid:96) · − n ( RN + ) · N = ( F α n N + ˜ I ( · , n ))( (cid:96) ) , where ˜ I = E N + , N ◦ F ◦ E RN + , N ◦ I ∈ C ( RN + ) × ( N + ) and α n = − n ( RN + ) N / . Sincethe 1D FFT and 1D frFT require only O ( N log N ) operations for a vector of size N , the total complexity of this algorithm for computing the pseudo-polar Fouriertransform from Definition 1 is indeed O ( N log N ) for an image of size N × N . We would like to also remark that for a different choice of constant m , one cancompute the pseudo-polar Fourier transform also with complexity O ( N log N ) foran image of size N × N . This however requires application of the fractional Fouriertransform in both directions u and v of the image, which results in a larger constantfor the computational cost; see also [6]. Next we tackle Step (S2), which is more delicate than it might seem, since theweights will not be derivable from simple density compensation arguments.
For this, we now aim to choose weights w : Ω R → R + so that the extended PPFTfrom Definition 1 becomes an isometry, i.e., N / − ∑ u , v = − N / | I ( u , v ) | = ∑ ( ω , ω ) ∈ Ω R w ( ω , ω ) · | ˆ I ( ω , ω ) | . (8)Observing the symmetry of the pseudo-polar grid, it seems natural to select weightfunctions w which have full axis symmetry properties, i.e., for all ( ω , ω ) ∈ Ω R ,we require w ( ω , ω ) = w ( ω , ω ) , w ( ω , ω ) = w ( − ω , ω ) , w ( ω , ω ) = w ( ω , − ω ) . (9)Then the following ‘Plancherel theorem’ for the pseudo-polar Fourier transform on Ω R – similar to the one for the Fourier transform on the cartesian grid – can beproved. Theorem 1 ([20]).
Let N be even, and let w : Ω R → R + be a weight function satis-fying (9) . Then (8) holds, if and only if, the weight function w satisfies δ ( u , v ) = w ( , )+ · ∑ (cid:96) = , N / RN / ∑ n = w ( nR , nR · − (cid:96) N ) · cos ( π u · nm R ) · cos ( π v · nm R · (cid:96) N )+ · N / − ∑ (cid:96) = RN / ∑ n = w ( nR , nR · − (cid:96) N ) · cos ( π u · nm R ) · cos ( π v · nm R · (cid:96) N ) (10) for all − N + ≤ u , v ≤ N − .Proof. We start by computing the right hand side of (8): igital Shearlet Transforms 11 ∑ ( ω , ω ) ∈ Ω R w ( ω , ω ) · | ˆ I ( ω , ω ) | = ∑ ( ω , ω ) ∈ Ω R w ( ω , ω ) · (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) N / − ∑ u , v = − N / I ( u , v ) e − π im ( u ω + v ω ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) = ∑ ( ω , ω ) ∈ Ω R w ( ω , ω ) · (cid:34) N / − ∑ u , v = − N / N / − ∑ u (cid:48) , v (cid:48) = − N / I ( u , v ) I ( u (cid:48) , v (cid:48) ) e − π im (( u − u (cid:48) ) ω +( v − v (cid:48) ) ω ) (cid:35) = ∑ ( ω , ω ) ∈ Ω R w ( ω , ω ) · N / − ∑ u , v = − N / | I ( u , v ) | + N / − ∑ u , v , u (cid:48) , v (cid:48) = − N / ( u , v ) (cid:54) =( u (cid:48) , v (cid:48) ) I ( u , v ) I ( u (cid:48) , v (cid:48) ) · (cid:34) ∑ ( ω , ω ) ∈ Ω R w ( ω , ω ) · e − π im (( u − u (cid:48) ) ω +( v − v (cid:48) ) ω ) (cid:35) . Choosing I = c u , v δ ( u − u , v − v ) + c u , v δ ( u − u , v − v ) for all − N / ≤ u , v , u , v ≤ N / − c u , v , c u , v ∈ C , we can conclude that (8) holds if andonly if ∑ ( ω , ω ) ∈ Ω R w ( ω , ω ) · e − π im ( u ω + v ω ) = δ ( u , v ) , − N + ≤ u , v ≤ N − . By the symmetry of the weights (9), this is equivalent to ∑ ( ω , ω ) ∈ Ω R w ( ω , ω ) · [ cos ( π m u ω ) cos ( π m v ω )] = δ ( u , v ) (11)for all − N + ≤ u , v ≤ N −
1. From this, we can deduce that (11) is equivalent to(10), which proves the theorem. (cid:117)(cid:116)
Notice that (10) is a linear system with RN / + RN / + ( N − ) equations, wherefore, in general, one needs the oversampling factor R to be atleast 16 to enforce solvability. The computation of the weights satisfying Theorem 1 by solving the full linearsystem of equations (10) is much too complex. Hence, we relax the requirementfor exact isometric weighting, and represent the weights in terms of undercompletebasis functions on the pseudo-polar grid.More precisely, we first choose a set of basis functions w , . . . , w n : Ω R → R + such that n ∑ j = w j ( ω , ω ) (cid:54) = ( ω , ω ) ∈ Ω R . We then represent weight functions w : Ω R → R + by w : = n ∑ j = c j w j , (12)with c , . . . , c n being nonnegative constants. This approach now enables solving(10) for the constants c , . . . , c n using the least squares method, thereby reducing thecomputational complexity significantly. The ‘full’ weight function w is then givenby (12).We next present two different choices of weights which were derived by thisrelaxed approach. Notice that ( ω , ω ) and ( n , (cid:96) ) will be used interchangeably. Choice 1 . The set of basis functions w , . . . , w is defined as follows: Center : w = ( , ) , Boundary : w = { ( ω , ω ) : | n | = NR / , ω = ω } and w = { ( ω , ω ) : | n | = NR / , ω (cid:54) = ω } , Seam lines : w = | n | · { ( ω , ω ) :1 ≤| n | < NR / , ω = ω } , Interior : w = | n | · { ( ω , ω ) :1 ≤| n | < NR / , ω (cid:54) = ω } . Choice 2 . The set of basis functions w , . . . , w N / + is defined as follows: Center : w = ( , ) , Radial Lines : w (cid:96) + = { ( ω , ω ) :1 < | n | < NR / , ω = (cid:96) N / ω } , (cid:96) = , , . . . , N / . The associated weight functions are displayed in Fig. 2. In general, suitableweight functions usually obey the pattern of linearly increasing values along theradial direction. Thus, this is a natural requirement for the basis functions.
A visual comparison shows that the patterns of the weight functions associated withChoices 1 and 2 are seemingly similar (see Fig. 2). However, carefully chosen mea-sures reveal that their performances can in fact be quite different.One essential criterion for the quality of a weight function is the degree to whichit allows a Plancherel theorem for the pseudo-polar Fourier transform as studied inTheorem 1. This can be measured in the following way – the reader might want tocompare this performance measure with the measures introduced in Subsection 4.2: igital Shearlet Transforms 13Choice 1 Choice 2
Fig. 2
Weight functions on the pseudo-polar grid for N =
256 and R = Let P and P (cid:63) denote the operators for the pseudo-polar Fourier transform and itsadjoint, respectively, and let w – by slightly abusing notation – denote the weightingoperator on the pseudo-polar grid Ω R . Letting R =
8, a sequence of 5 random images I , . . . , I of size N × N with standard normally distributed entries is generated tocompute M : = ∑ i = (cid:107) P (cid:63) wPI i − I i (cid:107) (cid:107) I i (cid:107) . The performance of the weight functions arising from Choices 1 and 2 with respectto this measure is presented in Table 1.
Table 1
Comparison of Choices 1 and 2 based on performance measure M . N
32 64 128 256 512Choice 1 4.2E-3 4.0E-3 1.8E-3 1.5E-3 8.8E-4Choice 2 9.8E-3 6.2E-3 3.4E-3 2.1E-3 N/A
Interestingly, a structured image, e.g., by using the measure M : = (cid:107) P (cid:63) wPI − I (cid:107) (cid:107) I (cid:107) , I the image ‘Barbara’,yields an even better performance and a better distinction, see Table 2. One couldreason that this behavior is due to the fact that the energy of most ‘real’ imagesis concentrated in the low frequency region, in which density compensation of thepseudo-polar grid is not as necessary as in the high frequency regions.These two tables show firstly, that with growing N , the weighted pseudo-polarFourier transform seems to converge to being an isometry on the testing imageclass. Secondly, judging from the relatively small deviation from being an isometry,it seems quite reasonable to choose a basis of weight functions forcing the weightsto linearly increase along the radial direction. And, thirdly, although Choice 2 con-tains many more basis functions than Choice 1, the performance results are worse, Table 2
Comparison of Choices 1 and 2 based on the performance measure M . N
32 64 128 256 512Choice 1 2.8E-3 1.2E-3 8.3E-4 3.9E-4 1.5E-4Choice 2 5.6E-3 2.8E-3 2.2E-3 9.1E-4 N/A which is very counterintuitive. The reason for this is the numerical instability whencomputing a minimizing set of coefficients for the basis of weight functions, whichcauses these effects.
For the FDST – as also in the implementation in
ShearLab – the coefficientsin the expansion (12) will be computed off-line, and then hardwired in the code.This enables the weighting of a function on the pseudo-polar grid to simply be apoint-wise multiplication in each sampling point. That is, letting J : = ˆ I : Ω R → C bethe pseudo-polar Fourier transform of an N × N image I and w : Ω R → R + be anysuitable weight function on Ω R , the values J w ( ω , ω ) = J ( ω , ω ) · (cid:112) w ( ω , ω ) for all ( ω , ω ) ∈ Ω R need to be computed.Let us comment on why the square root of the weight is utilized. If the weights w satisfy the condition in Theorem 1, we obtain P ∗ wP = Id, which can be written ina symmetric form as follows: ( √ wP ) ∗ √ wP = Id. This form shows that the operator √ wP can be inverted by taking the adjoint ( √ wP ) ∗ . In other words, each image canbe reconstructed from its weighted pseudo-polar Fourier transform by applying theadjoint of the weighted pseudo-polar Fourier transform. This issue will be discussedin further detail in Subsection 2.4.2. We next aim at deriving a faithful digitization of the shearlet transform associatedwith a band-limited cone-adapted discrete shearlet system to the pseudo-polar grid.This would settle Step (S3).
For this, let us recall the definition of the discrete shearlet transform associatedwith (2); taking the particular form (1) of the shearlet ψ ∈ L ( R ) into account. igital Shearlet Transforms 15 Restricting our attention to the cone C , we obtain f (cid:55)→ (cid:68) ˆ f , − j ˆ ψ ( S Tk A − j · ) χ C e π i (cid:104) A − j S k m , · (cid:105) (cid:69) = (cid:68) ˆ f , − j ˆ ψ ( − j ξ ) ˆ ψ ( k + j ξ ξ ) χ C e π i (cid:104) A − j S k m , · (cid:105) (cid:69) , for scale j , orientation k , position m , and cone ι .To approach a faithful digitization, we first have to partition Ω R according to thepartitioning of the plane into C , C , C , and C , as well as a centered rectangle R . The center C as defined in (6) will play the role of R . Thus it remains to partitionthe set Ω R beyond the already defined partitioning into Ω R and Ω R (cf. (4) and (5))by setting Ω R = Ω R ∪ C ∪ Ω R and Ω R = Ω R ∪ C ∪ Ω R , where Ω R = { ( − nR · (cid:96) N , nR ) : − N ≤ (cid:96) ≤ N , ≤ n ≤ RN } , Ω R = { ( − nR · (cid:96) N , nR ) : − N ≤ (cid:96) ≤ N , − RN ≤ n ≤ − } , Ω R = { ( nR , − nR · (cid:96) N ) : − N ≤ (cid:96) ≤ N , ≤ n ≤ RN } Ω R = { ( nR , − nR · (cid:96) N ) : − N ≤ (cid:96) ≤ N , − RN ≤ n ≤ − } . When restricting to the cone Ω R , say, the exact digitization of the coefficients ofthe discrete shearlet system is ∑ ω : =( ω , ω ) ∈ Ω R J ( ω , ω ) − j ˆ ψ ( S Tk A − j ω ) e − π i (cid:104) A − j S k m , ω (cid:105) = ∑ ( ω , ω ) ∈ Ω R J ( ω , ω ) − j W ( − j ω x ) V ( k + j ω ω ) e − π i (cid:104) A − j S k m , ω (cid:105) = RN ∑ n = N ∑ (cid:96) = − N J ( ω , ω ) − j W ( − j nR ) V ( k − j + (cid:96) N ) e − π i (cid:104) m , S Tk A − j ω (cid:105) , (13)where V and W as well as the ranges of j , k , and m are to be carefully chosen.Our main objective will be to achieve a digital shearlet transform, which is anisometry. This – as in the continuum domain situation – is equivalent to requiringthe associated shearlet system to form a tight frame for functions J : Ω R → C . Forthe convenience of the reader let us recall the notion of a Parseval frame in thisparticular situation. A sequence ( ϕ λ ) λ ∈ Λ – Λ being some indexing set – is a tightframe for all functions J : Ω R → C , if ∑ λ ∈ Λ (cid:12)(cid:12)(cid:12) ∑ ( ω , ω ) ∈ Ω R J ( ω , ω ) ϕ λ ( ω , ω ) (cid:12)(cid:12)(cid:12) = ∑ ( ω , ω ) ∈ Ω R | J ( ω , ω ) | . In the sequel we will define digital shearlets on Ω R and extend the definition tothe other cones by symmetry. We start by defining the scaling function, which will depend on two functions V and W , and the generating digital shearlet, which will depend on again two functions V and W . W and W will be chosen to be Fourier transforms of wavelets, and V and V will be chosen to be ‘bump’ functions, paralleling the construction of classicalshearlets.First, let W be the Fourier transform of the Meyer scaling function such thatsupp W ⊆ [ − , ] and W ( ± ) = , (14)and let V be a ‘bump’ function satisfyingsupp V ⊆ [ − / , / ] with V ( ξ ) ≡ | ξ | ≤ , ξ ∈ R . Then we define the scaling function φ for the digital shearlet system to beˆ φ ( ξ , ξ ) = W ( − j L ξ ) V ( − j L ξ ) , ( ξ , ξ ) ∈ R . For now, we define it in continuum domain, and will later restrict this function tothe pseudo-polar grid.Let next W be the Fourier transform of the Meyer wavelet function satisfying thesupport constraintssupp W ⊆ [ − , / ] ∪ [ / , ] and W ( ± / ) = W ( ± ) = , (15)as well as, choosing the lowest scale j L to be j L : = −(cid:100) log ( R / ) (cid:101) , | W ( − j L ξ ) | + (cid:100) log N (cid:101) ∑ j = j L | W ( − j ξ ) | = | ξ | ≤ N , ξ ∈ R . (16)We further choose V to be a ‘bump’ function satisfyingsupp V ⊆ [ − , ] and V ( ± ) = , (17)as well as | V ( ξ − ) | + | V ( ξ ) | + | V ( ξ + ) | = | ξ | ≤ , ξ ∈ R . (18)Then the generating shearlet ψ for the digital shearlet system on Ω R is defined asˆ ψ ( ξ , ξ ) = W ( ξ ) V ( ξ ξ ) , ( ξ , ξ ) ∈ R . (19) igital Shearlet Transforms 17 Notice that (18) implies j ∑ s = − j | V ( j ξ − s ) | = | ξ | ≤ , ξ ∈ R ; j ≥ , (20)which will become important for the analysis of frame properties. For the particularchoice of V , W , V , and W in ShearLab , we refer to Subsection 2.3.7.
We from now on assume that R and N are both positive even integers and that N = n for some integer n ∈ N . This poses no restrictions, since both parameters canbe enlarged to satisfy this condition.We start by analyzing the range of j . Recalling the definition of the shearlet ψ in (19) and the support properties of W and V in (15) and (17), respectively, weobserve that the digitized shearlet2 − j W ( − j nR ) V ( k − j + (cid:96) N ) e π i (cid:104) m , S Tk A − j ω (cid:105) (21)from (13) has radial support n = j − R + t , t = , . . . , j − · R (22)on the cone Ω R . To determine the appropriate range of j , we will analyze theprecise support in radial direction. If j < −(cid:100) log ( R / ) (cid:101) , then n <
1, which corre-sponds to only one point – the origin – and is dealt with by the scaling function. If j > (cid:100) log N (cid:101) , we have n ≥ RN . Hence the value W ( / ) = j ∈ { j L , . . . , j H } , where j L : = −(cid:100) log ( R / ) (cid:101) and j H : = (cid:100) log N (cid:101) . Next, we determine the appropriate range of k . Again recalling the definition ofthe shearlet ψ in (19), the digitized shearlet (21) has angular support (cid:96) = − j − N ( k − ) + t , t = , . . . , − j N (23)on the cone Ω R . To compute the range of k , we start by examining the case j ≥ k > j , we have (cid:96) ≥ N /
2. Hence the value V ( − ) = k ≥ − j . Thus the shearing parameter will be chosen to be k ∈ {− j , . . . , j } . We next compute the support as well as the support size of scaled and sheared ver-sion of digital shearlets. This will be used for the normalization of digital shearlets.As before, we first analyze the radial support. By (22), the radial supports of thewindows associated with scales j L < j < j H is n = j − R + t , t = , . . . , j − · R , (24)and the radial support of the windows associated with the scale j L = −(cid:100) log ( R / ) (cid:101) and j H = (cid:100) log N (cid:101) are n = t , t = , . . . , j L + R , for j = j L , n = j H − R + t , t = , . . . , RN − j H − R , for j = j H . (25)Turning to the angular direction, by (23), the angular support of the windows atscale j associated with shears − j < k < j is (cid:96) = − j − N ( k − ) + t , t = , . . . , − j N , (26)the angular support at scale j associated with the shear parameter k = − j is (cid:96) = − j − N ( − j − ) + t , t = − j N , . . . , − j N , and for k = j it is (cid:96) = − j − N ( j − ) + t , t = , . . . , − j N . (27)For the case j <
0, we simply let k = (cid:96) = − N / + t with t = , . . . , N . Also,for this lower frequency case, the window function W ( − j ω ) V ( k + j ω ω ) is slightlymodified to be W ( − j ω ) V ( k + j ω ω ) .These computations now allow us to determine the support size of the function W ( − j ω ) V ( k + j ω ω ) in terms of pairs ( n , (cid:96) ) , which for scale j and shear k , is L j = j + R : j = j L , j − · R + j L < j < j H , RN − j − R + j = j H , (28)and L j , k = − j N + − j < k < j with j ≥ , − j N + k ∈ {− j , j } with j ≥ , N + j < . (29) igital Shearlet Transforms 19 We next digitize the exponential term in (21), which can be rewritten as e − π i (cid:104) m , S Tk A − j ω (cid:105) = e − π i (cid:104) m , ( − j ω , − j k ω + − j ω ) (cid:105) = e − π i (cid:104) m , ( − j nR , − j k nR − − j (cid:96) nRN ) (cid:105) . We observe two obstacles: • The change of variables τ : = S Tk A − j ω possible in (13) can not be performedsimilarly in this situation due to the fact that the pseudo-polar grid is not invari-ant under the action of S Tk A − j . This is however the first step in the continuumdomain reasoning for tightness; see Chapter [1]. • The Fourier transform of a function defined on the pseudo-polar grid does not satisfy any Plancherel theorem.These problems require a slight adjustment of the exponential term, which will bethe only adaption we allow us to make when digitizing. This will circumvent thetwo obstacles and enable us to construct a Parseval frame as well as derive a directapplication of the inverse Fast Fourier transform in FDST.The adjustment will be made by using the mapping θ : R \ { } → R defined by θ ( x , y ) = ( x , yx ) . This yields the modified exponential term e − π i (cid:104) m , ( θ ◦ ( S Tk ) − )( − j nR , − j k nR − − j (cid:96) nRN ) (cid:105) = e − π i (cid:104) m , ( − j nR , − j + (cid:96) N ) (cid:105) , (30)which can be rewritten as e − π i (cid:104) m , ( − j nR , − j + (cid:96) N ) (cid:105) = e − π i ( m +( − k ) m ) e − π i (cid:68) m , ( − j t R , − j + t N ) (cid:69) , with t and t ranging over an appropriate set defined by (24), (25), and (26)–(27).Fig. 3 illustrates this adjustment. ( S Tk ) − θ Fig. 3
Adjustment of the exponential term through the map θ ◦ ( S Tk ) − . Now, taking into account of the support size of each W ( − j ω ) V ( k + j ω ω ) asgiven in (28) and (29), we obtain the following reformulation of (30):exp (cid:26) − π i (cid:28) m , (cid:18) L j − j ( / R ) L j t , − L j , k j + ( / N ) L j , k t (cid:19)(cid:29)(cid:27) , t , t . (31) This version shows that we might regard the exponential terms as characters of asuitable locally compact abelian group (see [18]) with associated annihilator identi-fied with the rectangle R j , k = (cid:40)(cid:32) j R · r L j , − N j + · r L j , k (cid:33) : r = , . . . , L j − , r = , . . . , L j , k − (cid:41) , where L j and L j , k were defined in (28) and (29), respectively. This viewpointwill be crucial to guarantee that the digital shearlet system defined in Subsection2.3.6 provides a Parseval frame on the pseudo-polar grid Ω R . In practice, (31) alsoensures that in Step (S3) on each windowed image on the pseudo-polar grid onlya 2D-iFFT – in contrast to a fractional Fourier transform – needs to be performed,thereby reducing the computational complexity.For the low frequency square, we further require the set R = { ( r , r ) : r = − , . . . , , r = − N , . . . , N } , which will be shown to be sufficient for guaranteeing that digital shearlet systemforms a Parseval frame. We are now ready to define digital shearlets, which we define as functions on thepseudo-polar grid Ω R . The spatial domain picture can thus be derived by the inversepseudo-polar Fourier transform. Definition 2.
Retaining the definitions and notations from Subsection 2.3, for all ( ω , ω ) ∈ Ω R , we define digital shearlets at scale j ∈ { j L , . . . , j H } , shear k =[ − j , j ] ∩ Z , and spatial position m ∈ R j , k by σ j , k , m ( ω , ω ) = C ( ω , ω ) √ | R j , k | W ( − j ω ) V j ( k + j ω ω ) χ Ω R ( ω , ω ) e π i (cid:68) m , ( − j ω , j ω ω ) (cid:69) , where V j = V for j ≥ V j = V for j <
0, and C ( ω , ω ) = ( ω , ω ) (cid:54)∈ S R ∪ S R , √ : ( ω , ω ) ∈ ( S R ∪ S R ) \ C , √ ( N + ) : ( ω , ω ) ∈ C . The shearlets σ j , k , m , σ j , k , m , σ j , k , m on the remaining cones are defined accordingly bysymmetry with equal indexing sets for scale j , shear k , and spatial location m . For ι = , ( ω , ω ) ∈ Ω ι R , and n ∈ R , we define the scaling function ϕ ι n ( ω , ω ) = C ( ω , ω ) √ | R | ˆ φ ( ω , ω ) χ Ω ι R ( ω , ω ) e π i (cid:104) n , ( n , (cid:96) N + ) (cid:105) . igital Shearlet Transforms 21 Then the digital shearlet system DSH is defined by
DSH = { ϕ ι n : ι = , , n ∈ R } ∪ { σ ι j , k , m : j ∈ { j L , . . . , j H } , k ∈ {− j , j } , m ∈ R j , s , ι = , , , } . As desired, the digital shearlet system
DSH , which we derived as a faithful digiti-zation of the continuum domain band-limited cone-adapted discrete shearlet system,forms a Parseval frame for J : Ω R → C . Theorem 2 ([20]).
The digital shearlet system DSH defined in Definition 2 forms aParseval frame for functions J : Ω R → C .Proof. Letting J : Ω R → C , we claim that (cid:104) J , J (cid:105) Ω R = ∑ ι , n |(cid:104) J , ϕ ι n (cid:105) Ω R | + ∑ ι , j , k , m |(cid:104) J , σ ι j , k , m (cid:105) Ω R | (32)which proves the result. Here (cid:104) J , J (cid:105) Ω R : = ∑ ( ω , ω ) ∈ Ω R J ( ω , ω ) J ( ω , ω ) for J , J : Ω R → C .We start by analyzing the first term on the RHS of (32). Let ι ∈ { , } and J C : Ω R → C be defined by J C ( ω , ω ) : = C ( ω , ω ) · J ( ω , ω ) for ( ω , ω ) ∈ Ω R .Using the support conditions of ˆ φ , ∑ n |(cid:104) J , ϕ ι n (cid:105) Ω R | = ∑ n (cid:12)(cid:12)(cid:12) ∑ ( ω , ω ) ∈ Ω ι R J ( ω , ω ) ϕ ι n ( ω , ω ) (cid:12)(cid:12)(cid:12) = | R | ∑ n (cid:12)(cid:12)(cid:12) ∑ ( ω , ω ) ∈ Ω ι R J C ( ω , ω ) · ˆ φ ( ω , ω ) · e − π i (cid:104) n , ( n , (cid:96) N + ) (cid:105) (cid:12)(cid:12)(cid:12) = | R | ∑ n (cid:12)(cid:12)(cid:12) ∑ n = − N / ∑ (cid:96) = − N / J C ( ω , ω ) · ˆ φ ( ω , ω ) · e − π i (cid:104) n , ( n , (cid:96) N + ) (cid:105) (cid:12)(cid:12)(cid:12) . The choice of R now allows us to use the Plancherel formula, see Subsection 2.3.5.Exploiting again support properties (see Subsection 2.3.5), we conclude that ∑ n |(cid:104) J , ϕ ι n (cid:105) Ω R | = ∑ ( ω , ω ) ∈ Ω ι R | C ( ω , ω ) · J ( ω , ω ) | · | ˆ φ ( ω , ω ) | . Combining ι = , ∑ ι ∑ n |(cid:104) J , ϕ ι n (cid:105) Ω R | = ∑ ( ω , ω ) ∈ Ω R | J ( ω , ω ) | · | W ( ω ) | . (33)Next we study the second term on the RHS in (32). By symmetry, it suffices toconsider the case ι =
21. By the support conditions on W and V (see (15) and (17)), ∑ j , k , m |(cid:104) J , σ j , k , m (cid:105) Ω R | = ∑ j , k ∑ m ∈ R j , k (cid:12)(cid:12)(cid:12) ∑ ( ω , ω ) ∈ Ω R J ( ω , ω ) σ j , k , m ( ω , ω ) (cid:12)(cid:12)(cid:12) = ∑ j , k | R j , k | ∑ m ∈ R j , k (cid:12)(cid:12)(cid:12) ∑ ( ω , ω ) ∈ Ω R J C ( ω , ω ) · W ( − j ω ) · V j ( k + j ω ω ) · e − π i (cid:68) m , ( − j ω , j ω ω ) (cid:69) (cid:12)(cid:12)(cid:12) = ∑ j , k | R j , k | ∑ m ∈ R j , k (cid:12)(cid:12)(cid:12) j + ( R / ) ∑ n = j − ( R / ) − j − N ( k + ) ∑ (cid:96) = − j − N ( k − ) J C ( ω , ω ) · W ( − j ω ) · V j ( k + j ω ω ) · e − π i (cid:104) m , ( − j nR , − j + (cid:96) N ) (cid:105) (cid:12)(cid:12)(cid:12) . Similarly as before, the choice of R j , k does allow us to use the Plancherel formula,see Subsection 2.3.5. Hence, ∑ j , k , m |(cid:104) J , σ j , k , m (cid:105) Ω R | = ∑ j , k ∑ ( ω , ω ) ∈ Ω R (cid:12)(cid:12)(cid:12) J C ( ω , ω ) · W ( − j ω ) V j ( k + j ω ω ) (cid:12)(cid:12)(cid:12) . Next we use (20) to obtain ∑ j , k ∑ ( ω , ω ) ∈ Ω R (cid:12)(cid:12)(cid:12) J C ( ω , ω ) · W ( − j ω ) · V j ( k + j ω ω ) (cid:12)(cid:12)(cid:12) = ∑ ( ω , ω ) ∈ Ω R | J C ( ω , ω ) | j H ∑ j = j L | W ( − j ω ) | · j ∑ k = − j | V j ( k + j ω ω ) | = ∑ ( ω , ω ) ∈ Ω R | J C ( ω , ω ) | j H ∑ j = j L | W ( − j ω ) | . Thus the second term on the RHS in (32) equals ∑ ι ∑ j , k , m |(cid:104) J , σ ι j , k , m (cid:105) Ω R | = ∑ ( ω , ω ) ∈ Ω R | J ( ω , ω ) | · j H ∑ j = j L . | W ( − j ω ) | . (34)Finally, our claim (32) follows from combining (33), (34), and (16). (cid:117)(cid:116) The final Step (S3) of the FDST then consists in decomposing the data on the pointsof the pseudo-polar grid given by the previously – in Steps (S1) and (S2) – computedweighted pseudo-polar image J w : Ω R → C into rectangular subband windows ac-cording to the digital shearlet system DSH defined in Definition 2, followed by a2D-iFFT. More precisely, given J w , the set of digital shearlet coefficients c ι n : = (cid:10) J w , ϕ ι n (cid:11) Ω R for all ι , n igital Shearlet Transforms 23 and c ι j , k , m : = (cid:68) J w , σ ι j , k , m (cid:69) Ω R for all j , k , m , ι is computed followed by application of the 2D-iFFT to each windowed image J w ϕ ι and J w σ ι j , k , restricted on the support of ϕ ι and σ ι j , k , , respectively.The definition of the digital shearlet system DSH in Definition 2 requires appro-priate choices of the functions φ , V , V , W , and W , and the required conditions arestated throughout Subsection 2.3.2. We now discuss one particular choice, which ischosen in ShearLab . We start selecting the ‘wavelets’ W and W . In Subsection2.3.2, these functions were defined to be Fourier transforms of the Meyer scalingfunction and wavelet function, respectively, i.e., W ( ξ ) = | ξ | ≤ , cos (cid:2) π ν ( | ξ | − ) (cid:3) : ≤ | ξ | ≤ , , and W ( ξ ) = sin (cid:2) π ν ( | ξ | − ) (cid:3) : ≤ | ξ | ≤ , cos (cid:2) π ν ( | ξ | − ) (cid:3) : 1 ≤ | ξ | ≤ , , where ν ≥ C k function or C ∞ function such that ν ( x ) + ν ( − x ) = ≤ x ≤
1. One possible choice for ν is the function ν ( x ) = x ( − x + x − x ) ,0 ≤ x ≤
1, which then automatically fixes W and W . Since | W ( ξ ) | + | W ( ξ ) | = | ξ | ≤
1, the required condition (16) is satisfied. The graphs of this choice offunctions W , W , and ν are illustrated in Fig. 4. ν W W Fig. 4
The graphs of ν , W , and w . The function ν can be also used to design the ‘bump’ function V as well,which needs to satisfy (18). One possible choice for V is to define it by V ( ξ ) = (cid:112) ν ( + ξ ) + ν ( − ξ ) , − ≤ ξ ≤ V can then simply be chosen as V ≡ φ is defined depending on V and W , wherefore fixingthese two functions determines φ uniquely. We have previously discussed all main ingredients of the fast digital shearlet trans-form (FDST) – Fast PPFT, Weighting, and Digital Shearlet Windowing –, and willnow summarize those findings. Depending on the application at hand, a fast inversetransform is required, which we will also detail in the sequel. In fact, we will presenttwo possibilities: the Adjoint FDST and the Inverse FDST depending on whether theweighting allows to use the adjoint for reconstruction or whether an iterative proce-dure is required for higher accuracy. Fig. 5 provides an overview of the main stepsof of the FDST and its inverse. For a more detailed description of FDST, AdjointFDST, and Inverse FDST in form of pseudo-code, we refer to [20].
Fig. 5
Flowcharts of the FDST (left) and its inverse (right).
For the sake of brevity, we now let P , w , and W denote the Fast PPFT fromSubsection 2.1.2, the weighting on the pseudo-polar grid described in Subsection2.2.4, and windowing operator consisting of the application of the shearlet windowsfollowed by 2D-iFFT to each array as detailed in Subsection 2.3.7, respectively. We can summarize the steps of the algorithm FDST as follows: • Step (S1):
For a given image I , apply the Fast PPFT as described in Subsection2.1.2 to obtain the function PI : Ω R → C . • Step (S2):
Apply the square root of an off-line computed weight function w : Ω R → C to PI as described in Subsection 2.2.4, yielding √ wPI : Ω R → C . • Step (S3):
Apply the shearlet windows to the function wPI , followed by a 2D-iFFT to each array to obtain the shearlet coefficients W √ wPI , which we denoteby c ι n , ι , n and c ι j , k , m , j , k , m , ι . igital Shearlet Transforms 25 Assuming that the weight function w used in Step (S2) satisfies the condition inTheorem 1, and using the Parseval frame property of the digital shearlet system(Theorem 2), we obtain ( W √ wP ) (cid:63) W √ wP = P (cid:63) √ w ( W (cid:63) W ) √ wP = P (cid:63) wP = Id . Hence in this case, the FDST, which is abbreviated by W √ wP can be inverted byapplying the adjoint FDST, which cascades the following steps: • Step 1:
For given shearlet coefficients C , i.e., c ι n , ι , n and c ι j , k , m , j , k , m , ι ,compute the linear combination of the shearlet windows with coefficients c ι n , ι , n and c ι j , k , m , j , k , m , ι . This gives the function W (cid:63) C : Ω R → C . • Step 2:
Apply the square root of an off-line computed weight function w : Ω R → C to W (cid:63) C , yielding the function √ wW (cid:63) C : Ω R → C . • Step 3:
Apply the Fast Adjoint PPFT by running the Fast PPFT ‘backwards’.For this, we just notice that the adjoint fractional Fourier transform of a vector c ∈ C N + with respect to a constant α ∈ C is given by F − α N + c . Also, for m > N , the adjoint padding operator E (cid:63) m , N applied to a vector c ∈ C m is given by ( E (cid:63) m , N c )( k ) = c ( k ) , k = − N / , . . . , N / −
1. The Adjoint PPFT gives an image P (cid:63) √ wW (cid:63) C . Normally – as also with the relaxed form of weights debated in Subsection 2.2.2 –the weights will not satisfy the conditions of Theorem 1 precisely. A measure forwhether application of the adjoint is still feasible was already discussed in Subsec-tion 2.2.3 (see also Subsection 4.2). If higher accuracy of the reconstruction is re-quired, one might use iterative methods, such as conjugate gradient methods. Sincethe digital shearlet system forms a Parseval frame, we always have W (cid:63) W √ wP = √ wP . Hence, iterative methods need to be ‘only’ applied to reconstruct an image I fromknowledge of J : = √ wPI , i.e., to solve the equation P (cid:63) wPI = P (cid:63) wJ for I . Since J might not be in the range of P , I is typically computed by solvingthe weighted least square problem min I (cid:107)√ wPI − √ wJ (cid:107) . Since the matrix cor-responding to P (cid:63) P is symmetric positive definite, iterative methods such as theconjugate gradient method are applicable. The conjugate gradient method is thenapplied to the equation Ax = b with A = P (cid:63) wP and b = P (cid:63) wJ . Its performancecan be measured by the condition number of the operator P (cid:63) wP : cond ( P (cid:63) wP ) = λ max ( P (cid:63) wP ) / λ min ( P (cid:63) wP ) , and it turns out that the weight function serves as a pre-conditioner. We remark that this measure is more closely studied in Subsection 4.2.To illustrate the behavior of the weights with respect to this measure, in Table 3we compute cond ( P (cid:63) wP ) for the weight functions arising from Choices 1 and 2 (cf.Subsection 2.2.2) with oversampling rate R =
8. Notice that the condition numbersof P (cid:63) wP are generally smaller than 2. Table 3
Comparison of Choices 1 and 2 based on the performance measure cond ( P (cid:63) wP ) . N
32 64 128 256 512Choice 1 1.379 1.503 1.621 1.731 1.833Choice 2 1.760 1.887 2.001 2.104 N/A
In this section, we will discuss two implementation strategies for computing shearletcoefficients associated with a cone-adapted discrete shearlet system now based on compactly supported shearlets, as introduced in Chapter [1]. Again, one main focuswill be on deriving a digitization which is faithful to the continuum setting.Recall that in the context of wavelet theory, faithful digitization is achieved bythe concept of multiresolution analysis, where scaling and translation are digitizedby discrete operations: Downsampling, upsampling and convolution. In the case ofdirectional transforms however, three types of operators: Scaling, translation anddirection, need to be digitized. In this section, we will pay particular attention toderiving a framework in which each of the three operators is faithfully interpretedas a digitized operation in digital domain. Both approaches will be based on thefollowing digitization strategies: • Scaling and translation: A multiresolution analysis associated with anisotropicscaling A j can be applied for each shear parameter k . • Directionality: A faithful digitization of shear operator S − j / k has to be achievedwith particular care.After stating and discussing the two main obstacles we are facing when consideringcompactly supported shearlets in Subsection 3.1, we present the digital separableshearlet transform (DSST), which is associated with a shearlet system generatedby a separable function alongside with discussions on its properties, e.g., its redun-dancy; see Subsection 3.2. Subsection 3.3 then presents the digital non-separableshearlet transform (DNST), whose shearlet elements are generated by non- separa-ble shearlet generator. igital Shearlet Transforms 27 Compactly supported shearlets have several advantages, and we exemplarily men-tion superior spatial localization and simplified boundary adaptation. However, wehave to face the following two problems:(P1) Compactly supported shearlets do not form a tight frame, which prevents uti-lization of the adjoint as inverse transform.(P2) There does not exist a natural hierarchical structure, mainly due to the appli-cation of a shear matrix, which – unlike for the wavelet transform – does notallow a multiresolution analysis without destroying a faithful adaption of thecontinuum setting.Let us now comment on these two obstacles, before delving into the details ofthe implementation in Subsection 3.2.
Let us first comment on the problem of non-tightness. Letting ( σ i ) i ∈ I denote a framefor L ( R ) – for example, a shearlet frame –, each function f ∈ L ( R ) can bereconstructed from its frame coefficients ( (cid:104) f , σ i (cid:105) ) i ∈ I by f = ∑ i ∈ I (cid:104) f , σ i (cid:105) S − ( σ i ) , where S = ∑ i ∈ I (cid:104)· , σ i (cid:105) σ i is the associated frame operator on L ( R ) , see Chapter [1].However, in case that ( σ i ) i ∈ I does not form a tight frame, it is in general difficult toexplicitly compute the dual frame elements S − ( σ i ) .Nevertheless, it is well known that the inverse frame operator S − can be effec-tively approximated using iterative schemes such as the Conjugate Gradient methodprovided that the frame ( σ i ) i ∈ I has ’good’ frame bounds in the sense of their ratiobeing ‘close’ to 1, see also [24]. Therefore, now focussing on the situation of shear-let frames, we may argue that input data f can be efficiently reconstructed from itsshearlet coefficients, if we use a compactly supported shearlet frame with ’good’frame bounds. In fact, the theoretical frame bounds of compactly supported shearletframes have been theoretically estimated as well as numerically computed in [19].These results were derived for the class of 2D separable shearlet generators ψ al-ready described in Chapter [1], which we briefly recall for the convenience of thereader:For positive integers K and L fixed, let the 1D lowpass filter m be defined by | m ( ξ ) | = ( cos ( πξ )) K L − ∑ n = (cid:18) K − + nn (cid:19) ( sin ( πξ )) n , for ξ ∈ R . Further, define the associated bandpass filter m by | m ( ξ ) | = | m ( ξ + / ) | , ξ ∈ R , and the 1D scaling function φ byˆ φ ( ξ ) = ∞ ∏ j = m ( − j ξ ) , ξ ∈ R . Using the filter m and scaling function φ , we now define the 2D scaling function φ and separable shearlet generator ψ byˆ φ ( ξ , ξ ) = ˆ φ ( ξ ) ˆ φ ( ξ ) and ˆ ψ ( ξ , ξ ) = m ( ξ ) ˆ φ ( ξ ) ˆ φ ( ξ ) . (35)In [19], it was shown that compactly supported shearlets ψ j , k , m generated by theshearlet generator ψ form a frame for L ( C ) ∨ with appropriately chosen parameters K and L , where C = { ξ ∈ R : | ξ / ξ | ≤ , | ξ | ≥ } . This construction is directly extended to construct a cone-adapted discrete shearletframe for L ( R ) (cf. also Chapters [1] and [4]).Table 4 provides some numerically estimated frame bounds in L ( C ) for certainchoice of K and L . It shows that indeed the ratio of the frame bounds of this classof compactly supported shearlet frames is sufficient small for utilizing an iterativescheme for efficient reconstruction; in this sense the frame bounds are ’good’. Table 4
Numerically estimated frame bounds for various choices of the parameters K and L . c and c are the sampling constants in the sampling matrix M c for translation (see Chapter [1]). K L c c B/A39 19 0.90 0.15 4.108439 19 0.90 0.20 4.108539 19 0.90 0.25 4.110439 19 0.90 0.30 4.132839 19 0.90 0.40 5.2495
The frequency covering by compactly supported shearlets ψ j , k , m , | ˆ φ ( ξ ) | + ∑ j ≥ ∑ k ∈ K j | ˆ ψ ( S Tk A j ξ ) | + | ˆ˜ ψ ( ˜ S Tk ˜ A j ξ ) | , is closely related to the ratio of frames bounds and, in particular, which areas infrequency domain cause a larger ratio. This function is illustrated in Fig. 6, whichshows that its upper and lower bounds are as expected well controlled. igital Shearlet Transforms 29 (a) Whole frequency plane. (b) Horizontal cone. (c) Vertical cone. Fig. 6
Frequency covering by shearlets | ˆ ψ j , k , m | : (a) Frequency covering of the entire frequencyplane. (b) Frequency covering of the horizontal cone. (c) Frequency covering of the vertical cone. Let us finally comment on the problem to achieve a hierarchical structuring. To al-low fast implementations, the data structure of the transform is essential. The hierar-chical structure of the wavelet transform associated with a multiresolution analysis,for instance, enables a fast implementation based on filterbanks. In addition, sucha hierarchical ordering provides a full tree structure across scales, which is of par-ticular importance for various applications such as image compression and adaptivePDE schemes. It is in fact mainly due to this property – and the unified treatment ofthe continuum and digital setting – that the wavelet transform became an extremelysuccessful methodology for many practical applications.From a certain viewpoint, shearlets ψ j , k , m can essentially be regarded as waveletsassociated with an anisotropic scale matrix A j , when the shear parameter k is fixed.This observation allows to apply the wavelet transform to compute the shearlet co-efficients, once the shear operation is computed for each shear parameter k . Thisapproach will be undertaken in the digital formulation of the compactly supportedshearlet transform, and, in fact, this approach implements a hierarchical structureinto the shearlet transform. The reader should note that this approach does not leadto a completely hierarchical structured shearlet transform – also compare our dis-cussion at the beginning of this section –, but it will be sufficient for deriving a fastimplementation while retaining a faithful digitization. We now describe a faithful digitization of the continuum domain shearlet transformbased on compactly supported shearlets as introduced in [22], which moreover ishighly computationally efficient.
We start by discussing those theoretical aspects which allow a faithful digitizationof the shearlet transform associated with the shearlet system generated by (35). Forthis, we will only consider shearlets ψ j , k , m for the horizontal cone, i.e., belonging to Ψ ( ψ , c ) . Notice that the same procedure can be applied to compute the shearlet co-efficients for the vertical cone, i.e., those belonging to ˜ Ψ ( ˜ ψ , c ) , except for switchingthe order of variables.To construct a separable shearlet generator ψ ∈ L ( R ) and an associated scalingfunction φ ∈ L ( R ) , let φ ∈ L ( R ) be a compactly supported 1D scaling functionsatisfying φ ( x ) = ∑ n ∈ Z h ( n ) √ φ ( x − n ) (36)for some ‘appropriately chosen’ filter h – we comment on the required condition be-low. An associated compactly supported 1D wavelet ψ ∈ L ( R ) can then be definedby ψ ( x ) = ∑ n ∈ Z g ( n ) √ φ ( x − n ) , (37)where again g is an ‘appropriately chosen’ filter. The selected shearlet generator isthen defined to be ψ ( x , x ) = ψ ( x ) φ ( x ) , (38)and the scaling function by φ ( x , x ) = φ ( x ) φ ( x ) . Let us comment on whether this is indeed a special case of the shearlet generatorsdefined in (35). The Fourier transform of ψ defined in (38) takes the formˆ ψ ( ξ , ξ ) = m ( ξ / ) ˆ φ ( ξ / ) ˆ φ ( ξ / ) , where m is a trigonometric polynomial whose Fourier coefficients are g ( n ) . Weneed to compare this expression with the Fourier transform of the shearlet generator ψ given in (35), which isˆ ψ ( ξ , ξ ) = m ( ξ ) ˆ φ ( ξ ) ˆ φ ( ξ ) , with 1D scaling function φ defined in (36). We remark that this later scaling func-tion is slightly different defined as in (35). This small adaption is for the sake ofpresenting a simpler version of the implementation; essentially the same implemen-tation strategy as the one we will describe can be applied to the shearlet generatorgiven in (35).The filter coefficients h and g are required to be chosen so that ψ satisfies acertain decay condition (cf. [19] of Chapter [1]) to guarantee a stable reconstructionfrom the shearlet coefficients. igital Shearlet Transforms 31 For the signal f ∈ L ( R ) to be analyzed, we now assume that, for J > f is of the form f ( x ) = ∑ n ∈ Z f J ( n ) J φ ( J x − n , J x − n ) . (39)Let us mention that this is a very natural assumption for a digital implementation inthe sense that the scaling coefficients can be viewed as sample values of f – in fact f J ( n ) = f ( − J n ) with appropriately chosen φ . Now aiming towards a faithful digi-tization of the shearlet coefficients (cid:104) f , ψ j , k , m (cid:105) for j = , . . . , J −
1, we first observethat (cid:104) f , ψ j , k , m (cid:105) = (cid:104) f ( S − j / k ( · )) , ψ j , , m ( · ) (cid:105) , (40)and, WLOG we will from now on assume that j / (cid:100) j / (cid:101) or (cid:98) j / (cid:99) would need to be taken. Our observation (40) shows us in fact preciselyhow to digitize the shearlet coefficients (cid:104) f , ψ j , k , m (cid:105) : By applying the discrete sepa-rable wavelet transform associated with the anisotropic sampling matrix A j to thesheared version of the data f ( S − j / k ( · )) . This however requires – compare the as-sumed form of f given in (39) – that f ( S − j / k ( · )) is contained in the scaling space V J = { J φ ( J · − n , J · − n ) : ( n , n ) ∈ Z } . It is easy to see that, for instance, if the shear parameter 2 − j / k is non-integer, thisis unfortunately not the case. The true reason for this failure is that the shear matrix S − j / k does not preserve the regular grid 2 − J Z in V J , i.e., S − j / k ( Z ) (cid:54) = Z . In order to resolve this issue, we consider the new scaling space V kJ + j / , J defined by V kJ + j / , J = { J + / j φ ( S k ( J + j / · − n , J · − n )) : ( n , n ) ∈ Z } . We remark that the scaling space V kJ + j / , J is obtained by refining the regular grid2 − J Z along the x -axis by a factor of 2 j / . With this modification, the new grid2 − J − j / Z × − J Z is now invariant under the shear operator S − j / k , since with Q = diag ( , ) , 2 − J − j / Z × − J Z = − J Q − j / ( Z ) = − J Q − j / ( S k ( Z ))= S − j / k ( − J − j / Z × − J Z ) . This allows us to rewrite f ( S − j / k ( · )) in (40) in the following way. Lemma 1.
Retaining the notations and definitions from this subsection, letting ↑ j / and ∗ denote the 1D upsampling operator by a factor of j / and the 1Dconvolution operator along the x -axis, respectively, and setting h j / ( n ) to be theFourier coefficients of the trigonometric polynomial H j / ( ξ ) = j / − ∏ k = ∑ n ∈ Z h ( n ) e − π i k n ξ , (41) we obtain f ( S − j / k ( x )) = ∑ n ∈ Z ˜ f J ( S k n ) J + j / φ k ( J + j / x − n , J x − n ) , where ˜ f J ( n ) = (( f J ) ↑ j / ∗ h j / )( n ) . The proof of this lemma requires the following result, which follows from thecascade algorithm in the theory of wavelet.
Proposition 1 ([22]).
Assume that φ and ψ ∈ L ( R ) satisfy equations (36) and (37) respectively. For positive integers j ≤ j , we then have j φ ( j x − n ) = ∑ d ∈ Z h j − j ( d − j − j n ) j φ ( j x − d ) (42) and j ψ ( j x − n ) = ∑ d ∈ Z g j − j ( d − j − j n ) j φ ( j x − d ) , (43) where h j and g j are the Fourier coefficients of the trigonometric polynomials H j defined in (41) and G j defined byG j ( ξ ) = (cid:16) j − ∏ k = ∑ n ∈ Z h ( n ) e − π i k n ξ (cid:17)(cid:16) ∑ n ∈ Z g ( n ) e − π i j − n ξ (cid:17) for j > fixed.Proof (Proof of Lemma 1). Equation (42) with j = J and j = J + j / J / φ ( J x − n ) = ∑ d ∈ Z h J − j / ( d − j / n ) J / + j / φ ( J + j / x − d ) . (44)Also, since φ is a 2D separable function of the form φ ( x , x ) = φ ( x ) φ ( x ) , wehave that f ( x ) = ∑ n ∈ Z (cid:16) ∑ n ∈ Z f J ( n , n ) J / φ ( J x − n ) (cid:17) J / φ ( J x − n ) . By (44), we obtain f ( x ) = ∑ n ∈ Z ˜ f J ( n ) J + j / φ ( J Q j / x − n ) , igital Shearlet Transforms 33 where Q = diag ( , ) . Using Q j / S − j / k = S k Q j / , this finally implies f ( S − j / k ( x )) = ∑ n ∈ Z ˜ f J ( n ) J + j / φ ( J Q j / S − j / k ( x ) − n )= ∑ n ∈ Z ˜ f J ( n ) J + j / φ ( S k ( J Q j / x − S − k n ))= ∑ n ∈ Z ˜ f J ( S k n ) J + j / φ ( S k ( J Q j / x − n )) . The lemma is proved. (cid:117)(cid:116)
The second term to be digitized in (40) is the shearlet ψ j , k , m itself. A direct corol-lary from Proposition 1 is the following result. Lemma 2.
Retaining the notations and definitions from this subsection, we obtain ψ j , k , m ( x ) = ∑ d ∈ Z g J − j ( d − J − j m ) h J − j / ( d − J − j / m ) J + j / φ ( J Q j / x − d ) . As already indicated before, we will make use of the discrete separable wavelettransform associated with an anisotropic scaling matrix, which, for j and j > c ∈ (cid:96) ( Z ) , we define by W j , j ( c )( n , n ) = ∑ m ∈ Z g j ( m − j n ) h j ( m − j n ) c ( m , m ) , ( n , n ) ∈ Z . (45)Finally, Lemmata 1 and 2 yield the following digitizable form of the shearletcoefficients (cid:104) f , ψ j , k , m (cid:105) . Theorem 3 ([22]).
Retaining the notations and definitions from this subsection, andletting ↓ j / be 1D downsampling by a factor of j / along the horizontal axis, weobtain (cid:104) f , ψ j , k , m (cid:105) = W J − j , J − j / (cid:16)(cid:16) ( ˜ f J ( S k · ) ∗ Φ k ) ∗ h j / (cid:17) ↓ j / (cid:17) ( m ) , where Φ k ( n ) = (cid:104) φ ( S k ( · )) , φ ( · − n ) (cid:105) for n ∈ Z , and h j / ( n ) = h j / ( − n ) . Computing the shearlet coefficients using Theorem 3 now restricts to applying thediscrete separable wavelet transform (45) associated with the sampling matrix A j to the scaling coefficients S d − j / k ( f J )( n ) : = (cid:16) ( ˜ f J ( S k · ) ∗ Φ k ) ∗ h j / (cid:17) ↓ j / ( n ) for f J ∈ (cid:96) ( Z ) . (46)Before we state the explicit steps necessary to achieve this, let us take a closer lookat the scaling coefficients S d − j / k ( f J ) , which can be regarded as a new sampling of Fig. 7
Illustration of appli-cation of the digital shearoperator S d − / : The dashedlines correspond to the refine-ment of the integer grid. Thenew sample values lie on theintersections of the shearedlines associated with S / withthis refined grid. the data f J on the integer grid Z by the digital shear operator S d − j / k . This procedureis illustrated in Fig. 7 in the case 2 − j / k = − / Φ k ( n ) in (46) can in fact be easilyprecomputed for each shear parameter k . For a practical implementation, one maysometimes even skip this additional convolution step assuming that Φ k = χ ( , ) .Concluding, the implementation strategy for the DSST cascades the followingsteps: • Step 1:
For given input data f J , apply the 1D upsampling operator by a factorof 2 j / at the finest scale j = J . • Step 2:
Apply 1D convolution to the upsampled input data f J with 1D lowpassfilter h j / at the finest scale j = J . This gives ˜ f J . • Step 3:
Resample ˜ f J to obtain ˜ f J ( S k ( n )) according to the shear sampling matrix S k at the finest scale j = J . Note that this resampling step is straightforward,since the integer grid is invariant under the shear matrix S k . • Step 4:
Apply 1D convolution to ˜ f J ( S k ( n )) with h j / followed by 1D downsam-pling by a factor of 2 j / at the finest scale j = J . • Step 5:
Apply the separable wavelet transform W J − j , J − j / across scales j = , , . . . , J − Since the digital realization of a shear matrix S − j / k by the digital shear operator S d − j / k is crucial for deriving a faithful digitization of the continuum domain shearlettransform, we will devote this subsection to a closer analysis.We start by remarking that in fact in the continuum domain, at least two opera-tors exist which naturally provide directionality: Rotation and shearing. Rotation isa very convenient tool to provide directionality in the sense that it preserves impor-tant geometric information such as length, angles, and parallelism. However, thisoperator does not preserve the integer lattice, which causes severe problems for dig-itization. In contrast to this, a shear matrix S k does not only provide directionality,but also preserves the integer lattice when the shear parameter k is integer. Thus, it igital Shearlet Transforms 35 is conceivable to assume that directionality can be naturally discretized by using ashear matrix S k .To start our analysis of the relation between a shear matrix S − j / k and the asso-ciated digital shear operator S d − j / k , let us consider the following simple example:Set f c = χ { x : x = } . Then digitize f c to obtain a function f d defined on Z by set-ting f d ( n ) = f c ( n ) for all n ∈ Z . For fixed shear parameter s ∈ R , apply the sheartransform S s to f c yielding the sheared function f c ( S s ( · )) . Next, digitize also thisfunction by considering f c ( S s ( · )) | Z . The functions f d and f c ( S s ( · )) | Z are illus-trated in Fig. 8 for s = − /
4. We now focus on the problem that the integer lattice is (a) (b)
Fig. 8 (a) Original image f d ( n ) . (b) Sheared image f c ( S − / n ) . not invariant under the shear matrix S / . This prevents the sampling points S / ( n ) , n ∈ Z from lying on the integer grid, which causes aliasing of the digitized im-age f c ( S − / ( · )) | Z as illustrated in Fig. 9(a). In order to avoid this aliasing effect,the grid needs to be refined by a factor of 4 along the horizontal axis followed bycomputing sample values on this refined grid.More generally, when the shear parameter is given by s = − − j / k , one canessentially avoid this directional aliasing effect by refining a grid by a factor of2 j / along the horizontal axis followed by computing interpolated sample values onthis refined grid. This ensures that the resulting grid contains the sampling points (( − j / k ) n , n ) for any n ∈ Z and is preserved by the shear matrix S − − j / k .This procedure precisely coincides with the application of the digital shear oper-ator S d − j / k , i.e., we just described Steps 1 – 4 from Subsection 3.2.2 in which thenew scaling coefficients S d − j / k ( f J )( n ) are computed.Let us come back to the exemplary situation of f c = χ { x : x = } and S − / westarted our excursion with and compare f c ( S − / ( · )) | Z with S d − / ( f d ) | Z obtainedby applying the digital shear operator S d − / to f d . And, in fact, the directional alias-ing effect on the digitized image f c ( S − / ( n )) in frequency illustrated in Fig. 9(a)is shown to be avoided in Fig. 9 (b)-(c) by considering S d − / ( f d ) | Z . Thus applica-tion of the digital shear operator S d − j / k allows a faithful digitization of the shearingoperator associated with the shear matrix S − j / k . Fig. 9 (a) Aliased image: DFT of f c ( S − / ( n )) . (b) De-aliased image: S d − / ( f d )( n ) . (c) De-aliasedimage: DFT of S d − / ( f d )( n ) . One of the main issues which practical applicability requires is controllable redun-dancy. To quantify the redundancy of the discrete shearlet transform, we assume thatthe input data f is a finite linear combination of translates of a 2D scaling function φ at scale J as follows: f ( x ) = J − ∑ n = J − ∑ n = d n φ ( J x − n ) as it was already the hypothesis in (39). The redundancy – as we view it in ouranalysis – is then given by the number of shearlet elements necessary to represent f . Furthermore, to state the result in more generality, we allow an arbitrary samplingmatrix M c = diag ( c , c ) for translation, i.e., consider shearlet elements of the form ψ j , k , m ( · ) = j ψ ( S k A j · − M c m ) . We then have the following result.
Proposition 2 ([22]).
The redundancy of the DSST is (cid:16) (cid:17)(cid:16) c c (cid:17) . Proof.
For this, we first consider shearlet elements for the horizontal cone for afixed scale j ∈ { , . . . , J − } . We observe that there exist 2 j / + shearing indices k and 2 j · j / · ( c c ) − translation indices associated with the scaling matrix A j andthe sampling matrix M c , respectively. Thus, 2 j + ( c c ) − shearlet elements fromthe horizontal cone are required for representing f . Due to symmetry reasons, werequire the same number of shearlet elements from the vertical cone. Finally, about c − translates of the scaling function φ are necessary at the coarsest scale j = (cid:16) c c (cid:17)(cid:16) J − ∑ j = j + (cid:17) = (cid:16) c c (cid:17)(cid:16) J + (cid:17) igital Shearlet Transforms 37 The redundancy of each shearlet frame can now be computed as the ratio of thenumber of coefficients d n and this number. Letting J → ∞ proves the claim. (cid:117)(cid:116) As an example, choose a translation grid with parameters c = c = . / A further essential characteristics is the computational complexity (see also Subsec-tion 4.6), which we now formally compute for the discrete shearlet transform.
Proposition 3 ([19]).
The computational complexity of the DSST isO ( log ( / ( L / − )) L · N ) . Proof.
Analyzing Steps 1 – 5 from Subsection 3.2.2, we observe that the most timeconsuming step is the computation of the scaling coefficients in Steps 1 – 4 for thefinest scale j = J . This step requires 1D upsampling by a factor of 2 j / followedby 1D convolution for each direction associated with the shear parameter k . Letting L denote the total number of directions at the finest scale j = J , and N the size of2D input data, the computational complexity for computing the scaling coefficientsin Steps 1 – 4 is O ( j / L · N ) . The complexity of the discrete separable wavelettransform associated with A j for Step 5 requires O ( N ) operations, wherefore it isnegligible. The claim follows from the fact that L = ( · j / + ) . (cid:117)(cid:116) It should be noted that the total computational cost depends on the number L ofshear parameters at the finest scale j = J , and this total cost grows approximately bya factor of L as L is increased. It should though be emphasized that L can be chosenin such a way that this shearlet transform is favorably comparable to other redun-dant directional transforms with respect to running time as well as performance. Areasonable number of directions at the finest scale is 6, in which case the constantfactor 2 log ( / ( L / − )) in Proposition 3 equals 1. Hence in this case the running timeof this shearlet transform is only about 6 times slower than the discrete orthogo-nal wavelet transform, thereby remains in the range of the running time of otherdirectional transforms. In Subsection 3.1.1, we already discussed that this transform is not an isometry,wherefore the adjoint cannot be used as an inverse transform. However, the ‘good’ratio of the frame bounds in the sense as detailed in Subsection 3.1.1 leads to a fastconvergence rate of iterative methods such as the conjugate gradient method. Let usmention that using the conjugate gradient method basically requires computing theforward DSST and its adjoint, and we refer to [24] and also Subsection 2.4.2 formore details.
In this section, we describe an alternative approach to derive a faithful digitaliza-tion of a discrete shearlet transform associated with compactly supported shearlets.This algorithmic realization, which was developed in [23], resolves the followingdrawbacks of the DSST: • Since this transform is not based on a tight frame, an additional computationaleffort is necessary to approximate the inverse of the shearlet transform by itera-tive methods. • Computing the interpolated sampling values in (46) requires additional compu-tational costs. • This shearlet transform is not shift-variant, even when downsampling associatedwith A j is omitted.We emphasize that although this alternative approach resolves these problems, thealgorithm DSST provides a much more faithful digitalization in the sense that theshearlet coefficients can be exactly computed in this framework.The main difference between DSST and DNST will be to exploit non-separable shearlet generators, which give more flexibility. We start by introducing the non-separable shearlet generators utilized in DNST.First, for each scale parameter j ≥
0, define the shearlet generator ψ non j byˆ ψ non j ( ξ ) = P J − j / ( ξ ) ˆ ψ ( ξ ) , where P (cid:96) ( ξ ) = P ( (cid:96) + ξ , ξ ) for (cid:96) ≥ P is a 2Dfan filter (c.f. [12]). For an illustration of P we refer to Fig. 10(a). This in turn definesshearlets ψ non j , k , m generated by non-separable generator functions ψ non j for each scaleindex j ≥ ψ non j , k , m ( x ) = j ψ non j ( S k A j x − M c j m ) , where M c j is a sampling matrix given by M c j = diag ( c j , c j ) and c j and c j are sam-pling constants for translation.One major advantage of these shearlets ψ non j , k , m is the fact that a fan filter en-ables refinement of the directional selectivity in frequency domain at each scale.Fig. 10(a)-(b) show the refined essential support of ˆ ψ non j , k , m as compared to shearlets ψ j , k , m arising from a separable generator as in Subsection 3.2.1. igital Shearlet Transforms 39 −1 −0.5 0 0.5 1−1−0.500.5100.51 F x F y M agn i t ude (a) (b) (c) Fig. 10 (a) Magnitude response of 2D fan filter. (b)Non-separable shearlet ψ non j , k , m . (c)Separableshearlet ψ j , k , m . Next, our aim is to derive a digital formulation of the shearlet coefficients (cid:104) f , ψ non j , k , m (cid:105) for a function f as given in (39). We will only discuss the case of shearlet coefficientsassociated with A j and S k ; the same procedure can be applied for ˜ A j and ˜ S k exceptfor switching the order of variables x and x .In Subsection 3.2.3, we discretized a sheared function f ( S − j / k · ) using the dig-ital shear operator S d − j / k as defined in (46). In this implementation, we walk adifferent path. We digitize the shearlets ψ non j , k , m ( · ) = ψ non j , , m ( S − j / k · ) by combiningmultiresolution analysis and digital shear operator S d − j / k to digitize the wavelet ψ non j , , m and the shear operator S − j / k , respectively. This yields digitized shearlet fil-ters of the form ψ d j , k ( n ) = S d − j / k (cid:16) p J − j / ∗ w j (cid:17) ( n ) , where w j is the 2D separable wavelet filter defined by w j ( n , n ) = g J − j ( n ) · h J − j / ( n ) and p J − j / ( n ) are the Fourier coefficients of the 2D fan filter P J − j / .The DNST associated with the non-separable shearlet generators ψ non j is then givenby DNST j , k ( f J )( n ) = ( f J ∗ ψ d j , k )( J − j c j n , J − j / c j n ) , for f J ∈ (cid:96) ( Z ) . We remark that the discrete shearlet filters ψ d j , k are computed by using a similarideas as in Subsection 3.2.1. As before, those filter coefficients can be precomputedto avoid additional computational effort.Further notice that by setting c j = j − J and c j = j / − J , the DNST simply be-comes a 2D convolution. Thus, in this case, DNST is shear invariant. In case that c j = j − J and c j = j / − J , the dual shearlet filters ˜ ψ d j , k can be easilycomputed by deconvolution, and we obtain the reconstruction formula f J = ∑ j , k ( f J ∗ ψ d j , k ) ∗ ˜ ψ d j , k . Thus, no iterative methods are required for the inverse DNST.The frequency response of a discrete shearlet filter ψ d j , k and its dual ˜ ψ d j , k is illus-trated in Fig. 11. We observe that primal and dual shearlet filters behave similarly inthe sense that both of filters are very well localized in frequency. Fig. 11
Magnitude response of shearlet filter ψ d j , k and its dual filter ˜ ψ d j , k . We next present the framework for quantifying performance of implementations ofdirectional transforms, which was originally introduced in [20, 14]. This set of testmeasures was designed to analyze particular well-understood properties of a givenalgorithm, which in this case are the desiderata proposed at the beginning of thischapter. This framework was moreover introduced to serve as a tool for tuning theparameters of an algorithm in a rational way and as an objective common groundfor comparison of different algorithms. The performance of the three algorithmsFDST, DSST, and DNST will then be tested with respect to those measures. Thiswill give us insight into their behavior with respect to the analyzed characteristics,and also allow a comparison. However, the test values of these three algorithms willalso show the delicateness of designing such a testing framework in a fair manner,since due to the complexity of algorithmic realizations it is highly difficile to doeach aspect of an algorithm justice. It should though be emphasized that – apartfrom being able to rationally tune parameters – such a framework of quantifyingperformance is essential for an objective judgement of algorithms. The codes of allmeasures are available in
ShearLab .In the following, S shall denote the transform under consideration, S (cid:63) its adjoint,and, if iterative reconstruction is tested, G A J shall stand for the solution of the matrixproblem AI = J using the conjugate gradient method with residual error set to be10 − . Some measures apply specifically to transforms utilizing the pseudo-polargrid, for which purpose we introduce the notation P for the pseudo-polar Fourier igital Shearlet Transforms 41 transform, w shall denote the weighting applied to the values on the pseudo-polargrid, and W shall be the windowing with additional 2D iFFT. We require the transform to be the precise implementation of a theory for digital dataon a pseudo-polar grid. In addition, to ensure numerical accuracy, we provide thefollowing measure. This measure is designed for transforms utilizing the pseudo-polar grid.
Measure 1
Generate a sequence of (of course, one can choose any rea-sonable integer other than 5) random images I , . . . , I on a pseudo-polargrid for N = and R = with standard normally distributed entries. Ourquality measure will then be the Monte Carlo estimate for the operator norm (cid:107) W (cid:63) W − Id (cid:107) op given byM alg = max i = ,..., (cid:107) W (cid:63) W I i − I i (cid:107) (cid:107) I i (cid:107) . This measure applies to the FDST – not to the DSST or DNST – , for which weobtain M alg = . E − . This confirms that the windowing in the FDST is indeed up to machine precision aParseval frame, which was already theoretically confirmed by Theorem 2.
We next test the pseudo-polar transform itself which might be used in the algorithmunder consideration. For this, we will provide three different measures, each beingdesigned to test a different aspect.
Measure 2 • Closeness to isometry . Generate a sequence of random images I , . . . , I of size × with standard uniformly distributed entries. Our qualitymeasure will then be the Monte Carlo estimate for the operator norm (cid:107) P (cid:63) wP − Id (cid:107) op given by M isom = max i = ,..., (cid:107) P (cid:63) wPI i − I i (cid:107) (cid:107) I i (cid:107) . • Quality of preconditioning . Our quality measure will be the spread of theeigenvalues of the Gram operator P (cid:63) wP given byM isom = λ max ( P (cid:63) wP ) λ min ( P (cid:63) wP ) . • Invertibility . Our quality measure will be the Monte Carlo estimate for theinvertibility of the operator √ wP using conjugate gradient method G √ wP (residual error is set to be − , here G A J means solving matrix problemAI = J using conjugate gradient method) given byM isom = max i = ,..., (cid:107) G √ wP √ wPI i − I i (cid:107) (cid:107) I i (cid:107) . This measure applies to the FDST – not to the DSST or DNST – , for which weobtain the following numerical results, see Table 5.
Table 5
The numerical results for the test on isometry of the pseudo-polar transform. M isom M isom M isom FDST 9.3E-4 1.834 3.3E-7
The slight isometry deficiency of M isom ≈ We now test the overall frame behavior of the system defined by the transform.These measures now apply to more than pseudo-polar based transforms, in particu-lar, to FDST, DSST, and DNST.
Measure 3
Generate a sequence of random images I , . . . , I of size × with standard uniformly distributed entries. Our quality measure will thentwo-fold: igital Shearlet Transforms 43 • Adjoint transform . The measure will be the Monte Carlo estimate for theoperator norm (cid:107) S (cid:63) S − Id (cid:107) op given byM tight = max i = ,..., (cid:107) S (cid:63) SI i − I i (cid:107) (cid:107) I i (cid:107) . • Iterative reconstruction . Using conjugate gradient G √ wP , our measurewill be given by M tight = max i = ,..., (cid:107) G √ wP W (cid:63) SI i − I i (cid:107) (cid:107) I i (cid:107) . The following table, Table 6, presents the performance of FDST and DNST withrespect to these quantitative measures.
Table 6
The numerical results for the test on Parseval property. M tight M tight FDST 9.9E-4 3.8E-7DSST 1.9920 1.2E-7DNST 0.1829 5.8E-16 (with dual filters)
The transform FDST is nearly tight as indicated by the measures M tight = . M tight = . M tight = . M tight = . M tight = . . The next measure is designed to test the degree to which the analyzing elements,here phrased in terms of shearlets but can be extended to other analyzing elements,are space-frequency localized.
Measure 4
Let I be a shearlet in a × image centered at the origin ( , ) with slope of scale , i.e., σ , , + σ , , . Our quality measurewill be four-fold: • Decay in spatial domain . We compute the decay rates d , . . . , d alonglines parallel to the y-axis starting from the line [ , : ] and the de-cay rates d , . . . , d with x and y interchanged. By decay rate, forinstance, for the line [
257 : 512 , ] , we first compute the smallest mono-tone majorant M ( x , ) , x = , . . . , – note that we could also choosean average amplitude here or a different ‘envelope’ – for the curve | I ( x , ) | , x = , . . . , . Then the decay rate is defined to be the aver-age slope of the line, which is a least square fit to the curve log ( M ( x , )) ,x = , . . . , . Based on these decay rates, we choose our measure tobe the average of the decay ratesM decay = ∑ i = ,..., d i . • Decay in frequency domain . Here we intend to check whether the Fouriertransform of I is compactly supported and also the decay. For this, let ˆ Ibe the 2D-FFT of I and compute the decay rates d i , i = , . . . , asbefore. Then we define the following two measures:– Compactly supportedness .M supp = max | u | , | v |≤ | ˆ I ( u , v ) | max u , v | ˆ I ( u , v ) | . – Decay rate . M decay = ∑ i = ,..., d i . • Smoothness in spatial domain . We will measure smoothness by the aver-age of local H¨older regularity. For each ( u , v ) , we compute M ( u , v ) = | I ( u , v ) − I ( u , v ) | , < max {| u − u | , | v − v |} ≤ . Then the local H¨olderregularity α u , v is the least square fit to the curve log ( | M ( u , v ) | ) . Thenour smoothness measure is given byM smooth = ∑ u , v α u , v . • Smoothness in frequency domain . We compute the smoothness now for ˆ I,the 2D-FFT of I to obtain the new α u , v and define our measure to beM smooth = ∑ u , v α u , v . igital Shearlet Transforms 45 Let us now analyze the space-frequency localization of the shearlets utilized inFDST, DSST and DNST by these measures. The numerical results are presented inTable 7.
Table 7
The numerical results for the test on space-frequency localization. M decay M supp M decay M smooth M smooth FDST -1.920 5.5E-5 -3.257 1.319 0.734DSST − ∞ − ∞ The shearlet elements associated with FDST are band-limited and those asso-ciated with DSST and DNST are compactly supported, which is clearly indicatedby the values derived for M decay , M supp , and M decay . It would be expected that M decay = − ∞ for FDST due to the band-limitedness of the associated shearlets.The shearlet elements are however defined by their Fourier transform on a pseudo-polar grid, whereas the measure M decay is taken after applying the 2D-FFT to theshearlets resulting in data on a cartesian grid, in particular, yielding a non-preciselycompactly supported function.The test values for M smooth and M smooth show that the associated shearlets aremore smooth in spatial domain for FDST than for DSST and DNST, with the re-versed situation in frequency domain. Shearing naturally occurs in digital imaging, and it can – in contrast to rotation – beprecisely realized in the digital domain. Moreover, for the shearlet transform, shearinvariance can be proven and the theory implies (cid:68) j / ψ ( S − k A j · − m ) , f ( S s · ) (cid:69) = (cid:68) j / ψ ( S − k + j s A j · − m ) , f (cid:69) . We therefore expect to see this or a to the specific directional transform adapted be-havior. The degree to which this goal is reached is tested by the following measure.
Measure 5
Let I be an × image with an edge through the origin ( , ) of slope . Given − ≤ s ≤ , generates an image I s : = I ( S s · ) andlet S j be the set of all possible scales j such that j s ∈ Z . Our quality measurewill then be the curve M shear , j = max − j < k , k + j s < j (cid:107) C j , k ( SI s ) − C j , k + j s ( SI ) (cid:107) (cid:107) I (cid:107) , scale j ∈ S j , where C j , k is the shearlet coefficients at scale j and shear k. We present our results in Table 8.
Table 8
The numerical results for the test on shear invariance. M shear , M shear , M shear , M shear , FDST 1.6E-5 1.8E-4 0.002 0.003
This table shows that the FDST is indeed almost shear invariant. A closer inspec-tion shows that M shear , and M shear , are relatively small compared to the measure-ments with respect to finer scales M shear , and M shear , . The reason for this is thealiasing effect which shifts some energy to the high frequency part near the bound-ary away from the edge in the frequency domain.We did not test DSST and DNST with respect to this measure, since these trans-forms show a different – not included in this Measure 5 – type of shear invariancebehavior. Speed is one of the most fundamental properties of each algorithm to analyze. Here,we test the speed up to a size of N =
512 which regard as sufficient to computingthe complexity.
Measure 6
Generate a sequence of random images I i , i = , . . . , of size i × i with standard normally distributed entries. Let s i be the speed of theshearlet transform S applied to I i . Our hypothesis is that the speed behaveslike s i = c · ( i ) d ; i being the size of the input. Let now ˜ d a be the averageslope of the line, which is a least square fit to the curve i (cid:55)→ log ( s i ) . Let alsof i be the 2fft applied to I i , i = , . . . , . Our quality measure will then be three-fold: • Complexity . M speed = ˜ d a . • The Constant . igital Shearlet Transforms 47 M speed = ∑ i = s i ( i ) M speed , . • Comparison with 2D-FFT . M speed = ∑ i = s i f i . Table 9 presents the results of testing FDST, DSST and DNST with respect tothese speed measures.
Table 9
The numerical results for the test on speed. M speed M speed M speed FDST 1.156 9.3E-6 280.560DSST 0.821 4.5E-3 88.700DNST 1.081 9.9E-8 40.519
To interpret these results correctly, we remark that the DNST was tested only withtest images I i for i = , . . . ,
9, since it can not be implemented for small size images.Interestingly, the results also show that the 2D FFT based convolution makes DNSTcomparable to DSST with respect to these speed measures, although it is muchmore redundant than DSST. Finally, the results show that FDST is comparable withboth DSST and DNST with respect to complexity measure M speed . From this, itis conceivable to assume that FDST is highly comparable with respect to speedfor large scale computations. The larger value M speed = .
560 appears due tothe fact that the FDST employs fractional Fourier transforms on an oversampledpseudo-polar grid of size.
One major advantage of directional transforms is their sensitivity with respect togeometric features alongside with their ability to sparsely approximate those (cf.Chapter [4]). This measure is designed to analyze this property.
Measure 7
Let I , . . . , I be × images of an edge through the origin ( , ) and of slope [ − , − . , , . , ] and the transpose of the middlethree, and let c i , j be the associated shearlet coefficients for image I i and scalej. Our quality measure will two-fold: • Decay of significant coefficients . Consider the curve ∑ i = max | c i , j (of analyzing elements aligned with the line) | , scale j , let d be the average slope of the line, which is a least square fit to log ofthis curve, and define M geo = d . • Decay of insignificant coefficients . Consider the curve ∑ i = max | c i , j (of all other analyzing elements) | , scale j , let d be the average slope of the line, which is a least square fit to log ofthis curve, and define M geo = d . Table 10 shows the numerical test results for FDST, DSST, and DNST.
Table 10
The numerical results for the test on geometric exactness. M geo M geo FDST -1.358 -2.032DSST -0.002 -0.030DNST -0.019 -0.342
As expected, the decay rate of the insignificant shearlet coefficients of FDST,i.e., the ones not aligned with the line singularity, measured by M geo ≈ -2.032 ismuch larger than the decay rate of the significant shearlet coefficients measured by M geo ≈ -1.358. Notice that this difference is even more significant in the case ofthe DSST and DNST. To analyze robustness of an algorithm, we choose thresholding as the most commonimpact on a sequence of transform coefficients.
Measure 8
Let I be the regular sampling of a Gaussian function with mean 0and variance on {− , } generating an × -image. igital Shearlet Transforms 49 • Thresholding 1 . Our first quality measure will be the curveM thres , p = (cid:107) G √ wP W (cid:63) thres , p SI − I (cid:107) (cid:107) I (cid:107) , where thres , p discards · ( − − p ) percent of the coefficients (p =[ ] ). • Thresholding 2 . Our second quality measure will be the curveM thres , p = (cid:107) G √ wP W (cid:63) thres , p SI − I (cid:107) (cid:107) I (cid:107) , where thres , p sets all those coefficients to zero with absolute values be-low the threshold m ( − − p ) with m being the maximal absolute valueof all coefficients. (p = [ .
001 : 0 .
01 : 0 . ] ) Table 11 shows that even if we discard 100 ( − − ) ∼ .
9% of the FDST co-efficients, the original image is still well approximated by the reconstructed image.Thus the number of the significant coefficients is relatively small compared to thetotal number of shearlet coefficients. From Table 12, we note that knowledge ofthe shearlet coefficients with absolute value greater than m ( − / . )( ∼ .
1% ofcoefficients) is sufficient for precise reconstruction.DNST shows a similar behavior with worse values for relatively large p . Itshould be however emphasized that firstly, the redundancy of DNST used in thistest is 25 and this is lower than the redundancy of FDST, which is about 71. Thiseffect can be more strongly seen by the test results of DSST whose redundancy with4 even much smaller. Secondly, a significant part of the low frequency coefficientsin both DSST and DNST will be removed by a relatively large threshold, since theratio between the number of the low frequency coefficients and the total number ofcoefficients is much higher than FDST. This prohibits a similarly good reconstruc-tion of a Gaussian function.This test in particular shows the delicateness of comparing different algorithmsby merely looking at the test values without a rational interpretation; in this case,without considering the redundancy and the ratio between the number of the lowfrequency coefficients and the total number of coefficients. Table 11
The numerical results for M thres , p . p Table 12
The numerical results for M thres , p . p Acknowledgements
The first author would like to thank David Donoho and Morteza Shahramfor many inspiring discussions on topics in this area. She also acknowledges partial supportby Deutsche Forschungsgemeinschaft (DFG) Grant SPP-1324 KU 1446/13 and DFG Grant KU1446/14. The second author was supported by DFG Grant SPP-1324 KU 1446/13, and the thirdauthor was supported by DFG Grant KU 1446/14.
References
1. Introduction2. Chapter on Applications.3. Chapter on ShearletMRA.4. Chapter on Sparse Approximation.5. A. Averbuch, R. R. Coifman, D. L. Donoho, M. Israeli, and Y. Shkolnisky,
A framework fordiscrete integral transformations I – the pseudo-polar Fourier transform , SIAM J. Sci. Com-put. (2008), 764–784.6. D. H. Bailey and P. N. Swarztrauber, The fractional Fourier transform and applications , SIAMReview, (1991), 389–404.7. E. J. Cand`es, L. Demanet, D. L. Donoho and L. Ying, Fast discrete curvelet transforms,
Mul-tiscale Model. Simul. (2006), 861–899.8. E. J. Cand`es and D. L. Donoho, Ridgelets: a key to higher-dimensional intermittency?,
Phil.Trans. R. Soc. Lond. A. (1999), 2495–2509.9. E. J. Cand`es and D. L. Donoho,
New tight frames of curvelets and optimal representations ofobjects with C singularities , Comm. Pure Appl. Math. (2004), 219–266.10. E. J. Cand`es and D. L. Donoho, Continuous curvelet transform: I. Resolution of the wavefrontset , Appl. Comput. Harmon. Anal. (2005), 162–197.11. E. J. Cand`es and D. L. Donoho, Continuous curvelet transform: II. Discretization of frames ,Appl. Comput. Harmon. Anal. (2005), 198–222.12. M. N. Do and M. Vetterli, The contourlet transform: an efficient directional multiresolutionimage representation , IEEE Trans. Image Process. (2005), 2091–2106.13. D. L. Donoho, Wedgelets: nearly minimax estimation of edges , Ann. Statist. (1999), 859–897.14. D. L. Donoho, G. Kutyniok, M. Shahram, and X. Zhuang, A rational design of a digital shear-let transform , Proceeding of the 9th International Conference on Sampling Theory and Appli-cations, Singapore, 2011.15. D. L. Donoho, A. Maleki, M. Shahram, V. Stodden, and I. Ur-Rahman,
Fifteen years of Re-producible Research in Computational Harmonic Analysis, preprint.16. G. Easley, D. Labate, and W.-Q Lim,
Sparse directional image representations using the dis-crete shearlet transform , Appl. Comput. Harmon. Anal. (2008), 25–46.17. B. Han, G. Kutyniok, and Z. Shen, A unitary extension principle for Shearlet Systems , preprint.18. E. Hewitt and K.A. Ross,
Abstract Harmonic Analysis I, II , Springer-Verlag, Berlin/ Heidel-berg/New York, 1963.igital Shearlet Transforms 5119. P. Kittipoom, G. Kutyniok, and W.-Q Lim,
Construction of Compactly Supported ShearletFrames . J. Fourier Anal. Appl., 2010, to appear.20. G. Kutyniok, M. Shahram, and X. Zhuang,
ShearLab: A Rational Design of a DigitalParabolic Scaling Algorithm , preprint.21. G. Kutyniok and T. Sauer,
Adaptive directional subdivision schemes and Shearlet Multireso-lution Analysis , preprint.22. W.-Q Lim,
The Discrete Shearlet Transform: A new directional transform and compactly sup-ported shearlet frames , IEEE Trans. Imag. Proc. (2010), 1166–1180.23. W.-Q Lim, Shift Invariant Shearlet Transform , preprint.24. S. Mallat,