# Benchmarking Graph Neural Networks on Link Prediction

Xing Wang Alexander Vinel
<EMAIL_ADDRESS><EMAIL_ADDRESS>
Auburn University, AL, 36849, USA

###### Abstract

In this paper, we benchmark several existing graph neural network (GNN) models
on different datasets for link predictions. In particular, the graph
convolutional network (GCN), GraphSAGE, graph attention network (GAT) as well
as variational graph auto-encoder (VGAE) are implemented dedicated to link
prediction tasks, in-depth analysis are performed, and results from several
different papers are replicated, also a more fair and systematic comparison
are provided. Our experiments show these GNN architectures perform similarly
on various benchmarks for link prediction tasks.

##  1 Introduction

Our primary focus in this paper is link prediction, which aims to infer
unobserved/missing links or predicting the future ones based on the
connections of currently observed partial networks. Link prediction is a
fundamental problem in network science with tremendous real-world applications
[12].

Over the years, there are many link prediction methods have been developed.
Local link prediction methods based on the assumption that two nodes are more
likely to be connected if they have many common neighbors [18], but generally
obtained low accuracy especially on large and sparse graphs. The global
structural similarity was also taken consideration [11], but is generally
computationally infeasible for large graphs. There are also probabilistic
approaches assuming a known prior structure of the graph, such as hierarchical
or circles structures [1], though the accuracy was not satisfactory. And there
has been a surge of algorithms that make link prediction through
representation learning that learns low dimensional embeddings for nodes, most
of these methods are based on skip-gram or matrix factorization, examples
include DeepWalk [15], node2vec [4], etc. These algorithms are basically
unsupervised learning with no supervised information during training, although
the learned representations can assist downstream tasks such as node
classification, graph classification, as well as link prediction to achieve
higher accuracy.

As the graph neural networks (GNNs) have been an emerging research area in
recent years, significant advances and various architectures are proposed and
developed. Lots of the new models are evaluated on the citation network
datasets (Cora, CiteSeer, PubMed). However, different training/test splits
were often deployed across studies [14], which is challenging to make fair
comparisons among different methods. Also, all three citation networks are
derived from the same domain, with similar structural properties. Furthermore,
when most of the models developed, they were primarily tested on tasks such as
node classification and graph classification, not all performed link
prediction. For instance, in [9], [5] and [17] that we will illustrate in
Section 3, only node classification and graph classification tasks were
benchmarked, but no link prediction task was evaluated. [10] experimented link
prediction in the citation networks, and compared only with baselines that are
not recently developed, such as spectral clustering [16] and DeepWalk [15].
[14] introduced an additional dataset for various GNN benchmarks, but for link
prediction task, only variational graph auto-encoder (VGAE) was evaluated and
compared with the multi-layer perceptron (MLP) baseline. In this paper, we aim
at providing a fair comparison on several most influential GNN architectures
that were developed in recent years, evaluating them on the same datasets,
also with consistent training splits of data among different models. Moreover,
we will benchmark the models on additional dataset other than the commonly
used citation networks.

The rest of this paper is organized as follows. Section 2 introduces the four
datasets we test on. Section 3 illustrates several recently developed GNN
architectures that we apply and evaluate. In Section 4, implementation details
are discussed. We benchmark the models and show the evaluation results in
Section 5. Finally, we conclude our findings and discuss future work
directions in Section 6.

##  2 Datasets

For this paper, we used the following datasets:

  * •

Cora: This dataset is a citation network that is used as a standard benchmark
dataset. The network has 2708 nodes (scientific publications) which are
classified into one of seven classes. The citation network contains 5429 links
(undirected edges) between nodes. Each node is described by a 0/1 bag-of-words
vector with a dictionary size of 1433 [3].

  * •

CiteSeer: This dataset is another citation network similar to Cora. It
contains 3327 nodes (scientific publications), 4732 edges (links), and 6
classes. Each node in this dataset is also described by a 0/1 bag-of-words
vector with a dictionary size of 3703 [3].

  * •

PubMed: This dataset is a citation network of scientific publications
pertaining to diabetes from the Pubmed database. It contains 19717 nodes,
44338 edges, and 3 classes. Each node is described by a TF/IDF weighted word
vector with a dictionary size of 500 [3].

  * •

WikiCS: This dataset is based on Wikipedia article references. It consists of
11701 nodes (Wikipedia articles on Computer Science), 216213 edges (hyperlinks
between articles), and 10 classes. Each node is described by a 300-dimensional
vector which was obtained as the average of pretrained GloVe embeddings [14].

Cora, CiteSeer, and PubMed are commonly used benchmark datasets for tasks
related to graph neural networks such as link prediction, node classification,
and graph classification. While these three datasets are all citation
networks, the Wiki-CS is a newly published dataset in a different domain. It
also has different structural properties, the primary one being that there are
significantly more edges per node. It has a mean degree of 36.94 while PubMed
has a mean degree of 4.5 in spite of having nearly 70% more nodes [14], thus
has a significantly higher connectivity rate and provides a different
distribution as a testbed. In this paper we use the abovementioned four
dataset to benchmark several recently developed GNN architectures on the link
prediction task.

##  3 Approach

Our primary goal is to benchmark several recently developed GNN models for
link prediction task. With consistent training splits of data among different
methods, our work provides a more fair comparison. In particular, the
following four models will be implemented and evaluated for link prediction
tasks on the datasets we discussed earlier in this paper.

###  3.1 Graph Convolutional Networks (GCN)

GCN [9] borrows the concept of convolution from the convolutional neural
network (CNN) and convolve the graph directly according to the connectivity
structure of the graph as the filter to perform neighborhood mixing. The
architecture can be summarized concisely as

$H^{l+1}=\sigma(\tilde{D}^{-1/2}\tilde{A}^{-1/2}\tilde{D}^{-1/2}H^{(l)}W^{(l)}),$

where $A$ and $D$ are the adjacency matrix and the degree matrix of the graph,
respectively. Two tricks were applied that help the GCN being successful: a
self-connection is added to each node in the adjacency matrix, and then the
adjacency matrix is normalized according to the degrees, so that we obtained
$\tilde{A}$ and $\tilde{D}$. $H^{(l)}$ contains the embedding of the nodes in
$l^{th}$ layer, $W^{(l)}$ denotes the weight matrix for this layer, and
$\sigma$ is the nonlinearity.

###  3.2 GraphSAGE

GraphSAGE [5] is an inductive learning algorithm for GNNs, meaning that
instead of applying the whole adjacency matrix information among all node, it
learns aggregator functions that can induce the embedding of a new node given
its features and neighborhood information without retraining of the entire
model. To learn the embeddings with aggregators, a neighborhood embedding
would be mapped with the aggregation function and concatenated with existing
embedding of the node. The concatecated vector would be passed through a GNN
layer to update the node embedding, and finally the embeddings are normalized
to unit norm. Note the aggregator can be of various form, [5] proposed the
options of the mean (as used in GCN), pooling, as well as LSTM. The sampling
from neighbors strategy in GraphSAGE yields impressive performance on node
labeling task over several large-scale networks.

To learn the weights of aggregators and embeddings, the loss function of
GraphSAGE is defined as

$\mathcal{L}=-\log\sigma(z_{u}^{T}z_{v})-Q\mathbb{E}_{v_{n}\sim
P_{n}(v)}\log\sigma(-z_{u}^{T}z_{v_{n}}),$

where $\sigma$ represents the sigmoid function, $u$ and $v$ are two neighbors,
while $v_{n}$ is a negative sample, and $Q$ is the number of negative samples.
The first term aims at maximizing the similarity between embeddings of $u$ and
$v$, while the second term tries to set apart embeddings of negative samples.

###  3.3 Graph Attention Networks (GAT)

In recent years, attention mechanisms have become state-of-the-art (SOTA) for
many tasks on sequential data, such as in natural language processing. GAT
[17] combined attention mechanisms with GNN, aiming at more effective learning
power to the neighborhoods’ features. As the building block of GAT, the graph
attention layer plays the role of aggregation function for the GNN. It first
applies a shared linear transformation to every node, parameterized by the
weight matrix $\mathbf{W}\in\mathbb{R}^{F^{\prime}\times F}$. Then self-
attention is performed on the nodes where a shared attention mechanism is used
to calculate the attention coefficients that capture the importance of node
$j$’s features to node $i$, i.e.,

$e_{ij}=a(\mathbf{W}h_{i},\mathbf{W}h_{j}).$

GAT normalized the coefficients across $j$ using a softmax function. As a
result, the attention mechanism is a single-layer network parameterized by a
weight vector $a\in\mathbb{R}^{2F^{\prime}}$, followed by a nonlinear
activation (e.g., LeakyReLU), and obtains the normalized coefficients as

$\alpha_{ij}=\dfrac{\exp\Big{(}\text{LeakyReLU}(a[\mathbf{W}h_{i},\mathbf{W}h_{j}])\Big{)}}{\sum_{k\in
N_{i}}\exp\Big{(}\text{LeakyReLU}(a[\mathbf{W}h_{i},\mathbf{W}h_{k}])\Big{)}}.$

These attention coefficients are used to compute a linear combination of the
neighbors’ features to obtain aggregated features for each node, i.e.,

$h^{\prime}_{i}=\sum_{k\in N_{i}}\alpha_{ij}\mathbf{W}h_{k}.$

And finally, multi-head attention (i.e., a set of independent attention
mechanism) is often employed to stabilize the learning process of self-
attention.

###  3.4 Variational Graph Auto-Encoders

Kipf and Welling introduced the Variational Graph Auto-Encoders (VGAE)
framework for unsupervised learning on graph-structured data [10]. This
framework is an elegant extension of the Variational Auto-Encoders framework
[8] to apply to graph-structured data. The VGAE uses a graph convolutional
network [9] as the encoder to map the input graph’s node features ($X$) into
the latent representation, followed by an inner product decoder to generate
the conditional probabilities of the adjacency matrix ($A$) given the latent
representation. The hidden layer of the encoder can be a GCN layer as we
described in 3.1:

$\displaystyle\bar{X}=GCN(X,A)=ReLU(\bar{A}XW_{0})$
$\displaystyle\tilde{A}=D^{-\frac{1}{2}}AD^{-\frac{1}{2}}$

where $\tilde{A}$ is the symmetrically normalized adjacency matrix [6]. By
applying reparametrization trick [8], the latent distribution produced by the
encoder can be written as

$\displaystyle\mu$ $\displaystyle=GCN_{\mu}(X,A)=\tilde{A}\bar{X}W_{1}$
$\displaystyle\log\sigma$
$\displaystyle=GCN_{\sigma}(X,A)=\tilde{A}\bar{X}W_{1}$ $\displaystyle Z$
$\displaystyle=\mu+\sigma*\epsilon$

where $\epsilon\sim\mathcal{N}(0,1)$.

Thus the encoder is given by [10] as

$q(z_{i}|X,A)=\mathcal{N}(z_{i}|\mu_{i},diag(\sigma_{i}^{2}))$

and the decoder is

$p(A_{ij}=1|z_{i},z_{j})=\sigma(z_{i}^{T}z_{i})$

Finally, the loss function is a combination of the reconstruction loss and the
KL-divergence:

$\mathcal{L}=\mathbb{E}_{q(Z|X,A)}[\log{p(A|Z)}]-KL[q(Z|X,A)\|p(Z)]$

##  4 Implementation Details

The four GNN models illustrated in Section 3 are benchmarked for link
prediction tasks on the four dataset mentioned in Section 2. In this section,
we describe the architectures for the four models.

We constructed GCN, GraphSAGE, GAT and VGAE models for link prediction also
with PyTorchGeometric [2]. We used the same default architectures as the
original implementations of corresponding papers. 111 The PyTorchGeometric
implementation of GCN is adapted from
https://github.com/rusty1s/pytorch_geometric/blob/master/examples/link_pred.py,
and the implementation of other models follows the paradigm with the change of
network architectures.  For GCN, we built a two level graph convolutional
layers, with 128 and 64 filters, respectively. The two graph convolutional
layers play a similar role as the encoder in VGAE. The outputs of the graph
convolutional layers are used for predict the label of a link between two
nodes similar to the decoder, as we will describe later. As for GraphSAGE, we
first implemented the SAGE module as described in Section 3.2, and we build
the encoder-like module for our GraphSAGE link prediction model as two such
SAGE layers, each with 64 hidden channels. We built two graph attention layers
for GAT, with a ReLU nonlineararity in between. The first layer is of multi-
headed attention using 8 heads, each consists of 8 filters, and with a dropout
rate of 0.6, while the second graph attention layer has a single attention
head consists of 16 filters. And finally, for the VGAE, we defined the encoder
as a hidden graph convolutional layer consists of 32 convolutional filters,
followed by two graph convolutional layers representing the mean and log
standard deviation of the latent normal distribution respectively (as per
[10]). The latent dimension was set as 16. And the loss function optimized was
set to be the sum of the reconstruction loss and KL-divergence loss.

Our implementation of the GCN, GraphSAGE, and GAT models use the same decoder-
like module to output the predicted probabilities of the existence for a link.
Similar to the decoder in VGAE, we use dot product to calculate the
similarities of linkage for corresponding nodes pairs. During training,
negative sampling strategy was deployed to generate negative samples to match
the positive samples (i.e., existing edges in the original graphs) within each
minibatch, and binary cross entropy loss was backpropagated for such a binary
classification problem, which is different from VGAE which is a generative
model.

##  5 Experiments and Results

Method | Cora | CiteSeer | PubMed | Wiki-CS
---|---|---|---|---
GCN | $0.9273\pm 0.0051$ | $0.8911\pm 0.0058$ | $0.9619\pm 0.0011$ | $0.9682\pm 0.0030$
GraphSAGE | $0.9091\pm 0.0038$ | $0.8716\pm 0.0063$ | $0.9018\pm 0.0032$ | $0.8758\pm 0.0423$
GAT | $0.9027\pm 0.0023$ | $0.9084\pm 0.0035$ | $0.9179\pm 0.0002$ | $0.9293\pm 0.0014$
VGAE | 0.9179 $\pm$ 0.0013 | 0.9091 $\pm$ 0.0021 | 0.9293 $\pm$ 0.0028 | 0.9578 $\pm$ 0.0008
Table 1: AUC comparison for GNN models on four datasets. Method | Cora | CiteSeer | PubMed | Wiki-CS
---|---|---|---|---
GCN | $0.9385\pm 0.0043$ | $0.9028\pm 0.0050$ | $0.9636\pm 0.0010$ | $0.9720\pm 0.0025$
GraphSAGE | $0.9183\pm 0.0035$ | $0.8863\pm 0.0052$ | $0.9057\pm 0.0025$ | $0.8778\pm 0.0425$
GAT | $0.9027\pm 0.0023$ | $0.9084\pm 0.0035$ | $0.9037\pm 0.0021$ | $0.9277\pm 0.0013$
VGAE | 0.9297 $\pm$ 0.0013 | 0.9234 $\pm$ 0.0017 | 0.9336 $\pm$ 0.0027 | 0.9623 $\pm$ 0.0008
Table 2: Average precision scores comparison for GNN models on four datasets.

We conducted experiments of link prediction with the GCN, GraphSage, GAT, as
well as VGAE models on the datasets we described in Section 2 (Cora, CiteSeer,
PubMed, Wiki-CS). All four datasets are available via Pytorch Geometric APIs
[2]. For each of the datasets, the edges in the graph were split into
training/validation/testing sets with a ratio of 85%-5%-10%. As we mentioned
above, negative edge samples were used for calculating the loss, by randomly
sampling node pairs without an linkage edge with the same number as the
positive edges within each mini-batch. All of our models were trained using
Adam optimizers [7] for stochastic gradient descent, while the initial
learning rates were set to 0.01, with an exception of the GAT task on WikiCS,
where the initial learning rate was set to 0.001. All the tasks were trained
for 200 epochs, and in order to replicate [14], we run each task for 50 runs
so that the results are reported as the average of 50 runs along with 95%
confidence intervals calculated using bootstapping.

We report the area under the ROC curve (AUC) and average precision scores for
each of link prediction tasks with the four GNN models. The AUC results are
showed in Table 1, while Table 2 reports the average precision scores. From
the two tables, we observe that the GCN model performs best on Cora, PubMed,
as well as the relatively larger Wiki-CS tasks, while VGAE achieves the best
performance on CiteSeer. Overall, the four models perform quite similar on
most of the tasks, with the AUC around 90%. The results obtained by GraphSage
may seem a little poorer for CiteSeer and Wiki-CS, not only in terms of
smaller mean AUC and average precision, also with larger standard deviation on
the metrics, indicating the model is less robust in these tasks. Considering
GraphSage models sample from neighbors for each node during propagation over
the graph, less robustness could be expected especially for the graph with
imbalanced distributions. The GCN and VGAE perform significantly better on
Wiki-CS, as the mean AUCs and the average precisions exceed 95%.

Furthermore, we include the visualizations of the embeddings obtained from
each model for the four datasets. The embeddings were obtained as the the
latent space encodings for the VGAE model, and as the output of the graph
convolutional or attention layers which serves as feature maps for the
downstream classifiers in the GCN, GraphSage and GAT models. T-SNE algorithm
[13] was applied to project the high dimensional embeddings into a 2D plot.
Figures 1, 2, 3, and 4 in Appendix show the t-SNE visualization for the
embeddings of nodes in Cora, CiteSeer, PubMed, and Wiki-CS datasets,
respectively. We observe that after training on the link prediction tasks, our
GNN models are able to generally project the nodes with the same class label
into closed neighboring areas of the latent space with fairly reasonable
separation between the classes. This demonstrates that the models are able to
adequately learn good latent representations to distinguish between the types
of labels, particularly for VGAE in spite of that the model being
unsupervised. We also observe that for those models perform well on a
particular task, the embedding visualization presents a clearer distinction
among different classes, such as in Figures 1(a), 2(b), and 4(a). Opposite
examples include that was shown in Figure 4(b), where the node embeddings
present some spiral patterns, indicating some correlation between the two
projected dimensions, which might be a further implication that the GraphSage
model is overfitting on this task as the prediction performance is worse that
that of other models, and we probably could use this information to monitor
the training, tune the hyperparameters, or apply additional regularization
such as early stopping.

##  6 Conclusions and Discussion

In this paper, we implemented several existing GNN models, and benchmarked on
different datasets for link predictions. We not only reproduced the results in
[10] and [14], also provided a more fair and systematic comparison. Our
experiments show these GNN architectures perform similarly on other benchmarks
for link prediction tasks. There are several interesting directions to pursue
for the future of this paper. First off, the datasets we benchmarked on are
still relatively small, in the future we could evaluate the models on much
larger graphs, especially with applications in real world. The second
interesting direction in which to take this paper would be to implement more
recently developed GNN models. Furthermore, we could try to design and develop
our own GNN architecture and benchmark on link prediction tasks.

## References

  * [1] Aaron Clauset, Cristopher Moore, and Mark EJ Newman.  Hierarchical structure and the prediction of missing links in networks.  Nature, 453(7191):98–101, 2008.
  * [2] Matthias Fey and Jan E. Lenssen.  Fast graph representation learning with PyTorch Geometric.  In ICLR Workshop on Representation Learning on Graphs and Manifolds, 2019.
  * [3] Lise Gatoor.  Datasets.  https://linqs.soe.ucsc.edu/data.
  * [4] Aditya Grover and Jure Leskovec.  node2vec: Scalable feature learning for networks.  In Proceedings of the 22nd ACM SIGKDD international conference on Knowledge discovery and data mining, pages 855–864, 2016.
  * [5] Will Hamilton, Zhitao Ying, and Jure Leskovec.  Inductive representation learning on large graphs.  In Advances in neural information processing systems, pages 1024–1034, 2017.
  * [6] Fanghao Han.  Tutorial on variational graph auto-encoders, 2019.  https://towardsdatascience.com/tutorial-on-variational-graph-auto-encoders-da9333281129.
  * [7] Diederik P Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  arXiv preprint arXiv:1412.6980, 2014.
  * [8] Diederik P Kingma and Max Welling.  Auto-encoding variational bayes, 2014.
  * [9] Thomas N Kipf and Max Welling.  Semi-supervised classification with graph convolutional networks.  arXiv preprint arXiv:1609.02907, 2016.
  * [10] Thomas N Kipf and Max Welling.  Variational graph auto-encoders.  arXiv preprint arXiv:1611.07308, 2016.
  * [11] Haifeng Liu, Zheng Hu, Hamed Haddadi, and Hui Tian.  Hidden link prediction based on node centrality and weak ties.  EPL (Europhysics Letters), 101(1):18004, 2013.
  * [12] Linyuan Lü and Tao Zhou.  Link prediction in complex networks: A survey.  Physica A: statistical mechanics and its applications, 390(6):1150–1170, 2011.
  * [13] Laurens van der Maaten and Geoffrey Hinton.  Visualizing data using t-sne.  Journal of machine learning research, 9(Nov):2579–2605, 2008.
  * [14] Péter Mernyei and Cătălina Cangea.  Wiki-cs: A wikipedia-based benchmark for graph neural networks.  arXiv preprint arXiv:2007.02901, 2020.
  * [15] Bryan Perozzi, Rami Al-Rfou, and Steven Skiena.  Deepwalk: Online learning of social representations.  In Proceedings of the 20th ACM SIGKDD international conference on Knowledge discovery and data mining, pages 701–710, 2014.
  * [16] Lei Tang and Huan Liu.  Leveraging social media networks for classification.  Data Mining and Knowledge Discovery, 23(3):447–478, 2011.
  * [17] Petar Veličković, Guillem Cucurull, Arantxa Casanova, Adriana Romero, Pietro Lio, and Yoshua Bengio.  Graph attention networks.  arXiv preprint arXiv:1710.10903, 2017.
  * [18] Tao Zhou, Linyuan Lü, and Yi-Cheng Zhang.  Predicting missing links via local information.  The European Physical Journal B, 71(4):623–630, 2009.

##  Appendix A T-SNE Embedding Visualization

(a) GCN

(b) GraphSage

(c) GAT

(d) VGAE

Figure 1: Cora Embeddings Visualization

(a) GCN

(b) GraphSage

(c) GAT

(d) VGAE

Figure 2: CiteSeer Embeddings Visualization

(a) GCN

(b) GraphSage

(c) GAT

(d) VGAE

Figure 3: PubMed Embeddings Visualization

(a) GCN

(b) GraphSage

(c) GAT

(d) VGAE

Figure 4: WikiCS Embeddings Visualization


# Improving reconstructions in nanotomography for homogeneous materials via
mathematical optimization†

Sebastian Kreuz, Benjamin Apeleo Zubiri∗, Silvan Englisch, Sung-Gyu Kang,
Rajaprakash Ramachandramoorthy, Erdmann Spiecker, Frauke Liers, and Jan
Rolfes∗ Department of Data Science
Friedrich-Alexander-Universität Erlangen-Nürnberg
Cauerstr. 11, 91058 Erlangen, Germany
Email<EMAIL_ADDRESS>
Institute of Micro- and Nanostructure Research (IMN)
$\&$ Center for Nanoanalysis and Electron Microscopy (CENEM)
Department of Materials Science and Engineering
Friedrich-Alexander-Universität Erlangen-Nürnberg
Cauerstr. 3, 91058 Erlangen, Germany
Email<EMAIL_ADDRESS>
Department Structure and Nano- / Micromechanics of Materials
Max-Planck-Institut für Eisenforschung GmbH
Max-Planck-Str. 1, 40237 Düsseldorf, Germany
Email<EMAIL_ADDRESS>
Department of Materials Engineering and Convergence Technology
Gyeongsang National University
Jinju-daero 501, 52828 Jinju, Republic of Korea

###### Abstract.

Compressed sensing is an image reconstruction technique to achieve high-
quality results from limited amount of data. In order to achieve this, it
utilizes prior knowledge about the samples that shall be reconstructed.
Focusing on image reconstruction in nanotomography, this work proposes
enhancements by including additional problem-specific knowledge. In more
detail, we propose further classes of algebraic inequalities that are added to
the compressed sensing model. The first consists in a valid upper bound on the
pixel brightness. It only exploits general information about the projections
and is thus applicable to a broad range of reconstruction problems. The second
class is applicable whenever the sample material is of roughly homogeneous
composition. The model favors a constant density and penalizes deviations from
it. The resulting mathematical optimization models are algorithmically
tractable and can be solved to global optimality by state-of-the-art available
implementations of interior point methods. In order to evaluate the novel
models, obtained results are compared to existing image reconstruction
methods, tested on simulated and experimental data sets. The experimental data
comprise one 360° electron tomography tilt series of a macroporous zeolite
particle and one absorption contrast nano X-ray computed tomography (nano-CT)
data set of a copper microlattice structure. The enriched models are optimized
quickly and show improved reconstruction quality, outperforming the existing
models. Promisingly, our approach yields superior reconstruction results,
particularly when information about the samples is available for a small
number of tilt angles only.

##  1\. Introduction

Projection-based nanotomography techniques, including electron tomography
(ET), nano X-ray computed tomography (nano-CT) and micro X-ray computed
tomography (micro-CT), are designed to gain three-dimensional (3D) information
from a series of two-dimensional (2D) projections of an object on a nm to mm
scale.[Ape+21, Wit+21, Bur+14] The object, called sample or specimen, is
placed on a sample holder between a typically stationary X-ray or electron
beam source and a detector array in a transmission electron microscope or
X-ray microscope. The projections encode information on interaction of the
initial beam while penetrating through the sample. To obtain easily
interpretable and reconstructable contrast, the measured intensities should
exhibit a monotonic relationship to specific sample properties such as local
density (i.e., mass attenuation) and thickness. In ET, this can be realized by
using the so-called high-angle annular dark-field (HAADF) scanning
transmission electron microscopy (STEM) imaging mode [Van+12], whereas in
nano- and micro-CT measurements, this is the case using the absorption
contrast mode [Whi+19, Tak+20]. In all cases, the sample holder or stage is
tilted to perform projections from different tilt angles of the object. A
visual representation of this process is depicted in Figure 1.

Figure 1. Principle of projection-based nanotomography: a tilt series of
projections from an object is acquired and, after an alignment with respect to
the common tilt axis, the reconstruction of the object is calculated.

Mathematically, the projections can be modeled as an idealized, noiseless
system of linear equations

(1) $\displaystyle Rf=p,$

where the vector $f\in\mathbb{R}^{n}$ denotes the pixel brightness,
$p\in\mathbb{R}^{m}$ the projection data and the matrix
$R\in\mathbb{R}^{m\times n}$ a discretized variant Radon transform. To be more
precise, each row of $R$ corresponds to a single projection ray. An entry of a
row displays how much of a pixel is passed by the ray in relation to its side
length. Multiple projection geometries including parallel, fan, and cone beams
can be modeled in this way [ZPB15]. However, as (1) describes an idealized
setting, introducing uncertainties, such as noise, potentially render (1)
infeasible. Therefore, image reconstruction is usually performed via error
minimization:

(2) $\displaystyle\min\limits_{f}\,\lVert Rf-p\rVert_{2}^{2},$

which, due to the convexity of the norm
$\|x\|_{2}:=\sqrt{\sum_{i=1}^{n}x_{i}^{2}}$, can be solved efficiently to
global optimality with, e.g., gradient-based methods. In this work, we focus
on the realistic and very natural setting where only few projections of the
sample can be acquired. This is true, e.g., if projections are expensive,
time-consuming or if beam-sensitive materials are studied, which can only
withstand a certain X-ray or electron dose. Under these conditions, it is
crucial to take only a small number of projections, the fewer the better. Few
projections lead to $m<<n$. As a result, the linear equation system (1) is
highly underdetermined and, consequently, the information is highly
undersampled.
Established iterative reconstruction techniques to solve (2) are, e.g., the
simultaneous algebraic reconstruction technique (SART) [SK88] and the
simultaneous iterative reconstruction technique (SIRT) [Gil72, KS01].
Compressed sensing (CS), also named compressive sensing or compressive
sampling, is a thoroughly studied signal processing technique to achieve
improved reconstructions from such undersampled information by using prior
knowledge about the reconstructed signal. In particular, considering the
potential sparsity of the underlying signal in a known transform domain can
crucially improve the reconstruction as was demonstrated in the seminal paper
of Candès et al. [CRT06]. CS has a wide range of applications, see Qaisar et
al. [Qai+13], one of them is nanotomography. In that field, compressed sensing
electron tomography (CS-ET) [Lea+13] and total variation regularized discrete
algebraic reconstruction technique (TVR-DART)[ZPB15] are two exemplary
established CS-based algorithms. The mathematical reconstruction scheme of CS
is described in Section 2. This work aims at improving this scheme by
incorporating prior knowledge about the sample as additional constraint
classes. The goal is to reduce artifacts and achieve better reconstruction
quality in case of undersampled information. Two such classes of constraints
are individually derived and explained in Section 3. On the one hand, they
consist in the algebraic formulation of strict upper bounds for pixel
intensities that eliminate artifacts outside of the reconstructed samples. On
the other hand, we also introduce soft upper bounds for pixel intensities with
which deviations above a certain value can be penalized. This latter
constraint class is applicable to homogeneous samples that are composed of
only one material density and exploits this problem-specific fact to increase
reconstruction quality. Hence, the CS algorithm including these additional
constraints will be referred to as compressed sensing for homogeneous
materials (CSHM). For an optimal exploitation of these constraints, the
parameter selection is crucial and is therefore illustrated in Section 4 along
with our benchmark algorithms. To evaluate the practical impact of the new
constraint classes, they are applied to one simulated and two experimental
data sets, one acquired with electron tomography, the other with nano-CT (see
Section 5). The results are shown in Section 6, where comparisons to existing
state-of-the-art algorithms are performed. It turns out that the new models
lead to high-quality reconstructions within short time when compared to
existing models, in particular, if only a few number of projections are
available. A discussion of the results and some future extensions and
applications are given in Section 7.

##  2\. Compressed Sensing Framework

The following description of the CS framework is based on Leary et al.
[Lea+13]. The first major assumption in CS is that there exists a known basis
transformation $\psi\in\mathbb{R}^{n\times n}$ of the vector $f$, that is to
be recovered, to a vector $c\in\mathbb{R}^{m}$ with $m<<n$ nonzero
coefficients $c_{i}$. In this case, we say that the image $f$ is
_compressible_ in $\psi$. Consequently, it suffices to only gather information
in the compressed form, i.e., instead of recording samples of $f$ directly,
one records well-chosen samples $b_{i}$ of linear combinations
$b_{i}=\sum_{j=1}^{n}\phi_{ij}f_{j}$. Here, $\phi\in\mathbb{R}^{m\times n}$
denotes known sampling coefficients, which are often inherent to the
respective application, e.g., in a tomographic experiment, the measurements
$b$ denote the value of line integrals with respect to projected lines through
the sample, whereas the lines itself are defined by $\phi$.

The second assumption in CS is that in order to recover $f$ efficiently, one
requires the matrices $\psi$ and $\phi$ to be rather dissimilar. Similar
matrices lead to redundant constraints $c=\psi f,b=\phi f$ and add little
value by the measurements $b$ as most of the coefficients $b_{i}$ are close to
zero. To this end, one measures the dissimilarity or _incoherence_ of $\psi$
and $\phi$. In the present paper, as well as in other nanotomography
applications, incoherency with respect to the sensing matrix $R$ and the
transforming matrix $\psi$ is controlled by the measuring process. Using,
e.g., HAADF-STEM imaging, incoherency is ensured by the independence of the
different scanning position measurements [HRD96]. With these two requirements,
a reconstruction scheme, enforcing both sparsity in the transform domain and
data consistency with the measurements, can be formulated. In particular,
given the $l_{p}$ norm for $x\in\mathbb{R}^{k},p\geq 1$ defined by

$\lVert x\rVert_{p}:=\left(\sum\limits_{i=1}^{k}\lvert
x_{i}\rvert^{p}\right)^{\frac{1}{p}},$

we focus on the _Basis pursuit denoising_ (BPD) model introduced by Chen et
al. [CDS01] for image reconstruction. For a broader overview on BPD as well as
other image reconstruction models, we refer to the seminal book of Foucart and
Rauhut [FR12]. In BPD, data consistency is enforced by penalizing deviations
in the objective function, where $\lambda\geq 0$ weighs the sparsity

(3) $\displaystyle\min\limits_{f\geq 0}\quad\lVert
Rf-p\rVert_{2}^{2}+\lambda\lVert\psi f\rVert_{1}.$

Measuring the sparsity of $\psi f$ directly via the non-convex $l_{0}$ norm,
defined by $\lVert x\rVert_{0}:=\lvert\\{i:x_{i}\neq 0\\}\rvert$, generally
leads to a more complex problem (3). Instead, the $l_{1}$ norm is applied, to
ensure sufficient sparsity, ideally restricting $\psi f$ to at most as many
non-zero components as the number of projections $m$ [FR13]. The BPD problem
formulation (3) is widely used in the literature, e.g., in Lustig et al.
[LDP07] or in Block et al [BUF07].
The main focus of this work lies in the reconstruction of homogeneous
materials. The latter are of frequent interest in nanotomography, such as in
the 3D investigation of porous supports for applications in heterogeneous
catalysis [Wir+21], particle chromatography [Joh+18], fuel cells [Ven+17] or
battery electrodes [Tro+14]. For homogeneous materials, reconstruction
algorithms based on combinatorial optimization models have been presented in
Liers and Pardella [LP11]. Since homogeneous materials consist of only one
approximately constant density value (i.e., composition or material phase),
they provide additional information that is known beforehand. In particular,
being composed of only one density value, these reconstructions are expected
to exhibit either areas of said density or empty space (typically filled with
vacuum or air), with only a few sharp edges in between. This promotes sparsity
in its spatial finite differences. Exploiting this sparsity with compressed
sensing is called total variation (TV) minimization and is often applied in
the literature, e.g., by Sidky and Pan [SP08] or Leary et al. [Lea+13] It uses
the a priori knowledge that most pixels are surrounded by pixels of the same
grey level.
For a quadratic picture with $n=l^{2}$ pixels, the sparsifying TV transform
$\psi^{TV}\in\\{-1,0,1\\}^{2(n-l)\times n}$ is indirectly defined like in
Block et al. [BUF07] by:

(4) $\displaystyle\lVert
f\rVert_{\text{TV}}:=\lVert\psi^{TV}f\rVert_{1}:=\sum\limits_{j=1}^{n}\lVert\nabla
f_{j}\rVert_{1}=\sum\limits_{j=1}^{n}\lvert\nabla^{x}f_{j}\rvert+\lvert\nabla^{y}f_{j}\rvert,$

where for every $j\in\\{1,\ldots,n\\}$, we define

(5a) $\displaystyle\nabla f_{j}:=(\nabla^{x}f_{j},\,\nabla^{y}f_{j}),$ (5b)
$\displaystyle\nabla^{x}f_{j}:=\begin{cases}f_{j+1}-f_{j}&\text{ if $j\mod
l\neq 0$},\\\ 0&\text{ if $j\mod l=0$},\end{cases}$ (5c)
$\displaystyle\nabla^{y}f_{j}:=\begin{cases}f_{j+l}-f_{j}&\text{ if $j\leq
n-l$},\\\ 0&\text{ if $j>n-l$}.\end{cases}$

For solving the CS optimization problems, several efficient algorithms are
known, such as interior point methods (IPM), [Kim+07] TVR-DART [ZPB15] and CS-
ET. [Lea+13]

###  2.1. TVR-DART

We revisit the TVR-DART algorithm, developed by Zhuge et al., [ZPB15] which is
a compressed sensing technique specifically developed for samples consisting
of only a few different density values. We use the image reconstructions
computed by this algorithm as a benchmark for solutions to (14). It aims at
combining compressed sensing and discrete tomography. Homogeneous samples in
nanotomography are one of many applications. It starts with estimating
parameters by means of a SIRT reconstruction. Then the following optimization
problem is solved

(6)
$\displaystyle\min\limits_{\rho_{1},\ldots,\rho_{G},\tau_{1},\ldots,\tau_{G-1},f}$
$\displaystyle\quad\lVert R\,S(f,\mathcall{P},\mathcall{T})-p\rVert_{2}^{2}$
$\displaystyle+\lambda\sum\limits_{1\leq j\leq n}M_{\delta}\left(\nabla
S(\nabla f,\mathcall{P},\mathcall{T})\right).$

Here, $\rho_{1},\ldots,\rho_{G}$ are the unknown different material density
values of the sample and $\mathcall{P}:=\\{\rho_{1},\ldots,\rho_{G}\\}$. The
number of density values, $G\geq 1$, is assumed to be known in advance.
$\tau_{1},\ldots,\tau_{G-1}$ are the thresholds between the density values
($\mathcall{T}=\\{\tau_{1},\ldots,\tau_{G-1}\\}$), setting the turning points
of when a pixel density is pushed to either the lower or the higher density
value. This is done by applying a differentiable soft segmentation function
$S$ to each individual pixel density $f_{j}$. The function smoothly drives
them towards the values in $\mathcall{P}$. The optimization problem (6) is now
solved for the smooth image $S(f,\mathcall{P},\mathcall{T})$. Zhuge et al.
[ZPB15] chose $S$ as a smooth approximation of a piecewise constant staircase
function with values in $\mathcall{P}$ and jumps in $\mathcall{T}$. Instead of
using the TV-norm for the sparsity term as in usual TV minimization, a
differentiable norm is suggested using a Huber loss function $M_{\delta}$ with
$\delta=10^{-4}$.

Although the formulation is natural, it comes with the computational
complexity that (6) is an unconstrained continuous however non-convex
optimization problem. Zhuge et al. use an alternating heuristic minimization
procedure to find solutions in a short amount of time, starting with an
initial reconstruction produced by SIRT. While TVR-DART often appears
effective in practical scenarios, its authors acknowledge its divergence under
specific conditions, implying the absence of a guaranteed high-quality
solution. This behavior is common when aiming to solve non-convex optimization
problems through local heuristics. In order to overcome these drawbacks due to
nonconvexities, this work introduces convex modelling extensions that can be
used for roughly homogeneous materials, enabling the utilization of efficient
interior point methods.

##  3\. Incorporating Problem Specific Information

In this section, we introduce model extensions to increase the quality of the
reconstructed images. They counteract different kinds of artifacts that
typically appear in image reconstructions. In the following, the terms hard
and soft upper bounds are used. Hard bounds mean bounds on variables that are
enforced by a constraint and cannot be violated in a solution. Soft bounds, on
the contrary, are bounds on variables that are not enforced but promoted by
penalizing violations in the objective function weighted by an appropriately
chosen parameter.

###  3.1. Hard Upper Bound on $f$ in Image Reconstruction

With (1), i.e., assuming perfect projections without noise, it holds that

(7) $\displaystyle p_{i}=\sum\limits_{j=1}^{n}R_{ij}f_{j}$
$\displaystyle\forall 1\leq i\leq m,$

with $p\geq 0$, $R\geq 0$ and $f\geq 0$. This implies that if a pixel $j$ of
image $f$ gets passed by a projection ray $i$, meaning $R_{ij}>0$, its density
value $f_{j}$ satisfies $R_{ij}f_{j}\leq\sum_{j=1}^{n}R_{ij}f_{j}=p_{i}$ and
thus the following variable bound emerges:

(8) $\displaystyle f_{j}\leq\min\limits_{1\leq i\leq m:\,R_{ij}\neq
0}\,\,\frac{p_{i}}{R_{ij}}$ $\displaystyle\forall 1\leq j\leq n,$

We note that on the one hand, the bound $\frac{p_{i}}{R_{ij}}$ is loose for
rays $i$ that pass through multiple sample-containing pixels $j$. This then
causes the projection value $p_{i}=\sum_{j=1}^{n}R_{ij}f_{j}$ to be large.
However, on the other hand, it is advantageous for rays that hit no sample
pixels. Indeed, they have a rather small projection value $p_{i}$. As a
consequence, all pixels passed by such a ray are restricted by a very tight
bound (8). This nudges the algorithm towards reconstructions that correctly
consider little to no material at these pixels, even if that would be
beneficial for the overall error minimization. Hence, the algorithm provides
reduced undersampling artifacts (e.g., streaking artifacts) outside of the
reconstructed sample, assuming the sample does not fill out the whole field of
view of the microscope or detector and that the sample is surrounded only by
air or vacuum. This effect is observed in the results later in Section 6. We
note that the impact of (8) is dependent on the average noise level of the
background – the higher the noise level, the looser the bounds. Therefore, a
mean background noise subtraction preprocessing step should be applied before
reconstruction.

As already stated, equation (7) holds for many projection geometries, making
the proposed bounds (8) applicable to a broad range of different tomography
image reconstruction problems.

###  3.2. Soft Upper Bound on $f$ for Homogeneous Materials

Contrary to the hard upper bound (8), the following class of constraints
describes suitable bounds for sample pixels, assuming a constant material
density over the whole sample. In materials science, such homogeneous samples
occur frequently. For the following bounds, we assume such a homogeneous
sample with the single density value $\omega\in\mathbb{R}_{\geq 0}$. If
$\omega$ is not known in advance, it can be estimated. This can be done
efficiently by reconstructing the same sample at a lower resolution and
calculating a mean of all (bigger) positive values.

The idea of the constraints is to quadratically penalize the density values
that exceed $\omega$. Therefore, we define an auxiliary variable vector
$d\in\mathbb{R}^{n}$ and add the constraints

(9) $\displaystyle d_{j}\geq f_{j}-\omega$ $\displaystyle\forall 1\leq j\leq
n,$ (10) $\displaystyle d_{j}\geq 0$ $\displaystyle\forall 1\leq j\leq n,$

to the model. Additionally, we add the quadratic term

(11) $\displaystyle\mu\lVert
d\rVert_{2}^{2}=\mu\sum\limits_{j=1}^{n}d_{j}^{2}$

to the objective function of (3). The parameter $\mu$ has to be chosen
appropriately, which is further analyzed in Section 4.1.
As a result of the quadratic formulation, the higher the reconstructed
intensity of a pixel is above the density $\omega$, thus being a potential
artifact, the more it is penalized. This leads to a reconstructed image with
fewer large deviations above $\omega$, so that the density values tend to be
comparable over the whole sample. A possible downside of this model is that
projection artifacts, e.g., non-linear contrast effects like remaining Bragg
contrast contributions in HAADF-STEM imaging of crystalline specimen, which
otherwise would be reconstructed as a brighter spot in the material, could now
be wrongly distributed onto other pixels. Especially feature edges, e.g., pore
edges, are prone to receive such material, since in undersampled data, the
placement of edges is not always clear. The here proposed soft constraints can
thus potentially lead to decreasing pore sizes in porous specimen. However,
the quadratic nature of the penalty keeps this effect small. These effects can
be observed later in Section 6.
Next, we go one step further and discuss the usage of more advanced modelling
techniques from mathematical optimization. One could argue that, if
homogeneous samples contain only two pixel density values, $0$ and $\omega$,
density values of $f$ in the space between $0$ and $\omega$ have to be
penalized as well. In principle, this can be accomplished by the class of
constraints

(12) $\displaystyle d_{j}\geq\lvert f_{j}-\omega b_{j}\rvert$
$\displaystyle\forall 1\leq j\leq n,$ (13) $\displaystyle b_{j}\in\\{0,1\\}$
$\displaystyle\forall 1\leq j\leq n.$

Consequently, all pixel values are moved towards either $0$ or $\omega$, with
small differences allowed if they lead to better data consistency or smaller
spatial finite differences. Edge pixels, where only a fraction of a pixel
contains material and the rest is background, which are usually present in
experimental data, contradict the idea of (12) and (13). However, the higher
the resolution, the less proportion of such pixels exist, so this drawback can
be considered as neglectable for high-resolution images. Despite their
modelling strengths, including binary variables significantly increases the
problem’s complexity. These variable turn the problem into a mixed-integer
quadratic program (MIQP), which is $\mathrm{NP}$-hard in general, i.e., it is
widely expected to not be efficiently solvable. In particular, this makes the
problem considerably more difficult than the original convex quadratic
program.
Thus, in the remainder of this article, we focus on (9) and (10) and postpone
algorithms involving binary variables to future research. In the following, we
thus consider model

(14a) $\displaystyle\min\limits_{f}$ $\displaystyle\lVert
Rf-p\rVert_{2}^{2}+\lambda\lVert f\rVert_{TV}+\mu\lVert d\rVert_{2}^{2}$ (14b)
$\displaystyle d_{j}\geq f_{j}-\omega$ $\displaystyle\forall 1\leq j\leq n,$
(14c) $\displaystyle d_{j}\geq 0$ $\displaystyle\forall 1\leq j\leq n,$ (14d)
$\displaystyle f_{j}\leq\min\limits_{i\in A_{j}}\frac{p_{i}}{R_{ij}}$
$\displaystyle\forall 1\leq j\leq n,$ (14e) $\displaystyle f_{j}\geq 0$
$\displaystyle\forall 1\leq j\leq n.$

This problem will be referred to as compressed sensing for homogeneous
materials (CSHM) because it exploits the homogeneity of the sample.

##  4\. Algorithmic Details

In this work, we focus on two-dimensional (2D) image reconstructions. For a
parallel beam geometry and single tilt axis tomography experiments, 3D
problems can be split into independent 2D problems, making the proposed
optimization problem applicable by consecutively running it for 2D data.

The new proposed reconstruction scheme was implemented in Python. It utilizes
the ASTRA Toolbox, developed by Van Aarle et al. [Van+15] and the TVR-DART
Toolbox, developed by Zhuge et al. [ZPB15]. For solving the CS problem (14)
with hard constraints, a state-of-the-art primal-dual interior point method
(IPM), or barrier method [Gur23] (version 10.0.0) is used. A main advantage of
IPM consists in the fact that it can very efficiently solve large-sized,
linearly constrained optimization problems with convex quadratic objective
functions which is exactly the form of (14).

###  4.1. Selection of Parameters

The objective function parameters $\lambda$, the weight of the total
variation, and $\mu$, the penalty for exceeding the mean density $\omega$,
have to be chosen before solving (14).

Several proposals in the literature for choosing $\lambda$, e.g., by Chen et
al. [CDS01] or Jin and Rao [JR10], did not lead to satisfying results. For
reasons of simplicity, $\lambda$ was empirically chosen in the implementation
of this work. The projection data sets show that $\lambda$ can be chosen
independently of pixel size or number of projection angles and yet perform
well for all 2D slices of a data set. The fact that $\lambda$ does not have to
be increased with the number of projection angles or image pixels, despite the
increasing error minimization term, can be interpreted in the following way.
First, with more image pixels and thus higher resolution, or more information
about the picture through more projection angles, the quality of images
returned by only error minimization of (2) increases, and the undersampling
artifacts decrease. A regularizing term $\lambda\lVert f\rVert_{TV}$, which is
counteracting the undersampling artifacts, is needed less and less.
Furthermore, the two new constraints introduce an additional smoothing effect.
The soft upper bounds penalize bright spots that show bigger densities than
material density $\omega$ and push the exceeding material to other pixels.
Those other pixels are lying inside the sample because the hard upper bounds
enforce pixels outside of the sample to be $0$ (or nearly $0$). This smoothes
the reconstructed sample. In this implementation, the parameter was
empirically chosen as $\lambda=4000$ for the ET data set of the zeolite
particle and the simulated data set, and $\lambda=3$ for the nano-CT
acquisition of the copper microlattice.
The parameter $\mu$ in (14a) specifies the weight of the penalty for exceeding
$\omega$ with pixel densities $f_{j}$. Since more projection angles or a
higher pixel count generally does not lead to an increasing error term
$\|d\|_{2}^{2}$, but $\omega$-exceeding artifacts remain (in contrast to
undersampling artifacts), $\mu$ has to be scaled correspondingly.
In X-ray and electron microscopy imaging, $\omega$-exceeding artifacts may
occur due to non-linear contrast relationships between measured intensity and
the specimen’s mass and thickness (i.e., local mass attenuation coefficients)
caused by, e.g., non-linear mass-thickness contrast [Van+12] or remaining
Bragg scattering effects [Ven+13, Ven+14] in HAADF-STEM imaging.
Thus, in order to keep our penalty $\mu\|d\|_{2}^{2}$ relevant with increasing
problem size and corresponding increasing reconstruction error term
$\|Rf-p\|_{2}^{2}$, we scale $\mu$ with the problem size as follows:

(15) $\displaystyle\mu=5a\frac{l}{256},$

where $a$ specifies the number of projection angles and $l$ the number of
pixels per row of the reconstructed image. The constant factor of
$\frac{5}{256}$ was empirically chosen as it produced good results in
preliminary tests.

###  4.2. Reconstruction Benchmark and Hardware

To evaluate the quality of the developed constraints in CSHM, it is compared
to three existing reconstruction techniques, namely SIRT, CS and TVR-DART.
Here, CS refers to the solution of (3) with TV minimization of (4) and no
additional constraints, but using the same IPM as CSHM. The IPM is run with a
solution tolerance of $10^{-6}$. The number of iterations that the comparison
algorithms SIRT and TVR-DART are supposed to perform when reconstructing
images are set to $1000$ for SIRT and $250$ for TVR-DART. The initial SIRT
reconstruction, which TVR-DART uses as a starting point, also runs for $1000$
iterations. Note that TVR-DART can also terminate earlier by fulfilling a
convergence stopping criterion before reaching the maximum number of
iterations.
All following image reconstructions presented in Section 6 were performed on a
laptop computer with 16 GB RAM, an AMD Ryzen 9 5900HX with Radeon Graphics,
3301 MHz, 8 CPU cores, and 16 logical processors. For a GPU-accelerated
version of SIRT, a NVIDIA GeForce RTX 3050 Laptop GPU was used.

##  5\. Experimental Details

Image reconstructions considered in this article are performed on three data
sets: one simulated phantom object slice (Figure 4a) with added random Poisson
noise and two experimental data sets, one acquired by ET and one by nano-CT.
The ET data in Figure 2 comprises a macroporous MFI-type zeolite particle
synthesized by Machoke et al. [Mac+15]. A full tilt series, as shown in SI
Video 1, in a tilt-angle range of 180° with 1° tilt increment (in total 180
projections) of the particle on top of the plateau of a tomography tip (so-
called 360°-ET or on-axis ET) was acquired using a FEI Titan3 80-300
transmission electron microscope operated at an acceleration voltage of 200 kV
in HAADF-STEM imaging mode (image size 1024 pixels x 1024 pixels; pixel size
3.55 nm) and a Fischione Model 2050 On-Axis Rotation Tomography sample holder
(E.A. Fischione Instruments, Inc.). HAADF-STEM imaging assures an
approximately parallel beam geometry and a monotonous relationship between
measured intensities, sample mass (density), and thickness, so that the Radon
transform can be applied for reconstruction. For more experimental details,
please refer to Przybilla et al. [Prz+18].
The nano-CT data set in Figure 3 shows a copper microlattice. These copper
microlattices were manufactured using an additive micromanufacturing technique
based on localized electrodeposition in liquid (CERES system – Exaddon AG,
Switzerland) [Kan+23, Ram+22]. The electrochemical ink required for the
electrodeposition is supplied with a hollow atomic force microscope (AFM)
cantilever with a small orifice ($\textasciitilde$300nm) and the metal ions
get reduced on the conductive substrate that acts as the working electrode.
Copper microlattices were built in a voxel-by-voxel manner using precise piezo
positioners which move the AFM cantilever after each voxel deposition to a new
location based on the deflection of the AFM cantilever sensed using a laser
based optical tracking system.
The nano-CT experiment was performed using a ZEISS Xradia 810 Ultra
laboratory-scale X-ray microscope equipped with a 5.4 keV rotating anode Cr
source in absorption contrast and large-field-of-view (LFOV) mode. The LFOV
mode covers a field of view (FOV) of 65 µm x 65 µm with a spatial resolution
down to 150 nm and a pixel size of 63.89 nm. A tilt series with $40$
projections in a tilt-angle range of 180° (tilt increment 4.5°) and an
exposure time of 300 s/frame was acquired, as shown in SI Video 2. We assume
parallel beam geometry and a beam attenuation according to Beer-Lambert law,
so that the Radon transform can be applied for reconstruction using the
projection data $p=-\ln{I/I_{0}}$, where $I$ is the measured intensity and
$I_{0}$ is the unattenuated incident beam intensity. The nano-CT projection
images in Figure 3 are displayed with inverted contrast. Please refer to SI
Movie 2 for the animated video of the original absorption contrast nano-CT
tilt series.
The 360°-ET tilt series was aligned by cross correlation in FEI Inspect 3D
version 3.0. The nano-CT tilt series was aligned using the adaptive motion
compensation (AMC; based on Wang et al. [Wan+12]) procedure implemented in the
native ZEISS software (XMController).

Figure 2. Exemplary projections of experimental ET data from different
projection angles of a macroporous zeolite particle on the plateau of a
tomography tip: (a) 20°, (b) 90°, (c) 160°. See SI Video 1 for an animation of
the full tilt series.

Figure 3. Exemplary projections of experimental absorption-contrast nano-CT
data set from different projection angles of a copper microlattice structure:
(a) 0°, (b) 10°, (c) 30°. See SI Video 2 for an animation of the full tilt
series.

Figure 4. "Ground truth" images for following reconstructions: (a) simulated
phantom structure, (b) SIRT reconstruction with 1000 iterations using 180
projections in a 180° tilt-angle range of a macroporous zeolite particle, (c)
SIRT reconstruction with 1000 iterations using 40 projections in a 180° tilt-
angle range of a copper microlattice structure.

##  6\. Computational Results

As an objective quality measure for the results of the simulated data set
reconstructions, the relative mean error (RME)

(16) $\displaystyle\text{RME}:=\frac{\sum\limits_{j=1}^{n}\lvert
f_{j}-\tilde{f}_{j}\rvert}{\sum\limits_{j=1}^{n}\lvert\tilde{f}_{j}\rvert}$

can be calculated, like in Zhuge et al. [ZPB15], comparing the reconstructed
image $f$ and the original image without noise ($\tilde{f}$). Since there
exist no ground truth images for the experimental data, the RME cannot be
calculated there. Instead, one can calculate the raw data coverage (RDC)

(17)
$\displaystyle\text{RDC}:=\frac{\sum\limits_{i=1}^{\bar{m}}\lvert\bar{R}f^{*}-\bar{p}\rvert}{\sum\limits_{i=1}^{\bar{m}}\lvert\bar{p}\rvert}$

that displays how well the projection of a reconstruction $f^{*}$, attained
with (possibly limited) projection data $p$, aligns with all projection data
$\bar{p}$ of all $\bar{m}$ available angles. $\bar{R}$ is the corresponding
radon matrix. Note that the RDC is not as meaningful or reliable as the RME
concerning reconstruction quality. This is due to the fact that the noise in
the original projection data $\bar{p}$ is not known and therefore part of the
comparison. For a visual comparison, Figure 4 shows "ground truth" images of
the 2D slices that are reconstructed later in this section. Since there are no
ground truths for the experimental data, the standard image reconstruction
method SIRT is applied to the whole available projection data ($1024^{2}$
pixels each, $180$ angles in Figure 4b, $40$ angles in Figure 4c) to attain
comparison images.

###  6.1. Results

As compressed sensing is designed to solve image reconstruction problems with
undersampled data, meaning only a few projection angles are available or the
projection range contains a missing wedge, primarily those problem settings
are studied. The projection angles in the following reconstructions are always
equidistantly spread over the available angular range. That is, 0° to 179° if
there is no missing wedge, and $\frac{x}{2}\degree$ to
$(180-\frac{x}{2})\degree$ if there is a missing wedge of $x\degree$.

The CSHM technique is compared to the existing image reconstruction techniques
CS, SIRT, and TVR-DART. Except for SIRT, all algorithms are specifically
designed to solve problems with undersampled projection data. The evaluation
is separately performed for simulated, experimental, and missing wedge
projection data.

#### Simulated Data

First, we compare reconstructions of the simulated data set with $512^{2}$
pixels and only $5$ projection angles (no missing wedge). The results are
displayed in Figure 5, where corresponding running times and RMEs are noted
below the reconstructions.

(a) SIRT: $3$s, RME: $0.7220$

(b) CS: $150$s, RME: $0.0755$

(c) TVR-DART: $184$s, RME: $0.1618$

(d) CSHM: $56$s, RME: $0.0413$

Figure 5. Results for simulated data with $512^{2}$ pixels and $5$ projection
angles.

(a) SIRT: $5$s, RME: $0.3209$

(b) CS: $307$s, RME: $0.0525$

(c) TVR-DART: $20$s, RME: $0.0420$

(d) CSHM: $111$s, RME: $0.0301$

Figure 6. Results for simulated data with $512^{2}$ pixels and $20$ projection
angles.

The small number of only $5$ angles was chosen because the shape of the object
is very simple and $5$ angles are already enough for good reconstructions in
CS approaches. For SIRT, however, this number of projection angles is too
small. Its reconstruction in Figure 5(a) is far from accurate, which is
emphasized by an RME of $0.7220$. Both CS and CSHM in Figures 5(b) and 5(d)
produce good results with a qualitative advantage of CSHM. This is also
reflected in the relative mean errors, the RME of CSHM is about half the one
of CS (0.0413 vs. 0.0755). The TVR-DART result in Figure 5(c) is better than
the one from SIRT, but worse than CS and CSHM reconstructions. Concerning the
running times, SIRT is the fastest approach. CSHM does not only produce
qualitatively better results than CS, but is also about three times faster.
This is owed to the hard upper bounds eliminating a great number of variables.
The TVR-DART runs the longest, which, in combination with its reconstruction
quality, indicates that it did not converge and iterated until the given
maximum number of iterations.

In Figure 6, the same reconstructions are performed with $20$ projection
angles. The SIRT reconstruction in Figure 6(a) is still very blurry with a
large RME of $0.3209$. CSHM is again about three times faster than CS while
producing an about two times better RME. Comparing these two reconstructions
shows that the CS reconstruction in Figure 6(b) exhibits noisy spots, both
dark and bright, while the CSHM reconstruction in Figure 6(d) only shows the
dark spotted noise. The brighter spots are weakened or even eliminated by the
soft upper bounds by penalization. The TVR-DART reconstruction in Figure 6(c)
also only exhibits darker spots but it also shows a lot of spotted artifacts
around the sample. Those are less for CS and nearly absent for CSHM due to the
hard upper bounds for pixels outside of the sample. The RME of TVR-DART,
$0.0420$, is better than the $0.0525$ for CS, but still worse than the
$0.0301$ of CSHM. The much faster running time of TVR-DART compared with the
reconstruction with only $5$ projection angles suggests that it converged.

Table 1. Results for simulated projection data with $256^{2}$ pixels and different amount of projection angles. The best RME value for each number of projection angles is highlighted in bold. |  | Running Time (in seconds) |  | Relative Mean Error (RME)
---|---|---|---|---
Proj. Angles |  | SIRT | CS | TVR-DART | CSHM |  | SIRT | CS | TVR-DART | CSHM
5 |  | 3 | 20 | 61 | 16 |  | 0.5605 | 0.0601 | 0.1537 | 0.0274
10 |  | 3 | 30 | 8 | 19 |  | 0.3382 | 0.0411 | 0.0375 | 0.0262
15 |  | 3 | 49 | 8 | 20 |  | 0.2511 | 0.0414 | 0.0365 | 0.0269
20 |  | 3 | 57 | 4 | 22 |  | 0.2124 | 0.0421 | 0.0314 | 0.0247
30 |  | 3 | 114 | 4 | 24 |  | 0.1686 | 0.0468 | 0.0282 | 0.0243
45 |  | 3 | 227 | 7 | 54 |  | 0.1390 | 0.0539 | 0.0278 | 0.0240
60 |  | 3 | 478 | 6 | 92 |  | 0.1279 | 0.0578 | 0.0236 | 0.0232
90 |  | 3 | 850 | 6 | 428 |  | 0.1200 | 0.0628 | 0.0221 | 0.0223
180 |  | 3 | - | 13 | 2897 |  | 0.1244 | - | 0.0225 | 0.0202

To get some more insights, a comparison in RME and running time is displayed
in Table 1 for a smaller pixel count of $256^{2}$ and different numbers of
projection angles. The best RME of every row is displayed in bold. SIRT solved
every problem within $3$ seconds. Both CS and CSHM exhibit strictly
monotonically increasing running times, which is expected for increasing
problem sizes solved with an IPM. CSHM is consistently faster than CS. The
biggest problem with $180$ projection angles could not be solved anymore with
CS because the computer ran out of memory. The running times of TVR-DART are
very low for all reconstructions except for the one using projections from $5$
angles, indicating that the algorithm failed to converge with only $5$ angles
of information but managed to do so for every other number. This claim is
supported by comparing the relative RME differences between $5$ and $10$
projection angles for all algorithms, showing by far the largest for TVR-DART.

Despite a monotonically decreasing RME for SIRT (except for the last value),
the lowest RME of SIRT, $0.1200$, is bigger than the RME that all the other
algorithms achieved with any number of projection angles, except TVR-DART with
$5$ angles. The RME of CS decreases from $5$ to $10$ angles and then slightly
increases in every step. This might be caused by the chosen $\lambda$ (see
Section 4.1) not scaling well for an increasing number of projections without
the additional smoothing effect by the new constraints. The RME of TVR-DART
decreases in every step except for the last one where it increases slightly.
Despite two small exceptions, the RME value for CSHM overall decreases with an
increasing number of projection angles and reaches a minimum of $0.0202$ at
$180$ angles. This is the lowest RME any of the tested algorithms could
achieve for any number of projection angles.

In eight out of nine cases, CSHM reached the lowest RME out of the four
algorithms, only one times slightly outperformed by TVR-DART for $90$
projection angles. It produced a better reconstruction with only $5$ angles
than CS or SIRT did for any given number of projections while maintaining a
running time $\leq 60$ seconds for problems with $45$ or fewer angles. Note
that the RME of CSHM acquired with $5$ projection angles is already very small
and does not improve by much anymore with more projection angles. This proves
the strength and high efficiency of CSHM for sparse data sets. Indeed, it
yields very good and high-quality results with relatively low running times in
particular in the situation when reconstructing from only a few number of
projections. Having only a low number of projections at hand is a typical
situation in many settings. This advantage is further fostered in the
reconstructions of the following experimental data sets.

#### Experimental Data

Since the experimental data sets do not have a ground truth image, the
following reconstructions are only compared visually and by the RDC, a quality
measure that is not robust against artifacts appearing in the projection data
before reconstructing.

(a) SIRT: 4s, RDC: $18.0740$

(b) CS: 262s, RDC: $8.0577$

(c) TVR-DART: 65s, RDC: $11.3882$

(d) CSHM: 94s, RDC: $8.4071$

Figure 7. Results for experimental data set of a zeolite with $512^{2}$ pixels
and $20$ projection angles.

(a) SIRT: $4$s, RDC: $10.8665$

(b) CS: $588$s, RDC: $7.0235$

(c) TVR-DART: $75$s, RDC: $9.2214$

(d) CSHM: $129$s, RDC: $7.8802$

Figure 8. Results for experimental data set of a zeolite with $512^{2}$ pixels
and $30$ projection angles.

(a) SIRT: 2s, RDC: $30.7225$

(b) CS: 1215s, RDC: $31.2364$

(c) TVR-DART: 52s, RDC: $31.5497$

(d) CSHM: 640s, RDC: $30.5611$

Figure 9. Results for experimental data set of a Cu-microlattice with
$512^{2}$ pixels and $40$ projection angles.

In Figure 7, reconstructions of the experimental electron tomography
projection data of the zeolite particle sample with SIRT, CS, TVR-DART and
CSHM are displayed. The images contain $512^{2}$ pixels and $20$ projection
angles (no missing wedge) were used. Below the reconstructions, the
corresponding running times and RDCs are noted. The running times are
comparable to the ones from the simulated data set with $5$ projection angles
in Figure 5. SIRT is by far the fastest algorithm, CSHM is considerably faster
than CS, and TVR-DART takes relatively long, if not as long as for the
simulated data set. This indicates that TVR-DART converged this time, but not
very quickly. In terms of quality, the SIRT reconstruction in Figure 7(a) is
very blurry and exhibits a lot of artifacts. The TVR-DART reconstruction in
Figure 7(c) is less blurry but shows frequently alternating brightness inside
the material, which should not be the case for a homogeneous sample. The
reconstruction from CS in Figure 7(b) has considerably less bright spots than
TVR-DART but a few of those artifacts can be observed here, as well. CSHM
produces the best-looking reconstruction for the zeolite in Figure 7(d). It
displays almost no artifacts outside of the sample, due to the hard upper
bounds, and has a nearly constant material density over the sample, due to the
soft upper bounds. The shape of the pores in the CSHM reconstructions might be
a little more compressed than in reality, as, e.g., the indentation on the
upper left side of the particle in Figure 7(d) appears a little pressed
together compared to the $180$ angles SIRT reconstruction in Figure 4b. The
corresponding raw data coverage values support the visual evaluation that the
SIRT and TVR-DART results are worse than CS and CSHM. However, the RDC of the
CS result with $8.0577$ is slightly better than the RDC of $8.4071$ for CSHM,
despite CSHM being visually the better reconstruction. This may be owed to
noise during the projection data acquisition that is reconstructed with CS but
penalized and lessened by CSHM. The RDC is therefore only a rough quality
measure here.

In Figure 8, the same reconstructions are performed with $30$ projection
angles. The SIRT reconstruction in Figure 8(a) improved a lot, exhibiting less
blur and less artifacts outside of the sample. The CS and CSHM reconstructions
improved only by small amounts, the pore edges are a little sharper than
before. The TVR-DART reconstruction in Figure 8(c) is very different than for
$20$ angles in Figure 7(c). Now there seems to exist only one density value
over the whole sample, but a lot of small black spots, pepper noise, can be
found across the whole sample. The RDC values for all algorithms improved.
Again, the raw data coverage value of the CS reconstruction is better than the
one of the CSHM reconstruction, whereas visually, the one from CSHM is a more
realistic image of a homogeneous sample.
Figure 9 compares the different reconstructions of an exemplary slice from the
copper microlattice nano-CT tilt series in absorption contrast from Figure 3
using 40 projection angles. From a qualitative perspective, CSHM and CS
provide the best reconstructions with no streaking artifacts in the inside of
the lattice structure and especially also not in the outside regions due to
the hard upper bound. The bright main copper strut features appear more
homogeneous in intensity for CSHM due to the soft upper bound, with slightly
higher intensity variations in the CS slice. Although TV-DART provides the
better result, both SIRT and TVR-DART reconstructions show artifacts around
the outside of the specimen, and, in the sample interior, remaining streaking
artifacts and stronger contrast variations are observed. Due to the increased
amount of projection angles, also the running time of CS and CSHM are
increased, whereas SIRT is the fastest, and also TVR-DART seems to have
converged rather quickly. With respect to the RDC, all four reconstructions
exhibit similar values, with CSHM just showing the lowest RDC.

#### Missing Wedge

In some nanotomography applications, the angular tilt range of the projections
is limited due to, e.g., shadowing of the sample holder, leading to so-called
missing-wedge reconstruction artifacts. [Lea+13] In Figure 10, CSHM is
compared to SIRT, CS, and TVR-DART on the simulated and the zeolite particle
ET data sets with $512^{2}$ pixels and a missing wedge of $60\degree$. The
respective projection angles are equidistantly spread over the remaining
angular range of $120\degree$. Beneath the reconstructed slices, the
corresponding RME/RDC values and running times are displayed.

(a) SIRT, 4s, RME: $0.5775$

(b) CS, 229s, RME: $0.0850$

(c) TVR-DART, 193s, RME: $0.2164$

(d) CSHM, 69s, RME: $0.0495$

(e) SIRT, 3s, RDC: $21.5546$

(f) CS, 329s, RDC: $12.1703$

(g) TVR-DART, 62s, RDC: $15.9070$

(h) CSHM, 92s, RDC: $12.0015$

Figure 10. Results for a missing wedge of $60\degree$ with $512^{2}$ pixels
and $16$ projection angles for the simulated data set (a)-(d) and $21$
projection angles for the ET data set (e)-(h).

The magnitude of the running times of the algorithms are comparable to the
cases without missing wedge. CSHM outperforms the other algorithms for the
simulated data set (Figures 10(a)-10(d)), both visually and regarding to the
relative mean error. For the zeolite, CSHM also reaches the best RDC value,
while being visually the best reconstruction. The typical elongation of
reconstructions with missing wedge in vertical direction can be seen in both
CSHM reconstructions of Figure 10(d) and Figure 10(h), but the effect is
relatively small, especially when compared to the SIRT reconstruction. In
Figure 10(h), the pore at the lower part of the zeolite is most likely a
distorted and erroneously reconstructed indentation with merged borders, which
can be deduced from a comparison with the $180$ angles SIRT reconstruction
without missing wedge in Figure 4b. This merging of pore walls cannot be
observed in the full tilt-angle range SIRT reconstruction, nor in the missing-
wedge SIRT, CS, and TVR-DART reconstructions in Figures 10(e)-10(g). However,
in those reconstructions, the missing-wedge effect causes a much stronger
erroneous merging of pores (instead of pore walls), which are separate in
Figure 4b, compared to CSHM.

##  7\. Conclusion

The previous section provided results obtained with the two newly developed
constraints of this work. CS was combined with both hard and soft upper
bounds, resulting in the here-called compressed sensing for homogeneous
materials (CSHM) algorithm. CSHM was compared to SIRT, CS, and TVR-DART on
simulated and experimental two-dimensional projection data from different
nanotomography techniques, namely 360°-ET and nano-CT. Different problem
settings were studied, including problems with a missing wedge.

Regarding reconstruction quality, the results of CSHM are very convincing and
of high quality. In almost all cases for the simulated data set, CSHM achieved
the lowest relative mean error and was only one time slightly outperformed by
TVR-DART. CSHM also showed significantly decreased running times compared to
CS, when solved with the same IPM. Both artifacts outside and inside of the
reconstructed samples were reduced in comparison to CS, especially regarding
the intensity of unexpected brighter spots inside and streaking artifacts
outside of the samples. The artifacts outside of the samples were reduced by
the hard upper bounds. Bounding variables to (near) $0$ lets the solver
eliminate a lot of variables resulting in faster solving and less required
memory than without the constraint. This effect relies on a preprocessing step
of subtracting the mean background noise, to have a background intensity near
$0$. The preprocessing step can be applied, if a mean background noise value
can be estimated, e.g., if there are rays known that only pass through vacuum
or air and do not hit the sample. The hard bounds are not only applicable to
homogeneous samples, but to a much broader range of image reconstruction
problems. They are especially efficient, when the sample is surrounded by
vacuum or air. The artifacts inside the sample are reduced by the soft upper
bounds, which are only applicable to homogeneous samples. They penalize larger
deviations above the calculated homogeneous material density $\omega$, leading
to the density of the reconstructed sample being closer to homogeneity. A
slight drawback of CSHM was shown in form of a decrease in feature sizes,
e.g., pore sizes, by distributing exceeding material on the feature edges in
the CSHM reconstructions. However, this effect is very small.

An advantage of CSHM compared to TVR-DART is the consistency of good
solutions. From a qualitative perspective, TVR-DART performs very well for
some instances, possibly better than CSHM. However, for some other instances
it returns bad image reconstructions, as can be observed in Figure 7(c). CSHM
however, performs consistently well, even for very few projection angles.

Missing-wedge problem reconstructions of CSHM show typical missing-wedge
artifacts like elongation and the merging of features that are separate in
reconstructions without missing wedge. Similar effects are observed for
reconstructions of all the comparing algorithms. In comparison, CSHM exhibits
the fewest merging of pores and an elongation similar to CS but slightly worse
than TVR-DART.

Using more projection angles or a higher pixel count increases the running
time for CSHM. Then, TVR-DART might become attractive with shorter running
times, as the increment in quality of CSHM opposite TVR-DART might not be
worth the large increment in running time anymore. However, as Table 1
indicates, the CSHM reconstructions do not improve by much anymore for an
increasing number of projection angles, compared to a reconstruction with only
a few angles, which is already very accurate. This justifies reconstructing
images with only a few projection angles, keeping the running time within an
acceptable range. When reconstructing with only very little information (very
few projection angles), CSHM outperformed TVR-DART in both running time and
quality, as TVR-DART did not converge to a good reconstruction. Then, CSHM
represents the most efficient choice for image reconstruction. To decrease
high running times, existing work can be applied to the used IPM. Multiple
approaches were proposed in the literature, e.g., using a matrix-free IPM to
speed up the search step of the Newton method, developed by Fountoulakis et
al. [FGZ14], or accelerating the IPM by performing parallel computations on
the GPU, where promising results were obtained by Lee et al. [Lee+22].

In conclusion, CSHM represents a competitive alternative to existing
algorithms for smaller and medium-sized image reconstruction problems of
homogeneous samples. It provides high quality results for highly undersampled
data in nanotomography, including missing-wedge problems. Non-homogeneous
sample reconstructions can still benefit from the hard upper bounds that
reduce artifacts around the studied objects. Furthermore, an extension of the
soft upper bounds to a version that can handle samples composed of few, but
more than one, density values is thinkable. This is a question for future
work.

## Author Contributions

Sebastian Kreuz: Methodology, Software, Investigation, Writing - Original
Draft. Benjamin Apeleo Zubiri: Conceptualization, Methodology, Validation,
Writing - Review & Editing. Silvan Englisch: Validation, Data Curation. Sung-
Gyu Kang: Writing - Review & Editing. Rajaprakash Ramachandramoorthy: Writing
- Review & Editing. Erdmann Spiecker: Writing - Review & Editing. Frauke
Liers: Conceptualization, Methodology, Writing - Review & Editing. and Jan
Rolfes: Conceptualization, Methodology, Writing - Review & Editing.

## Data availability

The tilt series raw data and the Python scripts of the implemented algorithms
are publicly available on Zenodo DOI: 10.5281/zenodo.10283560.

## Conflicts of interest

There are no conflicts to declare.

## Acknowledgements

We acknowledge funding by the Deutsche Forschungsgemeinschaft (DFG, German
Research Foundation) - Project-ID 416229255 - SFB 1411. We thank Alexander
Götz for the help with the acquisition of the nano-CT data set.

## References

  * [Ape+21] Benjamin Apeleo Zubiri et al.  “Correlative Laboratory Nano-CT and 360 Electron Tomography of Macropore Structures in Hierarchical Zeolites”  In _Advanced Materials Interfaces_ 8.4, Wiley Online Library, 2021, pp. 2001154
  * [BUF07] Kai Tobias Block, Martin Uecker and Jens Frahm  “Undersampled radial MRI with multiple coils. Iterative image reconstruction using a total variation constraint”  In _Magnetic Resonance in Medicine: An Official Journal of the International Society for Magnetic Resonance in Medicine_ 57.6 Wiley Online Library, 2007, pp. 1086–1098
  * [Bur+14] TL Burnett et al.  “Correlative tomography”  In _Scientific reports_ 4.1 Nature Publishing Group UK London, 2014, pp. 4711
  * [CDS01] Scott Shaobing Chen, David L Donoho and Michael A Saunders  “Atomic decomposition by basis pursuit”  In _SIAM_ 43.1 SIAM, 2001, pp. 129–159
  * [CRT06] Emmanuel J Candès, Justin K Romberg and Terence Tao  “Stable signal recovery from incomplete and inaccurate measurements”  In _Communications on Pure and Applied Mathematics: A Journal Issued by the Courant Institute of Mathematical Sciences_ 59.8 Wiley Online Library, 2006, pp. 1207–1223
  * [FGZ14] Kimon Fountoulakis, Jacek Gondzio and Pavel Zhlobich  “Matrix-free interior point method for compressed sensing problems”  In _Mathematical Programming Computation_ 6 Springer, 2014, pp. 1–31
  * [FR12] Simon Foucart and Holger Rauhut  “A Mathematical Introduction to Compressive Sensing”  Springer, 2012, pp. 58
  * [FR13] Simon Foucart and Holger Rauhut  “A Mathematical Introduction to Compressive Sensing”  Birkhaeuser New York, 2013, pp. 56–57
  * [Gil72] Peter Gilbert  “Iterative methods for the three-dimensional reconstruction of an object from projections”  In _Journal of theoretical biology_ 36.1 Elsevier, 1972, pp. 105–117
  * [Gur23] Gurobi Optimization, LLC  “Gurobi Optimizer Reference Manual”, 2023  URL: https://www.gurobi.com
  * [HRD96] P Hartel, H Rose and C Dinges  “Conditions and reasons for incoherent imaging in STEM”  In _Ultramicroscopy_ 63.2 Elsevier, 1996, pp. 93–114
  * [Joh+18] TF Johnson et al.  “Three dimensional characterisation of chromatography bead internal structure using X-ray computed tomography and focused ion beam microscopy”  In _Journal of Chromatography A_ 1566 Elsevier, 2018, pp. 79–88
  * [JR10] Yuzhe Jin and Bhaskar D Rao  “Algorithms for robust linear regression by exploiting the connection to sparse signal recovery”  In _2010 IEEE International Conference on Acoustics, Speech and Signal Processing_ , 2010, pp. 3830–3833
  * [Kan+23] Sung-Gyu Kang et al.  “Fabrication and extreme micromechanics of additive metal microarchitectures”  In _arXiv preprint arXiv:2311.14018_ , 2023
  * [Kim+07] Seung-Jean Kim et al.  “An interior-point method for large-scale $l_{1}$ -regularized least squares”  In _IEEE Journal of Selected Topics in Signal Processing_ 1.4 IEEE, 2007, pp. 606–617
  * [KS01] Avinash C Kak and Malcolm Slaney  “Principles of computerized tomographic imaging”  SIAM, 2001
  * [LDP07] Michael Lustig, David Donoho and John M Pauly  “Sparse MRI: The application of compressed sensing for rapid MR imaging”  In _Magnetic Resonance in Medicine: An Official Journal of the International Society for Magnetic Resonance in Medicine_ 58.6 Wiley Online Library, 2007, pp. 1182–1195
  * [Lea+13] Rowan Leary, Zineb Saghi, Paul A Midgley and Daniel J Holland  “Compressed sensing electron tomography”  In _Ultramicroscopy_ 131 Elsevier, 2013, pp. 70–91
  * [Lee+22] Sanghyeon Lee et al.  “GPU-accelerated PD-IPM for real-time model predictive control in integrated missile guidance and control systems”  In _Sensors_ 22.12 MDPI, 2022, pp. 4512
  * [LP11] Frauke Liers and Gregor Pardella  “Simplifying maximum flow computations: The effect of shrinking and good initial flows”  In _Discret. Appl. Math._ 159.17, 2011, pp. 2187–2203  DOI: 10.1016/j.dam.2011.06.030
  * [Mac+15] Albert G Machoke et al.  “Micro/macroporous system: MFI-type zeolite crystals with embedded macropores”  In _Advanced Materials_ 27.6 Wiley Online Library, 2015, pp. 1066–1070
  * [Prz+18] Thomas Przybilla et al.  “Transfer of Individual Micro-and Nanoparticles for High-Precision 3D Analysis Using 360° Electron Tomography”  In _Small Methods_ 2.1 Wiley Online Library, 2018, pp. 1700276
  * [Qai+13] Saad Qaisar et al.  “Compressive sensing: From theory to applications, a survey”  In _Journal of Communications and Networks_ 15.5 KICS, 2013, pp. 443–456
  * [Ram+22] Rajaprakash Ramachandramoorthy et al.  “Anomalous high strain rate compressive behavior of additively manufactured copper micropillars”  In _Applied Materials Today_ 27 Elsevier, 2022, pp. 101415
  * [SK88] Malcolm Slaney and AC Kak  “Principles of computerized tomographic imaging”  IEEE press, 1988
  * [SP08] Emil Y Sidky and Xiaochuan Pan  “Image reconstruction in circular cone-beam computed tomography by constrained, total-variation minimization”  In _Physics in Medicine & Biology_ 53.17 IOP Publishing, 2008, pp. 4777–4807
  * [Tak+20] Satoshi Takeya et al.  “X-Ray attenuation and image contrast in the X-ray computed tomography of clathrate hydrates depending on guest species”  In _Physical Chemistry Chemical Physics_ 22.47 Royal Society of Chemistry, 2020, pp. 27658–27665
  * [Tro+14] Panagiotis Trogadas et al.  “X-ray micro-tomography as a diagnostic tool for the electrode degradation in vanadium redox flow batteries”  In _Electrochemistry Communications_ 48 Elsevier, 2014, pp. 155–159
  * [Van+12] W Van den Broek et al.  “Correction of non-linear thickness effects in HAADF STEM electron tomography”  In _Ultramicroscopy_ 116 Elsevier, 2012, pp. 8–12
  * [Van+15] Wim Van Aarle et al.  “The ASTRA Toolbox: A platform for advanced algorithm development in electron tomography”  In _Ultramicroscopy_ 157 Elsevier, 2015, pp. 35–47
  * [Ven+13] Singanallur V Venkatakrishnan et al.  “A model based iterative reconstruction algorithm for high angle annular dark field-scanning transmission electron microscope (HAADF-STEM) tomography”  In _IEEE Transactions on Image Processing_ 22.11 IEEE, 2013, pp. 4532–4544
  * [Ven+14] SV Venkatakrishnan et al.  “Model-based iterative reconstruction for bright-field electron tomography”  In _IEEE Transactions on Computational Imaging_ 1.1 IEEE, 2014, pp. 1–15
  * [Ven+17] Senthil Velan Venkatesan, Mohamed El Hannach, Steven Holdcroft and Erik Kjeang  “Probing nanoscale membrane degradation in fuel cells through electron tomography”  In _Journal of Membrane Science_ 539 Elsevier, 2017, pp. 138–143
  * [Wan+12] Jun Wang et al.  “Automated markerless full field hard x-ray microscopic tomography at sub-50 nm 3-dimension spatial resolution”  In _Applied Physics Letters_ 100.14, 2012, pp. 143107  DOI: 10.1063/1.3701579
  * [Whi+19] Robin T White et al.  “Correlative X-ray tomographic imaging of catalyst layer degradation in fuel cells”  In _Journal of The Electrochemical Society_ 166.13 IOP Publishing, 2019, pp. F914
  * [Wir+21] Janis Wirth et al.  “Unraveling Structural Details in Ga-Pd SCALMS Systems Using Correlative Nano-CT, 360° Electron Tomography and Analytical TEM”  In _Catalysts_ 11.7 MDPI, 2021, pp. 810
  * [Wit+21] Philip J Withers et al.  “X-ray computed tomography”  In _Nature Reviews Methods Primers_ 1.1 Nature Publishing Group UK London, 2021, pp. 18
  * [ZPB15] Xiaodong Zhuge, Willem Jan Palenstijn and Kees Joost Batenburg  “TVR-DART: A more robust algorithm for discrete tomography from limited projection data with automated gray value estimation”  In _IEEE Transactions on Image Processing_ 25.1 IEEE, 2015, pp. 455–468