ShearLab: A Rational Design of a Digital Parabolic Scaling Algorithm
Abstract
Multivariate problems are typically governed by anisotropic features such as edges in images. A common bracket of most of the various directional representation systems which have been proposed to deliver sparse approximations of such features is the utilization of parabolic scaling. One prominent example is the shearlet system. Our objective in this paper is threefold: We firstly develop a digital shearlet theory which is rationally designed in the sense that it is the digitization of the existing shearlet theory for continuous data. This implicates that shearlet theory provides a unified treatment of both the continuum and digital realm. Secondly, we analyze the utilization of pseudopolar grids and the pseudopolar Fourier transform for digital implementations of parabolic scaling algorithms. We derive an isometric pseudopolar Fourier transform by careful weighting of the pseudopolar grid, allowing exploitation of its adjoint for the inverse transform. This leads to a digital implementation of the shearlet transform; an accompanying Matlab toolbox called ShearLab is provided. And, thirdly, we introduce various quantitative measures for digital parabolic scaling algorithms in general, allowing one to tune parameters and objectively improve the implementation as well as compare different directional transform implementations. The usefulness of such measures is exemplarily demonstrated for the digital shearlet transform.
Key words. Curvelets, digital shearlet system, directional representation system, fast digital shearlet transform, parabolic scaling, performance measures, software package, tight frames
AMS subject classifications. Primary 42C40; Secondary 42C15, 65K99, 65T60, 65T99, 94A08
1 Introduction
In recent years, applied harmonic analysts have introduced several approaches for directional representations of image data, each one with the intent of efficiently representing highly anisotropic image features. Examples include curvelets [6, 7, 3], contourlets [9], and shearlets [14, 26]. These proposals are inspired by elegant results in theoretical harmonic analysis, which study functions defined on the continuum plane (i.e., not digital images) and address problems of efficiently representing certain types of functions and operators. One set of inspiring results concerns the possibility of highly compressed representations of ‘cartoon’ images, i.e., functions which are piecewise smooth with singularities along smooth curves. Another set of results concerns the possibility of highly compressed representations of wave propagation operators. In ‘continuum theory’, anisotropic directional transforms can significantly outperform wavelets in important ways.
Accordingly, one hopes that a digital implementation of such ideas would also deliver performance benefits over wavelet algorithms in realworld settings. Anticipated applications include [20], where missing sensors cause incomplete measurements, and the problem of texture/geometry separation in image processing – for example in astronomy when images of galaxies require separated analyses of stars, filaments, and sheets [29, 10].
In many cases, however, there are no publicly available implementations of such ideas, or the available implementations are only sketchily tested, or the available implementations are only vaguely related to the continuum transforms they are reputed to represent. Accordingly, we have not yet seen a serious exploration of the potential benefit of such transforms, carefully comparing the expected benefits with those delivered by specific implementations.
In this paper we aim at providing both:

A rationally designed shearlet transform implementation.

A comprehensive framework for quantifying performance of directional representations in general.
For (1), we developed an implementation of the fast digital shearlet transform (FDST) based on a digital shearlet theory which is a very natural digitization of the existing shearlet theory for continuous data. Other parabolicscaling transforms, for example, curvelets are inherently based on operations (rotation) which translate awkwardly into the digital realm. In contrast, when we consider shearlets, rotations are replaced by shearing, which has a natural digital realization, thus enabling a unified treatment for the continuum and digital realm similar to wavelets.
The framework in (2) has three benefits. First, it provides quantitative performance measures which can be used to tune the parameters of our implementation, which is publicly available at www.ShearLab.org. This allows us to specify ‘recommended choices’ for the parameters of our implementation. Second, the same ‘measure and tune’ approach may be useful to other implementers of directional transforms. Third, we show a way to improve the level of intellectual seriousness in applied mathematics which pretends to work in image processing. We believe that widespread adoption of this measure and tune framework can be very valuable, since many supposedly scientific presentations are now little more than vague, numbing ‘advertising’ or ‘marketing’ pitches. They could instead offer quantitative comparisons between algorithms, and thereby be far more informative. In fact the combination of quantitative evaluation with reproducible research [11] would be particularly effective at producing both intellectual seriousness and rapid progress.
1.1 Desiderata
We start by proposing the following desiderata for the fast digital shearlet transform FDST and its implementation:

Algebraic Exactness. The transform should be based on a shearlet theory for digital data on a pseudopolar grid, than merely being ‘somewhat close’ to the shearlet theory for continuous data.

Isometry of PseudoPolar Fourier Transform. We introduce oversampling and weights to obtain an isometric pseudopolar Fourier transform, which allows us to use the adjoint as inverse transform.

Tight Frame Property. The shearlet coefficients computed by the transform should be based on a tight frame decomposition, which ensures an isometric relation between the input image and the sequence of coefficients as well as allows us to use the adjoint as inverse transform. This property follows by combining [D1] and [D2], and allows the comparison with other transforms in contrast to those previous two tests.

TimeFrequencyLocalization. The spatial portrait of the analyzing elements should ‘look like’ shearlets in the sense that they are sufficiently smooth as well as timelocalized. Localization and smoothness in frequency domain is ensured by definition.

True Shear Invariance. Since the orientationrelated operator of shearlets is in fact the shear operator, we expect to see a shearing of the input image mirrored in a simple shift of the transform coefficients.

Speed. The transform should admit an algorithm of order flops, where is the number of digital points of the input image.

Geometric Exactness. The transform should preserve geometric properties parallel to those of the continuum theory, for example, edges should be mapped to edges in shearlet domain.

Robustness. The transform should be resilient against impacts such as (hard) thresholded and quantized coefficients.
1.2 Definition of the Shearlet Transform
The main idea for the construction of the shearlet transform with discrete parameters for functions in is the choice of a twoparameter dilation group, where one parameter ensures the multiscale property, whereas the second parameter provides a means to detect directions. The choice for a direction sensitive parameter is particularly important, since the most canonical choice, the rotation, would prohibit a unified treatment of the continuum and digital realm due to the fact that the integer grid is not invariant under rotation. Shearlets parameterize directions by slope rather than angles. And the shear matrix does preserve the structure of the integer grid, which is key to enabling an exact digitization of the continuum domain shearlets.
For each and , let denote the parabolic scaling matrix and denote the shear matrix of the form
respectively. To provide an equal treatment of the  and axis, the frequency plane is split into the four cones – (see Figure LABEL:fig:ShearletsCone), defined by
Let now be a wavelet with and supp , and let be a ‘bump’ function satisfying and supp . We define by
\hb@xt@.01(1.1) 
For cone , at scale , orientation , and spatial position , the associated shearlets are then defined by their Fourier transforms
\hb@xt@.01(1.2)  
where index scale, orientation, position, and cone. The shearlets for , , and are defined likewise by symmetry, as illustrated in Figure LABEL:fig:ShearletsCone, and we denote the resulting discrete shearlet system by
\hb@xt@.01(1.3) 
The definition shows that shearlets live on anisotropic regions of width and length at various orientations.
It should be mentioned that discrete shearlets – ‘discrete’ referring to the set of parameters and not to the domain – can also be defined with respect to the dilation matrix . However, in this case the odd scales have to be handled particularly carefully. The attentive reader will have also observed that recently introduced compactly supported shearlets [22, 26] do not require projecting the shearlets to the respective cones; however, despite other advantageous properties, they do not form a tight frame for . Finally, the generating window allows in fact more freedom than (LABEL:eq:psidef), but in this paper we restrict ourselves to this (customary) choice.
Setting , we have the following theorem from [14, Thm. 3] concerning the frame properties of the discrete shearlet system. For the definition of a tight (sometimes called Parseval) frame, we refer to [8].
Theorem 1.1 ([14])
The system (LABEL:eq:shearletsystem) is a tight frame for .
We remark that the low frequency part can be appropriately filled in to obtain a tight frame for .
The transform associated with this system is the discrete shearlet transform, which for a given function is defined to be the map
It is this transform, which we aim to exactly digitize.
1.3 Ingredients of the Fast Digital Shearlet Transform (FDST)
The shearlet transform for continuum domain data (see Figure LABEL:fig:ShearletsCone) implicitly induces a trapezoidal tiling of frequency space which is evidently not cartesian. By introducing a special set of coordinates on the continuum 2D frequency space, the discrete shearlet transform can be represented as a cascade of five operations:

Classical Fourier transformation.

Change of variables to pseudopolar coordinates.

Weighting by a radial ‘density compensation’ factor.

Decomposition into rectangular tiles.

Inverse Fourier transform of each tiles.
Surprisingly, this process admits a natural translation into the digital domain. The key observation is that the pseudopolar coordinates are naturally compatible with digital image processing (compare Figure LABEL:fig:PPgrid) and perfectly suited for a digitization of the discrete shearlet transform as a comparison with the frequency tiling generated by continuum domain shearlets in Figure LABEL:fig:ShearletsCone already visually evidences. Fortunately, in [1] a fast pseudopolar Fourier transform (PPFT) is already developed. This transform evaluates the Fourier transform of an image of size , say, on a pseudopolar grid of the form , where
Figure LABEL:fig:PPgrid shows an illustration of the case .
For an image , the pseudopolar Fourier transform of evaluated on the pseudopolar grid is then defined to be
where is an integer which for the PPFT is chosen to be for computational reasons.
The existence of PPFT suggests that we can easily and naturally get a faithful FDST using this algorithm. However, besides the delicateness of digitizing the continuum domain shearlets so that they form a tight frame on the pseudopolar grid, also the use of the PPFT is not at all straightforward. The PPFT as presented [1] is not an isometry. The main obstacle is the highly nonuniform arrangement of the points on the pseudopolar grid. This intuitively suggests to downweight points in regions of very high density by using weights which correspond roughly to the density compensation weights underlying the continuous change of variables. In fact, we will show that isometry is possible with sufficient radial oversampling of the pseudopolar grid; however, the weights will not be derivable from simple density compensation arguments.
Summarizing, the FDST of an image cascades the following steps:

Application of the PPFT with an oversampling factor in radial direction.

Weighting of the function values on the pseudopolar grid by ‘densitycompensationstyle’ weights.

Decomposing the pseudopolarindexed values by a scaled and sheared generating window into rectangular subbands followed by application of the 2D iFFT to each array.
This is an exact analogy of the discrete shearlet transform, in which the steps of Fourier transformation and pseudopolar coordinate change as well as the steps of decomposition into rectangular tiles and the inverse Fourier transform are collapsed into one step, respectively. With a careful choice of the weights and the windows, this transform is an isometry as we will show. Hence the inverse transform can be computed by merely taking the adjoint in each step.
1.4 Performance Measurement
The above sketch does not uniquely specify an implementation; there is freedom in choice of weights and windows. How can we decide if one choice is better than another one? It seems that currently researchers often use overall system performance on isolated tasks, such as denoising and compression of specific standard images like ‘Lena’, ‘Barbara’, etc. However, overall system performance for a system made up of a cascade of steps seems very opaque and at the same time very particular. It seems far better from an intellectual viewpoint to carefully decompose performance according to a more insightful array of tests, each one motivated by a particular wellunderstood property we are trying to obtain.
We have developed quantitative performance measures inspired by the desiderata we presented in Subsection LABEL:subsec:desiderata. Each performance measure produces a real value or a realvalued curve, thus providing a standardized framework for evaluation and, especially, comparison.
1.5 Relation with Previous Work
Since the introduction of directional representation systems by many pioneer researchers ([4, 5, 6, 7, 9, 14]), various numerical implementations of their directional representation systems have been proposed. The closest ones are the curvelet, contourlet, and previous shearlet algorithms, whose main features we now briefly survey.
Curvelets [3]. The discrete curvelet transform is implemented in the software package CurveLab, which comprises two different approaches. One is based on unequispaced FFTs, which are used to interpolate the function in the frequency domain on different tiles with respect to different orientations of curvelets. The other is based on frequency wrapping, which wraps each subband indexed by scale and angle into a fixed rectangle around the origin. Both approaches can be realized efficiently in flops with being the image size. The disadvantage of this approach is the lack of an associated continuum domain theory.
Contourlets [9]. The implementation of contourlets is based on a directional filter bank, which produces a directional frequency partitioning similar to the one generated by curvelets. The main advantage of this approach is that it allows a treestructured filter bank implementation, in which aliasing due to subsampling is allowed to exist. Consequently, one can achieve great efficiency in terms of redundancy and good spatial localization. A drawback of this approach is that various artifacts are introduced and that an associated continuum domain theory is missing.
Shearlets [12, 28]. In [12], Easley et. al. implemented the shearlet transform by applying the Laplacian pyramid scheme and directional filtering successionally. One drawback is the deviation from the continuum domain theory. Another drawback is that the associated code was not made publicly available. In contrast to this implementation which is based on bandlimited subband tiling – similar to the implementation of curvelets in [3] – in [28], Lim provided an implementation of the shearlet transform based on compactly supported shearlet systems (see also [22]). These compactly supported shearlets are separable and provide excellent spatial localization. The drawback is that they do not form a tight frame, hence, the synthesis process needs to be performed by iterative methods. We further wish to mention two novel approaches [27] and [19] for which however no implementation is yet available nor was their focus on deriving an exact digitization of the continuum domain transform.
Summarizing, all the above implementations of directional representation systems have their own advantages and disadvantages, one of the most common shortcomings being the lack of providing the unified treatment of the continuum and digital world. Our effort will now be put to provide a natural digitization of the shearlet theory (bandlimited shearlets) fulfilling the unified treatment requirement as well as a software package ShearLab quantifying performances of directional representation systems.
1.6 Contribution of this Paper
The contributions of this paper are twofold. Firstly, we introduce a fast digital shearlet transform (FDST) which is rationally designed based on a natural digitization of shearlet theory. Secondly, we provide a variety of quantitative performance measures for directional representations, which allow tuning and comparison of implementations. Our digital shearlet implementation was tuned utilizing this framework, so we can provide the user community with an optimized representation.
All presented algorithms and tests are provided at www.ShearLab.org in the spirit of reproducible research [11].
1.7 Contents
Section LABEL:sec:STfinitedata introduces the fast digital shearlet transform FDST and proves isometry. In Section LABEL:sec:ISTfinitedata, we then discuss two variants of an inverse digital shearlet transform, namely, a direct and an iterative approach. In Section LABEL:sec:props, we prove several mathematical properties of the FDST such as decay properties of digital shearlet coefficients. The following section, Section LABEL:sec:implementation, is concerned with details of the associated ShearLab implementation at www.ShearLab.org. The FDST is then analyzed in Section LABEL:sec:numerics according to the quantitative measures introduced in Section LABEL:sec:qualitymeasures.
2 FDST for Finite Data
We start by discussing the three steps in the FDST as described in Subsection LABEL:subsec:DST, which we for the convenience of the reader briefly repeat:

Application of the PPFT with an oversampling factor in radial direction.

Weighting of the function values on the pseudopolar grid by ‘densitycompensationstyle’ weights.

Decomposing the pseudopolarindexed values by a scaled and sheared generating window into rectangular subbands followed by application of the 2D iFFT to each array.
We will also show that careful selection of the oversampling factor, of the weights, and of the windows yields an isometric transform, which enables us to compute the inverse shearlet transform by its adjoint (see Section LABEL:sec:ISTfinitedata).
2.1 Weighted PseudoPolar Fourier Transform
Given an image , it is well known that the Fourier transform of evaluated on a rectangular grid is an isometry:
\hb@xt@.01(2.1) 
This is the Plancherel formula for a function defined on a finite group [21].
We now intend to obtain a similar formula for the Fourier transform of evaluated on the pseudopolar grid. For this, we first extend the definition of the pseudopolar grid slightly by introducing an oversampling parameter in radial direction. This new grid, which we will denote in the sequel by , is defined by
where
\hb@xt@.01(2.2)  
\hb@xt@.01(2.3) 
This grid is illustrated in Figure LABEL:fig:PPgridR.
Notice that the ‘original’ pseudopolar grid (see Figure LABEL:fig:PPgrid) as introduced in Subsection LABEL:subsec:DST is a special case of this definition when choosing . Also observe that the center
appears times in as well as , and the points on the seam lines
appear in both and . Later, we will also utilize a further partitioning of the sets and as
where
Now our goal is to choose weights so that, for any image ,
\hb@xt@.01(2.4) 
where here we modify the definition of the Fourier transform according to [1] and define it by
\hb@xt@.01(2.5) 
where . Also notice that the factor appearing in (LABEL:eq:Plancherel) will now be hidden in the weights .
We start by computing the right hand side of (LABEL:eq:ppPlancherel):
Choosing for all and for all , we can conclude that (LABEL:eq:ppPlancherel) holds if and only if
\hb@xt@.01(2.6) 
This is equivalent to the two conditions
\hb@xt@.01(2.7)  
\hb@xt@.01(2.8) 
for all . In view of the symmetry of the pseudopolar grid, it is natural to impose the following symmetry conditions on the weights:

, ,

, ,

, ,

, .
In this case, (LABEL:eq:isometry02) automatically holds. By the sum formula for trigonometric functions, (LABEL:eq:isometry01) is then equivalent to
for all . Again, by the symmetry of the weights, this is equivalent to
\hb@xt@.01(2.9) 
for all . This is a linear system of equations with unknows and equations, wherefore, in general, we need the oversampling factor to be at least to enforce solvability.
For symmetry reasons, we can now restrict our attention to one quarter of a cone, say . Using (LABEL:eq:OmegaR1), (LABEL:eq:OmegaR2), and (LABEL:eq:isometry1), we then obtain the following equivalent condition to (LABEL:eq:isometry0):
\hb@xt@.01(2.10)  
for all . Concluding, we have the following result.
Theorem 2.1
Let be even, let be the pseudopolar grid defined in (LABEL:eq:OmegaR1) and (LABEL:eq:OmegaR2), and let be a weight function satisfying the symmetry conditions [S1] – [S4]. Then
holds if and only if the weights satisfy condition (LABEL:eq:weightcondition). Moreover, in general, needs to be at least for such weights to exist.
2.2 Weight Functions
To avoid high complexity in the computation of the weights satisfying Theorem LABEL:theo:weightedisometry, we relax the requirement for exact isometric weighting. Instead of representing the weights as the solution of a large system of equations, they will be represented in terms of an undercomplete basis for functions on the pseudopolar grid. More precisely, we first design basis functions such that for all . We then define the weight function to be , with being nonnegative constants. These coefficients are determined by solving (LABEL:eq:weightcondition) with respect to this weight function using the least square method. We compute the coefficients in this expansion once for a given problem size; then hardwire them in the algorithm.
2.2.1 Recommended Choices of Weights
In what follows, we present several designs of basis functions, each one providing nearly isometric weighting. Notice that, slightly abusing notation, we will use and interchangeably. The weighting for the three choice we recommend is displayed in Figure LABEL:fig:001.
Choice 1:
Our first choice are the following seven functions :
Center: ,
Boundary: ,
Seam lines: ,
Interior: .
Choice 2:
This is a simplified version of ‘Choice 1’ using the 5 functions:
Center: ,
Boundary: ,
Seam lines: ,
Interior: .
Choice 3:
Finally, we suggest the following functions on the pseudopolar grid:
Center: ,
Radial Lines: .
2.2.2 Comparison of Weight Functions
The patterns of the weights are seemingly similar in view of Figure LABEL:fig:001. However, their performances can be quite different depending on the chosen performance measure. We mention that the measures chosen below are also part of our framework of performance measures for parabolic scaling algorithms discussed in Section LABEL:sec:qualitymeasures.
Letting , we generate a sequence of 5 random images , , of size with standard normally distributed entries. We use the following measure to compare the performance of different weights:
where denotes the PPFT from (LABEL:def:ppft) and – by abusing notation – denotes the ‘weighting operator’ , for a image and a weight function . Table LABEL:tab:1 displays the performance of the weights from Choice 1 – 3 with respect to this measure.
32  64  128  256  512  

Choice 1  4.3E3  2.6E3  2.2E3  1.4E3  9.3E4 
Choice 2  4.2E3  4.0E3  1.8E3  1.5E3  8.8E4 
Choice 3  9.8E3  6.2E3  3.4E3  2.1E3  N/A 
Next, we choose the real image ‘Barbara’, which we denote by , and the measure to compare the performance of the different choices of weights.
32  64  128  256  512  

Choice 1  2.4E3  1.6E3  1.1E3  5.2E4  2.2E4 
Choice 2  2.8E3  1.2E3  8.3E4  3.9E4  1.5E4 
Choice 3  5.6E3  2.8E3  2.2E3  9.1E4  N/A 
From Tables LABEL:tab:1 and LABEL:tab:2, it can be seen that the operator seems to converge to an identity operator as . The data also indicates that it is justifiable to define weights which are linearly increasing along the radial direction.
When using iterative methods, for instance, the conjugate gradient method, for computing the inverse, a weight plays the role of a preconditioner. Therefore its performance as such can be effectively measured by the condition number of the operator , i.e., . This measure is displayed in Table LABEL:tab:3 for our selected three choices of weights. Notice that, for each choice of a weight function with the exception of Choice 3, the condition numbers of are always smaller than .
32  64  128  256  512  

Choice 1  1.328  1.483  1.621  1.726  1.834 
Choice 2  1.379  1.503  1.621  1.731  1.833 
Choice 3  1.760  1.887  2.001  2.104  N/A 
2.3 Windowing
According to our discussion of the main steps of the FDST in Section LABEL:sec:STfinitedata, after performing a weighted pseudopolar Fourier transform, the data has to be windowed using scaled and sheared versions of a generating window function.
In this section, we now define such a set of window functions, which we will coin digital shearlets, and prove that these – similar to the continuum domain – form a tight frame for functions . We remark that this construction is a digitization of the continuum domain discrete shearlets introduced in [14] (compare also (LABEL:eq:shearlets) and (LABEL:eq:shearletsystem)). However, it is far from obvious that again a tight frame is derived, since we here consider a finite domain.
For the convenience of the reader, we first briefly recall the notion of a tight frame in this particular situation. Let be two functions defined on . Then, the inner product is defined to be . A sequence with being an indexing set is a tight frame for functions , if
It then follows from basic frame theory (see [8]) that this allows recovery of a function from its coefficients by computing
Despite the danger of repeating ourselves, let us mention that our fundamental goal is to introduce digital shearlets as the exact digitization of continuum domain shearlets. We now describe the construction step by step, which will give evidence to the fact that we achieved this goal. There will be one step though – when defining the modulation – where we have to slightly deviate from an exact digitization, and we will explain the reasons for this.
We start by defining the scaling function and the generating digital shearlet. For this, let , which will soon be shown to be the lowest possible scale. Let be the Fourier transform of the Meyer scaling function such that
\hb@xt@.01(2.11) 
and let be a ‘bump’ function satisfying
Then we define the scaling function for the digital shearlet system to be
We will later restrict this function to the pseudopolar grid.
Let next be the Fourier transform of the Meyer wavelet function with
\hb@xt@.01(2.12) 
as well as
\hb@xt@.01(2.13) 
We further choose to be a ‘bump’ function satisfying
\hb@xt@.01(2.14) 
and also
Notice that this implies
\hb@xt@.01(2.15) 
which will become important for the analysis of frame properties. For the choice of , , , and in our implementation, we refer to Section LABEL:sec:implementation. Then the generating shearlet is defined as
\hb@xt@.01(2.16) 
We will now define digital shearlets on and extend the definition to the other cones by symmetry. At this time, we assume and are both positive, even integers and for some integer .
For this, we first analyze the exact digitization of the coefficients of the discrete shearlet system from Subsection LABEL:subsec:defST for a function by using the shearlets defined in (LABEL:eq:digitalshearlet). This will lead to the appropriate range of scales and to the support of a scaled and sheared version of the shearlet .
When restricting to the cone , the exact digitization of the coefficients of the discrete shearlet system is
\hb@xt@.01(2.17) 
where , , and are to be determined. The choice of leads to the coefficients