Learning a microlocal prior for limited-angle tomography

Limited-angle tomography is a highly ill-posed linear inverse problem. It arises in many applications, such as digital breast tomosynthesis. Reconstructions from limited-angle data typically suffer from severe stretching of features along the central direction of projections, leading to poor separation between slices perpendicular to the central direction. A new method is introduced, based on machine learning and geometry, producing an estimate for interfaces between regions of different X-ray attenuation. The estimate can be presented on top of the reconstruction, indicating more reliably the true form and extent of features. The method uses directional edge detection, implemented using complex wavelets and enhanced with morphological operations. By using machine learning, the visible part of the wavefront set is first extracted and then extended to the full domain, filling in the parts of the wavefront set that would otherwise be hidden due to the lack of measurement directions.


Introduction
X-ray tomography aims to recover the attenuation coefficient inside a physical body from a collection of radiographs recorded along different directions of projection.After calibrating the data, one can model the inverse problem mathematically as recovering a non-negative function f : R n → R from a collection of line integrals of f , where n is 2 or 3.This study focuses on limited-angle tomography, where the directions of the lines in the dataset are severely restricted.Such cases appear in many important applications of X-ray tomography, for example, Digital Breast Tomosynthesis (DBT) and weld inspection.
DBT is an enhanced form of mammography where one collects several X-ray images of a compressed breast within a limited angle of view.See Figure 1(a) for an illustration of the DBT measurement geometry.This provides more information about the three-dimensional structure of tissue and potentially leads to improved diagnosis of breast cancer [9,33,21,31,26,23,18,19].
Long pipelines are used, for instance, in the process industry and the energy industry.They are often constructed by welding together several pieces of metal pipes.With welding, there is a risk of gas bubbles, cracks, or other defects forming inside the weld.One can use limited-angle X-ray tomography for non-destructive testing of the welds [14,30] with imaging geometry as shown in 1(b).
Limited-angle tomography is a very ill-posed inverse problem [20,8,29]; in other words, it is very sensitive to modelling errors and measurement noise.In particular, reconstructions typically contain stretching artefacts along the central measurement direction and streak artefacts towards the extreme angle projections [11].This makes it hard to recover the true form and extent of features.However, note that we know from [25] precisely which parts of the boundaries of objects are stably visible in limited-angle data.
We introduce a new, partly data-driven and partly geometric method for recovering interfaces between different areas of X-ray attenuation in the target.We show how complex wavelets and neural networks can be used to first detect the visible parts of the interfaces and then extend them to the invisible parts.Before going into details of our algorithm, we review the related methodology in the literature.
In recent years, deep learning has been used to improve the quality of computed tomography (CT) reconstructions in different ways [1,15,32,4].A common idea is to use a neural network as a postprocessing step to improve the quality of the CT reconstruction [22].While this approach can produce promising results in reducing noise or imaging artefacts, such post-processing neural networks are generally considered "black-box".Additionally, such networks typically require large amounts of training (a) (b) Figure 1: (a) Imaging geometry of digital breast tomosynthesis, based on collecting several radiographs of the compressed breast at various angles.Typically, the most extreme directions are ±20 • from the vertical, as shown in the image.While only three projections are depicted here, real DBT systems collect more (for example 15) images within the 40°angular range available.Limited-angle imaging geometry is forced by the compression of the breast against the detector plate.If the angle becomes too extreme, the rays passing through the tissue will not hit the detector.(b) Imaging geometry for X-ray inspection of weldings in pipelines used for example in process industry.Shown is a longitudinal cross-section of the pipe, with dark gray color indicating metal.The purple area is the weld connecting two pieces of pipe together.Limited-angle imaging geometry is necessary as both the X-ray source and the detector need to be outside the long pipeline.data, which can be impractical in certain applications, for instance in medical imaging.Also, the postprocessing solution might not be optimal for severely ill-posed problems, where the quality of the initial reconstruction is extremely flawed.
An alternative approach is to combine analytical inverse problems algorithms with deep learning.These kinds of reconstruction techniques take advantage of the strengths of both variational regularization and deep learning.Having the explicit forward operator improves robustness and generalizability of the network, compared to previously mentioned "black-box" neural networks.In addition, such model-based reconstruction reduces the amount of training data needed.This is due to the fact that the network has less parameters, and the forward operator encodes the imaging geometry so it does not need to be learned [1,2].
With modern machine learning methods, it has become possible to suppress streak artefacts with unprecedented efficiency.As shown in [5], it is possible to learn only the information which is not stably represented in the limited-angle tomographic data.The authors propose a hybrid framework with shearlets and a convolutional neural network for solving limited-angle CT problems.There, the network fills in only the missing shearlet coefficients, and the rest is done using rigorous inversion methods.Thus, the result can be seen as interpretable, as the network only affects a controlled part of the final reconstruction.In [2], the authors propose a deep learning framework for CT by unrolling a proximaldual optimization method.In their Learned Primal-Dual (LPD) algorithm, convolutional neural networks replace the proximal operators.In [3], the authors combine microlocal analysis to the LPD algorithm and propose a scheme for joint reconstruction and wavefront set inpainting, enabling significant reduction of streak artefacts.
Inherent in deep learning is the aforementioned "black box" problem; in other words, poor understanding of how did the algorithm arrive at the given reconstruction.Especially in medical imaging, it is often important to use algorithms that are as interpretable as possible.With this is mind, our present approach aims to only learn the unstable parts of the boundaries, which then can be displayed as extra information, on top of a classical (non-learned) reconstruction.
We introduce a new algorithm for improved reconstruction: a data-driven and geometry-guided hybrid.We first compute an initial reconstruction of the target using a suitable algorithm (in our case Primal-Dual Fixed-Point (PDFP) method [7] promoting frame sparsity).Due to the limited-angle imaging geometry, the reconstruction contains stretching artefacts.
Next, we aim to detect the stable parts of boundaries within the reconstruction.We transform the reconstruction using Dual-Tree Complex Wavelet Transform (DTCWT), which has directional selectivity oriented at ±15 • , ±45 • and ±75 • .DTCWT does not offer such a theoretically complete discretisation Figure 2: Idealized limited-angle imaging geometry we consider in our computational examples.We restrict to parallel-beam geometry.Shown here are the two extreme angles in our measurement model; we use a total of 50 directions distributed evenly in the angular variable between these two extremes.There is a significant difference in terms of imaging technology between our model and the two applications shown in Figure 1, due to the fan-beam geometry used there.Namely, it is cumbersome to realize parallelbeam measurements with an X-ray source.However, the mathematics of wavefront sets is essentially the same, so our findings should apply with straightforward modifications. of the wavefront (WF) set as curvelets and shearlets do.However, we combine DTCWT with nonlinear morphological operations and a neural network model in a novel way, resulting in a robust six-angle edge detector, capable of automatically excluding streak artefacts.Our learned edge detection algorithm is independently applicable to a wide range of image processing tasks.
Using a neural network model, we can then learn the missing parts of the wavefront set in the framework of six oriented subbands, extending it beyond the scope of the limited-angle measurements.Our approach manages to ignore the streak artefacts and estimates the unstable parts of the boundaries of objects with unprecedented efficiency.Then, we can use this extended information about the wavefront set to form a boundary estimate that can be added as an overlay to the PDFP reconstruction.
We illustrate our method with computational examples, working with parallel-beam geometry as shown in Figure 2.This simplifying assumption helps in computational experiments as we can divide the 3D tomography problem defined in xyz-space into a stack of independent 2D problems posed in xz-planes determined by a constant value of y.The results are viewed by browsing through the xy-planes of the reconstructed volume by varying the z-coordinate.Our algorithm leads to remarkable accuracy in separating objects located away from each other in the z-coordinate but having overlapping projections in the xy-plane.Such objects are typically merged together in traditional limited-angle reconstructions.We note that there is no obstruction in principle in extending our methodology to fan-beam or cone-beam geometries.
The structure of this paper is as follows.In section 2, we go through the main mathematical background information necessary to understand the proposed method.These topics include dual-tree complex wavelet transform (section 2.1), wavefront set, singular support and tomography (section 2.2), Primal-Dual Fixed-Point optimization (section 2.3), and finally, we explain the morphological operations used in the computations (section 2.4).In section 3, we describe our proposed method for learned wavefront set extraction for limited-angle tomography.Then, in section 4, we continue to explain the process of learning a microlocal prior.Next, in section 5, we show the results of our proposed method in the xz-plane (section 5.1) and in the xz-plane (section 5.2) for a simplistic three-dimensional computational phantom.Finally, in section 6, we end with some concluding remarks.

Mathematical background 2.1 Dual-tree complex wavelet transform
The dual-tree complex wavelet transform consists of a complex-valued scaling function φ and a complexvalued wavelet ψ.Consider a complex and analytic wavelet ψ(x), given by and associated with a high-pass filter H, and a complex and analytic scaling function φ(x), that is given by φ and associated with a low-pass filter L. Note that ψ h and φ h are the real parts of ψ(x) and φ(x), and ψ g and φ g are the imaginary parts of ψ(x) and φ(x), respectively.For computational purposes, we use finitely supported wavelets and thus approximately analytic wavelets.We denote the complex wavelet coefficients by d ν (j, n) ∈ C, where j = 0, . . ., J represents the number of scales used in the complex wavelet representation and ν ∈ I = {LH, HH, HL, HL, HH, LH} represents the index of the six oriented subbands in each scale j.Note that j = 0 corresponds to the mother wavelet and j = J corresponds to the finest scale.See Figure 3b for the frequency-domain content indexed by ν.
To build intuition, the expression for the 2D wavelet in subband HH is obtained by: Note that this wavelet is oriented at −45 • .To obtain the 2D wavelet oriented at +45 • , that is the subband HH, we compute ψ 2 (x 1 , x 2 ) = ψ(x 1 )ψ(x 2 ), where The approximate supports of the Fourier transforms of the frame functions corresponding to the six subbands in the dual-tree complex wavelet transform.(c) Shearlets offer very detailed analysis of orientations using parabolic scaling and increasing numbers of orientations with finer scales.However, the benefits come with a computational cost.In a similar manner, the six subbands of the 2D complex wavelets for each scale j are computed as follows: Here n = (n 1 , n 2 ) ∈ Z 2 .For discrete images of the size of power-of-two (which we consider throughout the paper), n 1 and n 2 range from 0 to 2 j − 1, depending on the scale j.In our computations, j will range from 0 to 6, 0 corresponding to the coarsest scale, and 6 corresponding to the finest scale.We remark that the digital wavelet transform starts by computing the finest scale by applying filters to a pixel image.The calculation progresses towards coarser scales by repeatedly applying the filters to the low-pass subbands.
For a more detailed look on the theoretical properties of the complex wavelet transform, please refer to [28].
The complex wavelet coefficients are defined by: Here, the function f can be recovered from the coefficients by where S is the frame operator where T is the analysis operator and T * is the synthesis operator corresponding to S. See Figure 4 for an example of the absolute values of complex wavelet detail coefficients for all six subbands of the finest scale.
The complex wavelet transform also gives rise to real-valued, final-level low-pass scaling coefficients, known as the approximation coefficients.They are given by: The motivation for using complex wavelets is that they offer a good compromise between directionality and computational cost.More precisely, real-valued wavelets offer only very limited information about the orientation of the edges in an image.On the other hand, curvelets [6] and shearlets [16,17] use parabolic scaling, along with rotation or shearing, to achieve remarkably accurate approximation properties for the wavefront set.However, for our purposes, the DTCW transform offers just the right balance between directionality and computational cost.See Figure 3 for an illustration of the frequency content of the building blocks of these three types of transforms.

Wavefront set, singular support and tomography
Microlocal analysis tells us that limited-angle tomography data specifies certain elements of the wavefront set of the X-ray attenuation coefficient, while information about the rest of the wavefront set is not present in any well-posed form.We recall briefly the definitions of singular support and wavefront set here for the convenience of the reader.
The wavefront set of a signal f describes both the location x 0 and the direction θ 0 of singularities.The signal f is smooth near x 0 if there is a cutoff function φ ∈ C ∞ c , φ(x 0 ) = 0, such that the Fourier transform of φf decays rapidly.That is, The signal f has a singularity in x 0 if for all cutoff functions φ, the Fourier transform of φf does not decay rapidly.The set of all singularities of f is called the singular support and is denoted by singsupp(f ).To define the orientation of the singularities x ∈ singsupp(f ), we look for directions along which the localized Fourier transform φf does not decay rapidly.
Define the wavefront set of a function f as the set W F (f )(x 0 , α 0 ) with location x 0 ∈ singsupp(f ) and direction α 0 , such that for all cutoff functions φ ∈ C ∞ c , φ(x 0 ) = 0 , the localized Fourier transform φf (ξ) does not decay rapidly in any polar wedge W δ = {(r, w) : |w − α 0 | < δ}, where (r, w) are the polar coordinates in the frequency domain.The direction α 0 of a singularity x 0 ∈ singsuppf can be considered as the direction of maximum change of f at x 0 .In particular, for a piecewise constant function having a jump along a Jordan curve γ, singsupp(f ) equals γ and W F (f ) consists of pairs (x 0 , N (x 0 )) with x 0 ∈ γ and N (x 0 ) the normal vector of γ at x 0 .

Primal-dual fixed-point (PDFP) optimization
In this study, we use the Primal-dual fixed-point algorithm, described by Chen, Huang and Zhang [7], to solve the following optimization problem.In the discrete setting, the inverse problem of reconstructing a tomographic image f ∈ R n based on X-ray measurements m ∈ R m is modeled by where A ∈ R m×n is a discretized linear forward operator and > 0 models additive Gaussian noise.We consider regularized solutions to the inverse problems, achieved by minimizing the following functional: where the term W C f 1 promotes solutions having sparse representation in the wavelet basis, and µ > 0 serves as a regularization parameter providing a balance between data fidelity and a priori information.
The component-wise inequality f ≥ 0 is based on the physical fact that X-radiation can only attenuate inside the target and not strengthen.Analogous regularization was used in [24], but with real wavelets instead of complex.
The PDFP algorithm can be used to iteratively solve the above minimization problem (equation 16).The algorithm is given by: where τ and λ are positive parameters, , and µ > 0 represents the regularization parameter.Parameters 0 < λ < , where λ max is the maximum eigenvalue, and 0 < τ < 2/τ lip , where τ lip is the Lipschitz constant for G(f ), need to be suitably chosen for convergence.The soft-thresholding operator T is defined radially for complex-valued inputs as In equation (17), the non-negative quadrant is denoted by C = R N 2 + and P C is the Euclidean projection, that is, the operator P C replaces any non-negative elements in the input vector by zeroes.

Morphological operations
In this work, we use morphological operations [12] to clean and pre-process greyscale and binary data.In mathematical morphology, we consider the original image as the object D ⊂ Z 2 .A structuring element S ⊂ Z 2 is a pre-defined binary image that is used to probe the object D, and conclude how it fits the shapes of the object.
The morphological dilation operation can be used to "grow" objects in a binary image.The shape and size of the structuring element S defines the extent of the dilation.The dilation operation is defined as where D s = {d + s|d ∈ D, s ∈ S} is the translation of D by the structuring element s.
On the other hand, morphological erosion thins out objects in a binary image.It is the dual operation of dilation.The erosion operation is defined as The morphological opening operation removes small details from the object but preserves the shape and size of larger objects.The opening of D by S is defined as where D S is the morphological erosion of D by S, and ⊕ is the dilation operation of the result by S. The opening of D by S can be seen as the union of all the translations of S so that it fits within the object D: where The morphological skeleton operation extracts the centerline of the objects in a binary image while preserving its topology.The operation transforms the object D to 1-pixel wide curved lines, while not changing the essential structure.It can be expressed in terms of morphological erosions and openings: where (D kS) expresses k successive erosion operations and K = max{k|(D kS) = ∅} denotes the last step before the object D erodes to an empty set.Morphological operations can be extended for greyscale images.In that case, D in formula ( 22) is a greyscale pixel image.In greyscale morphology, images d are functions mapping the Euclidean space E The greyscale structuring elements are also functions mapping Euclidean space into R ∪ {−∞, ∞}.
Denote an image by d(x) and a structuring function by b(x).Then, the dilation operation is defined as where sup denotes the supremum.The erosion of d by b is given by while the opening and closing are given by

Learned wavefront set extraction
In limited-angle tomography, the reconstructed image typically fails to reproduce certain parts of the wavefront set of the original target, as explained in section 2.2.This is often accompanied by severe stretching artefacts in the reconstruction along the central direction of projection, as illustrated in Figure 6.
Our goal is to fill in the missing parts of the wavefront set by machine-learning the geometric rules of how the known parts of the wavefront set extend into the unknown parts.See Figure 7 for the geometric idea behind this microlocal prior.After we have the full wavefront set available, we can use it to form a boundary estimate of features.
We compute the aforementioned boundary estimate in several steps, including the application of two distinct convolutional neural networks.(a) The first neural network learns to extract the known part of the wavefront set by transforming morphologically opened complex wavelet coefficients to a binary form.(b) The second neural network learns to predict the complete wavefront set from the incomplete wavefront set that has been morphologically dilated.We will discuss the second network, among with the operations related to it, in Section 4.

Complex wavelet coefficients of a limited-angle reconstruction
We simulate a parallel-beam limited-angle sinogram with 50 projection directions and a 40-degree opening angle.Further, we compute a reconstruction using variational optimization with complex wavelet regularization, as described in section 2.3.Because of the limited-angle imaging geometry, the features in the reconstruction are stretched along the central direction of projection, as can be seen in Figure 6.
First, we want to extract the wavefront set (see Section 2.2) from the reconstruction.This is accomplished by picking out a discrete approximation of the wavefront set by using the finest scale of dual-tree  complex wavelet coefficients d ν (j, n 1 , n 2 ), where j = 6, n = (n 1 , n 2 ) ∈ Z 2 , and ν ∈ I is as described in section 2.1.Then, we take the absolute value of the coefficients: A demonstration can be seen in Figure 8a.Due to the limited-angle data, we observe nonzero complex wavelet coefficients mainly in subbands d HL (6, n 1 , n 2 ) and d HL (6, n 1 , n 2 ), while subbands d LH (6, n 1 , n 2 ), d HH (6, n 1 , n 2 ), d HH (6, n 1 , n 2 ), and d LH (6, n 1 , n 2 ) contain no information beyond noise.This is because the reconstruction from a limitedangle sinogram can only encode information about certain directions of the wavefront set, as explained in Section 2.2.

Morphological opening of the complex wavelet coefficients
Due to stretching artefacts and imperfect measurements, the complex wavelet coefficients available in subbands d HL (6, n 1 , n 2 ) and d HL (6, n 1 , n 2 ) are noisy.Therefore, we improve this information of the wavefront set using two nonlinear steps: grayscale morphology followed by learning.
To extract useful information about the complex wavelet coefficients, we use the morphological opening operation on the absolute values of the complex wavelet coefficients.First, we set subbands d LH (6, n 1 , n 2 ), d HH (6, n 1 , n 2 ), d HH (6, n 1 , n 2 ), and d LH (6, n 1 , n 2 ) to zero.For the remaining two subbands d HL (6, n 1 , n 2 ) and d HL (6, n 1 , n 2 ), we perform morphological opening with oriented line structuring elements S ν,6 : where A ν,6 , defined in equation ( 29), is a grayscale image considered as a function A ν,6 : Z 2 → [0, 1].The line structuring elements S ν,6 are the binary images shown in Figure 9 and are considered as a function  ).The result is illustrated in Figure 8b.Now, one would think that hard-thresholding of the opened subbands is as follows: 1, A ν,6 (n 1 , n 2 ) ≥ , 0 otherwise, would give binary indicators for the discrete wavefront set.This would be the case if we had a good threshold > 0 available.However, in practice, imperfect knowledge of the noise amplitude and the absence of ground truth makes it difficult to design an automatic method for choosing the threshold.Therefore, we resort to learning the thresholding process.

Learning the binary mask
We use a convolutional neural network N 1 : R Z 2 ×I → [0, 1] Z 2 ×I to learn the binary indicators for the discrete wavefront set in the two subbands d HL (6, n 1 , n 2 ) and d HL (6, n 1 , n 2 ).The input of neural network N 1 is the tensor A = ( A ν,6 (n)) n∈Z 2 ,ν∈I obtained from equation (30) and the network maps it to The training of the neural network uses the morphologically opened complex wavelet coefficients A HL,6 (n 1 , n 2 ) and A HL,6 (n 1 , n 2 ).After the network N 1 performs the thresholding, it outputs a binary mask with ones indicating the location of the wavefront set.See Figure 8c for an example result.
Training inputs are generated as follows.We first generate 5000 two-dimensional phantoms of size 128 × 128 with ellipse shapes that have varying radius, tilt, and position.For each phantom, we consider 50 X-ray measurements over a 40-degree opening angle, spanning from 70 to 110 degrees.We add 5% noise to the measurement sinogram.Then, we compute the reconstruction of each phantom using the PDFP algorithm with complex wavelet regularization.Next, we compute the absolute values of the complex wavelet coefficients, A ν,6 (n 1 , n 2 ), of each reconstruction.We set subbands A LH,6 (n 1 , n 2 ), A HH,6 (n 1 , n 2 ), A HH,6 (n 1 , n 2 ), and A LH,6 (n 1 , n 2 ) to zero.We perform morphological opening of the subbands A HL,6 (n 1 , n 2 ) and A HL,6 (n 2 , n 2 ) by using formula (30).These morphologically opened sets of subbands A ν,6 serve as the training inputs.
The corresponding ground truth values are generated as follows.First, we compute the absolute value of the complex wavelet coefficients of the ground truth phantoms.We again set subbands d LH (6, n 1 , n 2 ), d HH (6, n 1 , n 2 ), d HH (6, n 1 , n 2 ), and d LH (1, n 1 , n 2 ) to zero.We perform morphological opening on subbands d HL (6, n 1 , n 2 ) and d HL (6, n 1 , n 2 ) using (30) to get A HL,6 (n 1 , n 2 ) and A HL,6 (n 1 , n 2 ).Then, we convert the data into a binary form by thresholding sufficiently large values of the morphologically opened coefficients to 1 and the rest are set to 0. That is, in training of neural network N 1 , the ground truth data is obtained using the formula using a suitably chosen > 0.
We implement neural network N 1 as a convolutional autoencoder with residual connections, as inspired by [27].For the encoder and decoder parts of the network, we use convolution kernels of size 3 × 3, stride Figure 11: Workflow of the proposed method: 1. Compute an initial reconstruction using the PDFP algorithm with complex wavelet regularization.2. Compute the finest scale complex wavelet coefficients for the six subbands and take their absolute value.3. Clean the coefficients using morphological opening with appropriately oriented line structuring elements.4. Give this to the first neural network, which thresholds the six subbands into a binary format.5. Compute an initial guess of the microlocal prior, by dilating the binary subbands with specific directed structuring elements.6. Give this initial guess for the second neural network, which outputs a prediction of the wavefront set in all six subbands.7. Combine the predicted wavefront set information in all six subbands to form the singular support.8. Compute the morphological skeleton of the singular support, estimating its boundary.9. Add this learned boundary estimate as an overlay on top of the PDFP reconstruction.

Learned microlocal prior
In the previous section, we explained how to extract the known part of the wavefront set from a limitedangle tomography reconstruction.In this section, we will go through the rest of the steps in order to learn a boundary estimate.See Figure 11 for a summary of the entire workflow of our proposed method.The basic idea behind our microlocal prior is based on the assumption of the smoothness and bounded curvature of the singular support of the X-ray attenuation coefficient.When a part of a wavefront set is detected in a subband, it is approximately a piece of straight line, following the singular support along a direction designated for that subband.Imagine traveling along the line.Eventually the singular support is bound to turn outside the orientation range of the subband you are on, as the interfaces between areas of different attenuation are closed curves.Then the line ends abruptly.But there is a continuation of it on one of the two neighboring subbands, depending on the direction of the turn.The representation of the singular support continues there as an almost straight line (along an orientation 30 degrees rotated from where you started), until another inevitable turn is forcing a jump to yet another subband.The structuring elements shown in Figure 12 are designed so as to capture the locations in a neighboring subband where the continuation can go.
We learn the boundary estimate with a second neural network N 2 , using the output of the first neural network N 1 .In order to do so, we first approximate the microlocal prior by morphological dilation.The details of the steps are explained in the sections below.

Approximation of the microlocal prior by morphological dilation
Now, we extend the wavefront set information from the two available subbands (see Figure 13a) to the remaining four subbands.We begin by using morphological dilation, as explained in section 2.4, to compute an initial guess for the location of nonzero coefficients.As the structuring elements for the operation, we have designed custom-made directional structuring elements that are able to estimate  the location of the wavefront set based on a neighbouring subband.The structuring elements S ν→ν , ν, ν ∈ I that we use for our specific limited-angle measurement geometry are shown in Figure 12.These structuring elements are designed so that the elements (a) and (b) in Figure 12 are the rotated and scaled images of the set and the elements (c) and (d) are the rotated and scaled images of the set The design is motivated by the assumption that the unknown object may have a discontinuity on a curve γ whose curvature is bounded.For example, the unknown object corresponding to function f may be a sum of the indicator function of a smooth domain D multiplied by a smooth function f 1 and a smooth function and Let us consider the case when the path γ is parametrized as γ(s) = (s, h(s)), s ∈ [−1, 1], where h : [−1, 1] → R is a function that satisfies h(0) = 0 and h (0) = 0.Moreover, assume that the normal vector of γ(x) turns clockwise when s > 0 and s grows.
Then for s > 0 the second derivative of h satisfies h (s) < 0. If we also have 0 > h (s) ≥ −1, then γ([0, 1]) ⊂ S − .When s < 0, we assume that the normal vector of γ(x) turns clockwise when −s grows and that 0 < h (s) ≤ 1, in which case γ([−1, 0]) ⊂ S − .Note that for example the structuring element S HL→HH corresponds to predicting the points x = (x 1 , x 2 ) on the discontinuity ∂D with a normal vector ν(x) of ∂D such that there exists a nearby point x = (x 1 , x2 ) ∈ ∂D with a normal vector ν(x) such that ν(x) is a vector close to ν(x) that has been turned in the clockwise direction (see Figure 7).The approximation of the microlocal prior is done as follows using the above structuring elements.We make no changes to the subbands HL and HL, since we have already extracted the wavefront set for these subbands in the previous steps.That is,  See Figure 13c for an example output.
To get the training inputs for network N 2 , we use the same phantoms as with network N 1 .We perform morphological dilation (as explained in section 4.1) on thresholded ground truth subbands A ν,6 (n 1 , n 2 ), see equation 31, to get the approximations for the wavefront set D ν,6 (n).These act as the training inputs for N 2 .
The corresponding ground truth values are computed from the ground truth phantoms, where we have the full wavefront set information available in all of the subbands.We simply take the absolute value, morphologically open, and threshold the data to a binary form.These serve as the ground truth values in the training.
For the neural network N 2 , we use a similar architecture as for network N 1 , however, here batch normalization is not used in the training (see Figure 14).Otherwise, the training scheme and use of loss function were identical to that of network N 1 .See section 3.3 for a more detailed explanation.Similarly, the training datasets were of size 5000 × 128 × 128 × 6, from which we used 500 samples for validation.

Boundary estimate
Now, we can combine the extended wavefront set D ν,6 (n 1 , n 2 ), that is the output of network N 2 , to form the singular support (see section 2.2).That is, we take the maximum of all six binary subbands: See Figure 15a for the full singular support.Then, we can compute the boundary estimate of the singular support using the morphological skeleton operation described in equation ( 23) of section 2.4.See Figure 15b.Finally, this learned boundary estimate can be used as an overlay for the PDFP reconstruction.See section 5 for final results.

Reconstructions in the xz-plane
For testing our proposed method, we created a three-dimensional phantom consisting of three L p -balls, where p = 1.5, in a constant background (see Figure 16).First, we simulated X-ray measurements of the phantom using a parallel-beam imaging geometry of a 40-degree opening angle with 50 projection images.We computed the PDFP reconstructions slice-by slice, separately for each xz-plane, treating each slice as an independent 2D tomographic reconstruction problem.The complex wavelet coefficients used for regularization were computed using Kingsbury Q-shift filters.Following the workflow explained in Figure 11, we learn the boundary estimate for each xz-slice.The learned boundary estimate is then overlayed on the PDFP reconstruction to indicate the extent of boundaries of the stretched features in the reconstruction.
In Figure 17, a reconstructed xz-slice with the learned boundary estimate is compared to the ground truth xz-slice of the test phantom.Also, the tomosynthesis reconstruction is shown for comparison.The results of additional xz-slices are shown in Figure 18.To measure how well the learned boundary estimate matches the true boundary, we have computed the dice similarity coefficient (DSC) of segmentations computed from the learned boundary estimate (X z ) and compared those to the ground truth slice (Y z ): where |X z |, |Y z | are the cardinalities, that is, the number of elements in the sets.The resulting DSC values are presented in table 1.

Conclusion
We introduced a new method for estimating the wavefront set of an X-ray coefficient function, given limited-angle tomography data.We found that detecting the wavefront set by thresholding complex wavelet coefficients was difficult to make automatic.However, we were able to train a neural network to find a good approximation.In addition, our idea of using geometric a priori information in learning to continue the visible part of the wavefront set to the invisible part was successful.We trained our method with simulated elliptical targets in 2D.Our data was severely limited with a 40-degree angle-of-view.We found that the method generalizes successfully to 2D problems arising from slices of simulated L p -ball targets in R 3 .Namely, the extent of objects in the direction of the z-axis was rather accurately recovered.However, the method did have difficulties near places where the boundaries of L p -balls had high curvature.We also noted that if some inclusions were much smaller than others (or of less contrast), the method might not detect them.
Our computational example suggests that our method can recover boundary curves with unprecedented accuracy for smooth inclusions.We expect that to be useful, for instance, in detecting bubbles when inspecting welds.However, for DBT, it remains unclear what happens with erratic boundaries of tumors.
The natural next step in this research direction is to fix an application and apply the method to measured data.The generalization to fan-beam geometry is straight-forward and in principle, there should be no issues in extending our method to 3D with cone-beam imaging geometry.

Figure 3 :
Figure 3: Frequency tilings of different multi-resolution transforms.(a) Real-valued wavelets are often understood to give horizontal, vertical, and diagonal details.This can be seen in the approximate supports of the basis functions in the frequency domain.(b)The approximate supports of the Fourier transforms of the frame functions corresponding to the six subbands in the dual-tree complex wavelet transform.(c) Shearlets offer very detailed analysis of orientations using parabolic scaling and increasing numbers of orientations with finer scales.However, the benefits come with a computational cost.

Figure 4 :
Figure 4: The first scale of complex wavelet detail coefficients of a phantom shown in (a).The orientations of coefficient subbands are denoted below each subband.Here, the absolute value of the complex-valued detail coefficients is shown.

Figure 5 :
Figure 5: Demonstration of the visible part of the wavefront set in limited-angle tomography.(a) Digital phantom, where attenuation coefficient f (x) = 0 in the black area and f (x) = 1 inside the white character "e".(b) Full singular support of f .(c) We simulated parallel-beam Radon transform data of f in the limited angular range −45 • ≤ θ ≤ 45 • , calculated the filtered back-projection reconstruction, and determined the singular support using dual-tree complex wavelets.Which parts of the jump curves are missing?Those whose tangents are not parallel to any X-rays in the data.(d) Same as (c) but with more limited angular range −20 • ≤ θ ≤ 20 • .
is the indicator function of the set B and d(x) = 1 D (x) is the indicator function of the set D, we have (d • b)(x) = 1 D•B (x).(28) Similar formulas hold for dilation d ⊕ b and erosion d b.

Figure 6 :
Figure 6: (a) Elliptical digital phantom with constant attenuation and zero background.(b) PDFP reconstruction with complex wavelet regularization from 50 X-ray projection images over a 40-degree opening angle.The central direction of projection is vertical.Note the vertical elongation or stretching of the ellipse in the reconstruction.This is a typical artefact in limited-angle tomography.

Figure 7 :
Figure7: The geometric idea behind our microlocal prior.Assume we know two elements of the wavefront set, located spatially at the two red dots and having the directions indicated by the two blue arrows.We work with piece-wise constant attenuation coefficients having reasonably regular boundary curves (shown as the black arc here) between the constant domains.Then we know that there should be another element in the wavefront set with a spatial location somewhere in the are shown as blurred red, and with direction within the range indicated with the blue triangle.

Figure 8 :
Figure 8: Going from (a) the imperfect and approximate wavefront set provided by DTCW coefficients into a robust wavefront set estimator.The technique involves two nonlinear steps: (b) morphology followed by (c) learning.The subband indexes from left to right are: LH, HH, HL, HL, HH, LH.

Figure 9 :
Figure9: Line structuring elements S ν,6 that are binary images, for each index ν of the oriented subbands, listed in(1).See (6)-(11) for the connection of the orientations of the line structuring elements to specific wavelet functions.

Figure 12 :
Figure 12: The custom-made directional structuring elements S ν→ν we use for our specific limited-angle measurement geometry.

Figure 13 :
Figure 13: Learning to fill in the incomplete wavefront set in six complex wavelet coefficient subbands, extending it to the entire domain.The subband indexes from left to right are: LH, HH, HL, HL, HH, LH.

D
HL,6 (n 1 , n 2 ), = A HL,6 (n 1 , n 2 ).(34) To get the wavefront set estimate for the adjacent subbands HH and HH, we dilate subbands HL and HL with the structuring elements S HL→HH and S HL→HH , respectively.That is, D HH,6 = D HL,6 ⊕ S HL→HH (35) D HH,6 = D HL,6 ⊕ S HL→HH .(36) Next, for subbands LH and LH, we treat the once dilated subbands HH and HH as the objects for the dilation with structuring elements S HH→LH and S HH→LH , respectively.That is, D LH,6 = D HH,6 ⊕ S HH→LH = D HL,6 ⊕ S HL→HH ⊕ S HH→LH (37) D LH,6 = D HH,6 ⊕ S HH→LH = D HL,6 ⊕ S HL→HH ⊕ S HH→LH .(38)Now,we have a coarse initial guess of the wavefront set for the network N 2 to refine.See Figure13bfor an example.

4. 2
Learning to extend the wavefront set into the entire domainWe want to learn to refine the wavefront set estimate with neural network N 2 , based on the approximation provided by the dilation operation.The network N 2 outputs a prediction of the full wavefront set D ν,6 (n 1 , n 2 ), extending it to all six subbands.We use a convolutional neural network N 2 : R Z 2 ×I → [0, 1] Z 2 ×I , where the input of N 2 is the tensor D = (D ν,6 (n)) n∈Z 2 ,ν∈I and the network maps it to D = ( D ν,6 (n)) n∈Z 2 ,ν∈I .That is, D = N 2 (D).

For a
three-dimensional phantom, after independently computing PDFP reconstructions and the boundary estimates for each of the xz-slices separately (see section 5.1), we can stack the reconstruction results back together to form a volume.Then, we can slice the reconstructed volume in the xy-direction.We present the results as xy-slices with a selection of different z-values.Reconstruction results of several xy-slices of the phantom are shown in Figure 19.xz-slice (y) DSC 22 0

Figure 17 :
Figure 17: Comparisons of the xz-slice, where y = 57, of (a) the test phantom, (b) the tomosynthesis reconstruction, and (c) the PDFP reconstruction with complex wavelet regularization.On the bottom row, the learned boundary estimate from the PDFP reconstruction is added to each image as an overlay for comparison purposes.Imaging geometry: parallel beam, 40-degree opening angle (70-110), 50 projections.

Figure 18 :
Figure 18: Different xz-slices through the three-dimensional test phantom shown in Figure 16.In the model, xz-slices are vertical (perpendicular to the detector surface).In our simplified assumption of parallel-beam geometry, there is an independent limited-angle tomography problem restricted to each xz-slice.(a) Ground truth slice.(b) Tomosynthesis reconstruction.(c) PDFP reconstruction with learned boundary estimate.(d) Boundary estimate curve.Compare the boundary curves to the true boundaries of features shown in column (a).

Table 1 :
Dice simila18.y coefficients for segmentations computed from the learned boundary estimate and compared to the ground truth phantom in the xz-plane.The chosen slices correspond to those shown inFigures 17,18.