Fast Finite Shearlet Transform: a tutorial

# Fast Finite Shearlet Transform: a tutorial

Sören Häuser Fachbereich für Mathematik, Technische Universität Kaiserslautern, Paul-Ehrlich-Str. 31, 67663 Kaiserslautern, Germany, {haeuser,steidl}@mathematik.uni-kl.de    Gabriele Steidl 11footnotemark: 1
\setkomafont

disposition

## 1 Introduction

Directional multiscale representation of images to address curved singularities has received much attention in harmonic analysis in the last 25 years. In particular, shearlets [13] and curvelets [1] provide an optimally sparse approximation of carton-like images, that is

 ∥f−fN∥2L2≤CN−2(logN)3as N→∞,

where is the nonlinear shearlet approximation of a function from this class obtained by taking the largest shearlet coefficients in absolute value. Shearlets possess a uniform construction for both the continuous and the discrete setting. They further stand out since they stem from a square-integrable group representation [4] and have the corresponding useful mathematical properties. Moreover, similarly as wavelets are related to Besov spaces via atomic decompositions, shearlets correspond to certain function spaces, the so-called shearlet coorbit spaces [5].

Figure 1 illustrates the directional information contained in the shearlet coefficients. Shearlets have been applied to a wide field of image processing tasks, e.g., denoising [10, 6], inversion of the Radon transform [3, 8], inverse halftoning [12], deconvolution [26], geometric separation [7], inpainting [17] and many more. A detailed summary can be found in [9]. In [16] the authors show how the directional information encoded by the shearlet transform can be used in image segmentation. To this end, we introduce a simple discrete shearlet transform which translates the shearlets over the full grid at each scale and for each direction. Using the FFT this transform can be still realized in a fast way.

This tutorial explains the details behind the Matlab-implementation of the transform and shows how to apply the transform. The software is available for free under the GPL-license at

http://www.mathematik.uni-kl.de/imagepro/software/

In analogy with other transforms we named the software FFSTFast Finite Shearlet Transform. The package provides a fast implementation of the finite (discrete) shearlet transform.

For shearlets there are currently three toolboxes available. They are

Local Shearlet Toolbox11footnotemark: 1

developed by Easley, Labate and Lim. This was the first shearlet implementation, for details see [11].

ShearLab33footnotemark: 3

developed by Donoho, Kutyniok, Lim, Shahram, Zhuang and Reisenhofer. This package consists of three different implementations: One is implemented on pseudo-polar grids, one on Cartesian grids and the newest one using compactly supported shearlets, for details see [19, 21, 20, 18].

Fast Finite Shearlet Transform (FFST)55footnotemark: 5

developed by the authors. The first fully finite and translation invariant shearlet implementation, described in [16, 15] and in this tutorial.

Recall that the Fourier transform and the inverse transform are defined by

 Ff(ω) =^f(ω) := ∫R2f(t)e−2πi⟨ω,t⟩dt, F−1^f(ω) =f(t) = ∫R2^f(ω)e2πi⟨ω,t⟩dω.

This tutorial is organized as follows: In Section 2 we introduce the continuous shearlet transform and prove the properties of the involved shearlets. We follow in Section 3 the path via the continuous shearlet transform, its counterpart on cones and finally its discretization on the full grid to obtain the translation invariant discrete shearlet transform. This is different to other implementations as, e.g., in ShearLab, see [20]. Our discrete shearlet transform can be efficiently computed by the fast Fourier transform (FFT). The discrete shearlets constitute a Parseval frame of the finite Euclidean space such that the inversion of the shearlet transform can be simply done by applying the adjoint transform. The second part of the section covers the implementation and installation details and provides some performance measures.

## 2 Closer Look at the Continuous Shearlet Transform in R2

In this section we combine some mostly well-known results from different authors. To make this tutorial self-contained and to obtain a complete documentation we also include the proofs. The functions are taken from [25, 24]. The construction of the shearlets is based on ideas from [11] and [22]. The shearlet transform and the concept of shearlets on the cone were introduced in [13].

### 2.1 Some Functions and their Properties

To define usable shearlets we need functions with special properties. We begin with defining these functions and prove their necessary properties. The results will be used later.

We start by defining an auxiliary function as

 v(x):=⎧⎨⎩0for x<0,35x4−84x5+70x6−20x7for 0≤x≤1,1for x>1. (1)

This function was proposed by Y. Meyer in [25, 24], see Remark 3.2 for more information on the construction of . Other choices of are possible, in [20] the simpler function

 ˜v(x)=⎧⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪⎩0for x<0,2x2for 0≤x≤12,1−2(1−x)2for 12≤x≤1,1for x>1,

was chosen.

As we will see, the useful properties of for our purposes are its symmetry around and the values at and with increase in between. A plot of is shown in Figure ((a))(a).

Next we define the function with

 b(ω):=⎧⎪ ⎪ ⎪⎨⎪ ⎪ ⎪⎩sin(π2v(|ω|−1))for 1≤|ω|≤2,cos(π2v(12|ω|−1))for 2<|ω|≤4,0otherwise. (2)

Note that is symmetric, positive, real-valued and . We further have that . A plot of is shown in Figure ((b))(b).

Because of the symmetry we restrict ourselves in the following analysis to the case . Let , , thus, and . Observe that is increasing for and decreasing for . Obviously all these properties carry over to . These facts are illustrated in the following diagram

 ω2j2j+12j+22j+3bj0↗1↘0bj+10↗1↘0

where stands for the increasing and for the decreasing function.

For the overlap between the support of and is empty except for . Thus, for and we have that . In this interval is decreasing with and is increasing with . Their sum in this interval is

 b2j(ω)+b2j+1(ω)=cos2(π2v(2−j−1|ω|−1))+sin2(π2v(2−j−1|ω|−1))=1.

Hence, we can summarize

 (b2j+b2j+1)(ω)=⎧⎪ ⎪⎨⎪ ⎪⎩b2jfor ω<2j+1,1for 2j+1≤ω≤2j+2,b2j+1for ω>2j+2.

Consequently, we have the following lemma.

###### Lemma 2.1.

For defined as above, the relations

 ∞∑j=−1b2j(ω)=∞∑j=−1b2(2−jω)=1for |ω|≥1

and

 ∞∑j=−1b2j(ω)=⎧⎪ ⎪ ⎪⎨⎪ ⎪ ⎪⎩0for |ω|≤12,sin2(π2v(2ω−1))for 12<|ω|<1,1for |ω|≥1 (3)

hold true.

###### Proof.

In each interval only and , , are not equal to zero. Thus, it is sufficient to prove that in this interval. We get that

 (b2j+b2j+1)(ω) =b2(2−jω\makebox[0.0pt][c]$∈2−j[2j+1,2j+2]=[2,4]$)+b2(2−j−1ω\makebox[0.0pt][c]$∈2−j−1[2j+1,2j+2]=[1,2]$) =cos2(π2v(12⋅2−jω−1))+sin2(π2v(2−j−1ω−1)) =cos2(π2v(2−j−1ω−1))+sin2(π2v(2−j−1ω−1)) =1.

The second relation follows by straightforward computation. ∎

We define the function via its Fourier transform as

 ^ψ1(ω):=√b2(2ω)+b2(ω). (4)

Figure ((a))(a) shows the function . The following theorem states an important property of .

###### Theorem 2.2.

The above defined function has and fulfills

 ∑j≥0|^ψ1(2−2jω)|2=1for |ω|>1.
###### Proof.

The assumption on the support follows from the definition of . For the sum we have

 ∑j≥0|^ψ1(2−2jω)|2=∞∑j=0b2(2⋅2−2jω)+b2(2−2jω)=∞∑j=0b2(2−2j+1ω)+b2(2−2jω)

where (odd) and (even). Thus, by Lemma 2.1, we get

 ∑j≥0|^ψ1(2−2jω)|2=∞∑j=−1b2(2−jω)=1.\qed

By (3) we have that

 ∑j≥0|^ψ1(2−2jω)|2=⎧⎪ ⎪ ⎪⎨⎪ ⎪ ⎪⎩0for |ω|≤12,sin2(π2v(2ω−1))for 12<|ω|<1,1for |ω|≥1. (5)

Next we define a second function —again in Fourier domain—by

 ^ψ2(ω):={√v(1+ω)for ω≤0,√v(1−ω)for ω>0. (6)

The function is shown in Figure ((b))(b). Before stating a theorem about the properties of we need the following two auxiliary lemmas.

Recall that a function is point symmetric with respect to if and only if

 f(a+x)−b=−f(a−x)+bfor all x∈R.

With the substitution this is equivalent to

 f(x)+f(2a−x)=2bfor all x∈R.

Thus, for a function symmetric to we have that .

###### Lemma 2.3.

The function in (1) is symmetric with respect to , in particular, for all .

###### Proof.

The symmetry is obvious for and . It remains to prove the symmetry for . We will see in Remark 3.2 that . With the fundamental theorem of calculus this implies that

 v(x)=−140∫x0t3(t−1)3dt.

Since we know that . Next, consider . Substituting yields

 v(1−x)=140∫x1(1−t3)(−t)3dt=−140∫1xt3(t−1)3dt.

Finally, we obtain for the sum

 v(x)+v(1−x)=−140∫x0t3(t−1)3dt−140∫1xt3(t−1)3dt=−140∫10t3(t−1)3dt=1.\qed

Note that is axially symmetric to the -axis.

###### Lemma 2.4.

The function fulfills

 ^ψ22(ω−1)+^ψ22(ω)+^ψ22(ω+1)=1for |ω|≤1.
###### Proof.

We have

 ^ψ22(ω)={v(1+ω)for ω≤0,v(1−ω)for ω>0.

Consequently, we get for that

 ^ψ22(ω−1)+^ψ22(ω)+^ψ22(ω+1) =v(1+ω−1)+v(1−ω)+v(1−ω−1) =v(ω)+v(1−ω)+v(−ω)=0=1,

and for we obtain similarly

 ^ψ22(ω−1)+^ψ22(ω)+^ψ22(ω+1) =v(1+ω−1)+v(1−ω)+v(1−ω−1) =v(−|ω|)=0+v(1−|ω|)+v(|ω|)=1.\qed

It can be seen in the proof that the sum reduces in both cases to two (different) summands, in particular

 1=^ψ22(ω−1)+^ψ22(ω)+^ψ22(ω+1)={^ψ22(ω−1)+^ψ22(ω)for 0≤ω≤1,^ψ22(ω)+^ψ22(ω+1)for −1≤ω<0.

With these lemmas we can prove the next theorem.

###### Theorem 2.5.

The function defined in (6) fulfills

 2j∑k=−2j|^ψ2(k+2jω)|2=1for |ω|≤1, j≥0. (7)
###### Proof.

With the assertion in (7) becomes

 2j∑k=−2j|^ψ2(k+~ω)|2=1for |~ω|≤2j, j≥0.

For a fixed (but arbitrary) we need for since . Thus, for , only the summands for do not vanish. But for we have and . In this case the entire sum reduces to one summand such that

 2j∑k=−2j|^ψ2(k+ω⋆)|2=|^ψ2(−ω⋆+ω⋆)|2=|^ψ2(0)|2=1.

If and the only non-zero summands appear for . Thus, , yields

 2j∑k=−2j|^ψ2(k+ω⋆)|2=|^ψ2(−⌊ω⋆⌋+ω⋆)|2+|^ψ2(−⌊ω⋆⌋−1+ω⋆)|2=|^ψ2(r+)|2+|^ψ2(1−r+)|2

which is equal to by Lemma 2.4. Analogously we obtain for , that the remaining non-zero summands are those for . With we get

 2j∑k=−2j|^ψ2(k+ω⋆)|2=|^ψ2(⌈ω⋆⌉+ω⋆)|2+|^ψ2(⌈ω⋆⌉+1+ω⋆)|2=|^ψ2(r−)|2+|^ψ2(1+r−)|2.

By Lemma 2.4 and since , we finally conclude

 |^ψ2(r−)|2+|^ψ2(1+r−)|2=|^ψ2(|r−|)|2+|^ψ2(1−|r−|)|2=1.\qed

### 2.2 The Continuous Shearlet Transform

For the shearlet transform we use the dilation matrix and shear matrix For and they read

 Aa=(a00√a),a∈R+,Ss=(1s01),s∈R. (8)

The shearlets emerge by dilation, shear and translation of a function as before

 ψa,s,t(x)=a−34ψ(A−1aS−1s(x−t))=a−34ψ⎛⎝⎛⎝1a−sa01√a⎞⎠(x−t)⎞⎠. (9)

We assume that can be written as

 ^ψ(ω1,ω2)=^ψ1(ω1)^ψ2(ω2ω1). (10)

Consequently, we obtain for the Fourier transform

 ^ψa,s,t(ω) =a−34F⎛⎝ψ⎛⎝⎛⎝1a−sa01√a⎞⎠(⋅−t)⎞⎠⎞⎠(ω) =a−34e−2πi⟨ω,t⟩F⎛⎝ψ⎛⎝⎛⎝1a−sa01√a⎞⎠⋅⎞⎠⎞⎠(ω) =a−34e−2πi⟨ω,t⟩(a−32)−1^ψ((a0s√a√a)ω) =a34e−2πi⟨ω,t⟩^ψ(aω1,√a(sω1+ω2)) =a34e−2πi⟨ω,t⟩^ψ1(aω1)^ψ2(a−12(ω2ω1+s)).

The shearlet transform of is given as

 SHψ(f)(a,s,t) =⟨f,ψa,s,t⟩ =⟨^f,^ψa,s,t⟩ =∫R2^f(ω)¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯^ψa,s,t(ω)dω =a34∫R2^f(ω)^ψ1(aω1)^ψ2(a−12(ω2ω1+s))e2πi⟨ω,t⟩dω =a34F−1(^f(ω)^ψ1(aω1)^ψ2(a−12(ω2ω1+s)))(t).

The same formula is derived by interpreting the shearlet transform as a convolution with the function and using the convolution theorem.

The shearlet transform is invertible if the function fulfills the admissibility property

 ∫R2|^ψ(ω1,ω2)|2|ω1|2dω1dω2<∞.

Easy calculations show that any shearlet of the form (10), where and are continuous and and , , is admissible. Figure 4 shows a dilated and sheared shearlet in Fourier and time domain.

### 2.3 Shearlets on the Cone

Up to now we have said nothing about the support of our shearlet . We use band-limited shearlets, thus, we have compact support in Fourier domain. In the previous section we assumed that , where we now define and as in (4) and (6), respectively. With the results shown for for and for , i.e., , it is natural to define the area

 Ch:={(ω1,ω2)∈R2:|ω1|≥12,|ω2|<|ω1|}.

We will refer to this set as the horizontal cone (see Figure 5). Analogously we define the vertical cone as

 Cv:={(ω1,ω2)∈R2:|ω2|≥12,|ω2|>|ω1|}.

To cover the entire we define two more sets

 C× :={(ω1,ω2)∈R2:|ω1|≥12,|ω2|≥12,|ω1|=|ω2|}, C0 :={(ω1,ω2)∈R2:|ω1|<1,|ω2|<1},

where is the intersection (or the seam lines) of the two cones and is the “low-frequency” part. Altogether we have with an overlapping domain

 C□:=(−1,1)2∖(−12,12)2. (11)

Obviously the shearlet defined above is perfectly suited for the horizontal cone. For each set , , we define a characteristic function which is equal to 1 for and 0 for . We need these characteristic functions as cut-off functions at the seam lines. We set

 ^ψh(ω1,ω2):=^ψ(ω1,ω2)=^ψ1(ω1)^ψ2(ω2ω1)χCh. (12)

For the non-dilated and non-sheared the cut-off function has no effect since the support of is completely contained in . But after the dilation and shear we have

 supp^ψa,s,0⊆{(ω1,ω2):12a≤|ω1|≤4a,∣∣∣s+ω2ω1∣∣∣≤√a}.

The question arises for which and this set remains a subset of the horizontal cone. For we have that is in but not in . Thus, we can restrict ourselves to .

With fixed the second condition for reads

 −√a ≤s+ω1ω2 ≤ √a, −√a−s ≤ω1ω2 ≤ √a−s. (13)

Since the right condition becomes and for the left condition , hence, we can conclude

 −1+√a≤s≤1−√a.

For such it holds that , in particular the indicator function is not needed for these (with respective ). One might ask for which (depending on ) the indicator function cuts off only parts of the function, i.e., . We take again (13) but now we do not use a condition to guarantee that but rather ask for a condition that allows . Thus, the right bound should be larger than and the left bound should be smaller than . Consequently, we obtain

 −1−√a≤s≤1+√a.

Summing up, we have for that . For parts of are also in , which are cut off. For the whole shearlet is set to zero by the characteristic function. If we go back to the definition of we see that the vertical range is determined by . By definition is axially symmetric with respect to the -axis, in other words the “center” of is taken for the argument equal to zero, i.e., . It follows that for the center of is at the seam-lines. Thus, for approximately one half of the shearlet is cut off whereas the other half remains. For larger larger parts would be cut. Consequently, we restrict ourselves to .

The shearlet for the vertical cone is defined analogously with the roles of and interchanged, i.e.,

 ^ψv(ω1,ω2):=^ψ(ω2,ω1)=^ψ1(ω2)^ψ2(ω1ω2)χCv. (14)

All the results from above apply to this setting. For , i.e., , both definitions coincide and we define

 ^ψ×(ω1,ω2):=^ψ(ω1,ω2)χC×. (15)

The shearlets , (and ) are called shearlets on the cone. This concept was introduced in [13].

We have functions to cover three of the four parts of . The remaining part will be handled with a scaling function which is presented in the next section.

### 2.4 Scaling Function

For the center part (also low-frequency part) we define another set of functions. To this end, we need the following scaling function

 φ(ω):=⎧⎪ ⎪ ⎪⎨⎪ ⎪ ⎪⎩1for |ω|≤12,cos(π2v(2|ω|−1))for 12<|ω|<1,0otherwise.

The full scaling function is then defined as

 ^ϕ(ω1,ω2):= {φ(ω1)for |ω2|≤|ω1|,φ(ω2)for |ω1|<|ω2| (16) = ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩1for |ω1|≤12,|ω2|≤12,cos(π2v(2|ω1|−1))for 12<|ω1|<1,|ω2|≤|ω1|,cos(π2v(2|ω2|−1))for 12<|ω2|<1,|ω1|<|ω2|,0otherwise.

The decay of the scaling function (respectively ) is chosen to match with the increase of . For we have by (5) that

 |^ψ1(ω)|2+|φ(ω)|2=sin2(π2v(2|ω|−1))+cos