Extraction of digital wavefront sets using applied harmonic analysis and deep neural networks

TL;DR

This paper presents a novel algorithm that combines shearlet transforms and deep neural networks to accurately extract wavefront sets from images, significantly improving edge and orientation detection in imaging sciences.

Contribution

It introduces the first algorithmic approach integrating data-driven and model-based methods for wavefront set extraction using shearlet transforms and deep learning.

Findings

Outperforms existing algorithms in edge-orientation detection

Achieves superior ramp-orientation detection accuracy

Effectively combines shearlet analysis with neural networks

Abstract

Microlocal analysis provides deep insight into singularity structures and is often crucial for solving inverse problems, predominately, in imaging sciences. Of particular importance is the analysis of wavefront sets and the correct extraction of those. In this paper, we introduce the first algorithmic approach to extract the wavefront set of images, which combines data-based and model-based methods. Based on a celebrated property of the shearlet transform to unravel information on the wavefront set, we extract the wavefront set of an image by first applying a discrete shearlet transform and then feeding local patches of this transform to a deep convolutional neural network trained on labeled data. The resulting algorithm outperforms all competing algorithms in edge-orientation and ramp-orientation detection.

Tables

Table 1. Table 1: Ellipses/parallelograms data set performance metrics in percentage.

Method	Test accuracy	MF-score
Logistic regression	45.7	48.9
Decision trees	75.2	75.8
Linear SVM	46.5	50.3
K-nearest neighbors	72.7	73.2
Random forest	86.0	86.7
DeNSE	96.2	97.1

Table 2. Table 2: Edge detection performances of edge detection algorithms on the Ellipses/parallelograms data set. The MF-Score is in percentage.

Method	MF-score
Canny [11]	49.1
Sobel [50]	40.0
BEL [13]	63.3
Yi-Labate-Easley-Krim [54]	70.3
CoShREM [44]	90.6
DeNSE	97.5

Table 3. Table 3: BSDS500 (Berkeley) data set performance metrics in percentage.

Method	MF-score
gPb-owt-ucm [2]	73.7
gPb [2]	71.5
Mean Shift [12]	64.0
Normalized Cuts [49]	64.2
Felzenszwalb, Huttenlocher [19]	61.0
Canny	60.3
CoShREM [44]	75.7
DeepEdge [4]	75.3
DeNSE	95.4

Table 4. Table 4: Performance on the SBD data set. All values are in percentage.

Method	MF-score
OBDC	62.5
CASENet	71.8
CASENet-S	75.8
CASENet-C	80.4
CoShREM	69.7
SEAL	81.1
DeNSE	96.8

Table 5. Table 5: Error of wavefront set estimation by different inversion techniques.

Inversion technique	Mean square error
Tikhonov	443.0
Total variation	380.9
Filtered backprojection	504.3
Canonical relations	168.1

Equations

$$\hat{f}(\xi ):={\int }_{{\mathbb{R}}^d}f(x)e^{-2\pi i⟨x,\xi ⟩}dx\text{ for }\xi \in {\mathbb{R}}^d.$$

$$PW_{\Lambda}\coloneqq\bigl{\{}f\in L^{2}(\mathbb{R}^{d})\colon\operatorname{supp\,}\hat{f}\in[-\Lambda,\Lambda]^{d}\bigr{\}}.$$

$$\bigl{|}\widehat{\phi f}(\xi)\bigr{|}\leq C_{k}\bigl{(}1+|\xi|\bigr{)}^{-k}\quad\text{for all $\xi\in\mathbb{R}^{2}\setminus\{0\}$ such that $\xi/|\xi|\in V_{\lambda}$}$$

$$\mathrm{WF}(f):=\bigcup _{k\in \mathbb{N}}{\mathrm{WF}}_k(f).$$

$$\bigl{\{}x\in\mathbb{R}^{2}:(x,\lambda)\in\operatorname{WF}(f)\text{ for some $\lambda\in\mathbb{S}^{1}$}\bigr{\}}.$$

$$\operatorname{WF}(f)=\bigl{\{}(x,\lambda)\in\mathbb{R}^{2}\times\mathbb{S}^{1}\colon x\in\partial D\text{ and }\lambda=n_{x}\text{ where $n_{x}$ is a normal on $\partial D$ at $x$}\bigr{\}}.$$

$$A_a:=\left(\begin{array}{cc}a& 0\\ 0& \sqrt{a}\end{array}\right),{\tilde{A}}_a:=\left(\begin{array}{cc}\sqrt{a}& 0\\ 0& a\end{array}\right),\text{and}{S}_s:=\left(\begin{array}{cc}1& s\\ 0& 1\end{array}\right)\text{for a>0 and s∈R.}$$

$${\psi }_{a,s,t,1}(x):=a^{-\frac{3}{4}}\psi \left(A_a^{-1}{S}_s^{-1}(x-t)\right)\text{and}{\psi }_{a,s,t,-1}(x):=a^{-\frac{3}{4}}\tilde{\psi }\left({\tilde{A}}_a^{-1}{S}_s^{-T}(x-t)\right),$$

$$\mathcal{SH}_{\psi}:L^{2}(\mathbb{R}^{2})\to L^{\infty}\bigl{(}\mathbb{R}^{+}\times\mathbb{R}\times\mathbb{R}^{2}\times\{-1,1\}\bigr{)}\quad\text{where}\quad\mathcal{SH}_{\psi}(f)(a,s,t,\iota)\coloneqq\langle f,\psi_{a,s,t,\iota}\rangle.$$

$$\int_{\mathbb{R}^{2}}\frac{\bigl{|}\widehat{\psi}(\xi_{1},\xi_{2})\bigr{|}^{2}}{|\xi_{1}|^{2m}}d\xi<\infty.$$

$$\bigl{|}\psi(x)\bigr{|}=\mathcal{O}\bigl{(}(1+|x|)^{-p}\bigr{)}\quad\text{for $|x|\to\infty$.}$$

$$\bigl{|}\mathcal{SH}_{\psi}(f)(a,s,x,\iota)\bigr{|}=\mathcal{O}\Bigl{(}a^{\frac{p}{2}-\frac{3}{4}}+a^{\frac{m}{4}}+a^{\frac{3k}{4}-\frac{3}{4}}+a^{\frac{3l}{4}}\Bigr{)}\quad\text{as $a\to 0$}$$

$$\lambda (s,\iota ):=\{\begin{array}{ll}\left(\frac{1}{\sqrt{s^2+1}},\frac{s}{\sqrt{s^2+1}}\right)& \text{if ι=1,}\\ \left(\frac{s}{\sqrt{s^2+1}},\frac{1}{\sqrt{s^2+1}}\right)& \text{if ι=−1.}\end{array}$$

$$\bigl{|}\mathcal{SH}_{\psi}(f)(a,s,x,\iota)\bigr{|}=\mathcal{O}(a^{n})\quad\text{as $a\to 0$},$$

$$\mathcal{PW}_{\Lambda}\coloneqq\left\{f\in L^{2}\left(\mathbb{R}^{d}\right)\colon\operatorname{supp}\big{(}\widehat{f}\,\big{)}\subset[-\Lambda,\Lambda]^{d}\right\}.$$

$${\mathrm{sinc}}_d(x):=\prod _{i=1}^d\frac{\sin (\pi x_i)}{\pi x_i},\text{where }x=(x_1,\dots ,x_d)\in {\mathbb{R}}^d.$$

$$f\in {\mathcal{PW}}_{\Lambda }⟺f(x)=\sum _{n\in {\mathbb{Z}}^d}f\left(\frac{n}{\Lambda }\right){\mathrm{sinc}}_d(\Lambda \cdot x-n)\text{ for all x∈Rd.}$$

$$f(x):=\sum _{k\in {\mathbb{Z}}^2}{y}_k{\mathrm{sinc}}_2(\Lambda \cdot x-k),\text{ for }x\in {\mathbb{R}}^d,$$

$${\mathrm{DWF}}_{\Lambda }:{\mathcal{PW}}_{\Lambda }\to P\left({\mathbb{R}}^2\times {\mathbb{S}}^1\right)\text{such that DWF(PΛf)=WF(f) for all f∈L2(R2).}$$

$$W_{x}\coloneqq\bigl{\{}\lambda\in\mathbb{S}^{1}:(x,\lambda)\in W\bigr{\}}.$$

$${\mathrm{DWF}}_j:{\mathcal{PW}}_j\to \mathcal{P}\left({\mathbb{R}}^2\times {\mathbb{S}}^1\right)\text{for j∈N}$$

$$d_H\left(\overline{{\mathrm{DWF}}_j(P_j(f){)}_x},\overline{\mathrm{WF}(f{)}_x}\right)\to 0.$$

$$\mathrm{WF}(f_1+f_2)=\mathrm{WF}(f_1)\text{for every f1∈L2(R2) and f2∈C∞(R2)∩L2(R2).}$$

$$q_0$$

$$q_n$$

$$\sum _{n\le N}\Vert q_n-q_{n-1}{\Vert }_{C^k}\lesssim \sum _{n\le N}{\Vert (1+|\cdot {|}^2{)}^{\frac{k}{2}}\left({\hat{q}}_n-{\hat{q}}_{n-1}\right)\Vert }_{L^1}\le {\Vert (1+|\cdot {|}^2{)}^{\frac{k}{2}}\hat{g}\Vert }_{L^1}<\infty .$$

$$(P_{n}g-P_{n-1}g)\perp (P_{m}g-P_{m-1}g)\text{for all n=m.}$$

$$\sum _{n\le N}\Vert P_{n}g-P_{n-1}g{\Vert }_2^2\le \Vert g{\Vert }_2^2<\infty .$$

$$\displaystyle d_{H}\bigl{(}\overline{\operatorname{DWF}_{j}(P_{j}q_{j})_{x}},\overline{\operatorname{WF}(q)_{x}}\bigr{)}<1\quad\text{for all $j\geq j_{0}$.}$$

$$q-P_{j_0}q=\sum _{j>j_0}(P_{j}g-P_{j-1}g)=g-P_{j_0}g.$$

$${\Vert f-P_{j_1}f\Vert }_2\le \frac{ϵ}{2}\text{and}{\Vert g-P_{j_1}g\Vert }_2\le \frac{ϵ}{2}.$$

Full text

Keywords: Wavefront set, deep learning, convolutional neural networks, shearlets.

Mathematics Subject Classification: 35A18, 65T60, 68T10.

1 Introduction

Many scientific and industrial real-world applications require a precise understanding of how a model parameter, represented by a function, is transformed under some process that is described by an operator. Such analysis quickly becomes very challenging, and one attempt at simplifying it is to treat the singular (non-smooth) and smooth parts of the function separately. In fact, a significant portion of the useful information is often contained in the singular part. For images, this singular part corresponds to edges, ridges, or ramps in the image.

Microlocal analysis is a powerful mathematical theory that aims to precisely describe how the singular part of a function, or more generally a distribution, is transformed when acted upon by an operator. Since its introduction in the early 1970s by Sato [47] and Hörmander [27], it has proven itself useful in both pure and applied mathematical research. The crucial underlying observation in microlocal analysis is that the information about the location of the singularities (singular support) needs to be complemented with specifying directional information of the singularities. This extra (“microlocal”) information is key in elucidating the propagation of singularities by a certain class of operators. This class includes Fourier integral operators such as most differential and pseudo-differential operators as well as many integral operators arising in integral geometry. Such operators are frequently encountered in analysis, scientific computing, and physical sciences [27, 9].

Microlocal analysis is also particularly useful in inverse problems, where the goal is to reliably recover a hidden model parameter (signal) from a noisy transformed version. The goal here is to recover the wavefront set of the signal (image) given the noisy realization of a transformed version of the signal. Such applications frequently arise when using imaging/sensing technologies where the transform is a pseudo-differential or Fourier integral operator [31].

A prime example exhibiting the utility of microlocal analysis in inverse problems appears in the context of imaging applications. Here, the image is represented by a function describing the interior structure of the object under study. Recovering the image is, however often not possible, either because the transformation relating the image to data is not invertible or because data is incomplete. On the other hand, full knowledge of the image is not always necessary. As an example, if one is looking for a tumor, then its location and shape are often sufficient for decision making whereas the exact values of the tumor density may be ignored or identified with another more specialized technique. The location and shape of the tumor can be readily determined from the singular part of the image.

For the above reason, there exist a plethora of applications of microlocal analysis to tomographic imaging. In these applications, the transformation relating the image to data is the ray transform which can be interpreted as a Fourier integral operator. Thereby, it is possible to explicitly describe the relation between the wavefront set of the function (image) and its transformed version (tomographic data). Such relationships are referred to as (microlocal) canonical relations.

The canonical relation can also be used to identify which singularities can be recovered from data without explicitly computing the inverse ray transform. This was done in [41] for the case when the 2D/3D ray transform is restricted to parallel lines and in [43] for an analysis in the region-of-interest limited angle setting. A related principle was derived in [17] for the 3D ray transform restricted to lines given by helical acquisition that is common in medical imaging. Similar principles hold for transforms integrating along other types of curves, for example, ellipses with foci on the x-axis and geodesics [53]. Another observation is that recovering a signal from its ray transform is less ill-posed if one knows the wavefront set a priori. This was demonstrated in [15] where the severely ill-posed reconstruction problem in limited angle tomography becomes mildly ill-posed if the wavefront set of the solution is provided as prior information, see [43] for an application of this principle to cryo-electron tomography.

Additionally, the concept of the wavefront set is of significant interest in mathematical image processing outside of inverse problems. Directed singularities, in particular edges, ridges, or ramps play an essential role in image processing to the extent that edge detection is one of the principal problems of this field. This is due to the fact, that edges in images present boundaries of objects and carry most of the information of the associated physical scene [52, 38, 6, 7]. Additionally, it has been argued in [38] that the human visual cortex performs multiple operations of image processing, the first of which is rough sketching involving edge detection.

1.1 Wavefront set extraction

A necessary first step before applying techniques from microlocal analysis in real-world applications in the way described above is to extract the wavefront set of functions. By the definition of the wavefront set, this amounts to estimating the asymptotic behavior of the localized Fourier transform of the function in question.

Extraction from finitely many samples.

In applications, where only finitely many point samples of the underlying function are known, estimating the asymptotic behavior of the Fourier transform is usually not possible. Indeed, we show in Section 3 that every method that seeks to approximately extract the wavefront set of functions from a function class $\mathcal{C}\subset L^{2}$ will fail on a dense subset of $\mathcal{C}$ as soon as $\mathcal{C}$ contains at least all $k$-times differentiable functions for an arbitrary $k\in\mathbb{N}$.

The problem of extracting the wavefront set of functions from a function class $\mathcal{C}$ certainly becomes easier the smaller $\mathcal{C}$ gets. Indeed, in applications, where, for example, $\mathcal{C}$ models a class of images, it is reasonable to assume that $\mathcal{C}$ is a very small subset of $L^{2}$ and does not contain all $k$-times differentiable functions for any $k\in\mathbb{N}$. Hence, in this situation, wavefront set extraction from point samples explicitly designed for the class of images could be feasible. In fact, if the set $\mathcal{C}$ is so small that every function in $\mathcal{C}$ is uniquely determined through its samples on the grid, then wavefront set detection could, in principle, be performed through a large data base, which could be learned from $\mathcal{C}$.

From the above, it is clear that a useful wavefront set extractor needs to closely adapted to the underlying function class $\mathcal{C}$. This strong adaptation to $\mathcal{C}$ will henceforth be referred to as our guiding principle.

Classification with applied harmonic analysis.

Certain transforms from applied harmonic analysis, like the curvelet and shearlet transform, offer an alternative possibility to identify the wavefront set. In particular, the connection between the behavior of these transforms, and the wavefront set has been analyzed in [10, 33], with a particular application for edge detection in [54]. We will recall these results in Section 2. These approaches characterize the wavefront set through the rate of decay of the respective transforms. In this way, one can transform the problem of extracting the wavefront set of a function to a classification problem on the decay of another function. While this point of view certainly makes the wavefront set more accessible, especially since it does not depend on an unspecified localization procedure, it is still subject to the fundamental restrictions of the previous subsection, i.e., it cannot produce a successful wavefront set extractor acting on sampled functions.

Data driven wavefront set extraction.

Bearing in mind the aforementioned guiding principle, a successful wavefront set extractor needs to be tailored to the function class of interest. The relevant function classes in applications are, however, difficult to characterize analytically. An alternative is, therefore, to adopt a data-driven model where the function class of interest is given empirically through examples.

Based on the above, we propose the following algorithm coined DeNSE: First, we assemble a supervised training set consisting of images with their associated wavefront set, or a suitable surrogate. Second, we train a classifier—in our case a deep neural network—to predict the wavefront set from the shearlet coefficients of the training data. Finally, we apply the resulting classifier on unseen data.

In this way, the algorithm combines two crucial elements: On the one hand, as mentioned previously, the interaction of the shearlet transform with the wavefront set is theoretically well understood and presents the microlocal information of a function in a more accessible way. On the other hand, the trained classifier allows a strong adaptation to the underlying function class, thereby complying with our guiding principle.

In Section 6, we present the construction of the algorithm alongside the training data that is used. In this context, we analyze the detection of edges and orientations in images as well as higher-order wavefront set detection in sinogram data. We shall see below, in Section 7, that this method outperforms all conventional edge-orientation estimators as well as alternative data-driven methods including the current state-of-the-art. Moreover, we are unaware of any wavefront set extractor in the literature that goes beyond edge or ramp detection, so that the following analysis can be seen as the first advance in this direction.

1.2 Expected impact

The DeNSE method has many applications, some of which are outlined below:

• •

Solving inverse problems: Many inverse problems involve a pseudo-differential or a Fourier integral operator. Moreover, the associated microlocal canonical relations can be explicitly derived for such operators. This gives a relation between the wavefront sets of measurement data and the unknown signal one seeks to recover. In particular, one may use DeNSE to extract the wavefront set of data and then use the relevant microlocal canonical relation to compute the wavefront set of the unknown signal without solving the inverse problem. Section 7.5 demonstrates this approach in the context of tomography.

• •

Regularization of inverse problems: A priori knowledge of the relation between the wavefront sets for data and signal (the solution) can be used to regularize an ill-posed inverse problem. An example of this, which used shearlets for the identification of the wavefront set, was demonstrated in [8] in the context of limited angle tomography. This approach is, however, specifically tailored for inverse problems involving the Radon transform. In contrast, our DeNSE based approach is more generic and applies readily to a wide range of inverse problems, see Section 5 for further details.

• •

Edge detection: Detecting edges, ridges, or points of higher-order non-smoothness is a sub-task of wavefront set extraction. As we will observe below, DeNSE outperforms state-of-the-art competing edge-orientation detectors on a wide range of test cases. Moreover, except for ridge and ramp detection, the detection of points of higher-order non-smoothness has—to the best of our knowledge—not been pursued in the literature, but is possible with DeNSE without adaptations.

1.3 Basic concepts and notation

Below, we collect the notation used throughout this manuscript. This notation is fairly standard in the literature, and hence this subsection can be skipped and only consulted if a notation is unclear.

$\mathbb{R}$, $\mathbb{N}$, and $\mathbb{Z}$ denote the set of real numbers, natural numbers, and integers, respectively. Next, given a point $x\in\mathbb{R}^{n}$ and $r>0$, we use $B_{r}(x)$ to denote the ball of radius $r$ in $\mathbb{R}^{n}$ with center at $x$. Likewise, $\mathbb{S}^{n-1}$ denotes the unit sphere in $\mathbb{R}^{n}$. Furthermore, the boundary of a domain $\Omega\subset\mathbb{R}^{d}$ for $d\in\mathbb{N}$ is denoted by $\partial\Omega$.

Let $\Omega\subset\mathbb{R}^{d}$ be a fixed domain for some $d\in\mathbb{N}$. Then, $L^{2}(\Omega)$ is the space of Lebesgue square-integrable functions on $\Omega$, $C^{n}(\Omega)$ is the space of $n$-times continuously differentiable functions defined on $\Omega$, and $H^{n}(\Omega)$ is the space of $n$-times weakly differentiable functions whose weak derivatives are in $L^{2}(\Omega)$. The support of a measurable function $f\colon\mathbb{R}^{d}\to\mathbb{R}$ is denoted by $\operatorname{supp\,}f$. Furthermore, the Fourier transform of a function $f\in L^{1}(\mathbb{R}^{d})$ is defined as

[TABLE]

The Fourier transform operator $f\mapsto\hat{f}$ can be extended to an isometry on $L^{2}(\mathbb{R}^{n})$ by Plancherel’s identity. The Paley-Wiener spaces are formed by functions the Fourier transform of which is compactly supported. More precisely, for $\Lambda\in\mathbb{R}_{+}$, we define

[TABLE]

Finally, we use the Landau symbol $\mathcal{O}$ to describe asymptotic behavior, i.e., for functions $f,g\colon\mathbb{R}\to\mathbb{R}$, we write $f(x)=\mathcal{O}(g(x))$ as $x\to a$ whenever there exists a constant $c>0$ such that $|f(x)|\leq c|g(x)|$ for all $x$ in a neighborhood of $a$. Similarly, we write $f(x)=\mathcal{O}(g(x))$ as $x\to\infty$ whenever there exists a constant $c>0$ such that $|f(x)|\leq c|g(x)|$ for $|x|$ sufficiently large.

2 Directional multiscale systems and the wavefront set

We start by formally introducing the notion of a wavefront set followed by the definition of the directional multiscale system of shearlets. We then show how shearlets can indeed be used to resolve the wavefront set of a distribution. Similar results also hold for other multiscale systems, like the curvelet transform [10] and the more general continuous parabolic molecules [22].

2.1 The wavefront set

The wavefront set can be defined for any distribution on a manifold, but since we only deal with $L^{2}(\mathbb{R}^{2})$ functions, we restrict the definition to this setting.

Definition 2.1** ([28, Section 8.1]).**

Let $f\in L^{2}(\mathbb{R}^{2})$ and $k\in\mathbb{N}$. A point $(x,\lambda)\in\mathbb{R}^{2}\times\mathbb{S}^{1}$ is a $k$-regular directed point of $f$ if there exist open neighborhoods $U_{x}$ and $V_{\lambda}$ of $x$ and $\lambda$, respectively and a smooth function $\phi\in C^{\infty}(\mathbb{R}^{2})$ with $\operatorname{supp\,}\phi\subset U_{x}$ and $\phi(x)=1$ such that

[TABLE]

holds for some $C_{k}>0$. The $k$-wavefront set $\operatorname{WF}_{k}(f)$ is the complement of the set of all $k$-regular directed points and the wavefront set $\operatorname{WF}(f)$ is defined as

[TABLE]

The definition of the wavefront set is based on the well-known characterization of smoothness of a function in terms of the decay of its Fourier transform. More precisely, a function is smooth at a point if its Fourier transform decays faster than any polynomial in any direction. As an example, the singular support of $f$, i.e., the smallest closed set $U$ such that $f_{|U^{c}}\in C^{\infty}(U^{c})$ can be characterized in terms of the wavefront set as

[TABLE]

The wavefront set is a refined notion of the singular support, since it not only indicates at which points a function is not smooth, but also contains the associated directions causing the non-smoothness. A common example considers $f$ with a jump singularity across the smooth boundary of a fixed domain $D\subset\mathbb{R}^{2}$. Then, one can show [51, Chapter VI, Exercise 1.1] that

[TABLE]

The wavefront set is a very powerful tool in mathematical analysis, but it is difficult to compute in practice. This is mainly due to the asymptotic criteria involved in its definition, which means computing the wavefront set requires computing the “full” Fourier transform at every point. Continuous transforms associated to certain directional multiscale systems offer a convenient remedy. As an example, the shearlet transform automatically performs the necessary time-frequency-orientation localization described above, thereby “resolving” the wavefront set in a sense described in the following subsection.

2.2 Shearlets

The shearlet transform, which was introduced in [23], is based on applying translation, anisotropic dilation, and shearing to generator functions. To dilate and shear a function, we define the following three matrices:

[TABLE]

Next, given $(a,s,t)\in\mathbb{R}_{+}\times\mathbb{R}\times\mathbb{R}^{2}$, $\psi\in L^{2}(\mathbb{R}^{2})$, and $x\in\mathbb{R}^{2}$, we define

[TABLE]

where $\widetilde{\psi}(x_{1},x_{2})\coloneqq\psi(x_{2},x_{1})$ for all $x=(x_{1},x_{2})\in\mathbb{R}^{2}$. Following [21], we define the continuous shearlet transform as follows:

Definition 2.2** (Continuous shearlet transform).**

Let $\psi\in L^{2}(\mathbb{R}^{2})$. Then the family of functions $\psi_{a,s,t,\iota}\colon\mathbb{R}^{2}\to\mathbb{R}$ parametrized by $(a,s,t,\iota)\in\mathbb{R}^{+}\times\mathbb{R}\times\mathbb{R}^{2}\times\{-1,1\}$ that are defined in (2.2) is called a shearlet system. The corresponding (continuous) shearlet transform is defined by

[TABLE]

As we shall see next, if the generator function $\psi$ has directional vanishing moments, then the asymptotic behavior as $a\to 0$ of the continuous shearlet transform of an $L^{2}$-function $f$ characterizes its wavefront set. The precise statement in [21] reads as follows.

Theorem 2.3**.**

Let $f\in L^{2}(\mathbb{R}^{2})$ and assume $(x_{0},\lambda_{0})\in\mathbb{R}^{2}\times\mathbb{S}^{1}$ is a $k$-regular directed point of $f$ for some $k\in\mathbb{N}$. Next, consider a continuous shearlet system with generator function $\psi\in H^{l}(\mathbb{R}^{2})$, $l\in\mathbb{N}$, with Fourier transform $\widehat{\psi}\in L^{1}(\mathbb{R}^{2})$ where $\psi$ has $m\in\mathbb{N}$ vanishing moments in $x_{1}$-direction, i.e.,

[TABLE]

Finally, assume $\psi$ displays the following asymptotic behavior for $p\in\mathbb{N}$:

[TABLE]

Then, there exist a neighborhood $U_{0}\subset\mathbb{R}^{2}$ of $x_{0}$ and a neighborhood $S_{0}\subset\mathbb{S}^{1}$ of $\lambda_{0}$ such that

[TABLE]

for all $x\in U_{0}$ and all $s\in\mathbb{R}$ and $\iota\in\{-1,1\}$ such that $\lambda(s,\iota)\in S_{0}$, where

[TABLE]

Remark 2.4**.**

Under suitable assumptions on the shearlet generator $\psi$, the converse of Theorem 2.3 holds as well. More precisely, following [21], assume that $\psi$ is sufficiently regular for any $k\in\mathbb{N}$. Next, let $(x_{0},\lambda_{0})\in\mathbb{R}^{2}\times\mathbb{S}^{1}$ and assume there exist a neighborhood $U_{0}\subset\mathbb{R}^{2}$ of $x_{0}$ and a neighborhood $S_{0}\subset\mathbb{S}^{1}$ of $\lambda_{0}$ such that

[TABLE]

holds for sufficiently large $n\in\mathbb{N}$ uniformly for $x\in U_{0}$ and all $s\in\mathbb{R}$, $\iota\in\{-1,1\}$ such that $\lambda(s,\iota)\in S_{0}$. Then, $(x_{0},\lambda_{0})\not\in\operatorname{WF}_{k}(f)$.

Theorem 2.3 and Remark 2.4 demonstrate that the wavefront set is completely determined by the decay properties of the shearlet transform. This implies that in the continuous setting, one can compute the wavefront set of a function by first computing its continuous shearlet transform, then analyzing the pairs of point and direction where this shearlet transform exhibits rapid decay as $a\to 0$.

Theorem 2.3 and Remark 2.4 were first reported in [33] in a setting restricted to a specific shearlet generator (called the “classical shearlet”). Moreover, results similar to Theorem 2.3 and Remark 2.4 were obtained in [10] for the curvelet transform and in [18] for transforms stemming from general group representations. As shown in [22], all transforms that belong to the category of continuous parabolic molecules admit a similar characterization of the wavefront set. Finally, one can also use the shearlet transform to classify certain geometric properties of the singularities of a function that goes beyond differentiating between rapid and non-rapid decay of the shearlet transform, see e.g. [54, 24, 36].

3 Wavefront set of sampled functions

Here we analyze whether it is possible to construct an operator that maps a finitely sampled function $f$ to an estimate of its wavefront set. This typically arises in practical applications, e.g., images are only given as pixels representing point samples of a real-valued function.

It is natural to use Shannon’s sampling theorem to make the connection between a sampled function and its wavefront set more precise. The theorem is stated in Subsection 3.1 and – based on it – Subsection 3.2 introduces the notion of an approximate wavefront set extractor. Finally, in Subsection 3.3 we show that any approximate wavefront set extractor on a function class $\mathcal{C}\subset L^{2}(\mathbb{R}^{2})$ that predicts the wavefront set of a function $f\in\mathcal{C}$ from a finite number of sample values will fail on a dense subset of $\mathcal{C}$ if $\mathcal{C}$ contains at least all $C^{k}(\mathbb{R}^{2})\cap L^{2}(\mathbb{R}^{2})$ functions for an arbitrary $k\in\mathbb{N}$. This result holds even if the sampling density is allowed to depend on the function $f$.

3.1 Sampling theorem and Paley-Wiener spaces

The sampling theorem states that every band-limited function $f$ can be written as a sum of shifted cardinal sine functions weighted by point samples of $f$. In other words, a band-limited function is fully determined by its values on a discrete grid. To give the precise statement, we introduce the Paley-Wiener spaces. Given $\Lambda>0$, the Paley-Wiener space $\mathcal{PW}_{\Lambda}\subset L^{2}(\mathbb{R}^{d})$ is defined as

[TABLE]

We define the $d$-dimensional sinc-function as

[TABLE]

Bearing in mind the above notation, we now state the sampling theorem, see, e.g., [37].

Theorem 3.1** (Sampling theorem).**

Let $d\in\mathbb{N}$, $f\in L^{2}(\mathbb{R}^{d})\cap C(\mathbb{R}^{d})$ and $\Lambda>0$. Then

[TABLE]

In particular, for every summable sequence $(y_{k})_{k\in\mathbb{Z}^{2}}$ we can define

[TABLE]

and, by Theorem 3.1, $f$ is band-limited. Furthermore, since $\operatorname{sinc}_{2}(m-n)$ vanishes for every $m,n\in\mathbb{Z}^{2}$ such that $m\neq n$, we observe that $f(m/\Lambda)=y_{m}$ for all $m\in\mathbb{Z}^{2}$. In other words, every sequence on a grid defines an associated interpolating band-limited function and conversely, every band-limited function is uniquely determined by its values on a discrete grid.

As a consequence, the problem of extracting the wavefront set of a function $f\in L^{2}(\mathbb{R}^{2})\cap C(\mathbb{R}^{2})$ from its discrete sampled values $(f(m/\Lambda))_{m\in\mathbb{Z}^{2}}$ can be re-stated as extracting the wavefront set of $f$ from its projection onto a Paley-Wiener space, i.e., $P_{\mathcal{PW}_{\Lambda}}(f)\eqqcolon P_{\Lambda}(f)$.

For functions $f\in L^{2}(\mathbb{R}^{2})$, which are only defined almost everywhere, we can even use the sampling theorem as a definition of a point evaluation. For $f\in L^{2}(\mathbb{R}^{2})$, we have that $P_{\Lambda}(f)\in C(\mathbb{R}^{2})$. Hence, we set $f(m/\Lambda)\coloneqq P_{\Lambda}(f)(m/\Lambda)$.

3.2 Wavefront set extractors

As already stated, the problem of extracting the wavefront set from samples on a grid is equivalent to extracting the wavefront set given the projection onto a Paley-Wiener space. There are multiple conceivable notions of a wavefront set extractor. First, for $\Lambda>0$, we could ask for a map

[TABLE]

Here $P(\mathbb{R}^{2}\times\mathbb{S}^{1})$ denotes the power set of $\mathbb{R}^{2}\times\mathbb{S}^{1}$. Essentially, this map requests extraction of the wavefront set of a function $f$ from knowledge of the samples of $f$ on a fixed grid. It is clear that such a map, $\operatorname{DWF}_{\Lambda}$, cannot exist, since it will fail for functions $f$ that have fine structures which cannot be detected by coarse sampling. For example, a function that vanishes on every grid point of $\mathbb{Z}^{2}/\Lambda$ while having a non-trivial wavefront set would be classified the same as the zero function.

A more reasonable model for a wavefront set extractor should give an approximate prediction of the wavefront set that eventually improves as the sampling density increases. To weaken this statement even further, we might only ask for approximate extraction of the wavefront set at one point and only for functions from a fixed function class $\mathcal{C}\subset L^{2}(\mathbb{R}^{2})$. For a fixed set $W\subset\mathbb{R}^{2}\times\mathbb{S}^{1}$ and a point $x\in\mathbb{R}^{2}$, we define

[TABLE]

We can now model the approximation described above by considering a sequence of wavefront set extractors given by

[TABLE]

such that, for fixed $x\in\mathbb{R}^{2}$, and all $f\in\mathcal{C}$,

[TABLE]

Here $d_{H}$ denotes the Hausdorff distance with the convention $d_{H}(X,\varnothing)=d_{H}(\varnothing,X):=1$ for any non-empty compact subset $X\subset\mathbb{S}^{1}$ and $d_{H}(\varnothing,\varnothing):=0$. Recall that with this definition $d_{H}$ is a metric on compact subsets of $\mathbb{S}^{1}$ (including the empty set).

A sequence as in (3.2) satisfying (3.3) yields an approximate extraction of the wavefront set of $f$ at $x$ from point samples of $f$ where the sampling density may depend on $f$. This observation motivates the following definition.

Definition 3.2**.**

Let $\mathcal{C}\subset L^{2}(\mathbb{R}^{2})$. A sequence $(\operatorname{DWF}_{j})_{j\in\mathbb{N}}$ of mappings as in (3.2) is called an approximate wavefront set extractor. We say that an approximate wavefront set extractor is

• •

clairvoyant* at $x\in\mathbb{R}^{2}$ if the sequence satisfies (3.3) at $x$ for all $f\in\mathcal{C}$, and*

• •

ignorant* to $f\in\mathcal{C}$ at $x\in\mathbb{R}^{2}$ if $d_{H}(\overline{\operatorname{DWF}_{j}(P_{j}(f))_{x}},\overline{\operatorname{WF}(f)_{x}})\not\to 0$ as $j\to\infty$.*

3.3 Non-existence of clairvoyant approximate wavefront set extractors

We will observe below that, for every $x\in\mathbb{R}^{2}$, there does not exist a clairvoyant approximate wavefront set extractor if $\mathcal{C}$ contains at least all $k$-times differentiable functions, for a $k\in\mathbb{N}$. Even more severely, in this situation, every approximate wavefront set extractor is ignorant to a dense subset of $\mathcal{C}$ at $x$.

Theorem 3.3**.**

Let $k\in\mathbb{N}$ and $\mathcal{C}\subset L^{2}(\mathbb{R}^{2})$ be such that $\mathcal{C}\supset C^{k}(\mathbb{R}^{2})\cap L^{2}(\mathbb{R}^{2})$. For every $x\in\mathbb{R}^{2}$ we have that, for every approximate wavefront extractor $(\operatorname{DWF}_{j})_{j\in\mathbb{N}}$, there exists a dense subset $\mathcal{M}\subset\mathcal{C}$ such that $(\operatorname{DWF}_{j})_{j\in\mathbb{N}}$ is ignorant to all $f\in\mathcal{M}$ at $x$. In particular, no approximate wavefront set extractor is clairvoyant at $x$.

Proof.

The proof proceeds in two steps. First, for a given approximate wavefront set extractor, $(\operatorname{DWF}_{j})_{j\in\mathbb{N}}$, and a point $x\in\mathbb{R}^{2}$, we construct a function $q\in\mathcal{C}$ such that $(\operatorname{DWF}_{j})_{j\in\mathbb{N}}$ is ignorant to $q$ at $x$. Second, we show that the set of such functions is dense in $\mathcal{C}$.

Step 1: Notice that Definition 2.1 implies

[TABLE]

For $x\in\mathbb{R}^{2}$, we now choose a function $g\in C^{k}\cap L^{2}(\mathbb{R}^{2})$ such that for a non-empty set $V$ we have that $\operatorname{WF}(g)\supset\{x\}\times V$. Since $k<\infty$, such a function always exists. We can also assume hat $\|(1+|\cdot|^{2})^{k/2}\widehat{g}\|_{L^{1}(\mathbb{R}^{2})}<\infty$.

Then, by (3.4), we can conclude that $\operatorname{WF}(g-P_{j}g)\supset\{x\}\times V$ holds for every $j\in\mathbb{N}$. Moreover, by construction, we have $d_{H}(\overline{\operatorname{WF}(g)_{x}},\varnothing)=1$. To define the desired function $q$, we first set

[TABLE]

for all $n\geq 1$. We have that by the Riemann-Lebesgue Lemma, for every $N\in\mathbb{N}$,

[TABLE]

Moreover, by definition,

[TABLE]

Hence, by the Pythagorean theorem, for every $N\in\mathbb{N}$,

[TABLE]

It now follows from (3.6), (3.7) and (3.3) that $q_{n}$ is a Cauchy sequence in $C^{k}(\mathbb{R}^{2})$ and $L^{2}(\mathbb{R}^{2})$. Therefore $q_{n}$ converges to a limit $q\in L^{2}(\mathbb{R}^{2})\cap C^{k}(\mathbb{R}^{2})\subset\mathcal{C}$. Furthermore, one of the following statements holds:

1. (1)

$\overline{\operatorname{DWF}_{j}(P_{j}q_{j})_{x}}$ does not converge for $j\to\infty$; 2. (2)

$\overline{\operatorname{DWF}_{j}(P_{j}q_{j})_{x}}$ converges to a limit $W$ such that $d_{H}({W},\overline{\operatorname{WF}(q)_{x}})\geq 1/2$; 3. (3)

$\overline{\operatorname{DWF}_{j}(P_{j}q_{j})_{x}}$ converges to a limit $W$ such that $d_{H}({W},\overline{\operatorname{WF}(q)_{x}})<1/2$.

In Cases (1) and (2), we directly obtain that $\operatorname{DWF}_{j}$ is ignorant to $q$ at $x$. In Cases (3), we obtain that there exists some $j_{0}$ such that

[TABLE]

We now consider the cases $\operatorname{WF}(q)_{x}=\varnothing$ and $\operatorname{WF}(q)_{x}\neq\varnothing$ separately. If $\operatorname{WF}(q)_{x}=\varnothing$, then (3.8) implies that $\operatorname{DWF}_{j}(P_{j}q_{j})_{x}=\varnothing$ for all $j\geq j_{0}$ since no subset of $P(\mathbb{S}^{1})\setminus\{\varnothing\}$ has a distance less than $1$ to the empty set. Therefore,

[TABLE]

We obtain from (3.9) that $\varnothing=\operatorname{WF}(q)_{x}=\operatorname{WF}(g)_{x}\neq\varnothing$ which is a contradiction.

If $\operatorname{WF}(q)_{x}\neq\varnothing$, then $d_{H}(\overline{\operatorname{WF}(q)_{x}},\varnothing)=1$. By the triangle inequality, this yields that there exists some $j_{0}$ such that $\operatorname{DWF}_{j}(P_{j}q_{j})_{x}\neq\varnothing$ for all $j\geq j_{0}$. Therefore, $q=q_{j_{0}}\in\mathcal{PW}_{j_{0}}$ by definition, which implies that $W(q)_{x}=\varnothing$. Hence, Case (3) does not occur, i.e., $(\operatorname{DWF}_{j})_{j\in\mathbb{N}}$ is ignorant to $q$ at $x$.

Step 2: For an arbitrary $f\in\mathcal{C}$, there exists $j_{1}\in\mathbb{N}$ such that

[TABLE]

Define $q_{j_{1}}=P_{j_{1}}f$ and, for every $n\geq j_{1}$, we define $q_{n}$ as in (3.3). It is clear that $q_{n}$ converges to a limit $q_{f}\in C^{k}(\mathbb{R}^{2})\cap L^{2}(\mathbb{R}^{2})$. Also, it is straightforward to show that $\|q_{f}-f\|_{2}\leq\epsilon$. Now using the same arguments as in Step 1, it follows that $(\operatorname{DWF}_{j})_{j\in\mathbb{N}}$ is ignorant to $q_{f}$. ∎

Remark 3.4**.**

• (1)

The Theorem 3.3 and its proof also hold when “wavefront set” is replaced by “$(k+1)$-wavefront set” or “singular support”.

• (2)

The arguments in the proof of Theorem 3.3 are independent from the domain $\mathbb{R}^{2}$. Indeed, the same result holds for functions defined in an open domain $\Omega\subset\mathbb{R}^{2}$ and $x\in\Omega$. Here we define the wavefront set of $f\in L^{2}(\Omega)$ as

[TABLE]

• (3)

Theorem 3.3 demonstrates that there does not exist any clairvoyant approximate wavefront set extractor for sufficiently large function classes. Even more severely, for every $k\in\mathbb{N}$, every approximate wavefront set extractor fails on a dense subset of $L^{2}(\mathbb{R}^{2})\cap C^{k}(\mathbb{R}^{2})$.

On the other hand, if the function class is so small that every function is uniquely determined by its values on the grid, then one can construct a wavefront set extractor by storing the wavefront set for each function in a data base.

This shows that for most classical function classes it is impossible to analytically derive a wavefront set extractor. For function classes in applications, such as sets of images, wavefront set extractors could exist, but are potentially highly sensitive to the choice of $\mathcal{C}$.

4 Wavefront set and edge detectors

The characterization of a wavefront set given in (2.1) implies that detecting the wavefront set of a piece-wise smooth function with singularities along a smooth curve is equivalent to detecting edges and their normal direction. Edge detection, which is one of the most well-studied problems in image processing, is therefore a sub-problem to wavefront set extraction. However, we are unaware of any wavefront set extractor in the literature that goes beyond edge or ramp detection.

We next recall some edge-orientation detectors from the literature. We start by approaches based on filters, which is followed by a review of methods based on directional systems from applied harmonic analysis. We then compare these methods to the guiding principle. Finally, we will comment on recent data-driven algorithms for edge detection.

4.1 Filter-based edge-orientation detectors

The traditional way of detecting edges in digital images is to convolve the image with suitable convolution kernels to enhance edge-like features. These features can then be extracted using simple rules. For example, convolution with local difference filters leads to the Roberts [45], Sobel [14, 50], and Prewitt [40] operators. In [38] the authors convolve an image with the Laplacian of a Gaussian function. Here, the zeros of the resulting image are taken as estimate for the positions of edges in the original. Potentially the most famous edge-orientation detection algorithm of this category is the Canny edge detector [11], which is based on convolution with a Gaussian kernel with standard deviation $\sigma>0$. The obvious drawback of this algorithm is that the choice of the standard deviation $\sigma$ of the Gaussian window strongly influences the performance of the algorithm to the extent that a high $\sigma$ improves the robustness against noise, but might also remove high-frequency components.

There have been many efforts in addressing the shortcomings of the Canny edge detector. In [39] one extracts the positions and orientations of not only jump singularities but also of composite edges like ramp- and hat-like singularities. Most recently, [2] develops algorithms that apply higher-level heuristics to a set of patches retrieved from applying oriented gradient operators.

4.2 Edge-orientation detectors based on directional systems

An alternative to the filter-based approaches is to base the detector on multiscale directional transforms. These transforms include shearlets, which were introduced in Subsection 2.2, as well as curvelets, ridgelets, wavelets, or bandlets, see [29] for a survey.

The idea is to first transform the given image using a directional system. Next, one applies certain heuristics to the result by using the theoretical information on the behavior of the underlying transform at directed edge points. In principle, the transforms given by directional systems mentioned above can be written as a series of convolutions with respect to directed filters. Hence, this approach can, in a sense, be understood as a special case of the edge detectors in the previous section.

An example of such an approach, which served as the main inspiration of the forthcoming algorithm in Section 6, is the shearlet-based algorithm of [54], which we will now explain in more detail. This algorithm seeks to detect and classify edges of piece-wise smooth functions $f$ with piece-wise smooth jump curves. The algorithm is based on computing a sampled version of the shearlet transform of $f$, i.e., the algorithm computes

[TABLE]

for given $j_{0},J\in\mathbb{N}$ and $c>0$. This algorithmic step can be performed using standard software libraries such as the ShearLab toolbox [34], see also Subsection 6.1. This computation is followed by a series of heuristics that are applied to the coefficients (4.1) and which lead to a classification of points as point singularities, directed edges, corners, or smooth points.

Additionally, we mention an algorithm in [44] that is based on the complex shearlet transform. This algorithm computes two shearlet transforms, one with a symmetric and one with an anti-symmetric generator. The relationship between both is then used to determine if a point is an edge, a ridge, or a smooth point.

4.3 Data-driven edge-orientation detectors

As we have already advocated previously, the only way to comply with the guiding principle in applications is to use a data-driven algorithm. For edge detection, this idea has already been followed in a series of papers. For example, [13] uses supervised learning of an edge classifier based on a technique called probabilistic boosting tree. Another approach is DeepEdge that uses a convolutional neural network taking as an input candidate edges produced by a Canny edge detector as well as patches of the original image [4]. To this input one then applies parts of the KNet [32] for feature extraction and a network with two fully-connected layers for classification. Finally, the algorithms SEAL (simultaneous edge alignment and learning) [56] and the CASENet (category-aware semantic edge detection network) [55] perform high-level edge analysis with highly complex deep neural networks. The underlying CASENet is a 101-layered network. These methods are state-of-the-art for segmentation and edge detection.

5 Higher-order wavefront set extraction and inverse problems

The discussion of the previous section was limited to edge-orientation detection. In the context of certain inverse problems, also the detection of higher-order wavefront sets is relevant. In particular, this is the case for inverse problems, whose forward operator is either pseudo differential or Fourier integral. We will focus on the so-called X-ray transform which is the basis of the computed tomography (CT) imaging technique.

The goal of computed tomography (CT) is to get an image of the internal structure of an object by striking X-rays through the object coming from a known source. The rays travel on a line $L$ from the source through the object to a detector, being attenuated by the material on the specific line. One can mathematically model this measuring procedure by defining the Radon transform:

Definition 5.1** (Radon transform, [42]).**

Let $f\in L^{1}(\mathbb{R}^{2})$, $\varphi\in[0,\pi)$ and $s\in\mathbb{R}^{2}$. Moreover, the unit vector in direction $\varphi$ denoted as $\theta$ and its orthogonal vector $\theta^{\perp}$ are given by:

[TABLE]

The Radon transform of $f$ is given by the line integral

[TABLE]

where the line $L(\varphi,s)$ is defined as

[TABLE]

We will refer to $g\coloneqq\mathcal{R}(f)$ as the sinogram of $f$. The task of CT reconstruction is to recover the function $f$ from noisy samples of its sinogram. This task is a classic example of an ill-posed inverse problem [15]. The standard approach to overcome the ill-posedness is to incorporate a priori information, an approach that is known as regularization. There are plenty of regularization methods for computed tomography ranging from established model-based methods [5, 16, 48], to more recent end-to-end data-driven approaches, for the latter see the survey in [3]. Additionally, the sparsity of the solution under the shearlet representation [8] or the adjoint operator [1] can be used for regularization. Below, we outline a general approach to implement a priori knowledge of the wavefront of the sinograms into a regularization procedure.

Assuming for the moment that we could successfully build a digital wavefront set extractor. In the context of the tomography problem, this algorithm could be applied to the sampled sinogram $g$. To compute from this knowledge the wavefront set of $f$, we can invoke a certain relationship between the wavefront set of $f$ and its sinogram $g$. This relation between the wavefront set of a distribution and the wavefront set of its transformation under a Fourier integral operator (FIO) is known as micro-canonical relation. The Radon transform $\mathcal{R}$ and its adjoint $\mathcal{R}^{*}$ are examples of FIOs (see e.g. [42]), and its micro-canonical relation is described in the next theorem.

Theorem 5.2** (Radon transform micro-canonical relation, [42]).**

Let $f$ be a distribution with compact support. Let $x_{0}\in L(s_{0},\varphi_{0})$, $\theta_{0}=\theta(\varphi_{0})$, $\lambda_{0}=(1,-x_{0}\cdot\theta_{0}^{\perp})/(1+(x_{0}\cdot\theta_{0}^{\perp})^{2})^{1/2}$ and $k\in\mathbb{N}_{0}$. The micro-canonical relation is the wavefront set correspondence

[TABLE]

in particular, for $\operatorname{WF}(f)=\bigcup_{k\in\mathbb{N}_{0}}\operatorname{WF}_{k}(f)$ it holds

[TABLE]

The relation (5.2) naturally induces an invertible map $\mathbf{can}:\mathbb{R}^{2}\times\mathcal{P}(\mathbb{S}^{1})\to\mathbb{R}^{2}\times\mathcal{P}(\mathbb{S}^{1})$ such that

[TABLE]

holds for all distributions $f$ with compact support.

After extracting the wavefront set of the sinogram one can use it in a regularization functional to regularize the inverse problem, i.e. minimize jointly the residual $\|\mathcal{R}f-g\|_{2}^{2}$ and the functionals

[TABLE]

where $d$ is a suitable distance function between wavefront sets.

The method that we propose below is able to compute $\operatorname{WF}(g)$ or a digitized version thereof very efficiently.

6 Computing the digital wavefront set with shearlets and deep learning

We propose an algorithm that replaces the heuristic approach of the shearlet-based edge detection and classification algorithm of [54] by a data-driven approach. Concretely, instead of hand-crafted heuristics, we train a deep neural network using a variety of training data, adapted to the classification procedure at hand. The neural network takes as input the shearlet coefficients of an image and produces a set of point-direction pairs that are classified as elements of the wavefront set. We will present the construction of the classifier below and then present the computational realization of our algorithm in Subsection 6.4 at the end of this section.

6.1 Digital shearlet transform

The classifier to be constructed below is based on the shearlet transform of a digital image. Therefore, we need to work with a digitized shearlet transform, defined on a digital domain of pixel images. The digital shearlet transform was introduced in [35] and is defined as follows:

Let $M\in\mathbb{N}$, $J\subset\mathbb{N}$ be finite, $k_{j}\subset\mathbb{N}$ for all $j\in J$ and $K_{j}\coloneqq[-k_{j},\dots,0,\dots,k_{j}]$. We pick $2\sum_{j\in J}K_{J}+1$ matrices in $\mathbb{R}^{M\times M}$. We denote these matrices by $\phi^{dig}$ and $\psi_{j,k,\iota}^{dig}\text{ for }j\in J,k\in K_{j},\iota\in\{-1,1\}$. To make the connection to the classical shearlet transform, we can think of $\psi_{j,k,\iota}^{dig}$ as a digitized version of $\psi_{2^{-j},2^{-j/2}k,0,\iota}$ and of $\phi^{dig}$ as a digitized version of a low frequency filter. A concrete construction of the matrices $\phi^{dig}$ and $\psi_{j,k,\iota}^{dig}$ can be found in [35]. Then, we define the digital shearlet transform of an image $I\in\mathbb{R}^{M\times M}$ by

[TABLE]

where $j\in J,k\in K_{j}$, $m\in\{1,\dots,M\}^{2}$, and $T_{m}:\mathbb{R}^{M\times M}\to\mathbb{R}^{M\times M}$ circularly shifts the entries of the elements of a matrix by $m$.

Thus the digital shearlet transform of an image $I\in\mathbb{R}^{M\times M}$ is a stack of $2\sum_{j\in J}(K_{j}-1)+1$ matrices of dimension $M\times M$. In all our numerical experiments, we fixed $J=4$ and $K_{j}=2^{\lceil j/2+1\rceil}+1$ so $2\sum_{j\in J}(K_{j}-1)+1=49$. The computation of the digital shearlet transform is performed by using the Julia implementation of ShearLab [34] (www.shearlab.org/software).

6.2 Network architecture

We train a neural network which, given a patch of the shearlet coefficients of a function, produces a prediction of which directions belong to the wavefront set of the function at the position associated with this patch. These patches are of size $21\times 21\times 49$.

The network architecture consists of four convolutional layers, with $2\times 2$ max pooling, ReLU activation and batch normalization, followed by a fully-connected layer with 1024 neurons, softmax activation function and a one-dimensional output. The network architecture is depicted in Figure 1. We chose this architecture since it performed well in a series of tests while being of moderate size. Here we focused on networks with only a few layers because we expect that the shearlet transform already acts as the correct feature extractor of the problem. Therefore, the classifier does not need to learn the correct data representation. Nonetheless, it is conceivable that a deeper and larger neural network architecture could potentially lead to improvements for the classification results below.

We pick 180 directions $(\theta_{i})_{i=1}^{180}$. For each $\theta_{i}$, we then train a network $\Phi_{i}$ with the described architecture by passing patches of shearlet coefficients of images $I\in\mathbb{R}^{M\times M}$ of the form

[TABLE]

where $m^{*}\in\{11,\dots,M-10\}^{2}$, to the network. The associated label to a batch of (6.1) is $1$ if $I$ has a singularity with direction $\theta_{i}$ at $m^{*}$ and [math] else. In total, this procedure yields $180$ digital classifiers. We train one more network with the same data, but the label is $1$ if $I$ has no singularity at $m^{*}$ and [math] else. This additional classifier is used in test cases where all competing algorithms only perform edge detection and not edge-orientation detection.

The final classifier is constructed by putting all of these 181 networks in parallel, producing one large network with 181 outputs. For every $21\times 21\times 49$ patch of shearlet coefficients, this classifier generates a vector of length 181 indicating if the underlying function is smooth at the center point of the patch and listing all directions of edges present at the center point.

6.3 Network training

We train the network as described in Subsection 6.2 using stochastic gradient descent to minimize the cross-entropy over a variety of training sets. We use five different data sets to train our classifier and test our algorithm:

The first data set consists of patches of the shearlet transform of images made of random sums of ellipses and parallelograms of different contrasts, sizes, and orientations. 2. 2)

The second data set contains random sums of ellipses and parallelograms convolved with a kernel to generate functions with a higher-order wavefront set. 3. 3)

The third data set is the BSDS500 (Berkeley Segmentation data set) provided by the Computer Vision Group of UC Berkeley. It comprises 503 natural images of different types. 4. 4)

The fourth data set is the Semantic Boundaries data set (SBD) with 11355 natural images, again provided by the Computer Vision Group of UC Berkeley. 5. 5)

The fifth data set is constructed by taking Radon transforms of phantoms made of ellipses. The associated wavefront sets are computed analytically through the canonical relations of Definition 6.1.

We depict examples of functions from each of the data sets in Figures 2, 3, 4, 5, and 6.

To make these data sets suitable for our purposes, we need to equip each image of the data sets with an associated set of labels indicating the associated wavefront set or the set of edges. For the first two and the last data set, standard theoretical results on the wavefront sets of characteristic functions allow us to compute the associated wavefront sets analytically. The segmentation and semantic boundaries data sets, on the other hand, are natural images where such an approach is not possible. These data sets are used to assess the quality of segmentation and contour detection applications, see [25] and [2]. Therefore, every image in these data sets was annotated and has a set of ground truth edges. However, we should point out that this annotated ground truth does not contain all edges of the images, but only those between semantically different parts of the images. We depict the annotated edges in Figures 2, 3, 4, 5, and 6.

In the following subsections, we describe the computation of the associated wavefront sets in detail.

6.3.1 Ellipses and parallelograms

The wavefront sets of characteristic functions of ellipses and parallelograms can be identified by (2.1) and the fact that if $x$ is a vertex of a parallelogram $P$ then $\{x\}\times\mathbb{S}^{1}\subset\operatorname{WF}(\chi_{P})$. For sums of these functions, we have, by basic properties of the Fourier transform that

[TABLE]

Note that in this relation we do not have equality in general. Indeed, if $\operatorname{WF}(\chi_{P_{1}})\cap\operatorname{WF}(\chi_{P_{2}})\neq\varnothing$ then cancellations can occur. We shall neglect this technicality as the probability of cancellations is sufficiently small and assume that the wavefront set of characteristic functions as described above is the union of the respective wavefront sets.

We build this data set by randomly choosing a number of parallelograms and ellipses with random positions and computing the associated ground truth of the wavefront set as described above.

6.3.2 Higher-order wavefront data set

The ellipses/parallelograms data set contains images with jump singularities only. To test our method on functions with higher-order singularities, such as ramp singularities, we computed the convolution of the elements of the ellipses/parallelograms data set with a filter $h$ the Fourier transform of which is given by:

[TABLE]

It is not hard to see that $P:f\mapsto h*f$ is an elliptic pseudo-differential operator and hence $\operatorname{WF}(h*g)=\operatorname{WF}(g)$ for all $g\in L^{2}(\mathbb{R}^{2})$, see [20, Chapter 8 G] for details. Thus, the convolutions of the elements of the ellipses/parallelograms data set with $h$ have the same wavefront set as the associated ellipses or parallelograms, but of a higher order.

6.3.3 Segmentation and semantic boundaries data sets

In the BSDS500 and the SBD data sets, the ground truth of the edges is given in form of binary images with 0’s at positions where the image is smooth and 1’s at locations associated to edges. This annotated edge set is depicted in Figure 4.

To compute the orientation of the edges, we used a five-point stencil derivative on the edges to approximate the normal vectors. To detect corners and assign the appropriate orientations we used the Harris corner detector [26]. From these images, we produce patches for the training of the network classifier. However, due to the fact that the annotated image does not contain all edges we only use patches that are close to these edges for training, validation and testing.

6.3.4 Tomographic data

We generate tomographic data simulated from 2D phantoms similar to the commonly used Shepp-Logan phantom. Each phantom is represented by an $N\times N$ matrix, where $N\times N$ is the resolution of the images. As in the case of the ellipses/prallelograms data set, these phantoms have an analytic wavefront set given by the union of the wavefront sets of all the ellipses. These are stored as binary matrices of dimensions $N\times N\times 180$ where a $1$ at position $(i,j,\theta)$ indicates that the direction with angle $\theta$ is in the wavefront set at point $(i,j)$. Next, we simulate tomographic data (sinograms) using the Radon transform implementation in the ODL (http://github.com/odlgroup/odl) python library for inverse problems.

To label each sinogram with the correct wavefront set, we map the known wavefront sets of the corresponding phantom using a digitized version of the microlocal canonical relation for the Radon transform (5.3) defined as follows:

Definition 6.1**.**

Let $N\in\mathbb{N}$, we define

[TABLE]

where $Y_{s,\varphi,\lambda}=1$ if and only if

[TABLE]

for $(i_{x},j_{x},\theta_{x})$ such that $X_{i_{x},j_{x},\theta_{x}}=1$. We call $\mathbf{can}^{N}$ the digital canonical map.

Moreover, we define an inverse digital canonical map by

[TABLE]

where $X_{i,j,\theta}=1$, if and only if

[TABLE]

for a $(s_{y},\varphi_{y},\lambda_{y})$ such that $Y_{s_{y},\varphi_{y},\lambda_{y}}=1$.

We take Definition 6.1 as an Ansatz for the definition of the digital wavefront set of the digital Radon transform of a function of which we know the wavefront set. We depict one example of the elements of this data set in Figure 6.

6.4 DeNSE: Deep Network Shearlet Edge Extractor

We present our algorithm extracting the wavefront set of a digital image below. For $M\in\mathbb{N}$, and a digital image $I\in\mathbb{R}^{M\times M}$, this algorithm produces, for every $m^{*}\in[11,M-10]^{2}$ a prediction of the wavefront set of $I$ at $m^{*}$. The algorithm proceeds along the following three steps:

Step 1:

Train the network classifier as of Subsection 6.2 on a set of labeled training data.

Step 2:

For a given test image $I\in\mathbb{R}^{M\times M}$, compute the digital shearlet transform of $I$ with $49$ shearlet generators: $(\mathrm{DSH}(I)(j,k,m,\iota))_{j\in J,k\in K_{j},\iota\in\{-1,0,1\},m\in[1,M]^{2}}$.

Step 3:

For every $m^{*}=(m^{*}_{1},m^{*}_{2})\in[11,M-10]^{2}$, pass the patch

[TABLE]

to the classifier of Step 1. If the classifier predicts that an edge with direction $\theta$ is present, then classify $(m^{*},\theta)$ as an element of the wavefront set of $I$.

Henceforth, we refer to the above algorithm as Deep Network Shearlet Edge Extractor (DeNSE).

7 Numerical results

We implemented the training as described in the previous section using the GPU version of Tensorflow. To evaluate the classification quality, we use two quality measures, a mini-batch test average taken over all mini-batches and the so-called MF-score. The MF-score is computed as the mean of the F-score defined as

[TABLE]

where $P$ is the precision, i.e., the number of true positives divided by the sum of true and false positives, and $R$ is the recall, i.e., the number of true positives divided by the sum of true positives and false negatives, [46]. The MF-score is often used for evaluating classification performance when the distribution of classes is uneven. This is, for example, the case in edge detection, since there usually are significantly fewer edge points than smooth points in an image. Moreover, these performance measures enable us to compare with the state-of-the-art [56] on the respective data sets. In addition, the required code to reproduce the results is publicly available in http://www.shearlab.org/applications.

The implementations of the other methods that were used to compare the performance of DeNSE were taken from the publicly available github repositories provided by the authors of the methods (all available in the topic: http://github.com/topics/edge-detection), with exception of the Canny, Sobel algorithms that were taken from the python library OpenCV (http://opencv.org/). The data-driven methods were trained using the hyper-parameters proposed by the authors for the given data set (Berkeley segmentation set and SBD) without further tuning.

7.1 Results for ellipses/parallelograms

We train each of the 181 subnetworks as of Subsection 6.2 using 10,000 images as training data, 1,000 images as validation data, and 2,000 images as test data. For each direction $\theta_{i}$ we trained the associated subnetworks using a mini-batch procedure with 86 examples per batch and 3,000 training steps for each. We obtained an average test accuracy of 96.2% (taking the average over all 181 classifiers) and an MF-score of 97.1%. We also notice that the test accuracy of the individual classifiers was higher when classifying angles aligned to the discrete orientations of the underlying shearlet system.

We compared our method on the data set of the shearlet coefficients patches to other classifiers commonly used in machine learning namely: Logistic regression, Decision trees, $K$-nearest neighbors, Linear SVM, and Random forest, using their implementations in the python library scikit-learn (http://scikit-learn.org/). We report the performances of these classifiers in Table 1.

By construction, the last of the 181 subnetworks corresponds to an edge-detector, where the achieved average test accuracy was 97.5%, and the MF-score was 97.9%, the performance benchmarks with other classical edge classifiers can be found in Table 2. Figure 2 shows the results on an example of the ellipses/parallelograms data set.

We depict the classification for one instance of the test set of the parallelograms/ellipses data set in Figure 2 and compare the results with the classification by the heuristic approach by Yi-Labate-Easley-Krim [54]. We observe that our method performs significantly better in low contrast regions. Moreover, our algorithm appears to be more precise when differentiating between corners and edges. Here, we classify a point as a corner point if the classifiers predict at least two different orientations that differ by more than 10 degrees. In Figures 2,  3, 4 and 5, we indicate corners by white dots.

7.2 Results for higher-order wavefront set data set

We performed wavefront set detection for the the higher-order wavefront set data set, using the same procedure as in the ellipses/parallelograms classification. In this case, we used 30,000 patches as training data, 3,000 patches as validation data, and 6,000 patches as test data. We trained on 86-sized mini-batches, with 200,000 training steps. We obtained an average test accuracy of 93.4% and an MF-score of 94.6%. We are not aware of any algorithms specifically build for higher-order wavefront set detection, which is why we do not provide a comprehensive list of results of alternative algorithms in this case.

For completeness, we added Figure 4 showing an example of the obtained results. We also add two predictions by the algorithm of Yi-Labate-Easley-Krim [54] and the method CoShREM [44].

It is important to mention that the algorithm of Yi-Labate-Easley-Krim is constructed to detect jump singularities and not higher-order singularities. Hence this algorithm is expected to fail on this data set. Indeed, the performance of the algorithm of achieves only an MF-score of 30.5%. CoShREM, on the other hand, is built to detect edges and ridges. The performance was significantly better than that of Yi-Labate-Easley-Krim and resulted in an MF-score of 65.4%.

7.3 Results for Berkeley segmentation set

In the Berkeley segmentation data set, the complexity of the images is considerably higher compared to the images from the ellipses/parallelogram data set. Therefore, we use a significantly larger training set to train the associated classifier. For the classification of each angle, we used 30,000 patches as training data (around 600 patches per image), 3,000 patches as validation data, and 6,000 patches as test data. As in the case of the ellipses/parallelograms, we train using a mini-batch procedure, with 86 examples per batch, but in this case, using 30,000 training steps for each. We obtained an average test accuracy of 93.1% and MF-score of 95.4%, which is lower than the one obtained in the ellipses/parallelogram. This is likely due to the higher complexity of the patches. One advantage of this and the SBD data set is the existence of several benchmarks including state-of-the-art deep learning based algorithms.

We compared our method using the available benchmarks on this data set provided by the UC Berkeley Computer Vision Group, we refer to [2] for a more detailed explanation of these methods. In [2], just the MF-score of the competing algorithms was reported. We give the results in Table 3.

We present one example of the results obtained on the BSDS500 data set in Figure 4.

7.4 Results for semantic boundary set (SBD)

The SBD data set contains significantly more images than the BSDS500 which, as we will observe below, improves the overall classification performance slightly. In this case, we used 100,000 patches as training data, 10,000 patches as validation data, and 20,000 patches as test data. We train on 86-sized mini-batches, with 100,000 training steps. We obtained average test accuracy of 95.3% and MF-score of 96.8%.

This data set has recently been widely used for image segmentation tasks, in particular, it was used on the two deep learning based image segmentation frameworks proposed by Z. Yu et al., namely the SEAL (Simultaneous Edge Alignment and Learning) [56] and the CASENet (Category-Aware Semantic Edge Detection Network) [55]. We also compared them with the deep learning image boundary detector and classifier proposed (OBDC) by J. Y. Koh et al. [30]. The results can be found on Table 4.

Figure 5 shows the results obtained by DeNSE on an example image of the SBD data set, as in the case of the BSDS500 data set, the obtained result admits more edges than the ground truth due to the batch-based approach. Nonetheless, the method outperforms even the specialized algorithms for segmentation over the given data sets.

7.5 Results for tomographic data

Learning a data-driven wavefront set extractor for tomographic data uses 10 000 images with associated wavefront sets as supervised training data, 1 000 images as validation data, and 2 000 images as test data. In this case, we obtained an average test accuracy of $94.5\%$ and $95.7\%$ of MF-score. This score is comparable to the performances of DeNSE on the other image classes. In this case, we could not perform any comparison with alternative methods, because we are not aware of any competing algorithm for the detection of wavefront sets of sinograms.

Extracting the wavefront set of the data and then applying canonical relations has the computational advantage that we do not need to solve the inverse problem. On top of this, it seems to also give improved wavefront set detection of the original signal over strategies that first solve the inverse problem and then detect the wavefront set. To substantiate this claim, we perform the following test. Using three standard inversion schemes, filtered backprojection, Tikhonov-regularized and total variation regularized inversion, we first compute an approximate inverse from a low dose CT image and then compute the associated wavefront set using DeNSE. We then note the average mean square error to the true wavefront set of the data point. We compare this with the error resulting from computing the wavefront set of the sinogram and transforming this to the wavefront set of the image through the canonical relations. We depict the wavefront sets extracted in this way in Figure 7.

The results of Table 5 clearly demonstrate the advantage of first extracting the wavefront set of the sinogram and then applying the canonical relations over any first-invert-then-extract strategy.

7.6 Conclusion

We observe that in all performed tests, our novel algorithm DeNSE significantly outperforms all competitors. In doing so, DeNSE outperforms not only traditional methods, but also other, deep learning-based algorithms. Comparing the complexity of the involved neural networks reveals that the classifier of DeNSE uses a comparably small neural network.

One natural explanation for the jump in performance already with simple networks is the fact that the shearlet representation transforms the data in a much more convenient form for training and evaluation purposes at least from the point of view of wavefront set extraction.

Moreover, we observed strong performance of DeNSE for extraction of non-edge-like wavefront sets. In particular, the successful detection of the wavefront set of sinograms suggest a possibility to incorporate this information into regularizing schemes of inverse problems or for analysis purposes.

We observed that through the canonical relations this approach allows the detection of the wavefront set of a phantom without solving the inverse problem. This is particularly interesting since this approach does not only circumvent an inversion step, but also improves the performance over algorithms that first invert the Radon transform and then extract the wavefront set.

Acknowledgements

H.A.-L. is supported by the Berlin Mathematical School. G. K. acknowledges partial support by the Bundesministerium für Bildung und Forschung (BMBF) through the Berliner Zentrum for Machine Learning (BZML), Project AP4, by the Deutsche Forschungsgemeinschaft (DFG) through grants CRC 1114 ”Scaling Cascades in Complex Systems”, Project B07, CRC/TR 109 ”Discretization in Geometry and Dynamics’, Projects C02 and C03, RTG DAEDALUS (RTG 2433), Projects P1 and P3, RTG BIOQIC (RTG 2260), Projects P4 and P9, and SPP 1798 ”Compressed Sensing in Information Processing”, Coordination Project and Project Massive MIMO-I/II, by the Berlin Mathematics Research Center MATH+ , Projects EF1-1 and EF1-4, and by the Einstein Foundation Berlin. The work of O.Ö. was supported by the Swedish Foundation of Strategic Research grant AM13-004. P.P is supported by a DFG Research Fellowship ”Shearlet-based energy functionals for anisotropic phase-field methods”.