# Deep Unfolding of the DBFB Algorithm with Application to ROI CT Imaging with
Limited Angular Density

Marion Savanier, Emilie Chouzenoux, Jean-Christophe Pesquet, and Cyril Riddell
This work was supported by the European Research Council Starting Grant
MAJORIS ERC-2019-STG-850925, the ANRT CIFRE Convention 2018/1587, and the ANR
Research and Teaching Chair in Artificial Intelligence BRIDGEABLE.E.
Chouzenoux, J.-C. Pesquet and M. Savanier are with Univ. Paris-Saclay,
CentraleSupélec, CVN, Inria, Gif-sur-Yvette, France (e-mail:
firstname.lastname@centralesupelec.fr).M. Savanier and C. Riddell are with GE
Healthcare, Buc, France (e-mail<EMAIL_ADDRESS>cyril.riddell@ge.com).

###### Abstract

This paper addresses the problem of image reconstruction for region-of-
interest (ROI) computed tomography (CT). While model-based iterative methods
can be used for such a problem, their practicability is often limited due to
tedious parameterization and slow convergence. In addition, inadequate
solutions can be obtained when the retained priors do not perfectly fit the
solution space. Deep learning methods offer an alternative approach that is
fast, leverages information from large data sets, and thus can reach high
reconstruction quality. However, these methods usually rely on black boxes not
accounting for the physics of the imaging system, and their lack of
interpretability is often deplored. At the crossroads of both methods,
unfolded deep learning techniques have been recently proposed. They
incorporate the physics of the model and iterative optimization algorithms
into a neural network design, leading to superior performance in various
applications. This paper introduces a novel, unfolded deep learning approach
called U-RDBFB designed for ROI CT reconstruction from limited data. Few-view
truncated data are efficiently handled thanks to a robust non-convex data
fidelity function combined with sparsity-inducing regularization functions.
Iterations of a block dual forward-backward (DBFB) algorithm, embedded in an
iterative reweighted scheme, are then unrolled over a neural network
architecture, allowing the learning of various parameters in a supervised
manner. Our experiments show an improvement over various state-of-the-art
methods, including model-based iterative schemes, deep learning architectures,
and deep unfolding methods.

###### Index Terms:

region-of-interest, computed tomography, angular sub-sampling, deep unfolding,
forward-backward, iterative reweighted scheme

##  I Introduction

CT imaging is commonly used for diagnostic purposes and image guidance in
interventional radiology and surgery. In most interventions, such as follow-up
examinations of deployed stents and needle biopsies, only a small region of
the patient is of interest. Irradiating only the ROI by X-rays involves
focusing the X-ray beam with collimation techniques using radio-opaque blades
before passing through the patient [1]. A focused irradiation results in a
substantial reduction in patient dose [2]. The measurements (or projections)
obtained from such acquisitions are truncated. The inverse problem of
reconstructing an ROI from a set of truncated projections is an ill-posed
problem often combined with angular sub-sampling for fast acquisitions.
Previous research demonstrated the ability of model-based iterative
reconstruction (MBIR) for CT reconstruction from few-view measurements [3, 4].
MBIR implements nonlinear iterative algorithms aiming at minimizing a
penalized cost function. In CT, most works rely on total variation (TV)
regularization [5]. Although they achieve a reduction in angular and cone-beam
sub-sampling artifacts compared to analytical reconstruction methods,
truncated projections still challenge MBIR [6].
To implement MBIR, one must choose the reconstruction grid, i.e., the support
of the reconstructed area. When the grid matches the support of the ROI, the
data will not agree with the reprojection of the ROI. When the reconstruction
grid includes the support of the entire object [7], the reconstruction becomes
computationally expensive and less stable due to the increase of unknowns for
the same amount of data. Truncated data only allows a rough estimation of the
exterior anatomical background, which can never be clinically acceptable. One,
therefore, often considers an intermediate smaller grid size, bigger than the
ROI but possibly less than the entire object [8]. This achieves a faster and
more stable ROI reconstruction. An alternative is to reconstruct the entire
object with large voxels and then subtract the reprojection of the exterior to
the data before reconstructing only the ROI from the subtracted data [9, 10].
Both approaches produce a low-frequency approximation of the exterior of the
ROI. When dense objects such as cables or needles are located outside of the
ROI and, more importantly, outside of the chosen reconstruction grid, the
approximation is poor, and reprojection of the objects increases the error in
the data fidelity term. This means that the data should not be trusted equally
but through a statistical analysis different from measurement noise. When
angularly sub-sampled, they still produce streak artifacts over the ROI [6].
Thus, designing MBIR methods for ROI imaging with a low computational cost is
still challenging. Preconditioning is often too restrictive for CT
reconstruction when proximal algorithms are used [11]; hence heuristics have
been proposed to accelerate TV-based methods. In [6], for example, substantial
acceleration of the forward-backward algorithm is reached thanks to two
modifications. The first modification uses the ramp filter to weigh their
least-squares data fidelity term. Each gradient step then performs an
approximate inversion of the forward projection, and each proximity step acts
as a post-processing filter. The use of an approximate inverse of the
measurement operator in the data fidelity term is also advocated in [12, 13,
14]. The second modification is the linear decrease of the regularization
parameter’s strength along with the iterations. Despite these modifications,
the proximity operator of TV regularization does not have a closed form, so
the method still requires many sub-iterations.
Convolutional neural networks (CNN) are an attractive alternative to MBIR due
to their increased expressivity and fast inference. CNNs, and in particular
the U-net [15], have already been used for removing sub-sampling streaks in CT
reconstructions obtained using analytical methods [16, 17]. However, there are
concerns about the lack of guarantees and capacity for generalization of post-
processing CNNs, because these networks do not ensure data consistency [18].
Deep unfolding methods [19, 20] circumvent this issue by offering a way to
include a priori information in a neural network. They have been applied to
many fields of imaging, such as denoising [21, 22], deblurring [23], MRI
reconstruction [24], and CT reconstruction from few-view data [25] Unfolding
consists of untying each iteration of an optimization algorithm for MBIR,
defining a set of learnable parameters, and training each iteration (or layer)
in an end-to-end manner. Some authors allow the unfolding network to learn the
optimization algorithm hyperparameters [23] as well as linear operators in the
regularization, such as convolution kernels in ISTA-net [26] and in [27].
Others use CNNs to replace proximity operators, as in PD-net [25] and ADMM-net
[24]. Deep unfolding networks automatically inherit from the feedback
mechanism of MBIR for data consistency.
In this paper, we exploit the framework of deep unfolding for addressing ROI
image reconstruction from few-view truncated measurements. We design an MBIR
cost function, propose an iterative algorithm suitable for minimizing our cost
function, and present an unfolding of this algorithm. Specifically:

  * •

we formulate a penalized cost function combining a data fidelity term arising
from the theory of M-estimators and a variant of total-variation
regularization,

  * •

we derive an iterative optimization algorithm combining an instance of the
dual block forward-backward algorithm (DBFB) [28], also known as dual
coordinate ascent technique, combined with an iterative reweighted scheme;

  * •

we perform the unfolding of the above algorithm allowing the supervised
learning of some parameters involved both in our cost function (e.g., adjoint
of the linear operators) and in our optimization algorithm (e.g., step sizes);

  * •

we show competitive performance against state-of-the-art for limiting sub-
sampling artifacts from objects inside and outside the reconstruction grid.

Section II of this work introduces our notation, and Section III presents a
mathematical formulation of ROI imaging. Section IV presents our choices for
the data fidelity and regularization terms and introduces a convergent
iterative algorithm to minimize the resulting cost function. We then explain
how this algorithm is unfolded into a deep learning architecture, called
U-RDBFB, in Section V. Section VI presents the training strategy for U-RDBFB
architecture. Finally, Section VII shows and discusses the results of our
experiments, which were conducted with both our unfolded networks and other
approaches proposed in the literature.

##  II Notation

Throughout the paper, we consider the Euclidean space $\mathbb{R}^{L}$
equipped with the standard scalar product $\langle\cdot,\cdot\rangle$, and the
norm $\|\cdot\|$. Moreover, $|||{\mathbf{L}}|||$ denotes the spectral norm of
squared matrix ${\mathbf{L}}\in\mathbb{R}^{L\times L}$. Let
$\mathcal{S}^{+}_{L}$ be the set of symmetric positive definite matrices in
$\mathbb{R}^{L\times L}$. For ${\mathbf{Q}}\in\mathcal{S}^{+}_{L}$,
$\|\cdot\|_{{\mathbf{Q}}}$ denotes the ${\mathbf{Q}}$-weighted norm, i.e., for
every ${\mathbf{x}}\in\mathbb{R}^{L}$,
$\|{\mathbf{x}}\|_{{\mathbf{Q}}}=\sqrt{{\left\langle{{{\mathbf{x}}}\mid{{\mathbf{Q}}{\mathbf{x}}}}\right\rangle}}$.
The diagonal matrix with diagonal entries equal to vector
${\mathbf{z}}=(z_{l})_{\ell=1}^{L}\in\mathbb{R}^{L}$ is denoted
$\operatorname{diag}((z_{\ell})_{\ell=1}^{L})$. The class of functions which
are proper, convex, lower-semicontinuous, defined on $\mathbb{R}^{L}$ and
taking values in $\mathbb{R}\cup\\{+\infty\\}$ is denoted by
$\Gamma_{0}(\mathbb{R}^{L})$. The proximity operator of
$g\in\Gamma_{0}(\mathbb{R}^{L})$ at ${\mathbf{x}}\in\mathbb{R}^{L}$ is
uniquely defined as [29] $\text{\rm
prox}_{g}({\mathbf{x}})=\underset{\begin{subarray}{c}{{\mathbf{z}}\in\mathbb{R}^{L}}\end{subarray}}{\text{\rm
argmin}}\;\;\left(g({\mathbf{z}})+\frac{1}{2}\|{\mathbf{x}}-{\mathbf{z}}\|^{2}\right)$.
The indicator function of a nonempty subset $C$ of $\mathbb{R}^{L}$ is the
function $\iota_{C}$ equal to 0 on $C$ and $+\infty$ outside of $C$. If $C$ is
closed and convex, the projection onto $C$ is denoted by $\text{\rm
proj}_{C}$.

##  III ROI imaging

Figure 1: ROI imaging in 2D. The exterior of the ROI contributes only to a
subset of the projections

We consider a 1D detector array of $B$ bins rotating around a patient. The
detector is too short to image the entire object; its size defines a circular
ROI we aim to reconstruct. Let $S$ be the number of projection angles. The
vector of truncated sub-sampled tomographic data is
${\mathbf{y}}\in\mathbb{R}^{T}$ with $T=B\times S$.
A reconstruction of the patient attenuation map in the ROI can be obtained by
considering a model of the form:

${\mathbf{H}}\overline{{\mathbf{x}}}_{\rm G}={\mathbf{y}}+{\mathbf{n}},$ (1)

where ${\mathbf{n}}\in\mathbb{R}^{M}$ accounts for some acquisition noise,
$\overline{{\mathbf{x}}}_{\rm G}\in\mathbb{R}^{L}$, and
${\mathbf{H}}\in\mathbb{R}^{L\times M}$ is the projector that models
projection over an intermediary grid $G$ that does not contain the support of
the entire object. Operator ${\mathbf{H}}$ contains a subset of the columns of
the projector on the entire space, or equivalently it corresponds to setting a
subset of the columns of the complete projector, matching the pixels outside
of the ROI, to zero.
When the grid corresponds to the ROI, (1) assumes that the image values are 0
outside of the ROI. As this assumption is not necessarily true, (1) is not a
good model for such a grid. A better idea is to extend the support of the
reconstruction beyond the ROI so that no assumption is made about the values
outside of the ROI. We now adopt this strategy.
We find an estimate of $\overline{{\mathbf{x}}}_{\rm G}$ by computing a
minimizer of a penalized cost function consisting of the sum of a data
fidelity term $f$ involving ${\mathbf{H}}$ and ${\mathbf{y}}$, and a
regularization term $r$, as

$\underset{\begin{subarray}{c}{{\mathbf{x}}\in\mathbb{R}^{L}}\end{subarray}}{\text{\rm
argmin}}\;\;f({\mathbf{x}})+r({\mathbf{x}}).$ (2)

The next section presents our choices for $f$ and $r$, as well as a suitable
algorithm to solve (2).

##  IV Iterative reconstruction

###  IV-A Cost function

####  IV-A1 Cauchy data fidelity

When objects with high gradients (metallic wires or needles) do not belong to
$G$, the error between the data and the reprojection of the estimate over $G$
contains outliers at the projections of those high gradients. The effect of
these outliers is amplified by the angular under-sampling; it leads to streaks
originating from these objects, i.e., from outside the grid. To avoid these
outliers, we propose to decrease the influence of the largest errors between
${\mathbf{y}}$ and the reprojection of
${\mathbf{H}}\overline{{\mathbf{x}}}_{\rm G}$ using M-estimators. M-estimators
have been used in ultrasound imaging [30] and in CT reconstruction [31], where
the authors investigated their effect on ring artifacts that appear due to
defective detector bins.
We focus on the Cauchy estimator $\phi$, one redescending M-estimator [32],
defined as

$(\forall\zeta\in\mathbb{R})\qquad\phi(\zeta)=\frac{\beta\kappa^{2}}{2}\ln\left(1+\left(\frac{\zeta}{\kappa}\right)^{2}\right),$
(3)

where $(\beta,\kappa)\in\left]0,+\infty\right[^{2}$.
We introduce the ramp filter [33] as ${\mathbf{F}}\in\mathcal{S}_{L}^{+}$,
which acts as a metric matrix and a derivative operator, and we apply the
Cauchy function to the ramp-filtered reprojection errors, i.e.,

$(\forall{\mathbf{x}}\in\mathbb{R}^{L})\quad
f({\mathbf{x}})=g({\mathbf{F}}({\mathbf{H}}{\mathbf{x}}-{\mathbf{y}}))$ (4)

with

$(\forall{\mathbf{z}}=(z_{t})_{1\leqslant t\leqslant
T}\in\mathbb{R}^{T})\qquad g({\mathbf{z}})=\sum^{T}_{t=1}\phi(z_{t}).$ (5)

Note that the resulting cost function is non-convex. With such $f$, the
reprojection error involving only the flat anatomical background is down-
weighted with respect to the error due to outliers; the effect of the Cauchy
term is more focused on outliers.

####  IV-A2 Semi-local total variation

Sub-sampling streaks due to a low angular density of the detector positions
can be effectively handled by total variation (TV) regularization [5, 3, 6].
Semi-local variants (STV) [34] extend TV in a neighborhood of pixels indexed
in $\Lambda_{J}=\\{-J,\ldots,J\\}\setminus\\{0\\}$:

$\displaystyle(\forall{\mathbf{x}}\in\mathbb{R}^{L})$ $\displaystyle r_{\rm
STV}({\mathbf{x}})$
$\displaystyle=\sum_{j=1}^{J}\sum_{\ell=1}^{L}\alpha_{j,\ell}\sqrt{({\mathbf{x}}-{\mathbf{V}}_{j}{\mathbf{x}})_{\ell}^{2}+({\mathbf{x}}-{\mathbf{V}}_{-j}{\mathbf{x}})_{\ell}^{2}}$
$\displaystyle=\sum_{j=1}^{J}r_{j}({\mathbf{L}}_{j}{\mathbf{x}}).$ (6)

Hereinabove $\ell\in\\{1,\ldots,L\\}$ is the spatial index and
${\mathbf{V}}_{j},{\mathbf{V}}_{-j}\in\mathbb{R}^{L\times L}$ are shift
operators as shown in Figure 2 for $j\in\\{1,\ldots,J\\}$ and $J=6$. Moreover,
for every $j\in\\{1,\cdots,J\\}$, we define
${\mathbf{L}}_{j}=\begin{bmatrix}{\mathbf{V}}_{j}^{\top}&{\mathbf{V}}_{-j}^{\top}\end{bmatrix}^{\top}\in\mathbb{R}^{2L\times
L}$ and, for every
${\mathbf{z}}=({\mathbf{z}}_{1},{\mathbf{z}}_{2})\in\mathbb{R}^{2L}$,
$h_{j}({\mathbf{z}})=\sum_{\ell=1}^{L}\alpha_{j,\ell}\sqrt{({\mathbf{z}}_{1})_{\ell}^{2}+({\mathbf{z}}_{2})_{\ell}^{2}}$.
Parameters $(\alpha_{j,\ell})_{1\leqslant j\leqslant J,1\leqslant\ell\leqslant
L}$ are nonnegative weights that can be chosen to vary spatially, so making
STV adaptive to the spatial contents [35]. We recover the standard TV penalty
for constant values of these parameters and $J=1$.

Figure 2: Shift operators $(V_{j})_{j\in\Lambda_{6}}$ applied to a given pixel
position $\ell$

We complement the above STV regularization function $r_{\rm STV}$ with a
nonnegativity constraint on the pixel values, and a quadratic regularization
$\frac{1}{2}\|{\mathbf{x}}\|^{2}_{{\mathbf{M}}}={\mathbf{x}}^{\top}{\mathbf{M}}{\mathbf{x}}$.
Here, matrix
${\mathbf{M}}=\operatorname{diag}((m_{\ell})^{L}_{\ell=1})\in\mathcal{S}_{L}^{+}$
is such that, for every $\ell\in\\{1,\ldots,L\\}$, $m_{\ell}=1$ if the
$\ell$-th entry $x_{\ell}$ of vector ${\mathbf{x}}$ belongs to the ROI, and
$m_{\ell}=\xi>1$ otherwise. In such a way,
${\mathbf{M}}\in\mathcal{S}_{L}^{+}$ acts as a mask allowing us to smooth the
reconstruction outside of the ROI.
Altogether, our regularization function in (2) reads

$(\forall{\mathbf{x}}\in\mathbb{R}^{L})\quad
r({\mathbf{x}})=\sum_{j=1}^{J}r_{j}({\mathbf{L}}_{j}{\mathbf{x}})+\frac{1}{2}\|{\mathbf{x}}\|^{2}_{{\mathbf{M}}}+\iota_{[0,+\infty[^{L}}({\mathbf{x}}).$
(7)

###  IV-B Minimization of the cost function

####  IV-B1 Reweighting for non-convex criterion

Given our choices for $r$ and $f$, the optimization problem (2) becomes

$\underset{\begin{subarray}{c}{{\mathbf{x}}\in\mathbb{R}^{L}}\end{subarray}}{\text{\rm
argmin}}\;\;g({\mathbf{F}}({\mathbf{H}}{\mathbf{x}}-{\mathbf{y}}))+\sum_{j=1}^{J}r_{j}({\mathbf{L}}_{j}{\mathbf{x}})+\frac{1}{2}\|{\mathbf{x}}\|^{2}_{{\mathbf{M}}}+\iota_{[0,+\infty[^{L}}({\mathbf{x}}).$
(8)

Because of the non-convexity of the data fidelity term, we adopt an iterative
reweighting strategy where Problem (8) is replaced by a sequence of surrogate
convex problems built following a majoration principle. As shown in Appendix
A, the associated Majorization-Minimization strategy is given by Algorithm 1.

Algorithm 1 Iterative reweighting strategy for Problem (8)

${\mathbf{x}}_{0}\in\mathbb{R}^{L}$

for $k=0$ to $K-1$ do

Define majorant function $\quad Q({\mathbf{x}};\overline{{\mathbf{x}}})$

${\mathbf{x}}_{k+1}=\underset{\begin{subarray}{c}{{\mathbf{x}}\in\mathbb{R}^{L}}\end{subarray}}{\text{\rm
argmin}}\;\;Q({\mathbf{x}},{\mathbf{x}}_{k})$

end for

According to (56), for every $\overline{{\mathbf{x}}}\in\mathbb{R}^{L}$, the
surrogate function $Q(\cdot,\overline{{\mathbf{x}}})$ takes the form

$Q({\mathbf{x}},\overline{{\mathbf{x}}})=\iota_{[0,+\infty[^{L}}({\mathbf{x}})+h({\mathbf{B}}{\mathbf{x}};{\mathbf{B}}\overline{{\mathbf{x}}})+\frac{1}{2}\|{\mathbf{x}}\|^{2}_{{\mathbf{M}}},$
(9)

where $\forall({\mathbf{x}},\overline{{\mathbf{x}}})\in(\mathbb{R}^{L})^{2}$,

$\displaystyle
h({\mathbf{B}}{\mathbf{x}};{\mathbf{B}}\overline{{\mathbf{x}}})$
$\displaystyle=h_{0}({\mathbf{B}}_{0}{\mathbf{x}};{\mathbf{B}}_{0}\overline{{\mathbf{x}}})+h_{1}({\mathbf{B}}_{1}{\mathbf{x}})$
$\displaystyle h_{1}({\mathbf{B}}_{1}{\mathbf{x}})$
$\displaystyle=\sum_{j=1}^{J}r_{j}({\mathbf{L}}_{j}{\mathbf{x}}),$ (10)

with

$\displaystyle{\mathbf{B}}$
$\displaystyle=\begin{bmatrix}{\mathbf{B}}_{0}^{\top}&{\mathbf{B}}_{1}^{\top}\end{bmatrix}^{\top},$
(11) $\displaystyle{\mathbf{B}}_{0}$
$\displaystyle={\mathbf{F}}{\mathbf{H}},\quad{\mathbf{B}}_{1}=\begin{bmatrix}{\mathbf{L}}_{1}^{\top}&\cdots&{\mathbf{L}}_{J}^{\top}\end{bmatrix}^{\top}$
(12)

and, for every $\overline{{\mathbf{x}}}\in\mathbb{R}^{L}$,
$h_{0}(\cdot;{\mathbf{B}}_{0}\overline{{\mathbf{x}}})=\tilde{g}(\cdot-{\mathbf{F}}{\mathbf{y}};{\mathbf{B}}_{0}\overline{{\mathbf{x}}}-{\mathbf{F}}{\mathbf{y}})$.

####  IV-B2 Dual block forward-backward algorithm

Since ${\mathbf{M}}\in\mathcal{S}^{+}_{L}$, function
$Q(\cdot,\overline{{\mathbf{x}}})$ is strongly convex. Its minimization at
each reweighted step is well-defined, with a unique solution that can be
conveniently obtained using the dual forward-backward algorithm [36, 22]. An
accelerated version of the dual forward-backward algorithm is its block
coordinate version (DBFB) [37], which allows for splitting the terms of
function (IV-B1) and accessing their proximity operators separately. We
exploit the splitting of $h$ into two blocks $h_{0}$ and $h_{1}$ to yield the
form of the DBFB algorithm provided in (2). The estimate $\boldsymbol{x}_{n}$
generated by the algorithm at iteration $n$ has been denoted in bold italic to
distinguish it from the reweighting point $\overline{{\mathbf{x}}}$.

Algorithm 2 DBFB algorithm to minimize (9) with
$\overline{{\mathbf{x}}}\in\mathbb{R}^{L}$
$\displaystyle{\mathbf{z}}_{0}^{0}\in\mathbb{R}^{T},(\forall
j\in\\{1,\ldots,J\\})\;{\mathbf{z}}^{j}_{0}\in\mathbb{R}^{2L},{\mathbf{s}}^{1}_{0}=({{\mathbf{z}}^{1}_{0}},\ldots,{{\mathbf{z}}^{J}_{0}})$
$\displaystyle(\sigma,\tau_{1},\ldots,\tau_{J})\in]0,+\infty[^{J+1}$
$\displaystyle\Sigma=\operatorname{diag}(\tau_{1},\ldots,\tau_{J})$
$\displaystyle{\mathbf{w}}_{0}=-{\mathbf{M}}^{-1}({\mathbf{B}}_{0}^{\top}{\mathbf{z}}_{0}^{0}+{\mathbf{B}}_{1}^{\top}{\mathbf{s}}_{0}^{1})$
(13) $\displaystyle\text{For }n=0,1\ldots$
$\displaystyle\quad\left\lfloor\begin{array}[]{l}({{\mathbf{z}}^{1}_{n}},\ldots,{{\mathbf{z}}^{J}_{n}})\equiv{\mathbf{s}}^{1}_{n}\\\
\boldsymbol{x}_{n}=\text{\rm proj}_{[0,+\infty[^{L}}({\mathbf{w}}_{n})\\\
\mbox{Select $\varepsilon_{n}\in\\{0,1\\}$ and $\gamma_{n}\in]0,+\infty[$}\\\
\text{If }\varepsilon_{n}=0\quad\text{{(D)}}\\\
\quad\left\lfloor\begin{array}[]{l}\tilde{{\mathbf{z}}}_{n}^{0}={\mathbf{z}}_{n}^{0}+\gamma_{n}\sigma^{-1}{\mathbf{B}}_{0}\boldsymbol{x}_{n}\\\
{\mathbf{z}}_{n+1}^{0}=\tilde{{\mathbf{z}}}_{n}^{0}-\gamma_{n}\sigma^{-1}\text{\rm
prox}_{\gamma_{n}^{-1}\sigma
h_{0}(\cdot;{\mathbf{B}}_{0}\overline{{\mathbf{x}}})}(\gamma_{n}^{-1}\sigma\tilde{{\mathbf{z}}}_{n}^{0})\\\
{\mathbf{w}}_{n+1}={\mathbf{w}}_{n}-{\mathbf{M}}^{-1}{\mathbf{B}}_{0}^{\top}({\mathbf{z}}_{n+1}^{0}-{\mathbf{z}}_{n}^{0})\\\
{\mathbf{s}}_{n+1}^{1}={\mathbf{s}}_{n}^{1}\end{array}\right.\\\ \text{If
}\varepsilon_{n}=1\quad\text{{(R)}}\\\
\quad\left\lfloor\begin{array}[]{l}\tilde{{\mathbf{s}}}_{n}^{j}={\mathbf{s}}_{n}^{j}+\gamma_{n}\Sigma_{j}^{-1}{\mathbf{B}}_{1}\boldsymbol{x}_{n}\\\
{\mathbf{s}}_{n+1}^{j}=\tilde{{\mathbf{s}}}_{n}^{j}-\gamma_{n}\Sigma_{j}^{-1}\text{\rm
prox}_{\gamma_{n}^{-1}\Sigma_{j}h_{1}}(\gamma_{n}^{-1}\Sigma_{j}\tilde{{\mathbf{s}}}_{n}^{j})\\\
{\mathbf{w}}_{n+1}={\mathbf{w}}_{n}-{\mathbf{M}}^{-1}{\mathbf{B}}_{1}^{\top}({\mathbf{s}}^{1}_{n+1}-{\mathbf{s}}^{1}_{n})\\\
{\mathbf{z}}_{n+1}^{0}={\mathbf{z}}_{n}^{0}\end{array}\right.\\\
\end{array}\right.$ (22)

At each iteration of the DBFB algorithm, either the data fidelity term or the
regularization term is processed, through the data (D) and regularization (R)
steps, respectively. The implementation of each step is discussed in Appendix
B.

##  V Unfolded reconstruction

We now present our unfolding version for the reweighted DBFB, designated as
U-RDBFB (Unfolded Reweighted DBFB).

###  V-A Layering and relaxation of the DBFB algorithm

The unfolding principle recasts every algorithm iteration as one neural
network layer to allow selected algorithm parameters to be learned in a
supervised way. Each layer is thus made of costly and ill-conditioned linear
operators (${\mathbf{B}}_{0}$, ${\mathbf{B}}_{1}$ and their adjoints). Several
works in iterative reconstruction [38, 39, 40, 41] have shown the potential to
replace adjoint operators with surrogates to accelerate convergence [42, 43,
44]. Therefore, in order to minimize the number of U-RDBFB layers, we propose
relaxed steps (D) and (R) with surrogates, and thus mismatched adjoints.
We first focus on step (D). To gain more flexibility, we replace the transpose
of ${\mathbf{F}}$ by an operator $\tilde{{\mathbf{F}}}\in\mathbb{R}^{L\times
L}$. By setting $\nu_{n,0}=\gamma_{n}\sigma^{-1}$,
${\mathbf{u}}_{n}=\tilde{{\mathbf{z}}}_{n}^{0}-\nu_{n,0}{\mathbf{F}}{\mathbf{y}}$,
and by using the relation

$\text{\rm
prox}_{\nu_{n,0}h_{0}(\cdot;{\mathbf{B}}_{0}\overline{{\mathbf{x}}})}=\text{\rm
prox}_{\nu_{n,0}\tilde{g}(\cdot;{\mathbf{H}}\overline{{\mathbf{x}}}-{\mathbf{y}})}(\cdot-{\mathbf{F}}{\mathbf{y}})+{\mathbf{F}}{\mathbf{y}},$

we modify step (D) to define layer $\mathcal{L}_{\rm D}$, mimicking the $n$-th
iteration of DBFB for $\varepsilon_{n}=0$ and a given
$\overline{{\mathbf{x}}}$:

Data layer ($\mathcal{L}_{\rm D}$):
$\displaystyle\quad\left\lfloor\begin{array}[]{l}\boldsymbol{x}_{n}=\text{\rm
proj}_{[0,+\infty[^{L}}({\mathbf{w}}_{n})\\\
{\mathbf{u}}_{n}={\mathbf{z}}_{n}^{0}+\nu_{n,0}{\mathbf{F}}({\mathbf{H}}\boldsymbol{x}_{n}-{\mathbf{y}})\\\
{\mathbf{z}}_{n+1}^{0}={\mathbf{u}}_{n}-\nu_{n,0}^{-1}\,\text{\rm
prox}_{\nu_{n,0}\tilde{g}(\cdot;{\mathbf{H}}\overline{{\mathbf{x}}}-{\mathbf{y}})}(\nu_{n,0}{\mathbf{u}}_{n})\\\
{\mathbf{w}}_{n+1}={\mathbf{w}}_{n}-{\mathbf{M}}^{-1}{\mathbf{H}}^{\top}\tilde{{\mathbf{F}}}({\mathbf{z}}_{n+1}^{0}-{\mathbf{z}}_{n}^{0})\\\
{\mathbf{z}}_{n+1}^{j}={\mathbf{z}}_{n}^{j}\quad(\forall
j\in\\{1,\ldots,J\\}).\end{array}\right.$ (28)

Similarly, for step (R), we unfold by replacing the TV-based operators
$({\mathbf{L}}^{\top}_{j})_{j=1}^{J}$ by surrogate learnable operators
$(\tilde{{\mathbf{L}}}_{j})_{j=1}^{J}$. Using (59) and setting, for every
$j\in\\{1,\ldots,K\\}$, $\nu_{n,j}=\gamma_{n}\tau_{j}^{-1}$ yield the
following regularization layer ($\mathcal{L}_{\rm R}$), associated with the
$n$-th DBFB iteration when $\varepsilon_{n}=1$:

Regularization layer ($\mathcal{L}_{\rm R}$) :
$\displaystyle\quad\left\lfloor\begin{array}[]{l}\boldsymbol{x}_{n}=\text{\rm
proj}_{[0,+\infty[^{L}}({\mathbf{w}}_{n})\\\ \mbox{For
$j\in\\{1,\ldots,J\\}$}\\\
\quad\left\lfloor\begin{array}[]{l}(\forall\ell\in\\{1,\ldots,L\\})\\\
\quad\displaystyle({\mathbf{z}}^{j}_{n+1})_{\ell}=\frac{\left({\mathbf{z}}^{j}_{n}+\nu_{n,j}{\mathbf{L}}_{j}\boldsymbol{x}_{n}\right)_{\ell}}{\max\big{\\{}1,\|\left({\mathbf{z}}^{j}_{n}+\nu_{n,j}{\mathbf{L}}_{j}\boldsymbol{x}_{n}\right)_{\ell}\|_{2}/\alpha_{j,\ell}\big{\\}}}\end{array}\right.\\\
{\mathbf{w}}_{n+1}={\mathbf{w}}_{n}-{\mathbf{M}}^{-1}\sum_{j=1}^{J}\tilde{{\mathbf{L}}}_{j}({\mathbf{z}}^{j}_{n+1}-{\mathbf{z}}^{j}_{n})\\\
{\mathbf{z}}_{n+1}^{0}={\mathbf{z}}_{n}^{0}.\end{array}\right.$ (36)

We refer to the above regularization layer as $\mathcal{L}_{R}^{*}$ when using
it without learning surrogates adjoints, i.e. when, for every
$j\in\\{1,\ldots,J\\}$, $\tilde{{\mathbf{L}}}_{j}={\mathbf{L}}^{\top}_{j}$.
Note that $\mathcal{L}_{R}^{*}$ does not involve $\overline{{\mathbf{x}}}$.

###  V-B Total architecture

The $K$ groups of layers of U-RDBFB architecture correspond to $K$ reweighted
steps, each consisting of $N_{k}$ unfolded layers of DBFB. This yields an
architecture with $\sum_{k=0}^{K-1}N_{k}$ layers in total. We investigate
three variants of this architecture:

  * •

Architecture 1 ($\mathcal{A}_{1}$): the step sizes and the regularization
parameters are learned and $\tilde{{\mathbf{F}}}={\mathbf{F}}^{\top}$.
The $k$-th group of layers with $k\in\\{0,\ldots,K-1\\}$ implements the
following update:

$({\mathbf{x}}_{k+1},{\mathbf{z}}_{0,k+1})=\mathcal{L}_{1}^{k}({\mathbf{x}}_{k},\Theta_{1,k},{\mathbf{z}}_{0,k}).$
(37)

Hereinabove, $\mathcal{L}_{1}^{k}$ is a composition of $N_{k}$ layers
$\mathcal{L}_{\rm D}$ or $\mathcal{L}_{\rm R}^{*}$. ${\mathbf{x}}_{k}$
(corresponding to $\boldsymbol{x}_{0}$ in DBFB) is the initial value of the
primal variable and ${\mathbf{x}}_{k+1}$ the generated one (corresponding to
$\boldsymbol{x}_{N_{k}-1}$ in DBFB).
${\mathbf{z}}_{0,k}\in\mathbb{R}^{T}\times(\mathbb{R}^{2L})^{J}$ is the
initial value of the dual variables (corresponding to
$({\mathbf{z}}_{0}^{j})_{0\leqslant j\leqslant J}$ in DBFB) and
${\mathbf{z}}_{0,k+1}$ is the final value of these variables (corresponding to
$({\mathbf{z}}_{N_{k-1}}^{j})_{0\leqslant j\leqslant J}$ in DBFB). The
variable ${\mathbf{w}}_{0}$ required in the unfolded version of DBFB is set as
in (13). $\Theta_{1,k}$ is the vector of trainable parameters. As we allow all
these parameters to be layer-dependent, we thus index them by $k$ and $n$. If
the $n$-th layer of $\mathcal{L}_{1}^{k}$ is a $\mathcal{L}_{\rm D}$ layer,
its parameters are those of the Cauchy function $(\beta_{k,n},\kappa_{k,n})$,
the one of the quadratic regularization $\xi_{k,n}$, and a single step size
$\nu_{k,n,0}$. If this is a $\mathcal{L}^{*}_{\rm R}$ layer, regularization
parameters $(\alpha_{k,n,j,\ell})_{1\leqslant j\leqslant
J,1\leqslant\ell\leqslant L}$, $\kappa_{k,n}$, $\xi_{k,n}$, and step sizes
$(\nu_{k,n,j})_{j=1}^{J}$ are learned.

  * •

Architecture 2 ($\mathcal{A}_{2}$): similar to $\mathcal{A}_{1}$ with
$\tilde{{\mathbf{F}}}=\operatorname{Id}\,$.

  * •

Architecture 3 ($\mathcal{A}_{3}$): in addition to algorithm parameters (step
sizes, Cauchy function and STV weights), the surrogates
$(\tilde{{\mathbf{L}}}_{j})_{j=1}^{J}$ to the adjoints of operators
$({\mathbf{L}}_{j}^{\top})_{j=1}^{J}$ are also learned so that the $k$-th
group of layers of Architecture 3 implements the update rule:

$({\mathbf{x}}_{k+1},{\mathbf{z}}_{0,k+1})=\mathcal{L}_{2}^{k}({\mathbf{x}}_{k},\Theta_{2,k},{\mathbf{z}}_{0,k}),$
(38)

with $\tilde{{\mathbf{F}}}=\operatorname{Id}\,$. $\mathcal{L}_{2}^{k}$ is then
a composition of $N_{k}$ layers $\mathcal{L}_{\rm D}$ or $\mathcal{L}_{\rm
R}$. Vector $\Theta_{2,k}$ contains the same parameters as $\Theta_{1,k}$ and
the elements of matrices $(\tilde{{\mathbf{L}}}_{j})_{j=1}^{J}$ for each
$\mathcal{L}_{\rm R}$ layer.

Figure 3 illustrates one instance of Architecture 3, where
$\mathcal{A}_{3}=\mathcal{L}^{1}_{2}\circ\mathcal{L}^{0}_{2}$ with
$\mathcal{L}^{0}_{2}=(\mathcal{L}_{\rm R})^{2}\circ\mathcal{L}_{\rm D}$ and
$\mathcal{L}^{1}_{2}=\mathcal{L}_{\rm D}$ (i.e., $K=2$, $N_{0}=3$ and
$N_{1}=1$).

###  V-C Parameter learning

To infer the parameters gathered in vector $\Theta_{i,k}$, in the $k$-th group
of layers $\mathcal{L}_{i}^{k}$ of architecture $i\in\\{1,2\\}$, we introduce
hidden layers. In the following,
$(\mathcal{L}_{k,n}^{(\theta)})_{n=0}^{N_{k}-1}$ designates the hidden layers
used to estimate a given parameter $\theta$. We now discuss our choices for
$(\mathcal{L}_{k,n}^{(\theta)})_{n=0}^{N_{k}-1}$.

  * •

step size for $\mathcal{L}_{\rm D}$:
$\nu_{k,n,0}=\mathcal{L}_{k,n}^{(\nu)}=\operatorname{softplus}(a_{k,n})$ where
$a_{k,n}$ is a learnable real-valued parameter.

  * •

Parameters of the Cauchy function for $\mathcal{L}_{\rm D}$:

    * $*$

$\kappa_{k,n}=\mathcal{L}_{k,n}^{(\kappa)}=W_{\kappa}\operatorname{softplus}(c_{k,n})$
where $c_{k,n}$ is inferred from a fully connected layer, whose weights are
shared across the U-RDBFB network, applied on a histogram of the absolute
value of the filtered reprojection error i.e.
${\mathbf{F}}({\mathbf{H}}^{\top}{\mathbf{x}}_{n}-{\mathbf{y}})$. We
implemented the learnable histogram layer proposed in [45], which is piecewise
differentiable. More precisely, we built a cumulated histogram using 100 bins
from 0 to the maximum value of the filtered reprojection error.

    * $*$

$\beta_{k,n}=\mathcal{L}_{k,n}^{(\beta)}=W_{\beta}\operatorname{softplus}(d_{k,n})$.

  * •

Diagonal elements of ${\mathbf{M}}$ involved in (7) corresponding to the
locations of pixels outside of the ROI for both $\mathcal{L}_{\rm D}$ and
$\mathcal{L}_{\rm R}$/$\mathcal{L}_{\rm R}^{*}$:
$\xi_{k,n}=\mathcal{L}_{k,n}^{(\xi)}=\operatorname{softplus}(e_{k,n})$ where
$e_{k,n}$ is learned.

  * •

step size for $\mathcal{L}_{\rm R}$/$\mathcal{L}_{\rm R}^{*}$: For every
$j\in\\{1,\ldots,J\\}$,
$\nu_{k,n,j}=\mathcal{L}_{k,n}^{(\nu_{j})}=W_{\nu}\operatorname{softplus}(b_{k,n,j})$.

  * •

Parameters of the STV regularization for $\mathcal{L}_{\rm
R}$/$\mathcal{L}_{\rm R}^{*}$: For every $j\in\\{1,\ldots,J\\}$,

$\displaystyle\alpha_{k,n,j}$
$\displaystyle=(\alpha_{k,n,j,l})_{\ell=1}^{L}=\mathcal{L}_{k,n}^{(\alpha_{j})}$
$\displaystyle=W_{\alpha}\operatorname{softplus}(A_{k,n}\circ\operatorname{ReLU}\circ
B_{k,n}({\mathbf{L}}_{j}{\mathbf{x}}_{k})),$

where $A_{k,n}$ and $B_{k,n}$ are two convolutional operators with size
$3\times 3$ kernels and $8$ channels.

Finally, when using layer $\mathcal{L}_{R}$, we also learn operators
$(\tilde{\mathbf{L}}^{k,n}_{j})_{j=1}^{J}$ which have the same support as
$(\mathbf{L}^{k,n}_{j})_{j=1}^{J}$.

Figure 3: Example of Architecture 3 of U-RDBFB in the case where
$\mathcal{A}_{3}=\mathcal{L}^{1}_{2}\circ\mathcal{L}^{0}_{2}$ where
$\mathcal{L}^{0}_{2}=(\mathcal{L}_{R})^{2}\circ\mathcal{L}_{D}$ and
$\mathcal{L}^{1}_{2}=\mathcal{L}_{D}$ (i.e., $K=2$, $N_{0}=3$, $N_{1}=1$). Red
blocks represent the hidden structures to infer all the parameters
$\theta\in\Theta_{2}$. When $k=1$, the dual variables of DBFB are initialized
with the values of the dual variables at the end of the previous $N$
iterations of DBFB for $k=0$.

In $\mathcal{L}_{D}$, initial values of $a_{k,n}$, $d_{k,n}$ are set to $1$.
In $\mathcal{L}_{R}$ and $\mathcal{L}_{R}^{*}$, initial values for
$b_{k,n,j}$, $e_{k,n}$ are $1$. Normalization scalars $W_{\kappa}$ and
$W_{\nu}$ are set to $10^{-5}$ and $40$ for $\mathcal{L}^{*}_{R}$ and to
$10^{-5}$ and 10 for $\mathcal{L}_{R}$, respectively. When
$\tilde{{\mathbf{F}}}=\operatorname{Id}\,$, $W_{\beta}$ and $W_{\alpha}$ are
set to $10$ and $0.05$ and, when $\tilde{{\mathbf{F}}}={\mathbf{F}}^{\top}$,
$W_{\beta}$ and $W_{\alpha}$ are set to $2\times 10^{5}$ and $5$.
A schematic view of the data layer can be found in Figure 4, whereas Figure 5
shows an example of the regularization layer according to these settings.

Figure 4: Schematic view of layer $\mathcal{L}_{D}$, mimicking (V-A), for
$\overline{{\mathbf{x}}}={\mathbf{x}}_{k}$, i.e. $k$-th iterate of the
reweighted algorithm. The layer takes as inputs ${\mathbf{w}}_{n}$,
${\mathbf{z}}_{n}^{0}$ from the previous layers. The projections
${\mathbf{y}}$ are also used as input. Figure 5: Schematic view of layer
$\mathcal{L}_{R}$ mimicking (36). The layer takes as inputs
${\mathbf{w}}_{n}$, $({\mathbf{z}}_{n}^{j})_{1}^{J}$ from the previous layer.
The update of parameters for $j\in\\{2,\ldots,5\\}$ is hidden in the orange
block for the sake of readability. When $\mathcal{L}_{R}$ reduces to
$\mathcal{L}_{R}^{*}$, the blue blocks disappear. Only parameters $\alpha_{j}$
depend on the input.

##  VI Experiments

###  VI-A Datasets

Two datasets of CT images were used in our experiments.

####  VI-A1 Abdomen dataset

Our first dataset consists of 2D images obtained from 70 CT volumes of size
$512\times 512\times 512$ from the lower lungs to the lower abdomen of 70
patients, which were extracted from the public dataset CT Lymph Nodes from
https://www.cancerimagingarchive.net/. These volumes correspond to fully
sampled CT reconstructions. They were made isotropic by interpolating the
axial slices. A total of 50 out of 512 slices were kept per volume. We
randomly added intense metallic wires between 3000 and 5000 Hounsfield units
(HU) of varying sizes on the axial slices. We shifted the HU values of the
images by 1000 so that air is 0 HU and water is 1000 HU, using
$\boldsymbol{a}_{\rm tissue}\mapsto(\boldsymbol{a}_{\rm tissue}-\mu_{\rm
water})\times(1000/\mu_{\rm water})$, where $\mu_{\rm water}$ is the value of
the attenuation coefficient of water equal to 0.017 mm-1 and
$\boldsymbol{a}_{\rm tissue}\in\mathbb{R}^{Q}$ is the initial vector of
attenuation values. Finally the $512\times 512$ slices
$\overline{{\mathbf{x}}}_{\rm P}$ were normalized between $[0,1]$ (a value of
$1$ corresponding to an object of HU intensity equal to $5000$).
Projections were simulated for each slice of each volume in a 2D parallel
geometry with a short detector of 300 bins (bin size equal to pixel size, i.e.
1 mm) and an angular density of 110 projections over $180^{\degree}$. Noisy
projections ${\mathbf{y}}=(y_{t})_{t=1}^{T}$ are computed as

$(\forall t\in\\{1,\ldots,T\\})\quad
y_{t}=\mu\log\left(\frac{I_{0}}{\mathcal{P}(I_{0}\exp{(-\mu({\mathbf{H}}\overline{{\mathbf{x}}}_{\rm
P})_{t})})}\right),$ (39)

where we set $\mu=\mu_{\rm water}/1000$, $I_{0}=10^{4}$, and, for some
$\delta>0$, $\mathcal{P}(\delta)$ denotes a realization of a Poisson law with
mean $\delta$. In this context, the ROI was a centered disk of diameter $300$.
The resulting pairs of axial slice/projections $(\overline{{\mathbf{x}}}_{\rm
P},{\mathbf{y}})$ were split into a training of 2500 pairs from a pool of 50
patients, a testing set of 500 pairs from 10 other patients, and a validation
of 500 pairs from the 10 remaining patients.

####  VI-A2 Brain dataset

We used a second dataset containing 2D images extracted from 10 CT high-dose
brain reconstructions. These volumes are from the public repository 2016 Low
Dose CT Grand Challenge from https://www.cancerimagingarchive.net/. We
performed the same post-processing as above (with a pixel size of 0.5 mm) for
generating a training set of 2500 pairs of axial slices/projections
$(\overline{{\mathbf{x}}}_{\rm P},{\mathbf{y}})$ from 30 patients, validation
and testing sets of 500 pairs from 10 different patients each.

###  VI-B Training strategy

Let $i=\\{1,\ldots,I\\}$ be the index covering all $I$ instances of the
training set. The reconstruction grid ($G$) is a disk of diameter $400$. Let
${\mathbf{x}}^{*}_{\rm G,i}\in\mathbb{R}^{L}$ be the output of U-RDBFB for a
given projection input ${\mathbf{y}}_{i}$. Our networks are designed to
minimize $\sum_{i=1}^{I}\ell({\mathbf{C}}_{\rm G}{\mathbf{x}}^{*}_{{\rm
G},i},{\mathbf{C}}_{\rm P}\overline{{\mathbf{x}}}_{{\rm P},i})$, where
${\mathbf{C}}_{\rm G}$ is a cropping operator which extracts the ROI from the
grid $G$, ${\mathbf{C}}_{\rm P}$ is a cropping operator which extracts the ROI
from the entire $512\times 512$ grid, and $\ell(\cdot)$ is the loss retained
for training the network.
We use the following similar compositions of layers for defining
$\mathcal{A}_{1}$, $\mathcal{A}_{2}$ and $\mathcal{A}_{3}$:

$\displaystyle\mathcal{A}_{1}=\mathcal{A}_{2}=\mathcal{L}^{5}_{1}\circ\mathcal{L}^{4}_{1}\circ\mathcal{L}^{3}_{1}\circ\mathcal{L}^{2}_{1}\circ\mathcal{L}^{1}_{1}\circ\mathcal{L}^{0}_{1},$
(40)

where

$\displaystyle\mathcal{L}^{0}_{1}$
$\displaystyle=(\mathcal{L}_{R}^{*})^{5}\circ\mathcal{L}_{D}$
$\displaystyle\mathcal{L}^{1}_{1}$
$\displaystyle=\mathcal{L}^{2}_{1}=(\mathcal{L}_{R}^{*})^{4}\circ\mathcal{L}_{D}$
$\displaystyle\mathcal{L}^{3}_{1}$
$\displaystyle=\mathcal{L}^{4}_{1}=\mathcal{L}^{5}_{1}=(\mathcal{L}_{R}^{*})^{3}\circ\mathcal{L}_{D}$
(41)

and

$\displaystyle\mathcal{A}_{3}=\mathcal{L}^{5}_{2}\circ\mathcal{L}^{4}_{2}\circ\mathcal{L}^{3}_{2}\circ\mathcal{L}^{2}_{2}\circ\mathcal{L}^{1}_{2}\circ\mathcal{L}^{0}_{2},$
(42)

where

$\displaystyle\mathcal{L}^{0}_{2}$
$\displaystyle=(\mathcal{L}_{R})^{5}\circ\mathcal{L}_{D}$
$\displaystyle\mathcal{L}^{1}_{2}$
$\displaystyle=\mathcal{L}^{2}_{2}=(\mathcal{L}_{R})^{4}\circ\mathcal{L}_{D}$
$\displaystyle\mathcal{L}^{3}_{2}$
$\displaystyle=\mathcal{L}^{4}_{2}=\mathcal{L}^{5}_{2}=(\mathcal{L}_{R})^{3}\circ\mathcal{L}_{D}.$
(43)

These architectures correspond to setting $K=6$ and using a decreasing value
of $N_{k}$ for each $k\in\\{0,\cdots,K-1\\}$ ($N_{0}=6$, $N_{1}=N_{2}=5$,
$N_{3}=N_{4}=N_{5}=4$), resulting in a total of 28 layers per network.
For the first reweighted iteration ($k=0$), the dual estimates
$({\mathbf{z}}^{j}_{0})_{j=1}^{J}$ are initialized to zero while the dual
estimate ${\mathbf{z}}^{0}_{0}$ is set to the negative filtered projections,
$-{\mathbf{F}}{\mathbf{y}}$, leading to
${\mathbf{x}}_{0}={\mathbf{H}}^{\top}{\mathbf{F}}{\mathbf{y}}$. We implemented
all networks in Pytorch, using a Tesla V100 32 Gb GPU. We employed the toolbox
TorchRadon [46] to include Pytorch-compatible parallel-beam tomographic
operators in all architectures. The training procedure of $\mathcal{A}_{1}$,
$\mathcal{A}_{2}$ and $\mathcal{A}_{3}$ takes about one day and a half.
Standard auto-differentiation tools can compute all necessary derivatives for
backpropagation.

####  VI-B1 From layer-wise to greedy training

Our proposed initialization is a rough estimation of the parameters used for
$K$ reweighted iterations and $N_{k}$ inner iterations of DBFB. We circumvent
the initialization issue by using a layer-wise sequential training strategy,
as sometimes advocated for when initializing the weights of recurrent neural
networks [47, 48]. The learning in each layer $(k,n)$ ($n\geqslant 1$) starts
by considering only its previous layers from $(0,0)$ to $(k,n-1)$. As more and
more layers are trained simultaneously, all 28 layers are eventually trained
from end to end. We used ten epochs for training each layer $\mathcal{L}_{\rm
R}$ and twenty epochs for training each layer $\mathcal{L}_{\rm D}$; the only
exception was the last layer, for which we used 50 epochs. The learning rate
is decreased with a step decay by a factor of 0.995 from $10^{-2}$ every 4
epochs. The batch size used for each epoch varied from 20 to 8 as the number
of trained layers increased.

####  VI-B2 Choice for $\ell$

For all instances of the training set, $\ell({\mathbf{C}}_{\rm
G}{\mathbf{x}}^{*}_{\rm G,i},{\mathbf{C}}_{\rm P}\overline{{\mathbf{x}}}_{\rm
P,i})=\frac{1}{I}\|{\mathbf{C}}_{\rm G}{\mathbf{x}}^{*}_{\rm
G,i}-{\mathbf{C}}_{\rm P}\overline{{\mathbf{x}}}_{\rm P,i})\|^{2}$,
corresponding to the MSE loss for $\mathcal{A}_{1}$ and $\mathcal{A}_{2}$. For
$\mathcal{A}_{3}$, we used a regularized MSE loss, which contains an
additional term
$-\lambda\sum_{k=0}^{K-1}\sum_{n=0}^{N_{k}-1}\sum_{j=1}^{J}\operatorname{trace}\left(\tilde{{\mathbf{L}}}^{k,n}_{j}{\mathbf{L}}^{k,n}_{j}\right)$
weighted by parameter $\lambda>0$. To set $\lambda$, $\mathcal{A}_{2}$ was
overfit to an example of our dataset. A grid search on $\lambda$ was performed
and we selected the value which minimized the MSE on this example.

###  VI-C Comparative strategy

The quantitative metric used to assess the reconstruction quality of
${\mathbf{C}}_{\rm G}{\mathbf{x}}^{*}_{\rm G,i}$ is the PSNR. We also evaluate
the reconstruction performance using the structural similarity index (SSIM),
the PieApp value [49], and the Mean Absolute Error (MAE) of the difference
between ${\mathbf{C}}_{\rm G}{\mathbf{x}}^{*}_{\rm G,i}$ and
${\mathbf{C}}_{\rm P}\overline{{\mathbf{x}}}_{\rm P,i}$.
We compare our approach with five reconstruction methods, that we describe
hereinafter.

####  VI-C1 Hierarchical reconstruction

The first competing method is the iterative hierarchical reconstruction
algorithm of [6]. Inspired by homotopy strategies for solving
$\ell_{1}$-regularized least-squares problems [50, 51], this method amounts to
solving, for $k\in\\{1,\ldots,K\\}$, the following optimization problem
parameterized by $\lambda_{k}\in\left[0,+\infty\right[$

$\underset{\begin{subarray}{c}{{\mathbf{x}}\in\mathbb{R}^{L}}\end{subarray}}{\text{\rm
argmin}}\;\;\frac{1}{2}\|{\mathbf{H}}{\mathbf{x}}-{\mathbf{y}}\|^{2}_{{\mathbf{F}}}+\lambda_{k}r_{1}({\mathbf{L}}_{1}{\mathbf{x}}),$
(44)

where $r_{1}({\mathbf{L}}_{1}\cdot)$ is a special case of the DTV, reducing to
the isotropic TV.
The authors perform $N=1$ iteration of the forward-backward algorithm for each
problem (44) with $K=20$. The initial value for the image with regularizing
constant $\lambda_{k}$ is chosen as the previous iterate obtained with
$\lambda_{k-1}$, and the values $(\lambda_{k})_{k=1}^{K}$ are linearly
decreased from 100 to 0. This strategy reconstructs intense objects first and
then anatomical details, hence its name “hierarchical reconstruction”. Note
that it does not address statistical noise but subsampling only. For our
comparison, we use 100 inner iterations of the dual forward-backward algorithm
to compute the proximity operator of TV regularization (with warm restart),
within the forward-backward update.

####  VI-C2 Post-processing U-net

The second competing method is the CNN proposed in [17, 16]. This method is a
post-processing of FBP, which consists of computing ${\mathbf{H}}_{\rm
ROI}^{\top}{\mathbf{F}}{\mathbf{y}}$. It relies on a trained residual U-net,
with a depth of 4 levels, filters of size 32, and a batch size of 15. Batch
normalization is useful to improve the stability of training.

####  VI-C3 Preconditioned Neumann network (PNN)

Our third competing method is a preconditioned Neumann network (PNN) initially
introduced in [52] for MRI reconstruction. It builds on a method for solving
problem (2) with
$f({\mathbf{x}})=\frac{1}{2}\|{\mathbf{H}}{\mathbf{x}}-{\mathbf{y}}\|^{2}$.
For a differentiable function $r$, the resulting minimizer reads

$({\mathbf{H}}^{\top}{\mathbf{H}}+\nabla
r){{\mathbf{x}}}={\mathbf{H}}^{\top}{\mathbf{y}},$ (45)

which can be rewritten as

$({\mathbf{H}}^{\top}{\mathbf{H}}+\lambda\operatorname{Id}\,){{\mathbf{x}}}+(\nabla
r-\lambda\operatorname{Id}\,){{\mathbf{x}}}={\mathbf{H}}^{\top}{\mathbf{y}}.$
(46)

Setting
${\mathbf{T}}_{\lambda}=({\mathbf{H}}^{\top}{\mathbf{H}}+\lambda\operatorname{Id}\,)^{-1}$
yields

$(\operatorname{Id}\,-\lambda{\mathbf{T}}_{\lambda}+{\mathbf{T}}_{\lambda}\nabla
r){{\mathbf{x}}}={\mathbf{T}}_{\lambda}{\mathbf{H}}^{\top}{\mathbf{y}}.$ (47)

Using the Neumann identity
${\mathbf{B}}^{-1}=\sum_{n=0}^{\infty}(\operatorname{Id}\,-{\mathbf{B}})^{n}$,
the authors derive the architecture of PNN with $N\in\mathbb{N}^{*}$ layers
(see Figure 6)

$(\lambda{\mathbf{T}}_{\lambda}-{\mathbf{T}}_{\lambda}\nabla
r)^{N}\circ{\mathbf{T}}_{\lambda}({\mathbf{H}}^{\top}{\mathbf{y}}).$ (48)

All instances of ${\mathbf{T}}_{\lambda}$ are applied approximately using an
unrolling of 10 iterations of the conjugate gradient algorithm.
The operator ${\mathbf{T}}_{\lambda}\nabla r$ is replaced by a U-net, denoted
by $\Psi$ which has the same architecture as the aforementioned U-net without
the residual connection. The weights of the U-net are shared for all layers.
Following [52], no batch normalization is used.
We choose $N=3$, the batch size is 6, the initial value for $\lambda$ is
$0.01$, the inner U-net has a depth of 4, and the learning rate is set to
$10^{-4}$.

Figure 6: Architecture of PNN [52]: The network maps a linear function of the
measurements ${\mathbf{T}}_{\lambda}{\mathbf{H}}^{\top}{\mathbf{y}}$ to a
reconstruction ${\mathbf{x}}_{n}$ by successive applications of an operator of
the form $\lambda{\mathbf{T}}_{\lambda}-\Psi$, while summing the intermediate
outputs of each block. All instances of ${\mathbf{T}}_{\lambda}$ are replaced
by an unrolling of 10 iterations of the conjugate gradient algorithm. $\Psi$
is a trained network and the scale parameter $\lambda$ is also trained.

####  VI-C4 ISTA-net

Our fourth competing method is ISTA-net, which is derived from the work of
[26]. ISTA-net is designed to solve Problem (2) for
$f({\mathbf{x}})=\frac{1}{2}\|{\mathbf{H}}{\mathbf{x}}-{\mathbf{y}}\|^{2}$,
and $r({\mathbf{x}})=\lambda\|{\mathbf{W}}{\mathbf{x}}\|_{1}$ ($\lambda>0$),
where operator ${\mathbf{W}}$ is not known a priori but learned.
${\mathbf{W}}$ is an orthogonal linear operator in the initial ISTA algorithm,
whose iteration reads

${\mathbf{x}}_{n+1}={\mathbf{W}}^{\top}\operatorname{soft}\left({\mathbf{W}}({\mathbf{x}}_{n}-\tau{\mathbf{H}}^{\top}({\mathbf{H}}{\mathbf{x}}_{n}-{\mathbf{y}}),\lambda\tau\right),$
(49)

where $\operatorname{soft}$ is the soft-thresholding operation and $\tau>0$ is
the gradient step size. In ISTA-net, the authors replace ${\mathbf{W}}$ and
${\mathbf{W}}^{\top}$ by two decoupled nonlinear operators namely
${\mathbf{A}}_{n}\circ\operatorname{ReLU}\circ{\mathbf{B}}_{n}$ and
${\mathbf{C}}_{n}\circ\operatorname{ReLU}\circ{\mathbf{D}}_{n}$ (see Figure
7). The property of orthogonality of ${\mathbf{W}}$ is not imposed but favored
during training by adding a term, weighted by $\chi\in\left]0,+\infty\right[$,
penalizing the difference between
$({\mathbf{C}}_{n}\circ\operatorname{ReLU}\circ{\mathbf{D}}_{n})\circ({\mathbf{A}}_{n}\circ\operatorname{ReLU}\circ{\mathbf{B}}_{n}){\mathbf{x}}_{n}$
and ${\mathbf{x}}_{n}$ in the loss function. Each ${\mathbf{A}}_{n}$,
${\mathbf{B}}_{n}$, ${\mathbf{C}}_{n}$ and ${\mathbf{D}}_{n}$ is a 2D
convolutional operator. ${\mathbf{B}}_{n}$ and ${\mathbf{C}}_{n}$ are
associated with a kernel of size $3\times 3$ and 16 input and output channels;
${\mathbf{A}}_{n}$ has 1 input channel and 16 output channels and vice-versa
for ${\mathbf{D}}_{n}$. As suggested by the authors, we learn these
convolutional operators as well as $\lambda$ and $\tau$, which are allowed to
vary at each iteration.
Experiments are carried out with 10 layers, $\chi=0.1$, a batch size of 14 and
the learning rate is initialized to $10^{-3}$ (step-wise decreased by a factor
$0.9$ every $20$ epochs). The gradient norm is clipped to 1, and
${\mathbf{x}}_{0}$ is set at zero.

Figure 7: Architecture of ISTA-net [26]: Each layer is composed of a gradient
step followed by the application of a nonlinear operator which is the
combination of two learnable linear convolutional operators
(${\mathbf{A}}_{n}$, ${\mathbf{B}}_{n}$) separated by a ReLU, a soft-
thresholding operation and then two other learnable linear convolutional
operators (${\mathbf{C}}_{n}$, ${\mathbf{D}}_{n}$) separated by a ReLU. The
property
$({\mathbf{C}}_{n}\circ\operatorname{ReLU}\circ{\mathbf{D}}_{n})\circ({\mathbf{A}}_{n}\circ\operatorname{ReLU}\circ{\mathbf{B}}_{n})=\operatorname{Id}\,$
is favored during training.

####  VI-C5 PD-net

The last competing method is the learned Primal-Dual (PD-net) introduced in
[53] by unrolling the Primal-Dual Hybrid Gradient (PDHG) optimization
algorithm[54]. The authors consider problem (2) with a more generic data
fidelity term $f({\mathbf{x}})=G({\mathbf{H}}{\mathbf{x}};{\mathbf{y}})$. They
replace both the proximity operators of $G$ and $r$ in PDHG by residual CNN so
that one layer $n\in\mathbb{N}$ of their network reads

$\displaystyle{\mathbf{z}}_{n+1}$
$\displaystyle=\operatorname{CNN}({\mathbf{z}}_{n}+\sigma{\mathbf{H}}\tilde{{\mathbf{x}}}_{n};{\mathbf{y}})$
(50) $\displaystyle{\mathbf{x}}_{n+1}$
$\displaystyle=\operatorname{CNN}({\mathbf{x}}_{n}-\tau{\mathbf{H}}^{\top}{\mathbf{z}}_{n+1})$
(51) $\displaystyle\tilde{{\mathbf{x}}}_{n+1}$
$\displaystyle={\mathbf{x}}_{n+1}+\gamma({\mathbf{x}}_{n+1}-{\mathbf{x}}_{n}).$
(52)

The CNNs act both in the image and projection domains. Furthermore, buffers of
previous iterates of size $N_{p}\in\mathbb{N}$ in the primal domain (image)
and of size $N_{d}\in\mathbb{N}$ in the dual domain (projection) are kept to
enable the network to learn an acceleration. For our implementation, we use a
ReLU non-linearity instead of a PReLU. As advocated by the authors, the
parameters used here are $N=5$ layers, $N_{d}=N_{p}=3$, and 32 filters in the
convolutional layers. This network is illustrated in Figure 8. The batch size
is 20, and the learning rate is initialized to $10^{-3}$ (step-wise decreased
by a factor $0.9$ every $20$ epochs). The gradient norm is clipped to 1.

As advocated in their respective papers, the competing networks were also
trained with the MSE loss for 100 epochs in a standard end-to-end manner. Note
that codes are publicly available for these networks. We re-implemented them
in Pytorch and re-tuned some of the hyper-parameters starting from the values
proposed in the papers so that the number of parameters stays in the same
order of magnitude for all networks (see table I).

Figure 8: Architecture of PD-net [53]: The red and blue boxes represent the primal and dual networks, respectively. Buffers of 3 primal $({\mathbf{x}}_{n}^{1},{\mathbf{x}}_{n}^{2},{\mathbf{x}}_{n}^{3})$ and dual $({\mathbf{z}}_{n}^{1},{\mathbf{z}}_{n}^{2},{\mathbf{z}}_{n}^{3})$ estimates are used at each iteration. The initial primal estimates are set to the FBP reconstruction given by ${\mathbf{H}}^{\top}{\mathbf{F}}{\mathbf{y}}$ and the initial dual estimates are set to zero. | $\mathcal{A}_{1}$ | $\mathcal{A}_{2}$ | U-net | PNN | PD-net | ISTA-net
---|---|---|---|---|---|---
$|\Theta|$ | 21 242 | 31 311 | 483 153 | 447 214 | 117 393 | 48 980
TABLE I: Number of learnable parameters ($\Theta$)

##  VII Results

###  VII-A Comparison with non-deep learning methods and use of an M-estimator

Figure 9: Full ground truth $\overline{{\mathbf{x}}}_{\rm P}$: Example A
(left) and Example B (right). Figure 10: Evolution of the PSNR across layers
at the end of the end-to-end training phase of $\mathcal{A}_{2}$ trained on
Example B.

(a) $\overline{{\mathbf{x}}}_{\rm ROI}$

(b) $\mathcal{A}_{2}$

(c) TV

(d) FBP

Figure 11: Reconstructed ROIs of Example A (top) and Example B (bottom)
obtained by two overfitting architectures of type $\mathcal{A}_{2}$, the
hierarchical TV method and FBP.

First, we compare architecture $\mathcal{A}_{2}$ with the hierarchical method
and FBP. In $\mathcal{A}_{2}$, the mismatched adjoint is not learned; it is
such that the ramp filter is applied only once between ${\mathbf{H}}$ and
${\mathbf{H}}^{\top}$, as is the case for the hierarchical method. However,
unlike the latter, $\mathcal{A}_{2}$ takes into account the location of the
ROI through its loss. For the comparison, we thus designed an overfitting
network by training $\mathcal{A}_{2}$ on one instance of Dataset 1 at a time
for two instances (Example A and Example B, see figure 9). These examples
illustrate two phenomena. In Example A, the exterior of the ROI, which
contributes to the truncated projections, is made of anatomical tissues only.
In Example B, the projections contain the contribution of intense objects
outside the reconstruction grid.
Figure 10 shows the evolution of the PSNR for Example B along the 28 layers of
the same architecture. These plots show that, for our choice of $K$ and $N$,
the network can learn an optimized set of parameters such that the PSNR
gradually increases across each layer, though less rapidly in the last layers.
Given that larger values for $K$ and $N$ imply a higher memory and training
time, increasing the total number of layers here appeared unnecessary.
Figure 12 shows the learned scalar parameters in the data and regularization
layers for Example B. Interestingly, the regularization step sizes
($b_{k,n,j}$) tend to decrease for the neighbors of $x_{\ell}$ (see Figure 2),
except for its two furthest neighbors. The weight of the data term ($d_{k,n}$)
consistently decreases after the second layer and vice-versa for the
associated convergence parameter ($a_{k,n}$). The step sizes associated with
the data layer thus follow an opposite trend to those associated with the
regularization. This indicates that the number of regularization layers may be
too low, forcing the network to accelerate the convergence of the variables
$({\mathbf{z}}^{j}_{n})_{j=1}^{J}$ (high values for
$\gamma_{n}\tau_{j}^{-1}$).
Figure 13 displays the learned weights of the fully connected layers, which
infer the Cauchy function threshold from the reprojection error cumulative
histogram. Figure 11 shows reconstructed ROIs obtained with $\mathcal{A}_{2}$
(overfitted to each image) (b), the hierarchical TV method (c), and FBP (d).
In Example A, the sub-sampling streaks that appear in the FBP reconstruction
are not present anymore in the images produced by $\mathcal{A}_{2}$ and the TV
method. More precisely, our architecture $\mathcal{A}_{2}$ leads to similar
results as the TV method in a very reduced number of iterations. If the
hierarchical method uses few projection-backprojection operations, it requires
2000 sub-iterations in total to compute the proximity operator of the TV
regularization. In contrast, U-RDBDB only uses 25 regularization layers.
In Example B, an intense wire exists in the data but not in the reconstruction
grid, which produces high-intensity sub-sampling streaks in the FBP
reconstruction. Most of the streaks remain in the TV solution, whereas our
architecture $\mathcal{A}_{2}$ strongly attenuates them. In addition, the
noise is reduced, but the resolution of the image is decreased. Because
$\mathcal{A}_{2}$ is designed to minimize the $\ell_{2}$ norm of the error
with the noiseless ground truth and a very low number of regularization
layers, it tends to selectively smooth some parts of the image to remove
remaining streaks. On the contrary, the hierarchical method is data-agnostic;
it is designed to reduce sub-sampling streaks only, regardless of the noise.
Hence it is interesting to check whether an adjoint mismatch can help reduce
this smoothing effect and compensate for the low number of regularization
layers.

(a) $\operatorname{softplus}(b_{k,n,1})$

(b) $\operatorname{softplus}(b_{k,n,2})$

(c) $\operatorname{softplus}(b_{k,n,3})$

(d) $\operatorname{softplus}(b_{k,n,4})$

(e) $\operatorname{softplus}(b_{k,n,5})$

(f) $\operatorname{softplus}(b_{k,n,6})$

(g) $\operatorname{softplus}(e_{k,n})$

(h) $\operatorname{softplus}(d_{k,n})$

(i) $\operatorname{softplus}(a_{k,n})$

Figure 12: Evolution of learned scalar parameters of $\mathcal{A}_{2}$ along
the layers, when the network is trained to overfit Example B. Figure 13:
Learned weights of the fully connected layer $c_{k,n}$, shared across the
layers of $\mathcal{A}_{2}$, when the network is trained to overfit Example B.

###  VII-B Comparison with deep learning methods and impact of an adjoint
mismatch

Metrics | $\mathcal{A}_{1}$ | $\mathcal{A}_{2}$ | $\mathcal{A}_{3}$ | U-net | PPN | PD-net | ISTA-net | TV
---|---|---|---|---|---|---|---|---
PSNR | 38.5 | 40.5 | 42.1 | 40.0 | 40.3 | 39.7 | 35.3 | 33.9
SSIM | 0.959 | 0.970 | 0.974 | 0.960 | 0.936 | 0.952 | 0.959 | 0.903
MAE ($\times 10^{-3}$) | 6.65 | 4.32 | 3.12 | 5.34 | 7.80 | 5.68 | 5.62 | 12.5
PieApp | 0.592 | 0.421 | 0.389 | 0.622 | 0.791 | 0.442 | 0.651 | 0.653
TABLE II: Quantitative assessment of the reconstructed ROIs. Mean values computed over the test set of Dataset 1. Metrics | $\mathcal{A}_{1}$ | $\mathcal{A}_{2}$ | $\mathcal{A}_{3}$ | U-net | PPN | PD-net | ISTA-net | TV
---|---|---|---|---|---|---|---|---
PSNR | 36.2 | 41.0 | 42.8 | 40.8 | 35.0 | 36.6 | 33.5 | 28.1
SSIM | 0.967 | 0.971 | 0.979 | 0.978 | 0.972 | 0.977 | 0.958 | 0.948
MAE ($\times 10^{-3}$) | 4.05 | 3.81 | 3.12 | 3.74 | 7.98 | 6.87 | 11.6 | 16.5
PieApp | 0.362 | 0.289 | 0.267 | 0.310 | 0.403 | 0.269 | 0.438 | 0.476
TABLE III: Quantitative assessment of the reconstructed ROIs. Mean values
computed over the test set of Dataset 2.

We now compare all deep learning methods trained on the union of the training
sets of Dataset 1 and Dataset 2 ($I=5000$). Qualitative and quantitative
comparisons are made on test sets. Table II reports the performance of our
three architectures, U-net, and other deep unfolding networks on the first
dataset. The metrics obtained with the iterative hierarchical TV algorithm are
also added as a reference. First, we see that all supervised methods can
perform at least as well as the hierarchical method in a reduced number of
iterations (i.e., layers). Then, in both cases, $\mathcal{A}_{3}$ outperforms
$\mathcal{A}_{1}$ and $\mathcal{A}_{2}$. Moreover, both $\mathcal{A}_{2}$ and
$\mathcal{A}_{3}$ perform on average better than the other unfolding networks
(PNN, PD-net, and ISTA-net) and U-net for all considered metrics. The results
obtained with $\mathcal{A}_{1}$ are characterized by lower metrics in all
cases. This thus provides empirical evidence that performance can be optimized
by introducing mismatched adjoints even without learning.
Figure 14 illustrates the reconstruction results for four examples from the
test set of Dataset 1. The figure confirms that $\mathcal{A}_{2}$ and
$\mathcal{A}_{3}$ reduce streaks more effectively than the other unfolding
networks. In the last-row example, even though $\mathcal{A}_{2}$ removes
streaks over the colon, it is done at the expense of over-smoothing, which
does not appear with $\mathcal{A}_{3}$. Some images obtained with
$\mathcal{A}_{1}$ contain artificial background patterns, translating into
degradation for all metrics. At first sight, the images produced by U-net have
fewer artifacts than most deep unfolding networks. However, in the second- and
third-row images, U-net introduces an artificial background to the images; in
the second image, the background appears darker than in the ground truth. This
observation highlights that U-net can hallucinate structures under the sub-
sampling streaks of the FBP input. Unfolding networks avoid these
hallucinations; by simply alternating between U-net and several consistency
layers, PNN already minimizes this effect.

In Table III, we report the performance of the trained networks when tested on
Dataset 2. Both our architectures outperform the competing unfolding networks
as soon as we include an adjoint mismatch in the data layer. U-net ranks
third-best, while PD-net is fourth.
The four reconstructed images displayed in Figure 15 confirm this trend. U-net
effectively removes streaks in the first, second, and last images, but it
over-smoothes soft tissues, and an offset remains in some images. PNN manages
to correct these offsets but at the price of additional streaks artifacts.
Regarding our architectures, as for Dataset 1, the anatomical background is
better recovered when an adjoint mismatch is considered in the data layer
($\mathcal{A}_{2}$ and $\mathcal{A}_{3}$). Still, the images of the first and
third rows obtained with $\mathcal{A}_{2}$ have a patchy look, often
characteristic of TV regularization. This shows that $\mathcal{A}_{2}$ retains
the characteristics of the original optimization problem, avoiding the
generation of unexpected content as is possible with U-net. Learning the
adjoint of the regularization operators in $\mathcal{A}_{3}$ allows for more
flexibility. It modifies the piece-wise constant prior towards a prior that is
learned from the database with a potential acceleration of the iterative
reconstruction: the images do not display patchy artifacts. A visual
inspection of the images in the last row shows that none of the methods can
accurately recover the low contrasts in the brain. This suggests that the
number of projections is still too low to enable the visualization of such
tissues.

(a) $\overline{{\mathbf{x}}}_{\rm ROI}$

(b) U-net

(c) ISTA-net

(d) PD-net

(e) PNN

(f) $\mathcal{A}_{1}$

(g) $\mathcal{A}_{2}$

(h) $\mathcal{A}_{3}$

Figure 14: Reconstructed ROIs using different deep learning methods on four
examples in the test set of Dataset 1.

(a) $\overline{{\mathbf{x}}}_{\rm ROI}$

(b) U-net

(c) ISTA-net

(d) PD-net

(e) PNN

(f) $\mathcal{A}_{1}$

(g) $\mathcal{A}_{2}$

(h) $\mathcal{A}_{3}$

Figure 15: Reconstructed ROIs using different deep learning methods on four
examples in the test set of Dataset 2.

##  VIII Conclusion

In this paper, we benefited from the structure of the dual block forward-
backward algorithm to define two types of deep learning layers. We proposed
U-RDBFB, an architecture built from a nonuniform distribution of these two
layers, for reconstructing a region of interest from a limited number of
truncated CT projections. Based on the Cauchy function and the ramp filter, a
data consistency term and a non-local TV regularization are embedded to limit
sub-sampling streaks from intense objects outside the reconstruction grid. Our
data term allowed for fitting most data except for the projections of objects
with high gradients outside the reconstruction grid. The patchy artifacts
yielded by sparse priors on the non-sparse CT images for a very low number of
iterations were shown to be mitigated by learning layer-wise adjoints of the
regularization operators. Our experiments demonstrated that despite having
fewer learnable parameters, U-RDBFB performs favorably compared to state-of-
the-art deep learning methods. Our results suggest that including more a
priori knowledge can bring more robustness and reliability to the current
deep-learning methods. Future work will investigate the robustness of U-RDBFB
against changes in the acquisition settings and the trade-off between the
image quality in terms of resolution and the total number of layers of
U-RDBFB.

##  Appendix A Reweighting strategy

Let $\phi$ be given by (3). It was shown in [55] that, for every
$\overline{\zeta}\in\mathbb{R}$, the following convex quadratic function
$\tilde{\phi}(\cdot,\overline{\zeta})$, defined for every $\zeta\in\mathbb{R}$
as

$\tilde{\phi}(\zeta,\overline{\zeta})=\phi(\overline{\zeta})+\beta\frac{(\zeta-\overline{\zeta})\overline{\zeta}}{1+(\overline{\zeta}/\kappa)^{2}}+\frac{\beta}{2}\frac{(\zeta-\overline{\zeta})^{2}}{(1+(\overline{\zeta}/\kappa)^{2})},$
(53)

is a tangent majorant approximation to $\phi$ at $\overline{\zeta}$, that is

$(\forall\zeta\in\mathbb{R})\;\;\tilde{\phi}(\zeta,\overline{\zeta})\geqslant\phi(\zeta)\quad\text{
and
}\quad\tilde{\phi}(\overline{\zeta},\overline{\zeta})=\phi(\overline{\zeta}).$
(54)

This allows us to deduce a tangent majorant function $\tilde{g}$ of function
$g$ at any point $\overline{{\mathbf{z}}}\in\mathbb{R}^{T}$:
$(\forall{\mathbf{z}}\in\mathbb{R}^{T})$:

$\displaystyle\tilde{g}({\mathbf{z}},\overline{{\mathbf{z}}})=\sum_{t=1}^{T}\phi({\mathbf{z}}_{t},\overline{{\mathbf{z}}}_{t})$
$\displaystyle=g(\overline{{\mathbf{z}}})+\beta\operatorname{diag}\left(\left(\frac{\overline{{\mathbf{z}}}_{t}}{1+(\overline{{\mathbf{z}}}_{t}/\kappa)^{2}}\right)_{t=1}^{T}\right)({\mathbf{z}}-\overline{{\mathbf{z}}})$
$\displaystyle+\frac{\beta}{2}({\mathbf{z}}-\overline{{\mathbf{z}}})^{\top}\operatorname{diag}\left(\left(\frac{1}{1+(\overline{{\mathbf{z}}}_{t}/\kappa)^{2}}\right)_{t=1}^{T}\right)({\mathbf{z}}-\overline{{\mathbf{z}}})\geqslant
g({\mathbf{z}}).$

Finally, for every $\overline{{\mathbf{x}}}\in\mathbb{R}^{L}$, we set

$(\forall{\mathbf{x}}\in\mathbb{R}^{L})\quad\tilde{f}({\mathbf{x}},\overline{{\mathbf{x}}})=\tilde{g}({\mathbf{F}}({\mathbf{H}}{\mathbf{x}}-{\mathbf{y}});{\mathbf{F}}({\mathbf{H}}\overline{{\mathbf{x}}}-{\mathbf{y}})),$
(55)

that satisfies $\tilde{f}({\mathbf{x}},\overline{{\mathbf{x}}})\geqslant
g({\mathbf{F}}({\mathbf{H}}{\mathbf{x}}-{\mathbf{y}}))=f({\mathbf{x}})$. Given
this majoration, the iterative reweighting strategy approximates the solution
to (8) by the estimate produced by Algorithm 1, where

$\displaystyle(\forall({\mathbf{x}},\overline{{\mathbf{x}}})\in(\mathbb{R}^{L})^{2})\quad
Q({\mathbf{x}},\overline{{\mathbf{x}}})=$
$\displaystyle\tilde{f}({\mathbf{x}},\overline{{\mathbf{x}}})+r({\mathbf{x}}).$
(56)

##  Appendix B Implementation of DBFB

The DBFB iterates $({\mathbf{x}}_{n})_{n\in\mathbb{N}}$ converge to the
minimizer of (9) under the following assumptions on the algorithm parameters
[37]:

$\begin{cases}\sigma\geqslant|||{\mathbf{B}}_{0}{\mathbf{M}}^{-1}{\mathbf{B}}_{0}^{\top}|||,\\\
(\forall
j\in\\{1,\ldots,J\\})\quad\tau_{j}\geqslant|||{\mathbf{L}}_{j}{\mathbf{M}}^{-1}{\mathbf{L}}_{j}^{\top}|||,\\\
\gamma_{n}\in[\epsilon,2-\epsilon]\text{ with }\epsilon\in]0,1]\\\ (\exists
N\in\mathbb{N}\setminus\\{0,1\\})(\forall n\in\mathbb{N})\quad
0<\sum_{n^{\prime}=n}^{n+N-1}\varepsilon_{n^{\prime}}<N.\end{cases}$ (57)

The first three assumptions are step size range conditions. The last one means
that each step (D) and (R) is performed at least once, every $N$ iterations.
Step (D) involves the calculation of the proximity operator $\text{\rm
prox}_{\gamma_{n}^{-1}\sigma
h_{0}(\cdot,{\mathbf{A}}_{0}\overline{{\mathbf{x}}})}$, which has a closed-
form [29, Example 24.2], for
$(\forall({\mathbf{z}},\overline{{\mathbf{z}}})\in(\mathbb{R}^{T})^{2})$,

$\displaystyle\text{\rm prox}_{\gamma_{n}^{-1}\sigma
h_{0}(\cdot,{\mathbf{A}}_{0}\overline{{\mathbf{z}}})}({\mathbf{z}})=\text{\rm
prox}_{\gamma_{n}^{-1}\sigma\tilde{g}(\cdot-{\mathbf{F}}{\mathbf{y}},\overline{{\mathbf{z}}}-{\mathbf{F}}{\mathbf{y}})}({\mathbf{z}}),$
$\displaystyle={\mathbf{F}}{\mathbf{y}}+\text{\rm
prox}_{\gamma_{n}^{-1}\sigma\tilde{g}(\cdot,\overline{{\mathbf{z}}}-{\mathbf{F}}{\mathbf{y}})}({\mathbf{z}}-{\mathbf{F}}{\mathbf{y}})$
$\displaystyle=\left(({\mathbf{F}}{\mathbf{y}})_{t}+\frac{{\mathbf{z}}_{t}-({\mathbf{F}}{\mathbf{y}})_{t}}{{1+}\beta\gamma_{n}^{-1}\sigma\big{(}1+(\overline{{\mathbf{z}}}_{t}-({\mathbf{F}}{\mathbf{y}})_{t})^{2}/\kappa^{2}\big{)}^{-1}}\right)_{t=1}^{T}.$
(58)

Step (R) requires calculating the proximity operator of $h_{j}$ scaled by
parameter $\gamma\in\left]0,+\infty\right[$. It also has a closed form: for
${\mathbf{s}}=({\mathbf{s}}_{1},\ldots,{\mathbf{s}}_{J})\in\mathbb{R}^{2JL}$,
$\text{\rm prox}_{\gamma h_{1}}({\mathbf{s}})=\left(\text{\rm prox}_{\gamma
r_{j}}({\mathbf{s}}_{j})\right)_{j=1}^{J}$, where, for every
${\mathbf{z}}=({\mathbf{z}}_{1},{\mathbf{z}}_{2})\in\mathbb{R}^{2L}$,

$\text{\rm prox}_{\gamma
r_{j}}({\mathbf{z}})=\left(\max\Big{\\{}0,1-\frac{\gamma\alpha_{j,l}}{\|{\mathbf{z}}_{\ell}\|_{2}}\Big{\\}}{\mathbf{z}}_{\ell}\right)_{\ell=1}^{L},$
(59)

where, for every $\ell\in\\{1,\ldots,L\\}$,
${\mathbf{z}}_{\ell}=(({\mathbf{z}}_{1})_{\ell},({\mathbf{z}}_{2})_{\ell})\in\mathbb{R}^{2}$.

## References

  * [1] R. Chityala, K. Hoffmann, D. Bednarek, and S. Rudin.  Region of interest (ROI) computed tomography.  In Medical Imaging 2004: Physics of Medical Imaging, vol. 5368, pp. 534–541. SPIE, 2004.
  * [2] L. Yu, X. Liu, S. Leng, J. Kofler, J. Ramirez-Giraldo, M. Qu, J. Christner, J. Fletcher, and C. McCollough.  Radiation dose reduction in computed tomography: techniques and future perspective.  Imaging in Medicine, vol. 1, pp. 65–84, 2009.
  * [3] E. Sidky and X. Pan.  Image reconstruction in circular cone-beam computed tomography by constrained, total-variation minimization.  Physics in Medicine & Biology, vol. 53, pp. 4777–4807, 2008.
  * [4] J. Bian, J. Siewerdsen, X. Han, E. Sidky, J. Prince, C. Pelizzari, and X. Pan.  Evaluation of sparse-view reconstruction from flat-panel-detector cone-beam CT.  Physics in Medicine and Biology, vol. 55, pp. 6575–6599, 2010.
  * [5] L. Rudin, S. Osher, and E. Fatemi.  Nonlinear total variation based noise removal algorithms.  Physica D: Nonlinear Phenomena, vol. 60, pp. 259–268, 1992.
  * [6] H. Langet, C. Riddell, A. Reshef, Y. Trousset, A. Tenenhaus, E. Lahalle, G. Fleury, and N. Paragios.  Compressed-sensing-based content-driven hierarchical reconstruction: Theory and application to C-arm cone-beam tomography.  Medical Physics, vol. 42, pp. 5222–5237, 2015.
  * [7] H. Yu and G. Wang.  Compressed sensing based interior tomography.  Physics in Medicine & Biology, vol. 54, pp. 2791, 2009.
  * [8] P. Paleo and A. Mirone.  Efficient implementation of a local tomography reconstruction algorithm.  Advanced Structural and Chemical Imaging, vol. 3, pp. 5, 2017.
  * [9] B. Hamelin, Y. Goussard, J.-P. Dussault, G. Cloutier, G. Beaudoin, and G. Soulez.  Design of iterative ROI transmission tomography reconstruction procedures and image quality analysis.  Medical Physics, vol. 37, pp. 4577–4589, 2010.
  * [10] A. Ziegler, T. Nielsen, and M. Grass.  Iterative reconstruction of a region of interest for transmission tomography.  Medical Physics, vol. 35, pp. 1317–1327, 2008.
  * [11] M. Savanier, E. Chouzenoux, J.-C. Pesquet, and C. Riddell.  Unmatched preconditioning of the proximal gradient algorithm.  IEEE Signal Processing Letters, vol. 29, pp. 1122–1126, 2022.
  * [12] R. Tirer and R. Giryes.  On the convergence rate of projected gradient descent for a back-projection based objective.  SIAM Journal on Imaging Sciences, vol. 14, pp. 1504–1531, 2021.
  * [13] T. Tirer and R. Giryes.  Back-projection based fidelity term for ill-posed linear inverse problems.  IEEE Transactions on Image Processing, vol. 29, pp. 6164–6179, 2020.
  * [14] M. Nilchian, C. Vonesch, P. Modregger, M. Stampanoni, and M. Unser.  Iterative FBP for improved reconstruction of X-ray differential phase-contrast tomograms.  In Proceedings of the 2013 IEEE 10th International Symposium on Biomedical Imaging, pp. 1260–1263, 7-11 April 2013, San Francisco, CA, USA.
  * [15] O. Ronneberger, P. Fischer, and T. Brox.  U-Net: Convolutional Networks for Biomedical Image Segmentation.  In Proceedings of the 18th International Conference on Medical Image Computing and Computer Assisted Intervention (MICCAI) 2015, pp. 234–241, 5-9 October 2015, Munich, Germany.
  * [16] Y. Han, J. J. Yoo, and J. C. Ye.  Deep residual learning for compressed sensing CT reconstruction via persistent homology analysis.  Preprint 2016. Available: https://arxiv.org/abs/1611.06391.
  * [17] K. Jin, M. T. McCann, E. Froustey, and M. Unser.  Deep convolutional neural network for inverse problems in imaging.  IEEE Transactions on Image Processing, vol. 26, pp. 4509–4522, 2017.
  * [18] E. Sidky, I. Lorente, J. G. Brankov, and X. Pan.  Do CNNs solve the CT inverse problem.  IEEE Transactions on Biomedical Engineering, vol. 68, pp. 1799–1810, 2021.
  * [19] V. Monga, Y. Li, and Y. Eldar.  Algorithm unrolling: Interpretable, efficient deep learning for signal and image processing.  IEEE Signal Processing Magazine, vol. 38, pp. 18–44, 2021.
  * [20] G. Wang, J. C. Ye, and B. De Man.  Deep learning for tomographic image reconstruction.  Nature Machine Intelligence, vol. 2, pp. 737–748, 2020.
  * [21] S. Wang, S. Fidler, and R. Urtasun.  Proximal Deep Structured Models.  In Proceedings of the 30th Annual Conference on Neural Information Processing Systems, NIPS 2016, pp. 865–873, 5-10 December 2016, Barcelona, Spain.
  * [22] H. T. V. Le, N. Pustelnik, and M. Foare.  The faster proximal algorithm, the better unfolded deep learning architecture ? The study case of image denoising.  In Proceedings of the 30th European Signal Processing Conference, EUSIPCO 2022, 29 August - 2 September 2022, Belgrade, Serbia.
  * [23] C. Bertocchi, E. Chouzenoux, M.-C. Corbineau, J.-C. Pesquet, and M. Prato.  Deep unfolding of a proximal interior point method for image restoration.  Inverse Problems, vol. 36, pp. 034005, 2020.
  * [24] Y. Yang, J. Sun, H. Li, and Z. Xu.  Deep ADMM-Net for Compressive Sensing MRI.  In Proceedings of the 30th International Conference on Neural Information Processing Systems, pp. 10–18, 5-10 December 2016, Barcelona, Spain.
  * [25] J. Adler and O. Oktem.  Solving ill-posed inverse problems using iterative deep neural networks.  Inverse Problems, vol. 33, pp. 124007, 2017.
  * [26] J. Zhang and B. Ghanem.  ISTA-Net: Interpretable optimization-inspired deep network for image compressive sensing.  Proceedings of the 2018 IEEE/CVF Conference on Computer Vision and Pattern Recognition, pp. 1828–1837, 18-23 June 2018, Salt Lake City, UT, USA.
  * [27] N. Pustelnik, A. Benazza-Benhayia, Y. Zheng, and J.-C. Pesquet.  Wavelet-Based Image Deconvolution and Reconstruction, pages 1–34.  American Cancer Society, 2016.
  * [28] F. Abboud, E. Chouzenoux, J.-C. Pesquet, J.-H. Chenot, and L. Laborelli.  Dual Block Coordinate Forward-Backward Algorithm with Application to Deconvolution and Deinterlacing of Video Sequences.  Journal of Mathematical Imaging and Vision, vol. 59, pp. 415–431, November 2017.
  * [29] H. H. Bauschke and P. L. Combettes.  Convex Analysis and Monotone Operator Theory in Hilbert Spaces.  2017\.
  * [30] N. Ouzir, A. Basarab, O. Lairez, and J.-Y. Tourneret.  Robust optical flow estimation in cardiac ultrasound images using a sparse representation.  IEEE Transactions on Medical Imaging, vol. 38, pp. 741–752, 2019.
  * [31] D. Kazantsev, F. Bleichrodt, T. van Leeuwen, A. Kaestner, P. J. Withers, K. J. Batenburg, and P. D. Lee.  A novel tomographic reconstruction method based on the robust student’s t function for suppressing data outliers.  IEEE Transactions on Computational Imaging, vol. 3, pp. 682–693, 2017\.
  * [32] D. Andrews and F. Hampel.  Robust estimates of location.  In Robust Estimates of Location. Princeton University Press, 1972\.
  * [33] G. L Zeng.  Revisit of the ramp filter.  In Proceedings of the 2014 IEEE Nuclear Science Symposium and Medical Imaging Conference, pp. 1–6, 8-15 November 2014, Seattle, WA USA.
  * [34] S. Kindermann, S. Osher, and P. W. Jones.  Deblurring and denoising of images by nonlocal functionals.  Multiscale Modeling & Simulation, vol. 4, pp. 1091–1115, 2006.
  * [35] G. Gilboa and S. Osher.  Nonlocal operators with applications to image processing.  Multiscale Modeling & Simulation, vol. 7, pp. 1005–1028, 2009.
  * [36] P. Combettes and J.-C. Pesquet.  Proximal splitting methods in signal processing.  In Fixed-point algorithms for inverse problems in science and engineering, pp. 185–212. Springer, 2011.
  * [37] F. Abboud, E. Chouzenoux, J.-C. Pesquet, J.-H. Chenot, and L. Laborelli.  Dual block-coordinate forward-backward algorithm with application to deconvolution and deinterlacing of video sequences.  Journal of Mathematical Imaging and Vision, vol. 59, pp. 415–431, 2017\.
  * [38] E. Chouzenoux, J.-C. Pesquet, C. Riddell, M. Savanier, and Y. Trousset.  Convergence of proximal gradient algorithm in the presence of adjoint mismatch.  Inverse Problems, vol. 37, pp. 065009, 2021.
  * [39] T. Elfving and P. C. Hansen.  Unmatched projector/backprojector pairs: perturbation and convergence analysis.  SIAM Journal on Scientific Computing, vol. 40, pp. A573–A591, 2018.
  * [40] D. Lorenz, S. Rose, and F. Schöpfer.  The randomized Kaczmarz method with mismatched adjoint.  BIT Numerical Mathematics, vol. 58, pp. 1079–1098, 2018.
  * [41] E. Chouzenoux, A. Contreras, J.-C. Pesquet, and M. Savanier.  Convergence Results for Primal-Dual Algorithms in the Presence of Adjoint Mismatch.  Preprint, 2022. Available: https://hal.archives-ouvertes.fr/hal-03654126/.
  * [42] G. Zeng.  Counter examples for unmatched projector/backprojector in an iterative algorithm.  Chinese Journal of Academic Radiology, 2019.
  * [43] G. Zeng and G. Gullberg.  Unmatched projector/backprojector pairs in an iterative reconstruction algorithm.  IEEE Transactions on Medical Imaging, 2000.
  * [44] E. Sidky, E. C. Hansen, J. Jorgensen, and X. Pan.  Iterative image reconstruction for CT with unmatched projection matrices using the generalized minimal residual algorithm.  Preprint, 2022. Available: https://arxiv.org/abs/2201.07408.
  * [45] Z. Wang, H. Li, W. Ouyang, and W. Wang.  Learnable histogram: Statistical context features for deep neural networks.  In Proceedings of the 14th European Conference on Computer Vision (ECCV), Amsterdam, The Netherlands, 2016.
  * [46] M. Ronchetti.  TorchRadon: Fast differentiable routines for computed tomography.  Preprint, 2020. Available: https://arxiv.org/abs/2009.14788.
  * [47] I. Safran and O. Shamir.  On the quality of the initial basin in overspecified neural networks.  In Proceedings of the 33rd International Conference on Machine Learning, pp. 774–782, June 19 – June 24, 2016, New York City, NY, USA.
  * [48] J. Martens and I. Sutskever.  Learning recurrent neural networks with hessian-free optimization.  In Proceedings of the 28th International Conference on Machine Learning, pp. , 28 June to 2 July 2011, Bellevue, Washington, USA.
  * [49] E. Prashnani, H. Cai, Y. Mostofi, and P. Sen.  Pieapp: Perceptual image-error assessment through pairwise preference.  In Proceedings of the 2018 IEEE Conference on Computer Vision and Pattern Recognition, pp. 1808–1817, 18-23 June 2018, Salt Lake City, UT, USA.
  * [50] D. Donoho and Y. Tsaig.  Fast solution of $\ell_{1}$ -norm minimization problems when the solution may be sparse.  IEEE Transactions on Information Theory, vol. 54, pp. 4789–4812, 2008\.
  * [51] L. Xiao and T. Zhang.  A proximal-gradient homotopy method for the sparse least-squares problem.  SIAM Journal on Optimization, vol. 23, pp. 1062–1091, 2013.
  * [52] D. Gilton, G. Ongie, and R. Willett.  Neumann networks for linear inverse problems in imaging.  IEEE Transactions on Computational Imaging, vol. 6, pp. 328–343, 2020.
  * [53] J. Adler and O. Öktem.  Learned primal-dual reconstruction.  IEEE Transactions on Medical Imaging, vol. 37, pp. 1322–1332, 2018.
  * [54] J. Cheng, H. Wang, L. Ying, and D. Liang.  Model learning: Primal dual networks for fast MR imaging.  Preprint, 2019. Available: https://arxiv.org/abs/1908.02426.
  * [55] E. Chouzenoux, A. Jezierska, J.-C. Pesquet, and H. Talbot.  A majorize-minimize subspace approach for $\ell_{2}$-$\ell_{0}$ image regularization.  SIAM Journal on Imaging Sciences, vol. 6, pp. 563–591, 2013.