# 4D-image reconstruction directly from limited-angular-range data in
continuous-wave electron paramagnetic resonance imaging

Zheng Zhang  Boris Epel  Buxin Chen  Dan Xia  Emil Y. Sidky  Zhiwei Qiao
Howard Halpern  Xiaochuan Pan111Corresponding author. Email:
<EMAIL_ADDRESS>Department of Radiology, The University of Chicago,
Chicago, IL, USA  Department of Radiation & Cellular Oncology, The University
of Chicago, Chicago, IL, USA  School of Computer and Information Technology,
Shanxi University, Taiyuan, Shanxi, China

###### Abstract

Objective: We investigate and develop optimization-based algorithms for
accurate reconstruction of four-dimensional (4D)-spectral-spatial (SS) images
directly from data collected over limited angular ranges (LARs) in continuous-
wave (CW) electron paramagnetic resonance imaging (EPRI). Methods: Basing on a
discrete-to-discrete data model devised in CW EPRI employing the Zeeman-
modulation (ZM) scheme for data acquisition, we first formulate the image
reconstruction problem as a convex, constrained optimization program that
includes a data fidelity term and also constraints on the individual
directional total variations (DTVs) of the 4D-SS image. Subsequently, we
develop a primal-dual-based DTV algorithm, simply referred to as the DTV
algorithm, to solve the constrained optimization program for achieving image
reconstruction from data collected in LAR scans in CW-ZM EPRI. Results: We
evaluate the DTV algorithm in simulated- and real-data studies for a variety
of LAR scans of interest in CW-ZM EPRI, and visual and quantitative results of
the studies reveal that 4D-SS images can be reconstructed directly from LAR
data, which are visually and quantitatively comparable to those obtained from
data acquired in the standard, full-angular-range (FAR) scan in CW-ZM EPRI.
Conclusion: An optimization-based DTV algorithm is developed for accurately
reconstructing 4D-SS images directly from LAR data in CW-ZM EPRI. Future work
includes the development and application of the optimization-based DTV
algorithm for reconstructions of 4D-SS images from FAR and LAR data acquired
in CW EPRI employing schemes other than the ZM scheme. Significance: The DTV
algorithm developed may be exploited potentially for enabling and optimizing
CW EPRI with minimized imaging time and artifacts by acquiring data in LAR
scans.

###### keywords:

electron paramagnetic resonance imaging (EPRI) , spectral-spatial imaging ,
continuous-wave (CW) EPRI , Zeeman modulation , limited angular range ,
optimization-based reconstruction , primal-dual algorithm

††journal: Journal of Magnetic Resonance

##  1 Introduction

Electron paramagnetic resonance imaging (EPRI) has found applications to
biological and preclinical imaging [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]. In EPRI,
one seeks to determine quantitatively the four-dimensional (4D)-spectral-
spatial (SS) distribution of paramagnetic electron-spins of contrast
materials, referred to as the 4D-SS image, from which physical quantities of
biological or clinical interest such as oxygen concentration within the imaged
subject may be determined [11, 12, 13, 14, 15, 16]. Algorithms have been
developed for reconstructing 4D-SS images from data acquired in the full-
angular-range (FAR) scan in CW EPRI [17, 18, 19, 20]. It is of theoretical and
practical interest in reconstruction of 4D-SS images from data acquired in
limited-angular-range (LAR) scans in CW EPRI, and we thus investigate and
develop algorithms for accurate reconstruction of 4D-SS images directly from
data collected in LAR scans in CW EPRI.

Among various scanning schemes, the Zeeman modulation (ZM) scheme [21, 3, 22,
15] has been used to acquire data in CW EPRI. The filtered-backprojection
(FBP) algorithm is used currently for reconstructing the partial derivative
(PD) of a 4D-SS image with respect to its spectral dimension [15, 23] from
data collected in the FAR scan in CW-ZM EPRI; and carrying numerically out an
accumulative integration of the 4D-PD-SS image over the spectral dimension,
one subsequently obtains a 4D-SS image. In the work, we focus on the
development and investigation of an optimization-based algorithm for accurate
reconstruction of 4D-SS images directly from data acquired in FAR and LAR
scans in CW-ZM EPRI. The approach and algorithm developed can, however,
readily be extended to reconstruct 4D-SS images from data collected in FAR and
LAR scans in CW EPRI that employs other schemes than the ZM scheme for data
acquisition. We refer to data collected, respectively, in FAR and LAR scans
simply as FAR data and LAR data hereinafter.

One of the motivations of our work is based upon the observation of the
physical constraint on the gradient magnetic field, which is used for
achieving spatial encoding in CW EPRI. The gradient strength applied specifies
the tangent of the angle between the directions of the main magnetic field and
a sampling hyperplane in the 4D-SS space. When the angle is in the vicinity of
$\pm 90^{\circ}$, its tangent approaches infinity, and thus no physical
measurement at $\pm 90^{\circ}$ can be made because it requires an infinitely
large gradient strength, which is physically impossible. Therefore,
reconstruction of a 4D-SS image in practical CW EPRI can thus be an inherently
LAR problem. In particular, in practical CW-ZM EPRI [15], noise and other
physical perturbation in data acquired in the vicinity of $\pm 90^{\circ}$ may
be amplified in reconstruction of a 4D-SS image with existing algorithms such
as the FBP algorithm [24, 23], yielding a 4D-SS image of sub-optimal quality.

In CW-ZM EPRI, the standard data model that relates a 4D-PD-SS image to data
can be interpreted as the 4D-Radon transform of the 4D-PD-SS image [15], as
discussed in A, and existing algorithms such as the FBP algorithm are used
thus for reconstructing the 4D-PD-SS image from FAR data. There exists,
however, strong interests in obtaining 4D-SS images directly from FAR and LAR
data in CW-ZM EPRI. It remains unclear, to the best of our knowledge, that
algorithms exist for reconstruction of a 4D-SS image directly from FAR or LAR
data in CW-ZM EPRI. We investigate and develop such algorithms in the work.

We first formulate the image-reconstruction problem as a convex, constrained
optimization program, which includes a data-fidelity term and also constraints
on the directional total variations (DTVs) of the 4D-SS image. The choice of
the DTV constraints is motivated by a recent work that demonstrates the
potential utility of DTV constraints for accurate image reconstruction from
LAR data in computed tomography (CT) [25, 26]. We then develop a primal-dual-
based DTV algorithm, which is referred to simply as the DTV algorithm, to
reconstruct a 4D-SS image directly from FAR and LAR data by solving the
convex, constrained optimization problem formulated in CW-ZM EPRI.

We perform simulated- and real-data studies to verify and evaluate the
effectiveness and stability of the DTV algorithm developed for reconstruction
of 4D-SS images directly from FAR and LAR data acquired in CW-ZM EPRI. Results
of the studies reveal that the DTV algorithm developed can reconstruct
accurately 4D-SS images directly from FAR data, and also from LAR data with
visual appearance and numerical quantification comparable to their
counterparts obtained from FAR data in CW-ZM EPRI.

Following the introduction in Sec. 1, we discuss in Sec. 2 a data model in
practical CW-ZM EPRI, a formulation of the reconstruction problem as an
optimization program, the development of the DTV algorithm solving the
program, and strategies for verifying and evaluating the DTV algorithm.
Results of simulated- and real-data studies, including 4D-SS images
reconstructed by use of the DTV algorithm directly from FAR and LAR data, are
presented in Sec. 3, followed by discussion and conclusion in Secs. 4 and 5.

##  2 Materials and methods

Let three-dimensional (3D) vector $\vec{r}_{s}=(r_{s}{\rm cos}\phi\,{\rm
sin}\theta,r_{s}{\rm sin}\phi\,{\rm sin}\theta,r_{s}{\rm cos}\theta)^{\top}$
depicts a point in the 3D-spatial-coordinate system, where $r_{s}$ denotes the
magnitude of $\vec{r}_{s}$, “${\top}$” the transpose, $\theta$ the angle
between $z$-axis and $\vec{r}_{s}$, and $\phi$ the angle between the $x$-axis
and the projection of $\vec{r}_{s}$ onto the $x$-$y$ plane, as shown in Fig.
1a. Also, letting $B$ denote the magnetic field applied, we introduce 4D
vector $\vec{r}=(x,y,z,B)^{\top}$ to denote a point in 4D-image space. We use
a 4D unit-vector $\hat{\alpha}=({\rm cos}\phi\,{\rm sin}\theta\,{\rm
sin}\gamma,{\rm sin}\phi\,{\rm sin}\theta\,{\rm sin}\gamma,$ ${\rm
cos}\theta\,{\rm sin}\gamma,{\rm cos}\gamma)^{\top}$ to define the orientation
of a hyperplane in the 4D space, where $\gamma$ depicts the angle between the
$B$-axis and $\hat{\alpha}$, as shown in Fig. 1b. We use $f(\vec{r})$ to
denote a continuous 4D-SS image as 4D variable $\vec{r}$ can continuously
vary. Also, we use $g(\xi,\hat{\alpha})$ to denote a continuous 4D-data
function in a 4D space formed by scalar variable $\xi$ and 4D-unit-vector
$\hat{\alpha}$, where $\xi$ depicts the distance of a hyperplane of
orientation $\hat{\alpha}$ to the origin in the 4D-image space, and $\xi$ and
$\hat{\alpha}$ can vary continuously. As shown in Sec. A, for CW-ZM EPRI,
$g(\xi,\hat{\alpha})$ is related to the 4D-Radon transform, i.e., the
integration, of $\frac{\partial f(\vec{r})}{\partial B}$ computed over the
hyperplane, where $\frac{\partial f(\vec{r})}{\partial B}$ indicates the
partial derivative of continuous 4D-SS image $f(\vec{r})$ with respect to $B$.
As discussed in Sec. A, a continuous-to-continuous (CC)-data model shown in
Eq. (8) can be established for CW-ZM EPRI [15].

###  2.1 Discrete-to-discrete (DD)-data model in CW-ZM EPRI

In realistic CW-ZM EPRI, both image and data are presented on discrete arrays,
and one thus discretizes the variables of $f(\vec{r})$ and
$g(\xi,\hat{\alpha})$ onto discrete 4D-image and 4D-data arrays. Specifically,
in the work, variable $\vec{r}$ of $f(\vec{r})$ is discretized onto a discrete
array of $I$ identical 4D voxels, where $I=I_{x}\times I_{y}\times I_{z}\times
I_{B}$, and $I_{x}$, $I_{y}$, $I_{z}$, and $I_{B}$ depict, respectively, the
array sizes along dimensions $x$, $y$, $z$, and $B$. Each of the 4D voxels is
of identical size
$\Delta_{x}\times\Delta_{y}\times\Delta_{z}\times\Delta_{B}$, respectively,
filling completely in the 4D support of $L\times L\times L\times B_{w}$, i.e.,
intervals $\Delta_{x}={L}/I_{x}$, $\Delta_{y}={L}/I_{y}$,
$\Delta_{z}={L}/I_{z}$, and $\Delta_{B}={B}_{w}/I_{B}$. We use vector
$\mathbf{f}$ of size $I$ to denote the 4D-SS image represented on the 4D-image
array in a concatenated form in the order of $x$, $y$, $z$, and $B$, and entry
$i$ of vector $\mathbf{f}$ denotes the value of the 4D-SS image in voxel $i$,
where $i=1,2,...,I$. For simplicity, we refer to discrete 4D-SS image
$\mathbf{f}$ as an SS image hereinafter.

Similarly, variables $\xi$ and $\hat{\alpha}$ of 4D-data function
$g(\xi,\hat{\alpha})$ can be discretized onto a discrete array of $J$ grids,
where $J=J_{\alpha}\times J_{\xi}$, $J_{\alpha}$ denotes the total number of
discrete samples each of which is specified by $\hat{\alpha}$ on the
hypersphere in the 4D-image space, and let $J_{\xi}$ depict, for each
$\hat{\alpha}$, the total number of samples over
$\xi\in[-\xi_{\tau},\xi_{\tau}]$, where $2\xi_{\tau}$ depicts the maximum
support of the data function in the $\xi$-dimension. Also, using
$\gamma_{\tau}\leq 90^{\circ}$, $\theta_{\tau}\leq 90^{\circ}$, and
$\phi_{\tau}\leq 90^{\circ}$ to denote the maximum angles in a scan, we have
$\gamma\in[-\gamma_{\tau},\gamma_{\tau}]$,
$\theta\in[-\theta_{\tau},\theta_{\tau}]$, and
$\phi\in[-\phi_{\tau},\phi_{\tau}]$. Assuming that sample grids are obtained
by varying $\hat{\alpha}$ in the order of varying $\phi$, $\theta$, and
$\gamma$ and then by varying $\xi$, we can then use vectors $\mathbf{g}$ and
$\mathbf{g}^{[\mathcal{M}]}$ of size $J$ to denote discrete model data and
measured data, respectively, in a concatenated form in the order of $\phi$,
$\theta$, $\gamma$, and $\xi$. For simplicity, we refer henceforward to
$\mathbf{g}$ and $\mathbf{g}^{[\mathcal{M}]}$, as the model data and measured
data, respectively.

Basing upon the CC-data model in Eq. (8) and discrete image and data arrays
devised above, we devise a discrete-to-discrete (DD)-data model, which relates
$\mathbf{f}$ and $\mathbf{g}$ defined above, as

$\displaystyle\mathbf{g}=\mathcal{H}\mathbf{f}\quad\quad{\rm
with}\quad\quad\mathcal{H}=\mathcal{C}\mathcal{R}\mathcal{D}_{B},$ (1)

where system matrix $\mathcal{H}$ is of size $J\times I$, matrix $\mathcal{R}$
of size $J\times I$ denotes the 4D-Radon transform in a discrete form in which
element $r_{ji}$ is chosen to be the intersection hyperarea of hyperplane $j$
with voxel $i$ in the 4D-spectral-spatial space; $\mathcal{C}$ a diagonal
matrix of size $J\times J$ in which diagonal element $c_{jj}={\rm
cos}\gamma_{j}$; $\gamma_{j}$ the angle between the $B$-axis and vector
$\hat{\alpha}_{j}$; $\hat{\alpha}_{j}$ the orientation of a hyperplane at the
$j$th projection; and matrix $\mathcal{D}_{B}$ of size $I\times I$ represents
two-point difference along spectral dimension $B$. We note that
$|\gamma_{j}|\leq\gamma_{\tau}$ and that $\gamma_{\tau}<90^{\circ}$ because it
is physically impossible to collect data at $\gamma=90^{\circ}$ in CW EPRI.

Figure 1: (a) a 3D-spatial coordinate system in which the orientation of
vector $\vec{r}_{s}$, pointing to a spatial point $s$, is specified by angles
$\theta$ and $\phi$; (b) an effective 2D space formed by two orthogonal
dimensions $\vec{r}_{s}$ and $B$ (i.e., representing 4D space $\\{x,y,z,B\\})$
in which 4D unit vector $\hat{\alpha}$ has angle $\gamma$ relative to the
$B$-axis; (c) and (d): uniform sampling grids over $\gamma$ in FAR and LAR
scans with $\gamma_{\tau}=90^{\circ}$ and $60^{\circ}$; and ESA-sampling grids
over the entire hemisphere (e) specified by $\theta_{\tau}=90^{\circ}$ and
$\phi_{\tau}=90^{\circ}$, and over the partial hemisphere (f) specified by
$\theta_{\tau}=60^{\circ}$ and $\phi_{\tau}=45^{\circ}$, in the 3D-spatial
space $\\{x,y,z\\}$ for the FAR scan and a LAR scan described in Sec. 2.2.

###  2.2 Equal-solid-angle (ESA)-sampling scheme

The explicit form of the DD-data model, i.e., $\mathcal{H}$, depends upon not
only ranges $\xi_{\tau}$, $\gamma_{\tau}$, $\theta_{\tau}$, and $\phi_{\tau}$
but also the specific sampling schemes (i.e., the grids) over the ranges. In
Fig. 1a, we use 3D vector $\vec{r}_{s}$ to denote a point in 3D-spatial space
$\\{x,y,z\\}$. For a FAR or LAR scan in CW-ZM EPRI, system matrix
$\mathcal{H}$ in the DD-data model in Eq. (1) depends upon how the 4D-data
space is sampled. In the work, we divide its $\xi$-axis into a total of
$J_{\xi}$ samples uniformly distributed over $\xi\in[-\xi_{\tau},\xi_{\tau}]$
separated by constant interval $\Delta_{\xi}={2\xi_{\tau}}/J_{\xi}$, as shown
in Fig. 1b; and we then consider below the equal-solid-angle (ESA)-sampling
scheme that has been used for data acquisition in CW-ZM EPRI [7].

In an ESA-sampling scheme, $J_{\gamma}$ samples are evenly distributed over
$[-\gamma_{\tau},\gamma_{\tau}]$ with constant angular interval
$\Delta_{\gamma}$, as shown in Figs. 1c and 1d, whereas $J_{\theta}$ samples
evenly distributed over $[-\theta_{\tau},\theta_{\tau}]$ with constant angular
interval $\Delta_{\theta}$ and $J_{\phi}(\theta_{n})$ samples evenly
distributed over $\phi$ with angular interval $\Delta_{\phi}(\theta_{n})={\rm
min}(|\frac{\Delta_{\gamma}}{{\rm sin}(\theta_{n})}|,\pi)$, where function
${\rm min}(\cdot)$ returns the smallest value, yielding the respective ESA-
sampling grids shown. For example, for the FAR scan,
$\gamma_{\tau}=90^{\circ}$, $\theta_{\tau}=90^{\circ}$, and
$\phi_{\tau}=90^{\circ}$, and its ESA sampling grids are displayed in Figs. 1c
and 1e, whereas for a LAR scan specified by $\gamma_{\tau}=60^{\circ}$,
$\theta_{\tau}=60^{\circ}$, and $\phi_{\tau}=45^{\circ}$, its ESA sampling
grids are depicted in Figs. 1d and 1f.

###  2.3 Constrained optimization program

Image reconstruction in CW-ZM EPRI is tantamount to invert $\mathcal{H}$ for
obtaining SS image $\mathbf{f}$ from measured data
$\mathbf{g}^{[\mathcal{M}]}$ in a scan. In an attempt to invert the DD-data
model in Eq. (1), we formulate its inversion, i.e., the reconstruction problem
of $\mathbf{f}$, as a convex, constrained optimization program [25, 26, 27]
given by

$\displaystyle\mathbf{f}^{\star}\\!$ $\displaystyle=$
$\displaystyle\\!\underset{\mathbf{f}}{\mathsf{argmin}}\\!\left\\{\frac{1}{2}\\!\parallel\mathcal{C}\mathcal{R}\mathcal{D}_{B}{\mathbf{f}}-\mathbf{g}^{[\mathcal{M}]}\parallel_{2}^{2}\right\\}$
$\displaystyle{\rm
s.t.}\quad||\mathcal{D}_{x}{\mathbf{f}}||_{1}\\!\leq\\!t_{x},\,\,||\mathcal{D}_{y}{\mathbf{f}}||_{1}\\!\leq\\!t_{y},\,\,||\mathcal{D}_{z}{\mathbf{f}}||_{1}\\!\leq\\!t_{z},\,\,||\mathcal{D}_{B}{\mathbf{f}}||_{1}\\!\leq\\!t_{B},\,\,f_{i}\\!\geq\\!0,$

where matrices $\mathcal{D}_{x}$, $\mathcal{D}_{y}$, $\mathcal{D}_{z}$, and
$\mathcal{D}_{B}$ of sizes $I\times I$ denote two-point differences along
spatial dimensions $x$, $y$, and $z$ and spectral dimension $B$, respectively;
$||\cdot||_{1}$ depicts the $\ell_{1}$-norm of the input vector;
$||\mathcal{D}_{x}\mathbf{f}||_{1}$, $||\mathcal{D}_{y}\mathbf{f}||_{1}$,
$||\mathcal{D}_{z}\mathbf{f}||_{1}$, and $||\mathcal{D}_{B}\mathbf{f}||_{1}$
are referred to, respectively, as the DTVs of SS image $\mathbf{f}$ along
spatial dimensions $x$, $y$, and $z$ and spectral dimension $B$; and
constraint parameters $t_{x}$, $t_{y}$, $t_{z}$, and $t_{B}$ are the upper
bounds on the respective DTVs of SS image $\mathbf{f}$. Because a physical SS
image is non-negative, an additional constraint, i.e., $f_{i}\geq 0$, is
included in the optimization program in Eq. (2.3).

###  2.4 DTV algorithm for reconstruction of SS image $\mathbf{f}$

Basing upon the primal-dual approach to solving convex optimization problems
[28, 29, 25], we derive a primal-dual instance algorithm, referred to as the
DTV algorithm, to reconstruct SS image $\mathbf{f}$ through solving the
optimization program in Eq. (2.3), and show the pseudo-codes of the DTV
algorithm in Algorithm 1.

Algorithm 1 Pseudocode of the DTV algorithm reconstructing SS image
$\mathbf{f}$ through solving Eq. (2.3)

1:INPUT: $g^{[\mathcal{M}]}$, ${\mathcal{H}}$, ${t}_{x}$, ${t}_{y}$,
${t}_{z}$, ${t}_{B}$, and $\rho\approx 10^{-2}$

2:$\nu_{1}\leftarrow\sqrt{e_{L}}/||\mathcal{D}_{x}||_{2}$,
$\nu_{2}\leftarrow\sqrt{e_{L}}/||\mathcal{D}_{y}||_{2}$,
$\nu_{3}\leftarrow\sqrt{e_{L}}/||\mathcal{D}_{z}||_{2}$,
$\nu_{4}\leftarrow\sqrt{e_{L}}/||\mathcal{D}_{B}||_{2}$,
$\mu\leftarrow\sqrt{e_{L}}||_{2}/||\mathcal{I}||_{2}$, $\Sigma$, $\mathcal{T}$

3:$n\leftarrow 0$

4:INITIALIZE: ${\mathbf{f}}^{(0)}$, $\mathbf{w}^{(0)}$, $\mathbf{p}^{(0)}$,
$\mathbf{q}^{(0)}$, $\mathbf{r}^{(0)}$, $\mathbf{b}^{(0)}$, and
$\mathbf{t}^{(0)}$ to zero

5:$\bar{\mathbf{f}}^{(0)}\leftarrow\mathbf{f}^{(0)}$

6:repeat

7:
$\mathbf{w}^{(n+1)}=(\mathbf{w}^{(n)}+\Sigma(\mathcal{H}\bar{\mathbf{f}}^{(n)}-\mathbf{g}^{[\mathcal{M}]}))/(\mathcal{I}+\Sigma)$

8:
$\bar{\mathbf{p}}^{(n)}=\mathbf{p}^{(n)}+\Sigma\nu_{1}\mathcal{D}_{x}\bar{\mathbf{f}}^{(n)}$

9:
$\bar{\mathbf{q}}^{(n)}=\mathbf{q}^{(n)}+\Sigma\nu_{2}\mathcal{D}_{y}\bar{\mathbf{f}}^{(n)}$

10: $\bar{\mathbf{r}}^{(n)}\hskip
2.84544pt=\mathbf{r}^{(n)}+\Sigma\nu_{3}\mathcal{D}_{z}\bar{\mathbf{f}}^{(n)}$

11:
$\bar{\mathbf{b}}^{(n)}=\mathbf{b}^{(n)}+\Sigma\nu_{4}\mathcal{D}_{B}\bar{\mathbf{f}}^{(n)}$

12:
$\mathbf{p}^{(n+1)}=\bar{\mathbf{p}}^{(n)}-\Sigma\frac{\bar{\mathbf{p}}^{(n)}}{|\bar{\mathbf{p}}^{(n)}|}\ell_{1}{\rm
ball}_{\nu_{1}{t}_{x}}(\frac{\bar{\mathbf{p}}^{(n)}}{\Sigma})$

13:
$\mathbf{q}^{(n+1)}=\bar{\mathbf{q}}^{(n)}-\Sigma\frac{\bar{\mathbf{q}}^{(n)}}{|\bar{\mathbf{q}}^{(n)}|}\ell_{1}{\rm
ball}_{\nu_{2}{t}_{y}}(\frac{\bar{\mathbf{q}}^{(n)}}{\Sigma})$

14:
$\mathbf{r}^{(n+1)}=\bar{\mathbf{r}}^{(n)}-\Sigma\frac{\bar{\mathbf{r}}^{(n)}}{|\bar{\mathbf{r}}^{(n)}|}\ell_{1}{\rm
ball}_{\nu_{3}{t}_{z}}(\frac{\bar{\mathbf{r}}^{(n)}}{\Sigma})$

15:
$\mathbf{b}^{(n+1)}=\bar{\mathbf{b}}^{(n)}-\Sigma\frac{\bar{\mathbf{b}}^{(n)}}{|\bar{\mathbf{b}}^{(n)}|}\ell_{1}{\rm
ball}_{\nu_{4}{t}_{B}}(\frac{\bar{\mathbf{b}}^{(n)}}{\Sigma})$

16: $\mathbf{t}^{(n+1)}={\rm
neg}({\mathbf{t}^{(n)}+\Sigma\mu\bar{\mathbf{f}}^{(n)}})$

17:
$\mathbf{f}^{(n+1)}=\mathbf{f}^{(n)}-\mathcal{T}(\mathcal{H}^{\top}\mathbf{w}^{(n+1)}+\nu_{1}\mathcal{D}_{x}^{\top}{\mathbf{p}}^{(n+1)}+\nu_{2}\mathcal{D}_{y}^{\top}{\mathbf{q}}^{(n+1)}$

18:
$+\,\,\nu_{3}\mathcal{D}_{z}^{\top}{\mathbf{r}}^{(n+1)}+\nu_{4}\mathcal{D}_{B}^{\top}{\mathbf{b}}^{(n+1)}+\mu\mathbf{t}^{(n+1)})$

19: $\bar{\mathbf{f}}^{(n+1)}=2\mathbf{f}^{(n+1)}-\mathbf{f}^{(n)}$

20: $n\leftarrow n+1$

21:until the necessary convergence conditions in Eq. (10)are satisfied

22:OUTPUT: final SS image $\mathbf{f}^{\star}=\mathbf{f}^{(n)}$

In Algorithm 1, input parameter $\rho\neq 0$ is chosen for empirically
optimizing the convergence rate and thus the path of the DTV algorithm leading
to the feasible solution set specified by the optimization program in Eq.
(2.3); ${\mathcal{H}}=\mathcal{C}\mathcal{R}\mathcal{D}_{B}$; $\mathcal{I}$ is
an identity matrix of size $I\times I$; $\mathcal{T}$ and $\Sigma$ are non-
diagonal and diagonal matrices, respectively, of size $I\times I$:

$\displaystyle\mathcal{T}$ $\displaystyle=$
$\displaystyle\frac{1}{\rho}\left(\frac{\mathcal{I}}{e_{L}}+\sum_{l=1}^{L-1}\mathbf{u}_{l}(\frac{1}{e_{l}}-\frac{1}{e_{L}})\mathbf{u}_{l}^{\top}\right)\quad{\rm
and}\quad\Sigma=\rho\varsigma\mathcal{I},$ (3)

$e_{l}$ and $\mathbf{u}_{l}$ are the $l$th eigenvalue and eigenvector of
matrix $[{\mathcal{H}}^{\top}{\mathcal{H}}]$, $L$ is the number of
eigenvectors in the expansion and $L=10$ is used in the work;
$\varsigma={||\mathcal{T}({\mathcal{K}}^{\top}{\mathcal{K}})||_{2}}^{-1}$;
$\mathcal{K}^{\top}=({\mathcal{H}}^{\top},\nu_{1}\mathcal{D}_{x}^{\top},\nu_{2}\mathcal{D}_{y}^{\top},\nu_{3}\mathcal{D}_{z}^{\top},\nu_{4}\mathcal{D}_{B}^{\top}$,
$\mu\mathcal{I}$); superscript “$\top$” depicts the transpose; $||\cdot||_{2}$
indicates the largest singular value of a matrix; vector $\mathbf{w}^{(n)}$ is
of size $J$, whereas vectors $\bar{\mathbf{p}}^{(n)}$,
$\bar{\mathbf{q}}^{(n)}$, $\bar{\mathbf{r}}^{(n)}$, $\bar{\mathbf{b}}^{(n)}$,
$\mathbf{p}^{(n)}$, $\mathbf{q}^{(n)}$, $\mathbf{r}^{(n)}$,
$\mathbf{b}^{(n)}$, and $\mathbf{t}^{(n)}$ are of size $I$; operator ${\rm
neg}(\cdot)$ enforces the non-positivity constraint; operator $\ell_{1}{\rm
ball}_{\omega}(\cdot)$ projects a vector onto the $\ell_{1}$-ball of scale
$\omega$; $|\bar{\mathbf{q}}|$ depicts a vector of size $I$ with entry $i$
given by $(|\bar{\mathbf{q}}^{(n)}|)_{i}=|\bar{q}^{(n)}_{i}|$; and
$\bar{q}^{(n)}_{i}$ indicates the $i$th entry of vector
$\bar{\mathbf{q}}^{(n)}$. The final reconstruction of SS image
$\mathbf{f}^{\star}$ is obtained when multiple necessary convergence
conditions in Eq. (10) of the DTV algorithm [25] are reached numerically, as
shown in B.

###  2.5 Reconstruction evaluation and quantitative metrics

We conduct visual inspection and quantitative analysis to evaluate the
performance of the DTV algorithm in reconstruction of SS images in simulated-
and real-data studies. In visual inspection, we examine the SS-image
appearance, spectral profiles, and possible LAR artifacts as compared to their
corresponding reference images and spectral profiles, whereas in quantitative
analysis, we compute the quantitative metrics described below from SS images
reconstructed and compare them with those computed from their reference
images. In the simulated-data study, the truth SS image, i.e., the digital
phantom, is available and thus is used as the reference image; whereas in the
real-data study, we use the SS image obtained from FAR data as the reference
image.

The quantitative metrics used in the work include the normalized root-mean-
square-error (nRMSE) and Pearson-correlation-coefficient (PCC) [30, 31, 32,
25], defined as

$\displaystyle{\rm nRMSE}(\mathbf{f}^{\star})$ $\displaystyle=$
$\displaystyle||\mathbf{f}^{\star}-\mathbf{f}^{\rm[ref]}||_{2}/||\mathbf{f}^{\rm[ref]}||_{2},$
(4) $\displaystyle{\rm PCC}(\mathbf{f}^{\star})$ $\displaystyle=$
$\displaystyle\frac{|\rm
Cov(\mathbf{f}^{\star},\mathbf{f}^{\rm[ref]})|}{\sigma(\mathbf{f}^{\star})\,\sigma(\mathbf{f}^{\rm[ref]})},$
(5)

where ${\rm Cov}(\mathbf{f}^{\star},\mathbf{f}^{\rm[ref]})$ denotes the
covariance between final SS image $\mathbf{f}^{\star}$ reconstructed and
reference SS image $\mathbf{f}^{\rm[ref]}$; $\sigma^{2}(\mathbf{f})={\rm
Cov}(\mathbf{f},\mathbf{f})$ the variances of $\mathbf{f}$. Metric nRMSE
measures the accuracy of a reconstruction relative to the reference, whereas
metric PCC yields information about the visual correlation between
reconstructed and reference images. We note that ${\rm
nRMSE}(\mathbf{f}^{\star})\rightarrow 0$ and ${\rm
PCC}(\mathbf{f}^{\star})\rightarrow 1$, if
$\mathbf{f}^{\star}\rightarrow\mathbf{f}^{\rm[ref]}$.

Furthermore, spectral profile $V_{k}(B)$ of contrast material type $k$ at a
spatial location in EPRI can be modeled as the convolution of a Gaussian and a
Lorentzian, which is the Voigt function [33, 34], i.e.,

$\displaystyle V_{k}(B)=A_{k}N_{\sigma}(B)\circledast L_{\tau_{k}}(B),$ (6)

where $N_{\sigma}(B)$ denotes a Gaussian with mean $0$ and standard deviation
$\sigma$, and $L_{\tau_{k}}(B)=\frac{\tau_{k}/2}{B^{2}+(\tau_{k}/2)^{2}}$ the
Lorentzian of width $\tau_{k}$, where $k=1,2,...,K$, indicating type $k$ of
contrast material, $K$ is the total number of types of contrast materials
within the subject imaged, $\circledast$ denotes convolution, and parameter
$A_{k}$ is determined such that $\int{\rm d}B\,V_{k}(B)=1$. Lorentzian width
$\tau_{k}$ can be related to oxygen concentration $\eta_{k}$ within the
subject imaged through a linear relationship specified by parameters $l_{0}$
and $\beta$ as

$\displaystyle\tau_{k}=l_{0}+\beta\times\eta_{k}.$ (7)

We note that parameters $\sigma$, $l_{0}$, and $\beta$ are independent of
oxygen concentration given a contrast material type. In our quantitative
evaluation study, we also estimate width $\tau_{k}$ and oxygen concentration
$\eta_{k}$ from final SS images reconstructed from FAR and LAR data.

##  3 Results

While the DTV algorithm in Sec. 2.4 can be shown [29, 25, 26] theoretically to
solve the convex optimization program in Eq. (2.3), it remains necessary to
verify the correctness of its computer implementation through a numerical
study in which measured and model data are consistent, i.e.,
$\mathbf{g}^{[\mathcal{M}]}=\mathbf{g}$. This can be realized in a simulated-
data study in which $\mathbf{g}^{[\mathcal{M}]}$ is generated by using a
digital phantom in Eq. (1). We first perform a verification study and include
the study results in B.2 to avoid distraction from the presentation flow of
the main study results of the work below.

With the DTV algorithm verified, we first carry out a simulated-data study to
evaluate the performance of the DTV reconstructions directly from LAR data.
Specifically, we use system matrix $\mathcal{H}$ in the DD-data model of Eq.
(1) to generate LAR data from the digital phantom and also use $\mathcal{H}$
in Eq. (2.3) for image reconstruction. The purpose of the simulated-data study
is to evaluate reconstruction accuracy directly from LAR data relative to the
truth SS image for exploring empirically the upper bound of the DTV
algorithm’s performance in reconstruction of SS images from LAR data. In the
real-data study, we apply the DTV algorithm to reconstructing SS images
directly from FAR and LAR data collected from a physical phantom, and evaluate
the DTV algorithm’s performance and robustness with real data containing noise
and other physical factors that are inconsistent with (i.e., not included in)
the DD-data model in Eq. (1).

Figure 2: The 4D digital phantom (i.e., the truth SS image), which is used in
the simulated-data study, within slices specified by $B=0$ G and $y=0$ mm (a)
and $B=0$ G and $x=-1.875$ mm (b), respectively, along with the color bars;
and (c) the truth spectral profiles at spatial centers of the tubes labeled as
$k=1$, 2, and 3 in (a). Display window of [0., 0.13] arbitrary unit (AU) is
used in (a) and (b). The image profiles over the white vertical line overlaid
in (a) are plotted in Fig. 6 below.

###  3.1 Study materials and conditions

####  3.1.1 Digital and physical phantoms

In the simulated-data study, we design a digital phantom, i.e., the truth SS
image, on a discrete array of $32\times 32\times 32\times 32$ identical voxels
each of which is of size $0.31\times 0.31\times 0.31\times 0.16$
mm${}^{3}\,$G, covering a 4D support of size
$10\\!\times\\!10\\!\times\\!10\\!\times\\!0.5$ mm3G in 4D-SS space
$\\{x,y,z,B\\}$. As shown in Fig. 2, the digital phantom comprises three
ellipsoidal tubes embedded in the 3D-spatial space of zero background. The
three tubes contain contrast materials with Voigt functions, specified by Eq.
(6), of different widths $\tau_{k}$. We select truth $\sigma=44.2$ mG ,
$l_{0}=10$ mG, and $\beta=0.5$, as well as truth $\tau_{k}=$ 10, 22, and 35 mG
for tubes $k=$1, 2, and 3. Using truth $\tau_{k}$s, $l_{0}$, and $\beta$ in
Eq. (7), we obtain truth $\eta_{k}$=0, 24, and 50 torr, respectively, for the
three tubes. In Figs. 2a $\&$ 2b, we display the 4D phantom (i.e., the truth
spectral-spatial image) within slices specified, respectively, by $B=0$ G &
$y=0$ mm and by $B=0$ G & $x=-1.875$ mm. Additionally, we plot its truth
spectral profile (i.e., the truth Voigt function) in Fig. 2c, at the center
point within each of the three tubes in Fig. 2a.

Figure 3: The physical phantom used in the real-data study consists of three
borosilicate glass cylinders with diameters of 9.5 mm, 5 mm, and 5 mm, which
are indexed with #1, #2, and #3, respectively, in (a). The physical phantom is
placed in a cylindric-shaped resonator of a diameter of 19 mm of our EPRI
scanner, and the effective scan region along $y$-axis is indicated between two
black lines in (b). A US one-dime coin is placed aside for reference.

In the real-data study, we have fabricated a physical phantom consisting of
three borosilicate glass tubes of cylindrical shape with diameters of 5.0 mm,
5.0 mm, and 9.5 mm, respectively, and identical length of 45 mm, as shown in
Fig. 3. The bottom glass tube of 9.5 mm diameter is sealed following being
filled with 1 mM deoxygenated OX071 trityl solution in PBS, while the top two
glass tubes of 5mm diameter are filled with 1 mM OX071 trityl solution in PBS
with the addition of different amounts of Cr Maltolate salt to simulate
different degrees of oxygen linewidth broadening. Independent of imaging
scans, we have carried out an experiment in which we measured, in addition to
$\sigma=31.8$ mG, width $\tau^{(\rm e)}_{k}$ and oxygen concentration
$\eta^{(\rm e)}_{k}$ within tube $k$, where $k=1,2,3$ and superscript “e”
indicates that they are experimental measurements. Specifically, we obtain the
width and oxygen-concentration pairs as ($\tau^{(\rm e)}_{1}=57.5$ mG,
$\eta^{(\rm e)}_{1}=0$ torr), ($\tau^{(\rm e)}_{2}=93.5$ mG, $\eta^{(\rm
e)}_{2}=62.24$ torr), and ($\tau^{(\rm e)}_{3}=135$ mG, $\eta^{(\rm
e)}_{3}=142.6$ torr), respectively, for tubes $k=$1, 2, and 3. Through fitting
the three pairs into Eq. (7), we subsequently obtain parameters $l_{0}=57.5$
mG and $\beta=0.543$ mG/torr. In our real-data study, we use $\sigma$,
$l_{0}$, and $\beta$ determined experimentally in our estimation of $\tau_{k}$
and $\eta_{k}$ from 4D images reconstructed from FAR or LAR data, as presented
below.

####  3.1.2 FAR and LAR scans

Table 1: FAR and LAR scans considered in the simulated-data study. $\gamma_{\tau}$ ($\theta_{\tau}$, $\phi_{\tau}$) | (85.5∘, 90∘) | (85.5∘, 45∘) | (67.5∘, 45∘)
---|---|---|---
85.5∘ | FAR | LAR1 | LAR2
67.5∘ | LAR3 | LAR4 | LAR5
58.5∘ | LAR6 | LAR7 | LAR8

In the simulated-data study, for the FAR scan specified by $\xi_{\tau}=250$
mG, $\gamma_{\tau}=90^{\circ}$, $\theta_{\tau}=90^{\circ}$, and
$\phi_{\tau}=90^{\circ}$, its ESA-sampling grids are determined by use of
intervals $\Delta_{\xi}=7.8$ mG, $\Delta_{\gamma}=9^{\circ}$,
$\Delta_{\theta}=9^{\circ}$, and $\Delta_{\phi}(\theta_{n})={\rm
min}(|\frac{\Delta_{\gamma}}{{\rm sin}(\theta_{n})}|,\pi)$, where $\theta_{n}$
denotes the $n$th sample on $\theta$, as discussed in Sec. 2.2. As discussed
above, it is physically impossible to acquire data near the vicinity of
$\gamma_{\tau}=90^{\circ}$. Therefore, an effective FAR scan is specified in
Table 1 that has $\gamma_{\tau}=85.5^{\circ}$. Using the digital phantom of
Fig. 2 in Eq. (1), we generate FAR data on the ESA-sampling grids. For a LAR
scan specified in Table 1, we obtain the LAR data by choosing a subset of the
FAR data in the portion of the FAR scan corresponding to the angular ranges of
the LAR scan.

Table 2: FAR and LAR scans considered in the real-data study. $\gamma_{\tau}$ ($\theta_{\tau}$, $\phi_{\tau}$) | (84.4∘, 90∘) | (84.4∘, 45∘) | (73.1∘, 45∘)
---|---|---|---
84.4∘ | FAR | LAR1 | LAR2
73.1∘ | LAR3 | LAR4 | LAR5
61.9∘ | LAR6 | LAR7 | LAR8

In the real-data study, we utilized a spectroscopic imager [35, 36] to scan
the phantom operating at radio frequency (RF) of 250 MHz with a non-saturated
level of RF power of 40 $\mu$W. During the scan, we applied a 5 kHz Zeeman
modulation with an amplitude of 60 mG. The strength of the magnetic field
gradient ranges from -350 mG/mm to 350 mG/mm. A total of 2624 sweeps were
obtained for the FAR scan, with each sweep acquired at 256 field points along
the $\xi$-axis and a scan time of 3 ms per point.

For the FAR scan specified by $\xi_{\tau}=724$ mG, $\gamma_{\tau}=90^{\circ}$,
$\theta_{\tau}=90^{\circ}$, and $\phi_{\tau}=90^{\circ}$, its ESA-sampling
grids are determined by use of $\Delta_{\xi}=5.66$ mG,
$\Delta_{\gamma}=11.25^{\circ}$, $\Delta_{\theta}=11.25^{\circ}$, and
$\Delta_{\phi}(\theta_{n})={\rm min}(|\frac{\Delta_{\gamma}}{{\rm
sin}(\theta_{n})}|,\pi)$, where $\theta_{n}$ denotes the $n$th sample on
$\theta$, as discussed in Sec. 2.2. For the same reason discussed in the
simulated-data study, an effective FAR scan is specified in Table 2 that has
$\gamma_{\tau}=84.4^{\circ}$. Using our EPRI system, we collect FAR data from
the physical phantom in Fig. 3. For a LAR scan specified in Table 2, we obtain
the LAR data by choosing a subset of the FAR data in the portion of the FAR
corresponding to the angular ranges of the LAR scan.

####  3.1.3 Reconstruction parameters

In the simulated-data study, SS images are reconstructed on image arrays
identical to that of the truth SS image. Values of constraint parameters
($t_{x}$, $t_{y}$, $t_{z}$, $t_{B}$) are chosen to be the DTVs computed,
respectively, along $x$, $y$, $z$, and $B$ dimensions of the truth SS image,
i.e., the known digital phantom in the simulated-data study. In an attempt to
streamline the presentation flow, we opt not to include results of a noisy
simulated-data study, because we also perform in the work a real-data study in
which data measured contain not only noise but also other physical factors
inconsistent with the DD-data model Eq. (1) upon which the DTV algorithm is
developed, and all of these physical factors may impact SS images
reconstructed.

In the real-data study, from each set of data collected in FAR and LAR scans
specified in Table 2, we perform reconstructions of SS images by using the DTV
algorithm for multiple sets of values of constraint parameters ($t_{x}$,
$t_{y}$, $t_{z}$, $t_{B}$) and then select the set that yields an SS image
with visually minimal artifacts. The SS images are reconstructed on image
arrays of $64\times 64\times 64\times 64$ identical voxels each of which is of
size  $0.66\times 0.66\times 0.66\times 22.6$ mm${}^{3}\,$mG, covering a 4D
support of size  $42\\!\times\\!42\\!\times\\!42\\!\times\\!1448$ mm3mG in
4D-SS space $\\{x,y,z,B\\}$.

###  3.2 SS images reconstructed in the simulated-data study

The DTV algorithm is applied to reconstructing SS images from data of the
digital phantom for the FAR and LAR scans summarized in Table 1 above.

####  3.2.1 Reconstruction results

Figure 4: SS images within slice $B=0$ G and $y=0$ mm reconstructed by use of
the DTV algorithm (row 1) and FBP algorithm (row 2) from simulated data of the
FAR scan (column 1) and LAR scans LAR1-LAR8 (columns 2-9) summarized in Table
1, respectively. Display window: [0., 0.13] AU. Figure 5: SS images within
slice $B=0$ G and $x=-1.875$ mm reconstructed by use of the DTV algorithm (row
1) and FBP algorithm (row 2) from simulated data of the FAR scan (column 1)
and LAR scans LAR1-LAR8 (columns 2-9) summarized in Table 1, respectively.
Display window: [0., 0.13] AU.

We first perform reconstruction from FAR data and display the SS images
reconstructed within the $x$-$y$ plane specified by $B=0$ G & $y=0$ mm and the
$y$-$z$ plane specified by $B=0$ G & $x=-1.875$ mm, respectively, in Figs. 4
and 5. The result verifies that the DTV algorithm and its computer
implementation can solve the optimization program in Eq. (2.3) and accurately
reconstruct SS images from FAR data.

We subsequently apply the DTV algorithm verified to reconstructing SS images
from data of LAR scans specified in Table 1. The reconstructed SS images
within the $x$-$y$ plane specified by $B=0$ G & $y=0$ mm and the $y$-$z$ plane
specified by $B=0$ G & $x=-1.875$ mm are shown, respectively, in Figs. 4 and
5. For reference, we also perform FBP reconstructions from FAR and LAR data
and display them in the bottom rows of Figs. 4 and 5.

It can be observed that the DTV algorithm reconstructs accurately SS images
from data of LAR scans depicted in Table 1, as the SS images appear to be
largely free of LAR artifacts and visually resemble the truth SS image. On the
other hand, the FBP images reconstructed from LAR data appear to suffer
considerably from artifacts that distort the anatomic structures in, and
introduce biases to, the SS images. Also, while the SS images reconstructed
from FAR data by use of the DTV and FBP algorithms appear somewhat comparable,
the former is visually free of minor artifacts observable in the latter.

For revealing reconstruction accuracy, we also plot in Figs. 6 and 7 the
spatial profiles (over the white line in Fig. 2a) and the spectral profile at
a spatial point of $x=0$ mm, $y=0$ mm, and $z=1.24$ mm (in tube 1) in SS
images reconstructed from FAR and LAR data. The spatial and spectral profiles
in the SS images reconstructed by use of the DTV algorithm appear to be in
good agreement with their counterparts in the truth SS image, whereas those in
the FBP reconstructions differ significantly from that in the truth SS image.

####  3.2.2 Quantitative evaluation

Using metrics nRMSE and PCC in Eqs. (4) and (5), we quantitatively evaluate
reconstruction accuracy of SS images relative to the reference image, i.e.,
the truth SS image. In Fig. 8, we show the nRMSEs and PCCs computed from the
SS images reconstructed by using the DTV algorithm directly from data of the
FAR and LAR scans specified in Table 1. It can be observed that, as the
angular ranges of LAR scans increase, the nRMSE decreases while the PCC
remains largely 1, corroborating the observation above that the SS images
reconstructed appear visually comparable to the truth SS image. We also
calculate metrics nRMSE and PCC from the FBP reconstructions, which are shown
in Fig. 8. We observe that the nRMSEs of the FBP reconstructions are
significantly higher than those of the DTV reconstructions, while the PCCs of
the FBP reconstructions are considerably lower than 1.

We also fit the Voigt function with the spectral profile at a selected
3D-spatial voxel within an SS image reconstructed by using the DTV or FBP
algorithm to estimate [33, 34] the Lorentzian line width $\tau_{k}$, which is
subsequently used to compute the oxygen concentration in the solution by use
of Eq. (7). Specifically, we select three 3D-spatial regions of interest
(ROIs), and each contains 90 3D-voxels, respectively, at the centers of tubes
1, 2, and 3. For each ROI, we first estimate $\tau_{k}$ at each of its voxels
by fitting the spectral profile to a Voigt curve and then calculate spatial
mean $\bar{\tau}_{k}$ and (spatial) standard deviation $\sigma_{k}^{\tau}$ by
averaging $\tau_{k}$s estimated over the voxels within the ROI, which are
shown in row 1 of Fig. 9. Using $\bar{\tau}_{k}$ in Eq. (7), along with truth
$l_{0}$ and $\beta$, we calculate spatial mean $\bar{\eta}_{k}$. Similarly,
(spatial) standard deviation $\sigma_{k}^{\eta}$ can readily be computed by
using $\sigma_{k}^{\tau}$ in Eq. (7). It can be observed that $\bar{\tau}_{k}$
and $\bar{\eta}_{k}$ estimated from reconstructions from LAR data by use of
the DTV algorithm agree well with the corresponding truth values. On the other
hand, $\bar{\tau}_{k}$ and $\bar{\eta}_{k}$ estimated from the FBP
reconstructions appear to deviate considerably from the truth values.

Figure 6: Spatial profiles (over the white line overlaid in the SS image in
Fig. 2a) of the reference (i.e., the truth) SS image (solid), and of the SS
images reconstructed by use of the DTV algorithm (dashed) and the FBP
algorithm (dotted) from simulated data of FAR scan (a) and LAR scans LAR1-LAR8
(b-i) depicted in Table 1. Similar results of spatial profiles are obtained
also over other lines in SS images reconstructed from data of LAR scans
depicted in Table 1, and thus are not shown. Figure 7: Spectral profiles (at a
spatial point in tube 1 specified by $x=0$ mm, $y=0$ mm, and $z=-1.24$ mm) of
the reference (i.e., the truth) SS image (solid), and of the SS images
reconstructed by use of the DTV algorithm (dashed) and the FBP algorithm
(dotted) from simulated data of FAR scan (a) and LAR scans LAR1-LAR8 (b-i)
depicted in Table 1. Similar results are obtained for spectra at other spatial
locations in SS images reconstructed from data of LAR scans depicted in Table
1 and thus are not shown. Figure 8: Quantitative metrics nRMSE (a) and PCC (b)
of the SS images reconstructed by use of the DTV (black, dashed) and FBP (red,
dotted) algorithms from simulated data of FAR and LAR scans specified in
Table. 1. A zoomed-in view of nRMSEs of DTV reconstructions are shown in (c).
Figure 9: Mean widths $\bar{\tau}_{k}$ (row 1) and mean oxygen-concentrations
$\bar{\eta}_{k}$ (row 2), along with their error bars (i.e.,
$\sigma_{k}^{\tau}$s and $\sigma_{k}^{\eta}$s), where $k=1,2,$ and 3 (columns
1-3) within the 3 tubes of the numerical phantom in Fig. 2 estimated,
respectively, from the DTV (black, ‘$\boldsymbol{\times}$’) and the FBP
(green, ‘$\blacksquare$’) reconstructions from simulated data of FAR and LAR
scans specified in Table. 1. Solid red line indicates the corresponding
reference (i.e., truth) value.

###  3.3 SS images reconstructed in the real-data study

####  3.3.1 Reconstruction results

Figure 10: SS images within slice $B=0$ G and  $y=-3.27$ mm reconstructed by
use of the DTV algorithm (row 1) and FBP algorithm (row 2) from real data of
the FAR scan (column 1) and LAR scans LAR1-LAR8 (columns 2-9) summarized in
Table 2, respectively. Display window: [0, 190] AU. The image profiles over
the white vertical line overlaid on the image (row 1 and column 1) are plotted
in Fig. 12 below. Figure 11: SS images within slice $B=0$ G and  $x=-3.93$ mm
reconstructed by use of the DTV algorithm (row 1) and FBP algorithm (row 2)
from real data of the FAR scan (column 1) and LAR scans LAR1-LAR8 (columns
2-9) summarized in Table 2, respectively. Display window: [0, 190] AU.

Using our EPRI system [7], we acquire data from the physical phantom, shown in
Fig. 3, for the FAR and LAR scans summarized in Table 2 above. We then apply
the DTV algorithm to reconstructing SS images from data acquired. For each of
the FAR and LAR scans in Table 2, we first carry out reconstructions with
multiple sets of constraint-parameter values chosen, and subsequently select
the optimal set that yields the SS image with visually minimal artifacts. In
Figs. 10 and 11, we display the SS images reconstructed within the $x$-$z$
plane specified by $B=0$ G & $y=-3.27$ mm and the $x$-$y$ plane specified by
$B=0$ G & $x=-3.93$ mm, respectively, each of which is obtained when its
corresponding optimal set of constraint-parameter values selected is achieved
numerically [25].

It can be observed that the DTV algorithm reconstructs SS images directly from
LAR data with minimal LAR artifacts and visually resemble the reference SS
image reconstructed from FAR data. Again, as the angular ranges of LAR scans
further decrease, the visual quality of DTV reconstructions of the SS image
appears to be reduced slightly especially for LAR scans with smallest angular
ranges in Table 2. However, the FBP images reconstructed from LAR data appear
to be with considerable artifacts that distort the anatomic structures in, and
introduce biases to, the SS images. Also, the SS image reconstructed from FAR
data by use of the DTV algorithm appears visually free of much of the
artifacts observable in the FBP image.

In addition, we plot in Figs. 12 and 13 the spatial profiles over the white
line in Fig. 10a and spectral profile for a spatial point specified by
$x=0.66$ mm, $y=-3.27$ mm, and $z=-3.27$ mm in SS images reconstructed from
FAR and LAR data. The spatial and spectral profiles in DTV reconstructions
appear to be in agreement with their counterparts in the reference SS image
reconstructed from FAR data by use of the DTV algorithm, whereas those in the
FBP reconstructions differ significantly from that in the reference SS image.

####  3.3.2 Quantitative Evaluation

Using metrics nRMSE and PCC in Eqs. (4) and (5), we also perform a
quantitative evaluation of reconstruction accuracy relative to the reference
image, i.e., the SS image reconstructed from FAR data by use of the DTV
algorithm. In Fig. 14, we show nRMSEs and PCCs computed from the SS images
reconstructed by using the DTV algorithm directly from data of the LAR scans
specified in Table 2. It can be observed that, as the angular ranges of a LAR
scan increase, the nRMSE decreases while PCC remains largely 1, corroborating
the observation above that the SS images reconstructed appear visually
comparable to the reference SS image. We also show in Fig. 14 nRMSE and PCC
calculated from the FBP reconstructions, and it can be observed that the nRMSE
and PCC are considerably higher and lower, respectively, than those of the DTV
reconstructions.

In the real-data study, we have $\sigma=31.8$ mG, $l_{0}=57.5$ mG, and
$\beta=0.543$ mG/torr, as discussed in Sec. 3.1.1 above. For a 4D spectral-
spatial image reconstructed from FAR or LAR data by using the DTV or FBP
algorithm, we choose three 3D-spatial ROIs containing 72, 48, and 48
3D-voxels, respectively, in glass tubes $k=1$, 2, and 3. For each ROI, we
first estimate $\tau_{k}$ at each of its voxels by fitting the spectral
profile to the Voigt function with $\sigma=31.8$ mG; and we then compute
spatial mean $\bar{\tau}_{k}$, and corresponding (spatial) standard deviation
$\sigma_{k}^{\tau}$, from $\tau_{k}$s estimated within the ROI. In particular,
for $\bar{\tau}_{k}$ of images reconstructed from FAR data by using each of
the DTV and FBP algorithms, we obtain a calibration coefficient by matching
$\bar{\tau}_{k}$ with the corresponding $\tau^{(\rm e)}_{k}$ determined in
Sec. 3.1.1 above. Subsequently, applying the calibration coefficient
determined for the DTV or FBP algorithm to its corresponding reconstructions
from LAR data, we obtain $\bar{\tau}_{k}$s calibrated for the LAR scans. Using
each of these $\bar{\tau}_{k}$s in Eq. (7), along with $l_{0}=57.5$ mG and
$\beta=0.543$ mG/torr, we then obtain $\bar{\eta}_{k}$ for FAR and LAR scans.
Similarly, we can readily obtain $\sigma_{k}^{\eta}$ by using its
corresponding $\sigma_{k}^{\tau}$ in Eq. (7).

In Fig. 15, we display $\bar{\tau}_{k}$s and $\bar{\eta}_{k}$s, along with
their standard deviations, estimated from DTV and FBP reconstructions of the
SS images from FAR and LAR data. It can be observed that $\bar{\tau}_{k}$s and
$\bar{\eta}_{k}$s from the DTV reconstructions for LAR scans agree reasonably
well with the corresponding $\tau^{(\rm e)}_{k}$ and $\eta^{(\rm e)}_{k}$.
Conversely, $\bar{\tau}_{k}$s and $\bar{\eta}_{k}$s obtained from the FBP
reconstructions appear to have significant biases relative to $\tau^{(\rm
e)}_{k}$ and $\eta^{(\rm e)}_{k}$ for LAR scans, especially for scans
LAR3-LAR8.

Figure 12: Spatial profiles (over the white line overlaid in the SS image in
the top left panel of Fig. 10) of the reference SS image (solid) reconstructed
from the FAR data by use of the DTV algorithm, and of the SS images
reconstructed by use of the DTV algorithm (dashed) and the FBP algorithm
(dotted) from real data of FAR scan (a) and LAR scans LAR1-LAR8 (b-i) depicted
in Table 2. Similar results are obtained for spatial profiles over other lines
in SS images reconstructed from data of the FAR and LAR scans depicted in
Table 2 and thus are not shown. Figure 13: Spectral profiles (at a spatial
point in glass tube 1 specified by  $x=0.66$ mm, $y=-3.27$ mm, and $z=-3.27$
mm) of the reference SS image (solid) reconstructed from the FAR data by use
of the DTV algorithm, and of the SS images reconstructed by use of the DTV
algorithm (dashed) and the FBP algorithm (dotted) from real data of FAR scan
(a) and LAR scans LAR1-LAR8 (b-i) depicted in Table 2. Similar results are
obtained for spectra at other spatial locations in SS images reconstructed
from data of LAR scans depicted in Table 2 and thus are not shown. Figure 14:
Quantitative metrics nRMSE (a) and PCC (b) of the SS images reconstructed by
use of the DTV (black, dashed) and FBP (red, dotted) algorithms from real data
of FAR and LAR scans specified in Table. 2. Figure 15: Mean widths
$\bar{\tau}_{k}$ (row 1) and oxygen concentrations $\bar{\eta}_{k}$ (row 2),
along with their error bars (i.e., $\sigma_{k}^{\tau}$ and
$\sigma_{k}^{\eta}$), where $k=1,2,$ and 3 (columns 1-3) within the 3 glass
tubes of the physical phantom in Fig. 3 estimated, respectively, from the DTV
(black, ‘$\boldsymbol{\times}$’) and the FBP (green, ‘$\blacksquare$’)
reconstructions from real data of FAR and LAR scans specified in Table. 2.
Solid red line indicates the corresponding reference value.

##  4 Discussion

In the work, we investigate and develop the DTV algorithm for reconstruction
of SS images directly from LAR data in CW-ZM EPRI. The work is motivated by
(1) the observation that, because the physical gradient strength of the
spatial encoding magnetic field applied can be only of finite value, data
cannot be acquired physically in the vicinity of $|\gamma|=90^{\circ}$, thus
resulting in an inherent LAR reconstruction problem; (2) the interest in
reconstruction of an SS image directly from data in CW-ZM EPRI; and (3) the
interest in the development of CW-ZM EPRI with data acquired only over LARs
for possibly minimizing the scan time. Based upon the DD-data model in Eq. (1)
for CW-ZM EPRI, we first formulate the reconstruction problem as a constrained
optimization program with unique constraints on the DTVs of an SS image and
then developed the DTV algorithm to reconstruct the SS image through solving
the optimization program.

We carry out simulated- and real-data studies to verify and evaluate, in terms
of visual inspection and quantitative analysis, the DTV algorithm in
reconstructions of SS images directly from data acquired in LAR scans of
varying angular ranges in CW-ZM EPRI. The study results reveal that the DTV
algorithm can accurately reconstruct SS images directly from data collected in
LAR scans of angular ranges considerably smaller than that of the FAR scan,
that it can stably reconstruct SS images directly from real LAR data collected
in CW-ZM EPRI, and that it can yield, from FAR data, SS images of quality
improved over those obtained with the FBP algorithm used in current CW-ZM
EPRI. The DTV constraints are demonstrated to be effective in suppression of
LAR and other artifacts that are observed otherwise in images reconstructed
with the FBP algorithm.

While the spectral shape described in Eq. (6) is considered in the studies, we
note that the optimization-based reconstruction approach and the associated
DTV algorithm proposed can readily accommodate spectral shapes of other forms
[37, 12, 38], including regular and distorted/broadened Lorentzians. In the
simulated-data study, we have chosen the DTVs of the truth SS image, i.e., the
digital phantom, computed as the values of the DTV-constraint parameters,
whereas in the real-data study, we have selected the constraint parameter
values that yield the SS images with visually minimal artifacts. However, in
an application with a specific goodness metric devised, the values of the
constraint parameters can be chosen empirically for maximizing the metric. In
particular, one may perform reconstructions with multiple sets of the
constraint-parameter values chosen, then compute the metric values from the
reconstructions, and finally select the set of constraint-parameter values
yielding the numerically highest value of the application metric. We also note
that the reconstruction approach and DTV algorithm developed can readily be
applied to sampling schemes other than the ESA-sampling scheme considered.

We have performed additional quantitative analysis with digital and physical
phantoms of anatomies different than those in Figs. 2 and 3 in simulated- and
real-data studies. Based upon the results obtained, observations can be made
similar to those discussed in Sec. 3. While the constrained optimization
program designed in Eq. (2.3) is demonstrated to reduce LAR artifacts in SS
images reconstructed from LAR data, it remains of interest to investigate
different designs of optimization programs for potentially further minimizing
artifacts for SS images reconstructed from data of LAR scans with angular
ranges even smaller than those considered. For example, the data-fidelity term
considered is of the $\ell_{2}$-norm form in the work, it is worthy to
investigate data-fidelity terms of additional forms [39] such as the Kullback-
Leibler-divergence form for possibly reducing artifacts in SS images
reconstructed from data collected in LAR scans of very small angular ranges.

It is worthy to note that the study design and DTV algorithm developed in the
work can naturally be extended to investigating SS image reconstruction from
FAR or LAR data acquired at sparse views [40, 41, 42, 43] in CW-ZM EPRI. It
remains a future task of interest to tailor, optimize, and evaluate image
reconstructions from sparse-view data by using the DTV algorithm developed in
CW-ZM EPRI.

##  5 Conclusion

We have investigated and developed the optimization-based DTV algorithm for
accurate reconstruction of SS images directly from FAR and LAR data in CW-ZM
EPRI. Results of simulated- and real-data studies reveal that from LAR data,
SS images may be reconstructed that are visually and quantitatively comparable
to those obtained from FAR data in CW-ZM EPRI. The DTV algorithm may be
exploited possibly for enabling CW-ZM EPRI with reduced imaging time and
artifacts by acquiring only LAR data instead of FAR data. The approach and DTV
algorithm can readily be tailored to reconstruct SS images directly from FAR
and LAR data in CW EPRI with scanning schemes other than the ZM scheme.

## Declaration of Competing Interest

Authors HH and BE declare US patents 8,664,955, 10,568,537, and 9,392,957 on
aspects of the pO2 imaging technology and memebership in a start up company
O2MThe o to market the pO2 imaging technology. Other authors declare that they
have no known competing financial interests or personal relationships that
could have appeared to influence the work reported in this paper.

## Data availability

The authors will provide the reconstruction code used in this study upon
request.

## Acknowledgments

This work was supported in part by NIH Grants R01EB026282, R01EB023968,
R21CA263660, P41EB02034, R01CA98575, and R01CA236385. The computation of the
work was performed in part on Computer Cluster Mel in the Department of
Radiology at The University of Chicago. The contents of this article are
solely the responsibility of the authors and do not necessarily represent the
official NIH views.  We would like to thank the anonymous reviewers for their
insightful and constructive comments.

##  Appendix A Continuous-to-continuous (CC)-data model in CW-ZM EPRI

We summarize the data model in a continuous form in CW-ZM EPRI, which allows
for an intuitive incorporation of underlying physics involved in CW-ZM EPRI.
While a data model in a continuous form is not directly applicable to data
acquired and an SS image reconstructed on discrete arrays in practical EPRI,
it provides a basis for devising a DD-data model that can adequately relate
data and an SS image on discrete arrays in practical CW-ZM EPRI.

We consider an SS image of continuous form in a 4D-image space
$\\{x^{\prime},y^{\prime},z^{\prime},B\\}$ formed by 3D-spatial dimensions
$x^{\prime}$, $y^{\prime}$, and $z^{\prime}$ and 1D-spectral coordinate $B$.
Assuming $L^{\prime}$ mm, and $B_{w}$ G, to be the support along each of the
3D-spatial dimensions, and along the spectral dimension, respectively, within
the 4D-image space, we form a 4D-physical support of size $L^{\prime}\times
L^{\prime}\times L^{\prime}\times B_{w}$ mm3 G, with its center coinciding
with the origin of the 4D-image space, and assume that the SS image of
continuous form is completely enclosed in the support. It is noted that the
units of the 3D-spatial and 1D-spectral dimensions are different. For
discussion convenience, we introduce a new 4D space as
$\\{x,y,z,B\\}=\\{cx^{\prime},cy^{\prime},cz^{\prime},B\\}$, where
$c=B_{w}/L^{\prime}$ G mm-1, in which all of the dimensions are of identical
unit G. Subsequently, we have a 4D-image support of size $L\times L\times
L\times B_{w}$ G4, where $L=cL^{\prime}$, which completely encloses the SS
image of continuous form. For example, assuming a spatial support of $2\times
2\times 2$ mm3 and spectral support of 2 G, we have $c=1$ G/mm, and if
$x^{\prime}=0.41$ mm, $x=1\,{\rm G/mm}\times 0.41\,{\rm mm}=0.41$ G. However,
for clarity, we remain to use unit “mm” for each of spatial dimensions $x$,
$y$, and $z$, but with the understanding that it is converted automatically to
unit “G” in the mathematical analysis and algorithm development.

We also consider 4D-data space $\\{\xi,\vec{\alpha}\\}$, formed by scalar
variable $\xi$ and 4D unit-vector $\hat{\alpha}=({\rm cos}\phi\,{\rm
sin}\theta\,{\rm sin}\gamma,\,{\rm sin}\phi\,{\rm sin}\theta\,{\rm
sin}\gamma,\,{\rm cos}\theta\,{\rm sin}\gamma,\,{\rm cos}\gamma)^{\top}$,
where $\gamma$ depicts the angle between the $B$-axis and vector
$\hat{\alpha}$ and takes a positive or negative value when the inner product
of $\hat{\rm e}_{z}$ and $\hat{\alpha}$ is positive or negative. We also use
$\hat{\alpha}_{p}$ to denote the projection of $\hat{\alpha}$ onto the
hyperplane specified by $B=0$ in the 4D-image space. Moreover, in 3D spatial
space $\\{x,y,z\\}$ with unit vectors $\hat{\rm e}_{x}$, $\hat{\rm e}_{y}$,
and $\hat{\rm e}_{z}$, we use $\theta$ to depict the angle between $z$-axis
and $\hat{\alpha}_{p}$, which takes a positive or negative value when the
inner product of $\hat{\rm e}_{x}$ and $\hat{\alpha}_{p}$ is positive or
negative, and also $\phi$ to depict the angle between $x$-axis and the
projection of $\hat{\alpha}_{p}$ onto the $x$-$y$ plane, which takes a
positive or negative value when the inner product of $\hat{\rm e}_{y}$ and
$\hat{\alpha}_{p}$ is positive or negative. Considering a 3D hyperplane in the
4D-image space, we use scalar variable $\xi$, and unit vector $\hat{\alpha}$,
to specify its distance to the origin of, and its orientation in, the 4D image
space.

We use $g(\xi,\hat{\alpha})$ to denote the data function in CW-ZM EPRI, which
can be expressed as [15]

$\displaystyle g(\xi,\hat{\alpha})$ $\displaystyle=$ $\displaystyle{\rm
cos}\gamma\,\int d{\vec{r}}\,\frac{\partial f(\vec{r})}{\partial
B}\,\delta(\xi-\vec{r}\cdot\hat{\alpha})$ (8) $\displaystyle=$
$\displaystyle{\rm cos}\gamma\,\mathbb{R}\\!\left[\frac{\partial
f(\vec{r})}{\partial B}\right],$

where $\vec{r}$ is a point in the 4D space, $\mathbb{R}$ denotes the 4D-Radon
transform of $\frac{\partial f(\vec{r})}{\partial B}$, i.e., the partial
derivative of 4D-SS image $f(\vec{r})$ of continuous form, over a hyperplane
specified by $\xi$ and $\hat{\alpha}$, $\xi\in[-\xi_{\tau},\xi_{\tau}]$ with
$2\xi_{\tau}\geq{\rm max}(L,B_{w})$, and ${\rm max}(\cdot)$ returns the
largest value. In the FAR scan, $\gamma\in[-\frac{\pi}{2},\frac{\pi}{2}]$,
$\theta\in[-\frac{\pi}{2},\frac{\pi}{2}]$, and
$\phi\in[-\frac{\pi}{2},\frac{\pi}{2}]$. Because $\xi$, $\hat{\alpha}$, and
$\vec{r}$ can continuously vary, the data model in Eq. (8) is thus referred to
as the continuous-to-continuous (CC)-data model.

In CW-ZM EPRI, the magnetic field gradient is expressed as
$\mathbf{G}=G\hat{\mathbf{G}}$, where $G$ represents the magnitude of the
vector with units of G/mm, and unit vector $\hat{\mathbf{G}}$ is dimensionless
specifying the direction of $\mathbf{G}$. Angle $\gamma$ is defined as
$|\gamma|={\rm arctan}(G/c)$, while angles $\theta$ and $\phi$ specify
$\hat{\mathbf{G}}$. Algorithms have been developed for obtaining 4D-PD-SS
image $\frac{\partial f(\vec{r})}{\partial B}$ through inverting the 4D-Radon
transform in Eq. (8) from its full knowledge in
$\gamma\in[-\frac{\pi}{2},\frac{\pi}{2}]$,
$\theta\in[-\frac{\pi}{2},\frac{\pi}{2}]$, and
$\phi\in[-\frac{\pi}{2},\frac{\pi}{2}]$. In practical CW-ZM EPRI, the
physically achievable maximum gradient strength $G_{\tau}$ must be finite, and
thus knowledge of $g(\xi,\hat{\alpha})$ can be measured only over
$\gamma\in[-\gamma_{\tau},\gamma_{\tau}]$, where $\gamma_{\tau}={\rm
arctan}(G_{\tau}/c)<90^{\circ}$, instead of over full angular range
$\gamma\in[-\frac{\pi}{2},\frac{\pi}{2}]$. Therefore, image reconstruction in
practical CW-ZM EPRI is inherently a reconstruction problem from LAR data. In
Sec. 2.1, we use the CC-data model in Eq. (8) to devise a DD-data model in Eq.
(1) adopted in the work.

##  Appendix B Convergence of the DTV algorithm

We devise convergence conditions for the DTV algorithm and verify that the DTV
algorithm can reconstruct accurately an SS image from its FAR data through
solving the optimization program in Eq. (2.3).

###  B.1 Convergence conditions of the DTV algorithm

We devise below dimensionless convergence metrics as

$\begin{split}\widetilde{D}^{(n)}_{\mathbf{g}}&=|\sqrt{D_{\mathbf{g}}(\mathbf{f}^{(n)})}|/||\mathbf{g}^{[\mathcal{M}]}||_{2}\\\
\widetilde{D}^{(n)}_{\rm TV_{\it
x}}&=|(||\mathcal{D}_{x}\mathbf{f}^{(n)}||_{1}-t_{x})|/t_{x}\\\
\widetilde{D}^{(n)}_{\rm TV_{\it
y}}&=|(||\mathcal{D}_{y}\mathbf{f}^{(n)}||_{1}-t_{y})|/t_{y}\\\
\widetilde{D}^{(n)}_{\rm TV_{\it
z}}&=|(||\mathcal{D}_{z}\mathbf{f}^{(n)}||_{1}-t_{z})|/t_{z}\\\
\widetilde{D}^{(n)}_{\rm TV_{\it
B}}&=|(||\mathcal{D}_{B}\mathbf{f}^{(n)}||_{1}-t_{B})|/t_{B},\\\ \end{split}$
(9)

where
${D}_{\mathbf{g}}({\mathbf{f}^{(n)}})=\frac{1}{2}\\!\parallel\mathcal{C}\mathcal{R}\mathcal{D}_{B}{\mathbf{f}^{(n)}}-\mathbf{g}^{[\mathcal{M}]}\parallel_{2}^{2}$,
$\mathbf{f}^{(n)}$ denotes the SS image reconstructed at iteration $n$, and
$\mathbf{g}^{[\mathcal{M}]}$ depicts measured data. In the simulated-data
study considered in the work, $\mathbf{g}^{[\mathcal{M}]}=\mathbf{g}$ obtained
with the truth SS image in Eq. (1), whereas in the real-data study,
$\mathbf{g}^{[\mathcal{M}]}$ are acquired in the experiment, which contain
physical factors such as noise that are not included in (and thus are
inconsistent with) the DD-data model in Eq. (1). The necessary convergence
conditions for the DTV algorithm thus are given by

$\widetilde{D}^{(n)}_{\mathbf{g}}\rightarrow
C,\,\,\,\,\,\,\widetilde{D}^{(n)}_{\rm TV_{\it x}}\rightarrow
0,\,\,\,\,\,\,\widetilde{D}^{(n)}_{\rm TV_{\it y}}\rightarrow
0,\,\,\,\,\,\,\widetilde{D}^{(n)}_{\rm TV_{\it z}}\rightarrow
0,\,\,\,\,\,\,\widetilde{D}^{(n)}_{\rm TV_{\it B}}\rightarrow 0,$ (10)

as $n\rightarrow\infty$, where $C$ is a constant that is 0 for simulated data
or $>0$ for real data, respectively.

Figure B.1: Convergence metrics (a) $\widetilde{D}^{(n)}_{\mathbf{g}}$, (b)
$\widetilde{D}^{(n)}_{\rm TV_{\it x}}$, (c) $\widetilde{D}^{(n)}_{\rm TV_{\it
y}}$, (d) $\widetilde{D}^{(n)}_{\rm TV_{\it z}}$, (e)
$\widetilde{D}^{(n)}_{\rm TV_{\it B}}$, and (f) ${\rm
nRMSE}(\mathbf{f}^{(n)})$ of the DTV algorithm, as functions of iteration
number $n$ (in log-log scales).

###  B.2 Numerical verification of the DTV algorithm

It is critical to verify that the DTV algorithm and its computer
implementation can accurately solve the optimization program in Eq. (2.3) in
terms of numerically achieving the convergence conditions in Eq. (10). We have
conducted extensive verification studies by using FAR and LAR data in both
simulated- and real-data studies. Without loss of generality, we summarize
here the verification study involving FAR data in the simulated-data study.

In Figs. B.1a-B.1e, the convergence metrics of the DTV algorithm are shown as
functions of iteration number $n$ (in log-log scales); and these convergence
curves clearly demonstrate that the DTV algorithm (and its computer
implementation) can accurately solve the optimization program in Eq. (2.3) by
achieving the necessary convergence conditions in Eq. (10). As mentioned
above, convergence results similar to those in Fig. B.1 are obtained also for
LAR scans in the simulated-data study and for FAR and LAR scans in the real-
data study.

Furthermore, we also plot in Fig. B.1f metric ${\rm nRMSE}(\mathbf{f}^{(n)})$
of Eq. (4) in which the truth SS image is used as the reference image. It can
be observed that the ${\rm nRMSE}(\mathbf{f}^{(n)})$ curve approaches zero as
$n\rightarrow\infty$, suggesting that the final SS image reconstructed is
numerically identical to the truth SS image. The result reveals empirically
that the DTV algorithm (and its computer implementation) can reconstruct
accurately an SS image from FAR data consistent with the DD-data model in Eq.
(1). In other words, the DTV algorithm and its computer implementation are
verified empirically to invert system matrix $\mathcal{H}$ in the DD-data
model in Eq. (1) for the FAR scan in CW-ZM EPRI. Once verified, the DTV
algorithm (and its computer implementation) can be used for empirically
investigating reconstructions of SS images from LAR data and/or from FAR and
LAR data containing physical components such as noise inconsistent with the
DD-data model based upon which the DTV algorithm is developed.

## References

  * Berliner and Fujii [1985] Berliner LJ, Fujii H.  Magnetic resonance imaging of biological specimens by electron paramagnetic resonance of nitroxide spin labels.  _Science_ 1985;227(4686):517–9.
  * Kuppusamy et al. [1994] Kuppusamy P, Chzhan M, Vij K, Shteynbuk M, Lefer DJ, Giannella E, Zweier JL.  Three-dimensional spectral-spatial EPR imaging of free radicals in the heart: a technique for imaging tissue metabolism and oxygenation.  _Proc Nation Academ Scien_ 1994;91(8):3388–92.
  * Elas et al. [2003] Elas M, Williams BB, Parasca A, Mailer C, Pelizzari CA, Lewis MA, River JN, Karczmar GS, Barth ED, Halpern HJ.  Quantitative tumor oxymetric images from 4D electron paramagnetic resonance imaging (EPRI): methodology and comparison with blood oxygen level-dependent (BOLD) MRI.  _Magn Reson Med_ 2003;49(4):682–91.
  * Eaton et al. [2004] Eaton SS, Eaton GR, Berliner LJ.  Biomedical EPR-Part A: Free radicals, metals, medicine and physiology.  Springer Science & Business Media; 2004\.
  * Matsumoto et al. [2006] Matsumoto Ki, Subramanian S, Devasahayam N, Aravalluvan T, Murugesan R, Cook JA, Mitchell JB, Krishna MC.  Electron paramagnetic resonance imaging of tumor hypoxia: enhanced spatial and temporal resolution for in vivo pO2 determination.  _Magn Reson Med_ 2006;55(5):1157–63.
  * Elas et al. [2008] Elas M, Bell R, Hleihel D, Barth ED, McFaul C, Haney CR, Bielanska J, Pustelny K, Ahn KH, Pelizzari CA, Kocherginsky M, Halpern HJ.  Electron paramagnetic resonance oxygen image hypoxic fraction plus radiation dose strongly correlates with tumor cure in FSa fibrosarcomas.  _Int J Radiat Oncol Biol Phys_ 2008;71(2):542–9.
  * Epel et al. [2011] Epel B, Sundramoorthy SV, Barth ED, Mailer C, Halpern HJ.  Comparison of 250 MHz electron spin echo and continuous wave oxygen EPR imaging methods for in vivo applications.  _Med Phys_ 2011;38(4):2045–52.
  * Biller et al. [2016] Biller JR, Mitchell DG, Tseytlin M, Elajaili H, Rinard GA, Quine RW, Eaton SS, Eaton GR.  Rapid scan electron paramagnetic resonance opens new avenues for imaging physiologically important parameters in vivo.  _J Vis Exp_ 2016;(115):e54068.
  * Epel et al. [2017] Epel B, Kotecha M, Halpern HJ.  In vivo preclinical cancer and tissue engineering applications of absolute oxygen imaging using pulse EPR.  _J Magn Reson_ 2017;280:149–57.
  * Eaton et al. [2018] Eaton GR, Eaton SS, Ohno K, Press C.  EPR imaging and in vivo EPR.  CRC press Boca Raton, FL; 2018.
  * Maltempo et al. [1987] Maltempo MM, Eaton SS, Eaton GR.  Spectral-spatial two-dimensional EPR imaging.  _J Magn Reson_ 1987;72(3):449–55.
  * Kuppusamy et al. [1995] Kuppusamy P, Chzhan M, Samouilov A, Wang P, Zweier JL.  Mapping the spin-density and lineshape distribution of free radicals using 4D spectral-spatial EPR imaging.  _J Magn Reson_ 1995;107(2):116–25.
  * Mäder et al. [1997] Mäder K, Nitschke S, Stösser R, Borchert HH, Domb A.  Non-destructive and localized assessment of acidic microenvironments inside biodegradable polyanhydrides by spectral spatial electron paramagnetic resonance imaging.  _Polymer_ 1997;38(19):4785–94.
  * Kuppusamy et al. [1998] Kuppusamy P, Shankar RA, Zweier JL.  In vivo measurement of arterial and venous oxygenation in the rat using 3D spectral-spatial electron paramagnetic resonance imaging.  _Phys Med Biol_ 1998;43(7):1837–44.
  * Williams et al. [2005] Williams BB, Pan X, Halpern HJ.  EPR imaging: The relationship between CW spectra acquired from an extended sample subjected to fixed stepped gradients and the Radon transform of the resonance density.  _J Magn Reson_ 2005;174(1):88–96.
  * Pan et al. [2007] Pan X, Xia D, Halpern H.  Targeted-ROI imaging in electron paramagnetic resonance imaging.  _J Magn Reson_ 2007;187(1):66–77.
  * Tseitlin et al. [2007] Tseitlin M, Dhami A, Eaton SS, Eaton GR.  Comparison of maximum entropy and filtered back-projection methods to reconstruct rapid-scan EPR images.  _J Magn Reson_ 2007;184(1):157–68.
  * Tseitlin et al. [2008] Tseitlin M, Czechowski T, Eaton SS, Eaton GR.  Regularized optimization (RO) reconstruction for oximetric EPR imaging.  _J Magn Reson_ 2008;194(2):212–21.
  * Tseitlin et al. [2014] Tseitlin M, Biller JR, Elajaili H, Khramtsov VV, Dhimitruka I, Eaton GR, Eaton SS.  New spectral–spatial imaging algorithm for full EPR spectra of multiline nitroxides and pH sensitive trityl radicals.  _J Magn Reson_ 2014;245:150–5.
  * Komarov and Hirata [2017] Komarov DA, Hirata H.  Fast backprojection-based reconstruction of spectral-spatial EPR images from projections with the constant sweep of a magnetic field.  _J Magn Reson_ 2017;281:44–50.
  * Poole [1996] Poole CP.  Electron spin resonance: a comprehensive treatise on experimental techniques.  Courier Corporation; 1996.
  * Subramanian et al. [2004] Subramanian S, Matsumoto Ki, Mitchell JB, Krishna MC.  Radio frequency continuous-wave and time-domain EPR imaging and Overhauser-enhanced magnetic resonance imaging of small animals: instrumental developments and comparison of relative merits for functional imaging.  _NMR in Biomed_ 2004;17(5):263–94.
  * Som et al. [2007] Som S, Potter LC, Ahmad R, Kuppusamy P.  A parametric approach to spectral–spatial EPR imaging.  _J Magn Reson_ 2007;186(1):1–10.
  * Ahn and Halpern [2007a] Ahn KH, Halpern HJ.  Simulation of 4D spectral-spatial EPR images.  _J Magn Reson_ 2007a;187(1):1–9.
  * Zhang et al. [2021] Zhang Z, Chen B, Xia D, Sidky EY, Pan X.  Directional-TV algorithm for image reconstruction from limited-angular-range data.  _Med Imag Analy_ 2021;(70):102030.
  * Chen et al. [2021] Chen B, Zhang Z, Xia D, Sidky EY, Pan X.  Dual-energy CT imaging with limited-angular-range data.  _Phys Med Biol_ 2021;66(18):185020.
  * Zhang et al. [2022] Zhang Z, Chen B, Xia D, Sidky EY, Pan X.  Image reconstruction from data over two orthogonal arcs of limited-angular ranges.  _Med Phys_ 2022;49(3):1468–80.
  * Chambolle and Pock [2011] Chambolle A, Pock T.  A first-order primal-dual algorithm for convex problems with applications to imaging.  _J Math Imag Vis_ 2011;40(1):1 – 26.
  * Sidky et al. [2012] Sidky EY, Jørgensen JH, Pan X.  Convex optimization problem prototyping for image reconstruction in computed tomography with the Chambolle–Pock algorithm.  _Phys Med Biol_ 2012;57(10):3065–91.
  * Pearson [1895] Pearson K.  Notes on regression and inheritance in the case of two parents.  _Proc Roy Soc Lond_ 1895;58:240–2.
  * Pluim et al. [2003] Pluim JPW, Maintz JBA, Viergever MA.  Mutual-information-based registration of medical images: a survey.  _IEEE Trans Med Imag_ 2003;22(8):986–1004.
  * Bian et al. [2010] Bian J, Siewerdsen JH, Han X, Sidky EY, Prince JL, Pelizzari CA, Pan X.  Evaluation of sparse-view reconstruction from flat-panel-detector cone-beam CT.  _Phys Med Biol_ 2010;55(22):6575–99.
  * Robinson et al. [1999] Robinson B, Mailer C, Reese A.  Linewidth analysis of spin labels in liquids: I. Theory and data analysis.  _J Magn Reson_ 1999;138(2):199–209.
  * Mailer et al. [2003] Mailer C, Robinson BH, Williams BB, Halpern HJ.  Spectral fitting: The extraction of crucial information from a spectrum and a spectral image.  _Magn Reson Med_ 2003;49(6):1175–80.
  * Halpern et al. [1989] Halpern HJ, Spencer DP, van Polen J, Bowman MK, Nelson AC, Dowey EM, Teicher BA.  Imaging radio frequency electron-spin-resonance spectrometer with high resolution and sensitivity for in vivo measurements.  _Rev Sci Instrum_ 1989;60(6):1040–50.
  * Ahn and Halpern [2007b] Ahn KH, Halpern HJ.  Spatially uniform sampling in 4-D EPR spectral-spatial imaging.  _J Magn Reson_ 2007b;185(1):152–8.
  * Sotgiu et al. [1987] Sotgiu A, Gazzillo D, Momo F.  ESR imaging: spatial deconvolution in the presence of an asymmetric hyperfine structure.  _J Phy C_ 1987;20(36):6297–304.
  * Ardenkjaer-Larsen et al. [1998] Ardenkjaer-Larsen J, Laursen I, Leunbach I, Ehnholm G, Wistrand LG, Petersson J, Golman K.  EPR and DNP properties of certain novel single electron contrast agents intended for oximetric imaging.  _J Magn Reson_ 1998;133(1):1–12.
  * Zhang et al. [2016] Zhang Z, Han X, Pearson E, Pelizzari C, Sidky EY, Pan X.  Artifact reduction in short-scan CBCT by use of optimization-based reconstruction.  _Phys Med Biol_ 2016;61(9):3387–406.
  * Sidky and Pan [2008] Sidky EY, Pan X.  Image reconstruction in circular cone-beam computed tomography by constrained, total-variation minimization.  _Phys Med Biol_ 2008;53(17):4777–807.
  * Christodoulou et al. [2016] Christodoulou AG, Redler G, Clifford B, Liang ZP, Halpern HJ, Epel B.  Fast dynamic electron paramagnetic resonance (EPR) oxygen imaging using low-rank tensors.  _J Magn Reson_ 2016;270:176–82.
  * Johnson et al. [2014] Johnson DH, Ahmad R, He G, Samouilov A, Zweier JL.  Compressed sensing of spatial electron paramagnetic resonance imaging.  _Magn Reson Med_ 2014;72(3):893–901.
  * Qiao et al. [2018] Qiao Z, Zhang Z, Pan X, Epel B, Redler G, Xia D, Halpern H.  Optimization-based image reconstruction from sparsely sampled data in electron paramagnetic resonance imaging.  _J Magn Reson_ 2018;294:24–34.