# Exploiting the Transferability of Deep Learning Systems Across Multi-modal
Retinal Scans for Extracting Retinopathy Lesions
††thanks: This work is supported by a research fund from Khalifa University:
Ref: CIRA-2019-047.

Taimur Hassan† ⋄, Muhammad Usman Akram⋄, Naoufel Werghi† †Department of
Electrical Engineering and Computer Sciences, Khalifa University, Abu Dhabi,
United Arab Emirates.
⋄Department of Computer and Software Engineering, National University of
Sciences and Technology, Islamabad, Pakistan.

###### Abstract

Retinal lesions play a vital role in the accurate classification of retinal
abnormalities. Many researchers have proposed deep lesion-aware screening
systems that analyze and grade the progression of retinopathy. However, to the
best of our knowledge, no literature exploits the tendency of these systems to
generalize across multiple scanner specifications and multi-modal imagery.
Towards this end, this paper presents a detailed evaluation of semantic
segmentation, scene parsing and hybrid deep learning systems for extracting
the retinal lesions such as intra-retinal fluid, sub-retinal fluid, hard
exudates, drusen, and other chorioretinal anomalies from fused fundus and
optical coherence tomography (OCT) imagery. Furthermore, we present a novel
strategy exploiting the transferability of these models across multiple
retinal scanner specifications. A total of 363 fundus and 173,915 OCT scans
from seven publicly available datasets were used in this research (from which
297 fundus and 59,593 OCT scans were used for testing purposes). Overall, a
hybrid retinal analysis and grading network (RAGNet), backboned through
ResNet50, stood first for extracting the retinal lesions, achieving a mean
dice coefficient score of 0.822. Moreover, the complete source code and its
documentation are released at http://biomisa.org/index.php/downloads/.

###### Index Terms:

Retinal Lesions, Ophthalmology, Convolutional Neural Networks, Fundus
Photography, Optical Coherence Tomography.

##  I Introduction

Retinopathy or retinal diseases tend to damage the retina, which may result in
a non-recoverable loss of vision or even blindness if not timely treated. Most
of these diseases are associated with diabetes. However, they may also occur
due to aging, uveitis, and cataract surgeries. The two common retinal diseases
are macular edema (ME) and age-related macular degeneration (AMD). ME occurs
due to fluid accumulation within the macula mostly due to the associated
hyperglycemia, uveitis, and cataract surgeries. ME caused by diabetes is often
termed as diabetic macular edema (DME) which is identified by examining the
patient’s diabetic history and also by analyzing the retinal thickening
(caused due to retinal fluid) or hard exudates (HE) within the one-disc
diameter of the center of the macula (containing a small pit known as fovea)
[1].

Figure 1:  Retinal lesions in fundus and OCT scans of Ci-CSME (A, B) and dry
AMD pathology (C, D).

Early Treatment Diabetic Retinopathy Studies (ETDRS) classified clinically
significant ME as having: 1) either the thickening within 500$\mathrm{\mu}$m
of the macular center; 2) HE along with the adjacent thickening within the
macular center of 500$\mathrm{\mu}$m; or 3) retinal thickening regions of one
(or more) disc diameter in which some part of them are within the one-disc
diameter of the macular center [2]. But, with the advent of new imaging
techniques such as optical coherence tomography (OCT), the classification of
DME is redefined as centrally involved clinical significant macular edema (Ci-
CSME) if the presence of retinal thickening, due to retinal fluid or hard
exudates, is discovered within the central sub-field zone of the macula
(having 1mm or greater diameter). Otherwise, DME is classified as non-
centrally involved [3]. AMD is another retinal syndrome mostly found in elder
people. It is typically classified into two stages i.e. non-neovascular AMD
and the neovascular AMD. Non-neovascular AMD is the “dry” form of AMD in which
small, medium or large-sized drusen can be observed. With the increasing
disease severity, abnormal blood vessels intervene retina leading to
chorioretinal anomalies such as fibrotic scars and choroidal neovascular
membranes (CNVM). In such a case, AMD is classified as wet or neovascular AMD.
Fig. 1 shows some of the fundus and OCT scans showing retinal lesions at
different stages of AMD and DME.

##  II Related Work

In the literature, a large body of solutions assessing retinal regions
employed features extraction techniques coupled with classical machine
learning (ML) tools. The majority of these methods are validated on a limited
number of scans and thus exhibited feeble reproducibility. More recently, with
the advent of deep learning, a wide variety of end-to-end approaches,
operating on more massive datasets, have been proposed.

Traditional Approaches: Fundus imagery has been the modality of choice for
examining the retinal pathology for a while [5] and is still used as a
secondary examination technique in analyzing the complex retinal pathologies.
But, with the advent of OCT, most of the solutions for retinal image analysis
have migrated towards this new modality due to its ability to present
objective visualization of retinal abnormalities in early stages. Chiu et al.
[6] developed a kernel regression with graph theory and dynamic programming
(KR+GTDP) scheme to extract retinal layers and retinal fluid from DME affected
scans [6]. In [7], a Random Forest-based framework was proposed for the
automated extraction of retinal layers and fluid from scans affected by
central serous retinopathy (CSR). Wilkins et al. [8] presented an automated
method for the extraction of intra-retinal cystoid fluid using OCT images.
Vidal et al. [9] used a linear discriminant classifier, support vector
machines, and a Parzen window for the identification of intra-retinal fluid
(IRF). Apart from this, we have also proposed several methods for extracting
retinal layers, retinal fluid, and for classifying retinopathy using
traditional ML techniques [10, 11, 12, 13, 14].

Deep Learning Methods: Many researchers have applied deep learning for the
extraction of retinal layers [15] and retinal lesions such as IRF [16], sub-
retinal fluid (SRF) [17] and HE [21]. Seebock et al. [19] proposed a Bayesian
UNet based framework for recognizing different anomalies within the retinal
pathologies. Fang et al. [20] developed a lesion-aware convolutional neural
network (LACNN) model for the accurate classification of DME, choroidal
neovascularization, drusen (AMD) and normal pathologies. LACNN is composed of
a lesion detection network (LDN) and lesion-attention module where LDN first
generates a soft attention map to weight the lesion-aware features extracted
from the lesion-attention module and then these features are used for the
accurate classification of retinal pathologies. Apart from this, we have
recently proposed a hybrid retinal analysis and grading architecture (RAGNet)
[21] that utilizes a single feature extraction model for retinal lesions
segmentation, lesion-aware classification and severity grading of retinopathy
based on OCT images.

##  III Contributions

In this paper, we present a thorough evaluation of deep learning models for
the extraction of IRF, SRF, HE, drusen, and other chorioretinal anomalies such
as fibrotic scars and CNVM from multi-modal retinal images. Furthermore, we
exploited the transferability of these models for retinal lesions extraction
across multi-modal imagery. To the best of our knowledge, there is no
literature available to date providing a thorough transferability analysis of
encoder-decoder, fully convolutional, scene parsing, and hybrid deep learning
systems for extracting this multitude of lesions in one go from multi-modal
retinal imagery. Subsequently, the main contributions of this paper are:

  * •

A first comprehensive evaluation of semantic segmentation, scene parsing and
hybrid deep learning systems such as RAGNet [21], PSPNet [22], SegNet [23],
UNet [24], and FCN (8s and 32s) [25] for extracting multiple lesions from
multi-modal retinal imagery.

  * •

A comprehensive study encompassing seven publicly available datasets, and five
different retinal pathologies represented in a total of 363 fundus and 173,915
OCT scans from which 297 fundus and 59,593 OCT scans were used for testing
purposes.

  * •

A detailed exploration of the transferability of these models across multiple
scanner specifications.

##  IV Proposed Approach

We propose a novel study to analyze the transferability of the state-of-the-
art deep learning frameworks across fused fundus and OCT imagery for
extracting multiple retinal lesions in one go. The models which we considered
are as follows:

RAGNet: is a hybrid convolutional network that can perform pixel-level
segmentation and scan-level classification at the same time [21]. The
uniqueness in the RAGNet architecture is that it uses the same feature
extractor for the classification and segmentation purposes. So, if the problem
demands segmentation and classification from the same image based upon similar
features, then RAGNet would be an ideal choice rather than using two separate
models [21]. Here, we have only used RAGNet segmentation unit since we are
focusing on the retinal lesions segmentation.

PSPNet: is a state-of-the-art scene parsing network that contains a pyramid
pooling module to generate four pyramids of feature maps representing coarser
to finer details to minimize the loss of global scene context while generating
the latent representations [22]. The pooled outputs are then concatenated with
the original feature maps to generate the final segmentation results.

SegNet: is an encoder-decoder network for semantic segmentation. The
uniqueness in the SegNet model is that it uses pooling indices from the
corresponding encoder block to up-sample the feature maps at the decoder end
in a non-linear fashion. Afterward, the feature maps are convolved with
trainable filters to remove their sparsity. Moreover, SegNet has a smaller
number of trainable parameters due to which it is computationally more
efficient.

UNet: is an auto-encoder inspired by FCN for semantic segmentation. The key
feature of UNet is that it is fast and can generate good segmentation results
with a small number of training samples because of its in-built data
augmentation strategy [24]. UNet uses up-sampling instead of pooling
operations and generates a large number of feature maps to mitigate the
contextual information to the higher resolution layers [24].

FCN: is an end-to-end model proposed for semantic segmentation. FCN uses
learned representation from the pre-trained models, fine-tune them, and
generates finer pixel-level predictions in one go based upon up-sampling lower
network layers with finer strides. In this study, we have utilized FCN-8 and
FCN-32 (i.e. the finest and the coarsest version of FCN) for retinal lesions
extraction.

We have applied these models for extracting retinal lesions from both fundus
and OCT imagery. Here, we note that our study covers some of the most complex
and commonly occurring retinal pathologies including non-neovascular AMD,
neovascular AMD, Ci-CSME, and non-Ci-CSME. We also note that the related scans
were collected using machines from different manufacturers and exhibit varying
scan quality. To make the comparison objective and highly reproducible, we
have used publicly available datasets in our investigations. Furthermore, we
have tested the transferability of these models through an extensive cross-
dataset validation. The series of experiments we conducted in this work
provide a reliable benchmark for assessing the robustness and generalization
capacity of each model.

##  V Experimental Setup

This section reports the detailed description of the datasets which have been
used in this research. Furthermore, it contains implementation details and the
performance metrics on which the models are evaluated:

###  V-A Datasets

We have evaluated all the models on seven publicly available retinal image
datasets where the ground truths for retinal lesions were acquired through the
Armed Forces Institute of Ophthalmology, Rawalpindi Pakistan. The detailed
summary of each dataset is presented below:

Rabbani-I [26] is one of the few datasets which contains both OCT and fundus
images of each subject reflecting AMD, DME, and normal pathologies. The
dataset is acquired at Noor Eye Hospital, Tehran Iran and contains a total of
4,241 OCT and 148 fundus scans from 50 normal, 48 dry AMD, and 50 DME affected
subjects. In this paper, we considered 37 fundus scans and 1,061 OCT scans for
training and the rest for testing purposes.

Rabbani-II [27] contains 12,800 OCT and 100 color fundus scans from both eyes
of 50 healthy subjects. Since Rabbani-II only contain scans from the healthy
subjects so it served as an excellent benchmark to test the false positive
rate for all the models (indicating how many false positives each model
generates).

Duke-I [28] is one of the oldest retinal OCT datasets containing a total of
38,400 scans from which 26,900 scans reflect dry AMD and 11,500 scans show
controlled (healthy) pathology. In the proposed study, a total of 300 scans
have been used for training and 38,100 scans have been used for testing
purposes.

Duke-II [6] has a total of 610 OCT images from 10 severe DME affected
subjects. Moreover, the dataset contains highly detailed markings for retinal
layers and fluid from two clinicians. In this paper, a total of 305 scans were
used for training from the first five subjects, and the rest for testing.

Duke-III [29] is another dataset from Duke University which we used in our
research. The dataset contains 723 scans reflecting dry AMD pathologies, 1,101
scans reflecting DME pathology, and 1,407 scans showing normal retinal
pathology. For the experimentations, we considered 3,048 scans for training
and the rest 183 for testing.

BIOMISA dataset [30] contains a total of 5,324 OCT (657 dry AMD, 2,195 ME, 904
normal, 407 wet AMD, and 1,161 CSR) and 115 fundus scans from 99 subjects (17
healthy, 31 ME, 8 dry AMD, 19 wet AMD, 24 CSR). In this study, a total of
1,299 OCT and 29 fundus images from the BIOMISA dataset were used for training
and the rest for the evaluation purposes.

Zhang dataset [31] contains 109,309 scans representing wet AMD (CNV), dry AMD
(Drusen), DME, and healthy pathologies. The dataset also presents a clear
separation of 108,309 scans for training while 1,000 scans for testing
purposes which we followed in the experimentations as well.

###  V-B Implementation Details

All the models have been implemented using Keras, Python 3.7.4 on a machine
having Intel 8th generation Core i5, NVIDIA RTX 2080 GPU and 16 GB RAM where
ResNet50 was used as a backbone. Moreover, the optimizer used during the
training was an adaptive learning rate method (ADADELTA) with a default
learning and decay rate. The source code has been released at
http://biomisa.org/index.php/downloads/ for reproducibility.

###  V-C Evaluation Metrics

In the proposed study, all the segmentation models have been evaluated using
the following metrics:

Mean Dice Coefficient: Dice coefficient ($\mathrm{D_{C}}$) computes the degree
of similarity between the ground truth and the extracted results using
following relation: $(\mathrm{D_{C}=\frac{2T_{P}}{2T_{P}+F_{P}+F_{N}}})$,
where $\mathrm{T_{P}}$ indicates the true positives, $\mathrm{F_{P}}$
indicates the false positives and $\mathrm{F_{N}}$ indicates the false
negatives. After computing $\mathrm{D_{C}}$ for each lesion class. The mean
dice coefficient is computed for each network by taking an average of their
$\mathrm{D_{C}}$ scores.

Mean Intersection-over-Union: The mean intersection-over-union (IoU) is
computed by taking an average of IoU scores for each lesion class where the
IoU scores are computed through
$(\mathrm{IoU=\frac{T_{P}}{T_{P}+F_{P}+F_{N}}})$.

Recall, Precision and F-score: To further evaluate the models, we computed
pixel-level recall $\mathrm{(T_{PR}=\frac{T_{P}}{T_{P}+F_{N}}})$, precision
$\mathrm{(P_{PV}=\frac{T_{P}}{T_{P}+F_{P}}})$ and F-score
$\mathrm{(F_{1}=\frac{2xT_{PR}xP_{PV}}{T_{PR}+P_{PV}}})$.

Qualitative Evaluations: The performances of all the models for lesions
extraction have been also qualitatively evaluated through some visual
examples.

##  VI Results and Discussion

The evaluation of segmentation models has been conducted on the combination of
all seven datasets containing mixed OCT and fundus scans. In terms of
$\mathrm{T_{PR}}$ and $\mathrm{F_{1}}$ as shown in Table I, RAGNet achieves
9.48% and 3.36% improvements as compared to UNet and PSPNet, respectively.
However, in terms of precision, SegNet has a lead of 1.52% as compared to
PSPNet. This indicates that SegNet produces fewer false positives as compared
to the rest of the models. For pixel-level comparison, we have excluded
accuracy because it gives biased results towards a dominant-negative class
i.e. the background.

Table I: Performance evaluations in terms of pixel-level recall, precision and F-scores on combined dataset. Bold indicate the best performance. Network | $\mathrm{T_{PR}}$ | $\mathrm{P_{PV}}$ | $\mathrm{F_{1}}$
---|---|---|---
RAGNet | 0.8547 | 0.8606 | 0.8576
PSPNet | 0.7540 | 0.9200 | 0.8287
SegNet | 0.6388 | 0.9342 | 0.7587
UNet | 0.7736 | 0.8842 | 0.8252
FCN-8 | 0.6238 | 0.6165 | 0.6201
FCN-32 | 0.4755 | 0.5611 | 0.5147

Tables 2 and 3 reports the performance of all the models for extracting
retinal lesions in terms of mean $\mathrm{D_{C}}$ and mean IoU, respectively.
From Table 2, it can be observed that RAGNet achieves the best mean
$\mathrm{D_{C}}$ score of 0.822, leading PSPNet by 4.5% and FCN-32 by 51.33%.
Moreover, in terms of mean IoU, RAGNet also achieves the overall best
performance (mean IoU: 0.710) showing a neat gap over its competitors for
extracting IRF, SRF, and HE regions. In Table 3, the second-best performance
is achieved by PSPNet that lags from RAGNet by 6.9%. Also, we noticed that on
fundus images UNet achieves optimal lesion extraction results with an overall
performance comparable to that of PSPNet. Fig. 2 shows the qualitative results
of all the models when trained on multi-modal images from all seven datasets
at once, where we can notice the best overall performance of RAG-Net. It
should be noticed that extracting lesions accurately from both modalities at
once is quite challenging as their image features vary a lot.

Figure 2:  Comparison of retinal lesions extraction on combined dataset. (From left to right: original image, ground truth, RAGNet, PSPNet, UNet, FCN-8, FCN-32, SegNet). Blue, red, yellow, green, and pink indicate HE, IRF, SRF, CA, and drusen, respectively. Table II: Performance evaluations of deep segmentation models for retinal lesions extraction in terms of $\mathrm{D_{C}}$. Bold indicates the overall best performance. Network | IRF | SRF | CA | HE | Drusen | Mean
---|---|---|---|---|---|---
RAGNet | 0.846 | 0.850 | 0.941 | 0.633 | 0.840 | 0.822
SegNet | 0.810 | 0.610 | 0.886 | 0.373 | 0.695 | 0.675
PSPNet | 0.843 | 0.809 | 0.944 | 0.594 | 0.735 | 0.785
UNet | 0.816 | 0.757 | 0.878 | 0.581 | 0.864 | 0.779
FCN-8 | 0.681 | 0.568 | 0.761 | 0.124 | 0.410 | 0.509
FCN-32 | 0.651 | 0.434 | 0.638 | 0.032 | 0.243 | 0.400
Table III: Performance evaluations of deep segmentation models for retinal lesions extraction in terms of IoU. Bold indicates the overall best performance. Network | IRF | SRF | CA | HE | Drusen | Mean
---|---|---|---|---|---|---
RAGNet | 0.733 | 0.739 | 0.890 | 0.464 | 0.725 | 0.710
SegNet | 0.681 | 0.439 | 0.796 | 0.229 | 0.533 | 0.535
PSPNet | 0.728 | 0.680 | 0.895 | 0.423 | 0.581 | 0.661
UNet | 0.689 | 0.609 | 0.783 | 0.409 | 0.761 | 0.650
FCN-8 | 0.517 | 0.397 | 0.615 | 0.066 | 0.257 | 0.370
FCN-32 | 0.482 | 0.277 | 0.468 | 0.016 | 0.138 | 0.276
Table IV: Transferability analysis (Training $\rightarrow$ Testing) for all models in terms of mean IoU. Bold and blue indicates the first and second-best performance, respectively. (Datasets name are coded as follows: R: Rabbani, D: Duke, Z: Zhang and B: BIOMISA). The rest of the abbreviations are RN: RAGNet, PN: PSPNet, SN: SegNet, UN: UNet, F-8: FCN-8, F-32: FCN-32. | RN | PN | SN | UN | F-8 | F-32
---|---|---|---|---|---|---
R $\rightarrow$ D | 0.624 | 0.589 | 0.414 | 0.574 | 0.281 | 0.170
D $\rightarrow$ R | 0.649 | 0.601 | 0.426 | 0.612 | 0.301 | 0.194
R $\rightarrow$ Z | 0.657 | 0.615 | 0.468 | 0.604 | 0.322 | 0.225
Z $\rightarrow$ R | 0.663 | 0.632 | 0.472 | 0.629 | 0.329 | 0.245
B $\rightarrow$ R | 0.573 | 0.542 | 0.389 | 0.534 | 0.236 | 0.144
R $\rightarrow$ B | 0.554 | 0.535 | 0.374 | 0.521 | 0.213 | 0.115
Z $\rightarrow$ D | 0.809 | 0.752 | 0.613 | 0.734 | 0.476 | 0.342
D $\rightarrow$ Z | 0.794 | 0.741 | 0.598 | 0.721 | 0.445 | 0.335
D $\rightarrow$ B | 0.582 | 0.534 | 0.487 | 0.525 | 0.229 | 0.136
B $\rightarrow$ D | 0.571 | 0.529 | 0.479 | 0.517 | 0.214 | 0.127
B $\rightarrow$ Z | 0.564 | 0.513 | 0.465 | 0.524 | 0.221 | 0.152
Z $\rightarrow$ B | 0.552 | 0.507 | 0.457 | 0.512 | 0.235 | 0.178

In the second series of experiments, we have conducted a transferability
analysis to assess the generalization capabilities of all models. Here, we
combined Duke-I, II, and III as one dataset (i.e. Duke) and Rabbani-I and
Rabbani-II dataset as Rabbani to avoid redundant combinations as they have
similar image features. We report the results in Table 4 where it can be
observed that all the methods have shown good performance for Duke and Zhang
dataset pairs and this is natural because both datasets are acquired through
Spectralis, Heidelberg Inc. Moreover, RAGNet achieved the overall best
performance as evident from Table 4, whereas PSPNet stood 2nd best but its
performance is comparable with UNet. In another experiment, we have used
Rabbani-II dataset to test how many false positive does each model produce.
Since Rabbani-II contains only healthy scans, so there are no actual lesions
in this dataset. The best performance in this experiment is achieved by the
RAGNet with a true negative $\mathrm{(T_{N})}$ rate of 0.9999 indicating that
it produces a minimum number of false lesions. Apart from this, the worse
performance is achieved for FCN-32 ($\mathrm{T_{N}}$ rate: 0.9379). The worse
performance of FCN-32 is even above 90% because the ratio of $\mathrm{T_{N}}$
pixels and the $\mathrm{F_{P}}$ pixels is extremely high. The results for this
experiment are available in the codebase package for the readers at
http://biomisa.org/index.php/downloads/.

##  VII Conclusion and Future Research

In this paper, we presented a thorough evaluation of semantic segmentation,
scene parsing, and hybrid deep learning systems for extracting retinal lesions
from fused fundus and OCT imagery. We also assessed the generalization
capacity of each model through comprehensive cross-data validations where
RAGNet, due to its robustness in retaining lesion contextual information
during scan decomposition, produces superior results as compared to other
models. Furthermore, the benchmarking performed in this work will be of great
utility for both researchers and practitioners who want to employ deep
learning models for lesion-aware grading of the retina. In the future, we plan
to extend and exploit this study for the extraction of the optic disc and
retinal layers in the optic nerve head region for the glaucoma analysis.

## References

  * [1] G. M. Comers, “Cystoid macular edema,” in Kellog Eye Center.  Accessed: June 2019.
  * [2] “Diabetic macular edema,” in EyeWiki.  Accessed: November 4th, 2019.
  * [3] N. Relhan et al., “The early treatment diabetic retinopathy study historical review and relevance to today’s management of diabetic macular edema,” in Current Opinion in Ophthalmology.  Wolters Kluwer, May 2017.
  * [4] M. U. Akram et al., “An automated system for the grading of diabetic maculopathy in fundus images,” in 19th International Conference on Neural Information Processing.  November 12th-15th, 2012.
  * [5] T. Hassan et al., “Review of oct and fundus images for detection of macular edema,” in IEEE International Conference on Imaging Systems and Techniques (IST).  September, 2015.
  * [6] S. J. Chiu et al., “Kernel regression based segmentation of optical coherence tomography images with diabetic macular edema,” in Biomedical Optics Express.  Vol. 6, No. 4, April 2015.
  * [7] D. Xiang et al., “Automatic retinal layer segmentation of oct images with central serous retinopathy,” in IEEE Journal of Biomedical and Health Informatics.  Vol 23, No. 1, January 2019.
  * [8] G. R. Wilkins et al., “Automated segmentation of intraretinal cystoid fluid in optical coherence tomography,” in IEEE Transactions on Biomedical Engineering.  pp. 1109-1114, 2012.
  * [9] P. L. Vidal et al., “Intraretinal fluid identification via enhanced maps using optical coherence tomography images,” in Biomedical Optics Express.  October 2018.
  * [10] S. Khalid et al., “Automated segmentation and quantification of drusen in fundus and optical coherence tomography images for detection of armd,” in Journal of Digital Imaging.  December 2017.
  * [11] S. Khalid et al., “Fully automated robust system to detect retinal edema, central serous chorioretinopathy, and age related macular degeneration from optical coherence tomography images,” in BioMed Research International.  March 2017.
  * [12] T. Hassan et al., “Automated segmentation of subretinal layers for the detection of macular edema,” in Applied Optics.  55, 454-461, 2016.
  * [13] B. Hassan et al., “Structure tensor based automated detection of macular edema and central serous retinopathy using optical coherence tomography images,” in Journal of Optical Society of America A.  33, 455-463, 2016.
  * [14] A. M. Syed et al., “Automated diagnosis of macular edema and central serous retinopathy through robust reconstruction of 3D retinal surfaces,” in Computer Methods and Programs in Biomedicine.  137, 1-10, 2016.
  * [15] L. Fang et al., “Automatic segmentation of nine retinal layer boundaries in oct images of non-exudative amd patients using deep learning and graph search,” in Biomedical Optics Express.  Vol. 8, No. 5, May 2017.
  * [16] A. G. Roy et al., “ReLayNet: retinal layer and fluid segmentation of macular optical coherence tomography using fully convolutional networks,” in Biomedical Optics Express.  Vol. 8, No. 8, 1 August 2017.
  * [17] T. Schlegl et al., “Fully automated detection and quantification of macular fluid in oct using deep learning,” in Ophthalmology.  Vol. 125, No. 4, April 2018.
  * [18] B. Hassan et al., “Deep ensemble learning based objective grading of macular edema by extracting clinically significant findings from fused retinal imaging modalities,” in MDPI Sensors.  July 2019.
  * [19] P. Seebock et al., “Exploiting epistemic uncertainty of anatomy segmentation for anomaly detection in retinal oct,” in IEEE Transactions on Medical Imaging.  May 2019.
  * [20] L. Fang et al., “Attention to lesion: Lesion-aware convolutional neural network for retinal optical coherence tomography image classification,” in IEEE Transactions on Medical Imaging.  August 2019.
  * [21] T. Hassan et al., “RAG-FW: A hybrid convolutional framework for the automated extraction of retinal lesions and lesion-influenced grading of human retinal pathology,” in IEEE Journal of Biomedical and Health Informatics.  March 2020.
  * [22] H. Zhao et al., “Pyramid scene parsing network,” in IEEE CVPR.  2017\.
  * [23] V. Badrinarayanan et al., “Segnet: A deep convolutional encoder-decoder architecture for image segmentation,” in IEEE Transactions on Pattern Analysis and Machine Intelligence.  December 2017.
  * [24] O. Ronneberger et al., “U-net: Convolutional networks for biomedical image segmentation,” in MICCAI.  2015\.
  * [25] J. Long et al., “Fully convolutional networks for semantic segmentation,” in IEEE CVPR.  2015\.
  * [26] R. Rasti et al., “Macular oct classification using a multi-scale convolutional neural network ensemble,” in IEEE Transactions on Medical Imaging.  vol. 37, no. 4, pp. 1024-1034, April 2018.
  * [27] T. Mahmudi et al., “Comparison of macular octs in right andleft eyes of normal people,”  in Proc. SPIE, Medical Imaging, San Diego, California, United States Feb. 15-20, 2014.
  * [28] S. Farsiu et al., “Quantitative classification of eyes with and without intermediate age-related macular degeneration using optical coherence tomography,” in Ophthalmology.  121(1), 162-172 January 2014.
  * [29] P. P. Srinivasan et al., “Fully automated detection of diabetic macular edema and dry age-related macular degeneration from optical coherence tomography images,” in Biomedical Optics Express.  Vol. 5, No. 10 — DOI:10.1364/BOE.5.0035 68, 12 Sep 2014.
  * [30] T. Hassan et al., “BIOMISA Retinal Image Database for Macular and Ocular Syndromes,” in ICIAR-2018.  Portugal, June 2018.
  * [31] D. Kermany et al., “Identifying medical diagnoses and treatable diseases by image-based deep learning,” Cell.  172(5):1122-1131, 2018.


# Learning a microlocal prior
for limited-angle tomography

Siiri Rautio  Department of Mathematics and Statistics, University of Helsinki
Rashmi Murthy  Department of Mathematics and Statistics, University of
Helsinki  Tatiana A. Bubba  Department of Mathematical Sciences, University of
Bath
Matti Lassas  Department of Mathematics and Statistics, University of Helsinki
Samuli Siltanen  Department of Mathematics and Statistics, University of
Helsinki

###### Abstract

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.

##  1 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:\mathbb{R}^{n}\rightarrow\mathbb{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).

(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 $\pm 20^{\circ}$ 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.

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 post-processing 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 data, which can be
impractical in certain applications, for instance in medical imaging. Also,
the post-processing 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 proximal-dual 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.

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 parallel-beam
measurements with an X-ray source. However, the mathematics of wavefront sets
is essentially the same, so our findings should apply with straightforward
modifications.

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 $\pm
15^{\circ}$, $\pm 45^{\circ}$ and $\pm 75^{\circ}$. DTCWT does not offer such
a theoretically complete discretisation 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.

##  2 Mathematical background

###  2.1 Dual-tree complex wavelet transform

(a) Real-valued wavelets HLHLLHLHHHHHHHHH(b) Complex-valued
wavelets$\overline{H}L$$\overline{H}H$$\overline{L}H$$LH$$LH$$\overline{L}H$$HH$$\overline{H}H$$HH$$HL$$\overline{H}L$$HL$(c)
Shearlets 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.

The dual-tree complex wavelet transform consists of a complex-valued scaling
function $\phi$ and a complex-valued wavelet $\psi$. Consider a complex and
analytic wavelet $\psi(x)$, given by

$\psi(x)=\psi_{h}(x)+i\psi_{g}(x)$

and associated with a high-pass filter $H$, and a complex and analytic scaling
function $\phi(x)$, that is given by

$\phi(x)=\phi_{h}(x)+i\phi_{g}(x)$

and associated with a low-pass filter $L$. Note that $\psi_{h}$ and $\phi_{h}$
are the real parts of $\psi(x)$ and $\phi(x)$, and $\psi_{g}$ and $\phi_{g}$
are the imaginary parts of $\psi(x)$ and $\phi(x)$, respectively. For
computational purposes, we use finitely supported wavelets and thus
approximately analytic wavelets.

We denote the complex wavelet coefficients by $d_{\nu}(j,n)\in{\mathbb{C}}$,
where $j=0,\dots,J$ represents the number of scales used in the complex
wavelet representation and

$\nu\in\mathcal{I}=\\{\overline{L}H,\overline{H}H,\overline{H}L,HL,HH,LH\\}$
(1)

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 $\nu$.

To build intuition, the expression for the 2D wavelet in subband $HH$ is
obtained by:

$\displaystyle\psi(x_{1},x_{2})$ $\displaystyle=$
$\displaystyle[\psi_{h}(x_{1})+i\psi_{g}(x_{2})][\psi_{h}(x_{1})+i\psi_{g}(x_{2})]$
(2) $\displaystyle=$
$\displaystyle\psi_{h}(x_{1})\psi_{h}(x_{2})-\psi_{g}(x_{1})\psi_{g}(x_{2})+i[\psi_{g}(x_{1})\psi_{h}(x_{2})+\psi_{h}(x_{1})\psi_{g}(x_{2})].$
(3)

Note that this wavelet is oriented at $-45^{\circ}$. To obtain the 2D wavelet
oriented at $+45^{\circ}$, that is the subband $\overline{H}H$, we compute
$\psi_{2}(x_{1},x_{2})=\overline{\psi(x_{1})}\psi(x_{2})$, where

$\displaystyle\psi_{2}(x_{1},x_{2})$ $\displaystyle=$
$\displaystyle[\overline{\psi_{h}(x_{1})+i\psi_{g}(x_{1})}][\psi_{h}(x_{2})+i\psi_{g}(x_{2})]$
(4) $\displaystyle=$
$\displaystyle\psi_{h}(x_{1})\psi_{h}(x_{2})+\psi_{g}(x_{1})\psi_{g}(x_{2})+i[\psi_{h}(x_{1})\psi_{g}(x_{2})-\psi_{g}(x_{1})\psi_{h}(x_{2})].$
(5)

In a similar manner, the six subbands of the 2D complex wavelets for each
scale $j$ are computed as follows:

$\displaystyle\Psi_{LH,j,n}(x_{1},x_{2})$ $\displaystyle:=$
$\displaystyle\phi(2^{j}x_{1}-n_{1})\psi(2^{j}x_{2}-n_{2})$ (6)
$\displaystyle\Psi_{HL,j,n}(x_{1},x_{2})$ $\displaystyle:=$
$\displaystyle\psi(2^{j}x_{1}-n_{1})\phi(2^{j}x_{2}-n_{2})$ (7)
$\displaystyle\Psi_{HH,j,n}(x_{1},x_{2})$ $\displaystyle:=$
$\displaystyle\psi(2^{j}x_{1}-n_{1})\psi(2^{j}x_{2}-n_{2})$ (8)
$\displaystyle\Psi_{\overline{H}H,j,n}(x_{1},x_{2})$ $\displaystyle:=$
$\displaystyle\overline{\psi(2^{j}x_{1}-n_{1})}\psi(2^{j}x_{2}-n_{2})$ (9)
$\displaystyle\Psi_{\overline{L}H,j,n}(x_{1},x_{2})$ $\displaystyle:=$
$\displaystyle\overline{\phi(2^{j}x_{1}-n_{1})}\psi(2^{j}x_{2}-n_{2})$ (10)
$\displaystyle\Psi_{\overline{H}L,j,n}(x_{1},x_{2})$ $\displaystyle:=$
$\displaystyle\overline{\psi(2^{j}x_{1}-n_{1})}\phi(2^{j}x_{2}-n_{2}).$ (11)

Here $n=(n_{1},n_{2})\in{\mathbb{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:

$d_{\nu}(j,n):=\langle
f,\Psi_{\nu,j,n}\rangle=\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}f(x_{1},x_{2})\Psi_{\nu,j,n}(x_{1},x_{2})\,dx_{1}dx_{2}.$
(12)

Here, the function $f$ can be recovered from the coefficients by

$f(x_{1},x_{2})=\sum_{\nu\in
J}\sum_{j=1}^{J}\sum_{n_{1}=0}^{2^{j}-1}\sum_{n_{2}=0}^{2^{j}-1}d_{\nu}(j,n)S^{-1}\Psi_{\nu,j,n}(x_{1},x_{2}),$
(13)

where $S$ is the frame operator

$S=T^{\ast}T,$

where $T$ is the analysis operator and $T^{\ast}$ 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.

(a) Phantom

(b) $\overline{L}H$

(c) $\overline{H}H$

(d) $\overline{H}L$

(e) $HL$

(f) $HH$

(g) $LH$

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.

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:

$c(n)=\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}f(x_{1},x_{2})\phi(x_{1}-n_{1})\phi(x_{2}-n_{2})\,dx_{1}dx_{2},\quad
n=(n_{1},n_{2}).$ (14)

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.

###  2.2 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 $\theta_{0}$ of singularities. The signal $f$ is smooth near $x_{0}$
if there is a cutoff function $\phi\in C^{\infty}_{c}$, $\phi(x_{0})\neq 0$,
such that the Fourier transform of $\phi f$ decays rapidly. That is,

$\widehat{\phi
f}(\xi)=\mathcal{O}(|\xi|)^{-N}\quad|\xi|\rightarrow\infty\,\text{for
all}\,N>0.$

The signal $f$ has a singularity in $x_{0}$ if for all cutoff functions
$\phi$, the Fourier transform of $\phi f$ does not decay rapidly. The set of
all singularities of $f$ is called the singular support and is denoted by
$\texttt{singsupp}(f)$. To define the orientation of the singularities
$x\in\texttt{singsupp}(f)$, we look for directions along which the localized
Fourier transform $\widehat{\phi f}$ does not decay rapidly.

Define the wavefront set of a function $f$ as the set
$WF(f)(x_{0},\alpha_{0})$ with location $x_{0}\in\texttt{singsupp}(f)$ and
direction $\alpha_{0}$, such that for all cutoff functions $\phi\in
C^{\infty}_{c}$, $\phi(x_{0})\neq 0$ , the localized Fourier transform
$\widehat{\phi f}(\xi)$ does not decay rapidly in any polar wedge
$W_{\delta}=\\{(r,w):|w-\alpha_{0}|<\delta\\}$, where $(r,w)$ are the polar
coordinates in the frequency domain. The direction $\alpha_{0}$ of a
singularity $x_{0}\in\texttt{singsupp}{f}$ 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 $\gamma$, $\texttt{singsupp}(f)$
equals $\gamma$ and $WF(f)$ consists of pairs $(x_{0},N(x_{0}))$ with
$x_{0}\in\gamma$ and $N(x_{0})$ the normal vector of $\gamma$ at $x_{0}$.

Consider the continuous two-dimensional Radon transform $\mathcal{R}f$ of a
function $f$:

$\mathcal{R}f(\theta,s)=\int_{L(\theta,s)}f(x)dS(x),$

where
$L(\theta,s)=\\{x\in\mathbb{R}^{2}:x_{1}\text{cos}(\theta)+x_{2}\text{sin}(\theta)=s\\}$
is the line with normal direction $\theta$ and signed distance from the
orientation $s$. Given the Radon transform data $\mathcal{R}f(\theta,s)$ for
$(\theta,s)$ arbitrarily near $(\theta_{0},s_{0})$, one can reconstruct
singularities of $f$ stably with location $x\in L(\theta_{0},s_{0})$ and
direction $\theta_{0}$. In other words, a singularity $(x_{0},\theta_{0})$ of
$f$ is visible from a limited angle measurement if the line
$L(\theta_{0},x_{0}\cdot\theta_{0})$ is recorded in the measurements. See [13,
25, 11, 10] for references.

In (parallel-beam) limited-angle tomography, we know the Radon transform
$\mathcal{R}f(\theta,s)$ for all $s\in\mathbb{R}$ and for the angular range
$[\theta_{0}-\delta,\theta_{0}+\delta]$, for some $\theta_{0}\in\mathbb{R}$
and $0<\delta<\pi$. Then, the visible part of the wavefront set is
$WF(f)(x_{0},\alpha)$, for

$\pi/2+\theta_{0}-\delta\leq\alpha\leq\pi/2+\theta_{0}+\delta.$

In Figure 5, we show an example of the visible part of the wavefront set for
$\delta=\pi/4$, (45°), and $\delta=\pi/18$, (20°).

(a) Phantom(b) Full angle(c) $|\theta|\leq 45^{\circ}$(d) $|\theta|\leq
20^{\circ}$ 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^{\circ}\leq\theta\leq
45^{\circ}$, 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^{\circ}\leq\theta\leq 20^{\circ}$.

###  2.3 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
${\bf{f}}\in\mathbb{R}^{n}$ based on X-ray measurements
${\bf{m}}\in\mathbb{R}^{m}$ is modeled by

${\bf{m}}=A{\bf{f}}+\epsilon,$ (15)

where $A\in\mathbb{R}^{m\times n}$ is a discretized linear forward operator
and $\epsilon>0$ models additive Gaussian noise. We consider regularized
solutions to the inverse problems, achieved by minimizing the following
functional:

${\bf{f}}_{S}=\bigg{\\{}\mathop{\mathrm{arg\,min}}_{{\bf{f}}\in\mathbb{R}^{n},{\bf{f}}\geq
0}\frac{1}{2}\|A{\bf{f}}-m\|_{2}^{2}+\mu\|W_{\mathbb{C}}{\bf{f}}\|_{1}\bigg{\\}},$
(16)

where the term $\|W_{\mathbb{C}}{\bf{f}}\|_{1}$ promotes solutions having
sparse representation in the wavelet basis, and $\mu>0$ serves as a
regularization parameter providing a balance between data fidelity and _a
priori_ information. The component-wise inequality ${\bf{f}}\geq 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:

$\left\\{\begin{array}[]{ll}{\bf{y}}^{k+1}&=\mathbb{P}_{C}({\bf{f}}^{k}-\tau\bigtriangledown{\mathcal{G}}({\bf{f}}^{(k)})-\lambda(W_{\mathbb{C}})^{T}{\bf{v}}^{k}),\\\
{\bf{v}}^{k+1}&=(I-\mathcal{T}_{\mu\tau/\lambda})(W_{\mathbb{C}}{\bf{y}}^{k+1}+{\bf{v}}^{k}),\\\
{\bf{f}}^{k+1}&=\mathbb{P}_{C}({\bf{f}}^{k}-\tau\bigtriangledown{\mathcal{G}}({\bf{f}}^{(k)})-\lambda(W_{\mathbb{C}})^{T}{\bf{v}}^{(k+1)}),\end{array}\right.$
(17)

where $\tau$ and $\lambda$ are positive parameters,
${\mathcal{G}}({\bf{f}})=\frac{1}{2}\|A{\bf{f}}-{\bf{m}}\|_{2}^{2}$, and
$\mu>0$ represents the regularization parameter. Parameters
$0<\lambda<\frac{1}{\lambda_{\text{max}}(W_{\mathbb{C}}(W_{\mathbb{C}})^{T})}$,
where $\lambda_{\text{max}}$ is the maximum eigenvalue, and
$0<\tau<2/\tau_{\text{lip}}$, where $\tau_{\text{lip}}$ is the Lipschitz
constant for $\mathcal{G}({\bf{f}})$, need to be suitably chosen for
convergence. The soft-thresholding operator $\mathcal{T}$ is defined radially
for complex-valued inputs as

$\mathcal{T}_{\mu}(\theta)=(|\theta|-\mu)e^{i\theta}.$ (18)

In equation (17), the non-negative quadrant is denoted by
$C=\mathbb{R}^{N^{2}}_{+}$ and $\mathbb{P}_{C}$ is the Euclidean projection,
that is, the operator $\mathbb{P}_{C}$ replaces any non-negative elements in
the input vector by zeroes.

###  2.4 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 $\mathcal{D}\subset\mathbb{Z}^{2}$. A structuring
element $S\subset\mathbb{Z}^{2}$ is a pre-defined binary image that is used to
probe the object $\mathcal{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

$\mathcal{D}\oplus S=\bigcup_{s\in S}\mathcal{D}_{s},$ (19)

where $\mathcal{D}_{s}=\\{d+s|d\in\mathcal{D},s\in S\\}$ is the translation of
$\mathcal{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

$\mathcal{D}\ominus S=\bigcap_{s\in S}\mathcal{D}_{-s}.$ (20)

The morphological opening operation removes small details from the object but
preserves the shape and size of larger objects. The opening of $\mathcal{D}$
by $S$ is defined as

$\mathcal{D}\circ S=(\mathcal{D}\ominus S)\oplus S,$ (21)

where $\mathcal{D}\ominus S$ is the morphological erosion of $\mathcal{D}$ by
$S$, and $\oplus$ is the dilation operation of the result by $S$. The opening
of $\mathcal{D}$ by $S$ can be seen as the union of all the translations of
$S$ so that it fits within the object $\mathcal{D}$:

$\mathcal{D}\circ S=\bigcup_{z\in\mathcal{E}_{\mathcal{D}}}S_{z},$ (22)

where $\mathcal{E}_{\mathcal{D}}=\\{z\in\mathbb{Z}^{2}:\
S_{z}\subseteq\mathcal{D}\\}$.

The morphological skeleton operation extracts the centerline of the objects in
a binary image while preserving its topology. The operation transforms the
object $\mathcal{D}$ to 1-pixel wide curved lines, while not changing the
essential structure. It can be expressed in terms of morphological erosions
and openings:

$Skeleton(\mathcal{D})=\bigcup_{k=0}^{K}(\mathcal{D}\ominus
kS)-(\mathcal{D}\ominus kS)\circ S,$ (23)

where $(\mathcal{D}\ominus kS)$ expresses $k$ successive erosion operations
and $K=\max\\{k|(\mathcal{D}\ominus kS)\neq\emptyset\\}$ denotes the last step
before the object $\mathcal{D}$ erodes to an empty set.

Morphological operations can be extended for greyscale images. In that case,
$\mathcal{D}$ in formula (22) is a greyscale pixel image. In greyscale
morphology, images $d$ are functions mapping the Euclidean space $E$ into
$\mathbb{R}\cup\\{-\infty,\infty\\}$. That is,
$d:E\rightarrow\mathbb{R}\cup\\{-\infty,\infty\\}.$ The greyscale structuring
elements are also functions mapping Euclidean space into
$\mathbb{R}\cup\\{-\infty,\infty\\}$.

Denote an image by $d(x)$ and a structuring function by $b(x)$. Then, the
dilation operation is defined as

$(d\oplus b)(x)=\sup_{y\in E}[d(y)+b(x-y)],$ (24)

where $\sup$ denotes the supremum. The erosion of $d$ by $b$ is given by

$(d\ominus b)(x)=\inf_{y\in E}[d(y)-b(x-y)],$ (25)

while the opening and closing are given by

$d\circ b=(d\ominus b)\oplus b,$ (26) $d\circ b=(d\oplus b)\ominus b.$ (27)

If $b(x)={\bf 1}_{B}(x)$ is the indicator function of the set $B$ and
$d(x)={\bf 1}_{D}(x)$ is the indicator function of the set $D$, we have

$(d\circ b)(x)={\bf 1}_{D\circ B}(x).$ (28)

Similar formulas hold for dilation $d\oplus b$ and erosion $d\ominus b$.

##  3 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.

(a)

(b)

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: 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.

###  3.1 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_{\mathbb{\nu}}(j,n_{1},n_{2})$, where $j=6$,
$n=(n_{1},n_{2})\in{\mathbb{Z}}^{2}$, and $\nu\in\mathcal{I}$ is as described
in section 2.1. Then, we take the absolute value of the coefficients:

$A_{\nu,6}(n_{1},n_{2})=|{d}_{\nu}(6,n_{1},n_{2})|.$ (29)

A demonstration can be seen in Figure 8(a).

Due to the limited-angle data, we observe nonzero complex wavelet coefficients
mainly in subbands $d_{\overline{H}L}(6,n_{1},n_{2})$ and
$d_{HL}(6,n_{1},n_{2})$, while subbands $d_{\overline{L}H}(6,n_{1},n_{2})$,
$d_{\overline{H}H}(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 limited-angle sinogram can only encode information
about certain directions of the wavefront set, as explained in Section 2.2.

(a) Absolute values of the detail coefficients
$A_{\nu,6}(n_{1},n_{2})=|{d}_{\nu}(6,n_{1},n_{2})|$ of a limited-angle PDFP
reconstruction.

(b) Cleaned coefficients $\widetilde{A}_{\nu,6}$ obtained by morphological
opening with a line structuring element $S_{\nu,6}$. Input for the neural
network $\mathcal{N}_{1}$.

(c) Binary mask $\widehat{A}_{\nu,6}(n_{1},n_{2})$, output from neural network
$\mathcal{N}_{1}$ that performs binary thresholding.

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:
$\overline{L}H,\overline{H}H,\overline{H}L,HL,HH,LH$.

###  3.2 Morphological opening of the complex wavelet coefficients

(a) $\overline{L}H$

(b) $\overline{H}H$

(c) $\overline{H}L$

(d) $HL$

(e) $HH$

(f) $LH$

Figure 9: Line structuring elements $S_{\nu,6}$ that are binary images, for
each index $\nu$ 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.

Due to stretching artefacts and imperfect measurements, the complex wavelet
coefficients available in subbands $d_{\overline{H}L}(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_{\overline{L}H}(6,n_{1},n_{2}),d_{\overline{H}H}(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_{\overline{H}L}(6,n_{1},n_{2})$ and $d_{HL}(6,n_{1},n_{2})$, we perform
morphological opening with oriented line structuring elements $S_{\nu,6}$:

$\widetilde{A}_{\nu,6}=A_{\nu,6}\circ S_{\nu,6},$ (30)

where $A_{\nu,6}$, defined in equation (29), is a grayscale image considered
as a function $A_{\nu,6}:{\mathbb{Z}}^{2}\to[0,1]$. The line structuring
elements $S_{\nu,6}$ are the binary images shown in Figure 9 and are
considered as a function $S_{\nu,6}:{\mathbb{Z}}^{2}\to\\{0,1\\}$. For our
limited-angle problem, we use structuring elements $S_{\overline{H}L,6}$ and
$S_{HL,6}$, since we perform opening for subbands
$d_{\overline{H}L}(6,n_{1},n_{2})$ and $d_{HL}(6,n_{1},n_{2})$. The result is
illustrated in Figure 8(b).

Now, one would think that hard-thresholding of the opened subbands is as
follows:

$\left\\{\begin{array}[]{ll}1,&\widetilde{A}_{\nu,6}(n_{1},n_{2})\geq\epsilon,\\\
0&\textrm{otherwise,}\end{array}\right.$

would give binary indicators for the discrete wavefront set. This would be the
case if we had a good threshold $\epsilon>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.

###  3.3 Learning the binary mask

We use a convolutional neural network
$\mathcal{N}_{1}:\mathbb{R}^{{\mathbb{Z}}^{2}\times\mathcal{I}}\to[0,1]^{{\mathbb{Z}}^{2}\times\mathcal{I}}$
to learn the binary indicators for the discrete wavefront set in the two
subbands $d_{\overline{H}L}(6,n_{1},n_{2})$ and $d_{HL}(6,n_{1},n_{2})$. The
input of neural network $\mathcal{N}_{1}$ is the tensor
$\widetilde{A}=(\widetilde{A}_{\nu,6}(n))_{n\in{\mathbb{Z}}^{2},\nu\in\mathcal{I}}$
obtained from equation (30) and the network maps it to
$\widehat{A}=(\widehat{A}_{\nu,6}(n))_{n\in{\mathbb{Z}}^{2},\nu\in\mathcal{I}}$,
that is,

$\widehat{A}=\mathcal{N}_{1}(\widetilde{A}).$

The training of the neural network uses the morphologically opened complex
wavelet coefficients $\widetilde{A}_{\overline{H}L,6}(n_{1},n_{2})$ and
$\widetilde{A}_{HL,6}(n_{1},n_{2})$. After the network $\mathcal{N}_{1}$
performs the thresholding, it outputs a binary mask with ones indicating the
location of the wavefront set. See Figure 8(c) for an example result.

Training inputs are generated as follows. We first generate 5000 two-
dimensional phantoms of size $128\times 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_{\nu,6}(n_{1},n_{2})$, of each reconstruction. We set
subbands
$A_{\overline{L}H,6}(n_{1},n_{2}),A_{\overline{H}H,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_{\overline{H}L,6}(n_{1},n_{2})$ and $A_{HL,6}(n_{2},n_{2})$ by
using formula (30). These morphologically opened sets of subbands
$\widetilde{A}_{\nu,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_{\overline{L}H}(6,n_{1},n_{2}),d_{\overline{H}H}(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_{\overline{H}L}(6,n_{1},n_{2})$ and $d_{HL}(6,n_{1},n_{2})$ using
(30) to get ${\widetilde{A}}_{\overline{H}L,6}(n_{1},n_{2})$ and
${\widetilde{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 $\mathcal{N}_{1}$, the ground truth data is obtained using the
formula

$\overline{A}_{\nu,6}(n_{1},n_{2})=\left\\{\begin{array}[]{ll}1,&\widetilde{A}_{\nu,6}(n_{1},n_{2})\geq\epsilon,\\\
0&\textrm{otherwise}\end{array}\right.$ (31)

using a suitably chosen $\epsilon>0$.

128 x 128 x 6128 x 128 x 6464 x 64 x 6464 x 64 x 12832 x 32 x 12832 x 32 x
25616 x 16 x 25616 x 16 x 25632 x 32 x 25632 x 32 x 12864 x 64 x 12864 x 64 x
64128 x 128 x 64128 x 128 x 63 x 3 conv + ReLu1 x 1 conv + sigmoidBN +
down/upsampleconcatenate Figure 10: Network architecture of $\mathcal{N}_{1}$.

We implement neural network $\mathcal{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\times 3$, stride
of size 1, and rectified linear unit (ReLU) activation functions. After each
convolutional layer, we perform batch normalization, followed by down- or up-
sampling. We use skip connections to join the encoder and decoder parts of the
network to improve training. For the last layer, we use a $1\times 1$
convolution with a sigmoid activation. See the detailed architecture in Figure
10. The best results were achieved by using the dice coefficient loss function
given by:

$L(\overline{A},\widehat{A})=1-\frac{2\sum(\overline{A}\cdot\widehat{A})}{\sum(\overline{A}+\widehat{A})},$
(32)

where $\overline{A}$ denotes the ground truth values and $\widehat{A}$ the
network predictions.

The training input and output datasets were of size $5000\times 128\times
128\times 6$ each. During training, 500 training set examples were used as the
validation set.

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.

(a) $\widetilde{S}_{\overline{H}H\rightarrow\overline{L}H}$

(b) $\widetilde{S}_{\overline{H}L\rightarrow\overline{H}H}$

(c) $\widetilde{S}_{HL\rightarrow HH}$

(d) $\widetilde{S}_{HH\rightarrow LH}$

Figure 12: The custom-made directional structuring elements
$\widetilde{S}_{\nu\rightarrow\nu^{\prime}}$ we use for our specific limited-
angle measurement geometry.

##  4 Learned microlocal prior

In the previous section, we explained how to extract the known part of the
wavefront set from a limited-angle 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 $\mathcal{N}_{2}$,
using the output of the first neural network $\mathcal{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.

(a) Binary mask $\widehat{A}_{\nu,6}(n_{1},n_{2})$, starting point for the
dilation operation. Output of the first neural network $\mathcal{N}_{1}$.

(b) Initial guess of the wavefront set $D_{\nu,6}(n_{1},n_{2})$, computed
using morphological dilation. Input for the second network, $\mathcal{N}_{2}$.

(c) Full prediction of the wavefront set $\widehat{D}_{\nu,6}(n_{1},n_{2})$,
output of the second network $\mathcal{N}_{2}$.

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:
$\overline{L}H,\overline{H}H,\overline{H}L,HL,HH,LH$.

###  4.1 Approximation of the microlocal prior by morphological dilation

Now, we extend the wavefront set information from the two available subbands
(see Figure 13(a)) 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
$\widetilde{S}_{\nu\rightarrow\nu^{\prime}}$, $\nu,\nu^{\prime}\in\mathcal{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

$S^{-}=\\{(x,y)\in\mathbb{R}^{2}:\ x\in[-1,1],y=-a\,\hbox{sign}(x)x^{2}\hbox{
with some }0\leq a\leq 1\\}$

and the elements (c) and (d) are the rotated and scaled images of the set

$S^{+}=\\{(x,y)\in\mathbb{R}^{2}:\ x\in[-1,1],\
y=+a\,\hbox{sign}(x)x^{2}\hbox{ with some }0\leq a\leq 1\\}.$

The design is motivated by the assumption that the unknown object may have a
discontinuity on a curve $\gamma$ 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 $f_{2}$:

$f(x_{1},x_{2})=f_{1}(x_{1},x_{2}){\bf
1}_{D}(x_{1},x_{2})+f_{2}(x_{1},x_{2}),$

and $\gamma:[-1,1]\to\mathbb{R}^{2}$ is a path satisfying
$\gamma([-1,1])\subset\partial D$. Let us consider the case when the path
$\gamma$ is parametrized as $\gamma(s)=(s,h(s))$, $s\in[-1,1]$, where
$h:[-1,1]\to\mathbb{R}$ is a function that satisfies $h(0)=0$ and
$h^{\prime}(0)=0$. Moreover, assume that the normal vector of $\gamma(x)$
turns clockwise when $s>0$ and $s$ grows. Then for $s>0$ the second derivative
of $h$ satisfies $h^{\prime\prime}(s)<0$. If we also have
$0>h^{\prime\prime}(s)\geq-1$, then $\gamma([0,1])\subset S^{-}$. When $s<0$,
we assume that the normal vector of $\gamma(x)$ turns clockwise when $-s$
grows and that $0<h^{\prime\prime}(s)\leq 1$, in which case
$\gamma([-1,0])\subset S^{-}$. Note that for example the structuring element
$\widetilde{S}_{\overline{H}L\to\overline{H}H}$ corresponds to predicting the
points $x=(x_{1},x_{2})$ on the discontinuity $\partial D$ with a normal
vector $\nu(x)$ of $\partial D$ such that there exists a nearby point
$\tilde{x}=(\tilde{x}_{1},\tilde{x}_{2})\in\partial D$ with a normal vector
$\nu(\tilde{x})$ such that $\nu(x)$ is a vector close to $\nu(\tilde{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 $\overline{H}L$ and
$HL$, since we have already extracted the wavefront set for these subbands in
the previous steps. That is,

$\displaystyle D_{\overline{H}L,6}(n_{1},n_{2})$ $\displaystyle=$
$\displaystyle\widehat{A}_{\overline{H}L,6}(n_{1},n_{2})$ (33) $\displaystyle
D_{HL,6}(n_{1},n_{2}),$ $\displaystyle=$
$\displaystyle\widehat{A}_{HL,6}(n_{1},n_{2}).$ (34)

To get the wavefront set estimate for the adjacent subbands $\overline{H}H$
and $HH$, we dilate subbands $\overline{H}L$ and $HL$ with the structuring
elements $\widetilde{S}_{\overline{H}L\rightarrow\overline{H}H}$ and
$\widetilde{S}_{HL\rightarrow HH}$, respectively. That is,

$\displaystyle D_{\overline{H}H,6}$ $\displaystyle=$ $\displaystyle
D_{\overline{H}L,6}\oplus\widetilde{S}_{\overline{H}L\rightarrow\overline{H}H}$
(35) $\displaystyle D_{HH,6}$ $\displaystyle=$ $\displaystyle
D_{HL,6}\oplus\widetilde{S}_{HL\rightarrow HH}.$ (36)

Next, for subbands $\overline{L}H$ and $LH$, we treat the once dilated
subbands $\overline{H}H$ and $HH$ as the objects for the dilation with
structuring elements $\widetilde{S}_{\overline{H}H\rightarrow\overline{L}H}$
and $\widetilde{S}_{HH\rightarrow LH}$, respectively. That is,

$\displaystyle D_{\overline{L}H,6}$ $\displaystyle=$ $\displaystyle
D_{\overline{H}H,6}\oplus\widetilde{S}_{\overline{H}H\rightarrow\overline{L}H}=\big{(}D_{\overline{H}L,6}\oplus\widetilde{S}_{\overline{H}L\rightarrow\overline{H}H}\big{)}\oplus\widetilde{S}_{\overline{H}H\rightarrow\overline{L}H}$
(37) $\displaystyle D_{LH,6}$ $\displaystyle=$ $\displaystyle
D_{HH,6}\oplus\widetilde{S}_{HH\rightarrow
LH}=\big{(}D_{HL,6}\oplus\widetilde{S}_{HL\rightarrow
HH}\big{)}\oplus\widetilde{S}_{HH\rightarrow LH}.$ (38)

Now, we have a coarse initial guess of the wavefront set for the network
$\mathcal{N}_{2}$ to refine. See Figure 13(b) for an example.

###  4.2 Learning to extend the wavefront set into the entire domain

We want to learn to refine the wavefront set estimate with neural network
$\mathcal{N}_{2}$, based on the approximation provided by the dilation
operation. The network $\mathcal{N}_{2}$ outputs a prediction of the full
wavefront set $\widehat{D}_{\nu,6}(n_{1},n_{2})$, extending it to all six
subbands. We use a convolutional neural network
$\mathcal{N}_{2}:\mathbb{R}^{{\mathbb{Z}}^{2}\times\mathcal{I}}\to[0,1]^{{\mathbb{Z}}^{2}\times\mathcal{I}}$,
where the input of $\mathcal{N}_{2}$ is the tensor
$D=(D_{\nu,6}(n))_{n\in{\mathbb{Z}}^{2},\nu\in\mathcal{I}}$ and the network
maps it to
$\widehat{D}=(\widehat{D}_{\nu,6}(n))_{n\in{\mathbb{Z}}^{2},\nu\in\mathcal{I}}$.
That is,

$\widehat{D}=\mathcal{N}_{2}(D).$

See Figure 13(c) for an example output.

To get the training inputs for network $\mathcal{N}_{2}$, we use the same
phantoms as with network $\mathcal{N}_{1}$. We perform morphological dilation
(as explained in section 4.1) on thresholded ground truth subbands
$\overline{A}_{\nu,6}(n_{1},n_{2})$, see equation 31, to get the
approximations for the wavefront set $D_{\nu,6}(n)$. These act as the training
inputs for $\mathcal{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 $\mathcal{N}_{2}$, we use a similar architecture as for
network $\mathcal{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 $\mathcal{N}_{1}$. See section 3.3
for a more detailed explanation. Similarly, the training datasets were of size
$5000\times 128\times 128\times 6$, from which we used 500 samples for
validation.

128 x 128 x 6128 x 128 x 6464 x 64 x 6464 x 64 x 12832 x 32 x 12832 x 32 x
25616 x 16 x 25616 x 16 x 25632 x 32 x 25632 x 32 x 12864 x 64 x 12864 x 64 x
64128 x 128 x 64128 x 128 x 63 x 3 conv + ReLu1 x 1 conv + sigmoidBN +
down/upsampleconcatenate Figure 14: Network architecture of $\mathcal{N}_{2}$.

###  4.3 Boundary estimate

Now, we can combine the extended wavefront set
$\widehat{D}_{\nu,6}(n_{1},n_{2})$, that is the output of network
$\mathcal{N}_{2}$, to form the singular support (see section 2.2). That is, we
take the maximum of all six binary subbands:

$\texttt{singsupp}(\widehat{D}_{\nu,6})=\max\Big{(}\widehat{D}_{(\overline{L}H,6)},\widehat{D}_{(\overline{H}H,6)},\widehat{D}_{(\overline{H}L,6)},\widehat{D}_{(HL,6)},\widehat{D}_{(HH,6)},\widehat{D}_{(LH,6)}\Big{)}.$

See Figure 15(a) 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 15(b). Finally, this learned boundary estimate can be used as an
overlay for the PDFP reconstruction. See section 5 for final results.

(a)

(b)

Figure 15: (a) The singular support $\texttt{singsupp}(\widehat{D}_{\nu,6})$
computed from the predicted full wavefront set $\widehat{D}_{\nu,6}$ and (b)
its boundary estimate.

##  5 Results

###  5.1 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}$):

$\textrm{DSC}=\frac{2|X_{z}\cap Y_{z}|}{|X_{z}|+|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.

Figure 16: A computational 3D phantom of size $128\times 128\times 128$: three
$L^{p}$ -balls, where $p=1.5$, with a constant value of 1, in a 0-valued
background.

(a)

(b)

(c)

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.

###  5.2 Reconstructions in the $xy$-plane

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.68966
57 | 0.86774
61 | 0.85950
84 | 0.83652
94 | 0.79417
113 | 0.71338
Table 1: Dice similarity 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 in Figures 17, 18.

Slice
$y=22$

$y=61$

$y=84$

$y=94$

$y=113$

(a)

(b)

(c)

(d)

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).

Slice
$z=7$

$z=35$

$z=57$

$z=94$

$z=117$

(a)

(b)

(c)

(d)

Figure 19: Different $xy$-slices through the three-dimensional test phantom of
Figure 16 and through 3D-volumes achieved by stacking 2D-reconstructions from
$xz$-planes. In the model, $xy$-slices are parallel to the detector surface.
(a) Ground truth slice. (b) Tomosynthesis. (c) PDFP reconstruction with
learned boundary estimate. (d) Boundary estimate curve. Note that one can
deduce from the boundary curves that the gray features in the top and bottom
images of column (c) are artefacts.

##  6 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 $2$D. 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 $\mathbb{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.

## Acknowledgements

Siiri Rautio was supported by a University of Helsinki -funded doctoral
researcher position. Rashmi Murthy was funded by the Future Makers project
AIDMEI, funded by the Jane and Aatos Erkko Foundation and Technology
Industries of Finland Centennial Foundation. Tatiana A. Bubba was supported by
the Academy of Finland (310822) and by the Royal Society through the Newton
International Fellowship grant n. NIF\R1\201695. Matti Lassas and Samuli
Siltanen are supported by the Finnish Centre of Excellence in Inverse
Modelling and Imaging, 2018-2025, decision numbers 312339 and 336797, Academy
of Finland grants 284715 and 312110.

## References

  * [1] Jonas Adler and Ozan Öktem.  Solving ill-posed inverse problems using iterative deep neural networks.  Inverse Problems, 33(12):124007, 2017.
  * [2] Jonas Adler and Ozan Öktem.  Learned primal-dual reconstruction.  IEEE transactions on medical imaging, 37(6):1322–1332, 2018.
  * [3] Héctor Andrade-Loarca, Gitta Kutyniok, Ozan Öktem, and Philipp Petersen.  Deep microlocal reconstruction for limited-angle tomography.  arXiv preprint arXiv:2108.05732, 2021.
  * [4] Simon Arridge, Peter Maass, Ozan Öktem, and Carola-Bibiane Schönlieb.  Solving inverse problems using data-driven models.  Acta Numerica, 28:1–174, 2019.
  * [5] Tatiana A Bubba, Gitta Kutyniok, Matti Lassas, Maximilian März, Wojciech Samek, Samuli Siltanen, and Vignesh Srinivasan.  Learning the invisible: A hybrid deep learning-shearlet framework for limited angle computed tomography.  Inverse Problems, 35(6):064002, 2019.
  * [6] Emmanuel J Candès and David L Donoho.  New tight frames of curvelets and optimal representations of objects with piecewise c2 singularities.  Communications on Pure and Applied Mathematics: A Journal Issued by the Courant Institute of Mathematical Sciences, 57(2):219–266, 2004.
  * [7] Peijun Chen, Jianguo Huang, and Xiaoqun Zhang.  A primal–dual fixed point algorithm for convex separable minimization with applications to image restoration.  Inverse Problems, 29(2):025011, 2013.
  * [8] Mark E Davison.  The ill-conditioned nature of the limited angle tomography problem.  SIAM Journal on Applied Mathematics, 43(2):428–448, 1983.
  * [9] James T Dobbins III and Devon J Godfrey.  Digital x-ray tomosynthesis: current state of the art and clinical potential.  Physics in medicine & biology, 48(19):R65, 2003.
  * [10] Jürgen Frikel.  Sparse regularization in limited angle tomography.  Applied and Computational Harmonic Analysis, 34(1):117–141, 2013\.
  * [11] Jürgen Frikel and Eric Todd Quinto.  Characterization and reduction of artifacts in limited angle tomography.  Inverse Problems, 29(12):125007, 2013.
  * [12] Rafael C. Gonzalez and Richard E. Woods.  Digital image processing.  Pearson, 2018.
  * [13] Allan Greenleaf and Gunther Uhlmann.  Nonlocal inversion formulas for the x-ray transform.  Duke mathematical journal, 58(1):205–240, 1989.
  * [14] Misty I Haith, Peter Huthwaite, and Michael JS Lowe.  Defect characterisation from limited view pipeline radiography.  NDT & E International, 86:186–198, 2017.
  * [15] Kyong Hwan Jin, Michael T. McCann, Emmanuel Froustey, and Michael Unser.  Deep convolutional neural network for inverse problems in imaging.  IEEE Transactions on Image Processing, 26(9):4509–4522, 2017.
  * [16] Gitta Kutyniok and Demetrio Labate.  Shearlets: Multiscale analysis for multivariate data.  Springer Science & Business Media, 2012.
  * [17] Demetrio Labate, Wang-Q Lim, Gitta Kutyniok, and Guido Weiss.  Sparse multidimensional representation using shearlets.  In Wavelets XI, volume 5914, page 59140U. International Society for Optics and Photonics, 2005.
  * [18] G Landi, E Loli Piccolomini, and JG Nagy.  A limited memory bfgs method for a nonlinear inverse problem in digital breast tomosynthesis.  Inverse Problems, 33(9):095005, 2017.
  * [19] Germana Landi, E Loli Piccolomini, and J Nagy.  Nonlinear conjugate gradient method for spectral tomosynthesis.  Inverse Problems, 35(9):094003, 2019.
  * [20] Frank Natterer.  Computerized tomography.  In The Mathematics of Computerized Tomography, pages 1–8. Springer, 1986.
  * [21] Loren T Niklason, Bradley T Christian, Laura E Niklason, Daniel B Kopans, Donald E Castleberry, BH Opsahl-Ong, Cynthia E Landberg, Priscilla J Slanetz, Angela A Giardino, Richard Moore, et al.  Digital tomosynthesis in breast imaging.  Radiology, 205(2):399–406, 1997.
  * [22] Daniël M Pelt, Kees Joost Batenburg, and James A Sethian.  Improving tomographic reconstruction from limited data using mixed-scale dense convolutional neural networks.  Journal of Imaging, 4(11):128, 2018.
  * [23] E Loli Piccolomini and E Morotti.  A fast total variation-based iterative algorithm for digital breast tomosynthesis image reconstruction.  Journal of Algorithms & Computational Technology, 10(4):277–289, 2016.
  * [24] Zenith Purisha, Juho Rimpeläinen, Tatiana Bubba, and Samuli Siltanen.  Controlled wavelet domain sparsity for x-ray tomography.  Measurement Science and Technology, 29(1):014002, 2017.
  * [25] Eric Todd Quinto.  Singularities of the X-ray transform and limited data tomography in $\mathbb{R}^{2}$ and $\mathbb{R}^{3}$.  SIAM Journal on Mathematical Analysis, 24(5):1215–1225, 1993.
  * [26] Maaria Rantala, Simopekka Vanska, Seppo Jarvenpaa, Martti Kalke, Matti Lassas, Jan Moberg, and Samuli Siltanen.  Wavelet-based reconstruction for limited-angle x-ray tomography.  IEEE transactions on medical imaging, 25(2):210–217, 2006.
  * [27] Olaf Ronneberger, Philipp Fischer, and Thomas Brox.  U-net: Convolutional networks for biomedical image segmentation.  In International Conference on Medical image computing and computer-assisted intervention, pages 234–241. Springer, 2015.
  * [28] Ivan W Selesnick, Richard G Baraniuk, and Nick C Kingsbury.  The dual-tree complex wavelet transform.  IEEE signal processing magazine, 22(6):123–151, 2005.
  * [29] Samuli Siltanen, Ville Kolehmainen, Seppo Järvenpää, JP Kaipio, P Koistinen, Matti Lassas, J Pirttilä, and Erkki Somersalo.  Statistical inversion for medical x-ray tomography with few radiographs: I. general theory.  Physics in Medicine & Biology, 48(10):1437, 2003.
  * [30] Weslley Silva, Ricardo Lopes, Uwe Zscherpel, Dietmar Meinel, and Uwe Ewert.  X-ray imaging techniques for inspection of composite pipelines.  Micron, 145:103033, 2021.
  * [31] Srinivasan Vedantham, Andrew Karellas, Gopal R Vijayaraghavan, and Daniel B Kopans.  Digital breast tomosynthesis: state of the art.  Radiology, 277(3):663–684, 2015.
  * [32] Ge Wang, Jong Chul Ye, and Bruno De Man.  Deep learning for tomographic image reconstruction.  Nature Machine Intelligence, 2(12):737–748, 2020.
  * [33] Tao Wu, Richard H Moore, Elizabeth A Rafferty, and Daniel B Kopans.  A comparison of reconstruction algorithms for breast tomosynthesis.  Medical physics, 31(9):2636–2647, 2004.