Training End-to-End Unrolled Iterative Neural Networks for SPECT Image
  Reconstruction
Training end-to-end unrolled iterative neural networks for SPECT image
reconstruction requires a memory-efficient forward-backward projector for
efficient backpropagation. This paper describes an open-source, high
performance Julia implementation of a SPECT forward-backward projector that
supports memory-efficient backpropagation with an exact adjoint. Our Julia
projector uses only ~5% of the memory of an existing Matlab-based projector. We
compare unrolling a CNN-regularized expectation-maximization (EM) algorithm
with end-to-end training using our Julia projector with other training methods
such as gradient truncation (ignoring gradients involving the projector) and
sequential training, using XCAT phantoms and virtual patient (VP) phantoms
generated from SIMIND Monte Carlo (MC) simulations. Simulation results with two
different radionuclides (90Y and 177Lu) show that: 1) For 177Lu XCAT phantoms
and 90Y VP phantoms, training unrolled EM algorithm in end-to-end fashion with
our Julia projector yields the best reconstruction quality compared to other
training methods and OSEM, both qualitatively and quantitatively. For VP
phantoms with 177Lu radionuclide, the reconstructed images using end-to-end
training are in higher quality than using sequential training and OSEM, but are
comparable with using gradient truncation. We also find there exists a
trade-off between computational cost and reconstruction accuracy for different
training methods. End-to-end training has the highest accuracy because the
correct gradient is used in backpropagation; sequential training yields worse
reconstruction accuracy, but is significantly faster and uses much less memory.
Zongyu Li, Yuni K. Dewaraja, Jeffrey A. Fessler
# Training End-to-End Unrolled Iterative Neural Networks for SPECT Image Reconstruction

###### Abstract

Training end-to-end unrolled iterative neural networks for SPECT image reconstruction requires a memory-efficient forward-backward projector for efficient backpropagation. This paper describes an open-source, high performance Julia implementation of a SPECT forward-backward projector that supports memory-efficient backpropagation with an exact adjoint. Our Julia projector uses only \(\sim\)5% of the memory of an existing Matlab-based projector. We compare unrolling a CNN-regularized expectation-maximization (EM) algorithm with end-to-end training using our Julia projector with other training methods such as gradient truncation (ignoring gradients involving the projector) and sequential training, using XCAT phantoms and virtual patient (VP) phantoms generated from SIMIND Monte Carlo (MC) simulations. Simulation results with two different radionuclides (\({}^{50}\)Y and \({}^{177}\)Lu) show that: 1) For \({}^{177}\)Lu XCAT phantoms and \({}^{50}\)Y VP phantoms, training unrolled EM algorithm in end-to-end fashion with our Julia projector yields the best reconstruction quality compared to other training methods and OSEM, both qualitatively and quantitatively. For VP phantoms with \({}^{177}\)Lu radionuclide, the reconstructed images using end-to-end training are in higher quality than using sequential training and OSEM, but are comparable with using gradient truncation. We also find there exists a trade-off between computational cost and reconstruction accuracy for different training methods. End-to-end training has the highest accuracy because the correct gradient is used in backpropagation; sequential training yields worse reconstruction accuracy, but is significantly faster and uses much less memory.

 End-to-end learning, regularized model-based image reconstruction, backpropagatable forward-backward projector, quantitative SPECT.

## I Introduction

Single photon emission computerized tomography (SPECT) is a nuclear medicine technique that images spatial distributions of radioisotopes and plays a pivotal role in clinical diagnosis, and in estimating radiation-absorbed doses in nuclear medicine therapies [1, 2]. For example, quantitative SPECT imaging with Lutetium-177 (\({}^{177}\)Lu) in targeted radionuclide therapy (such as \({}^{177}\)Lu DOTATATE) is important in determining dose-response relationships in tumors and holds great potential for dosimetry-based individualized treatment. Additionally, quantitative Ytrium-90 (\({}^{90}\)Y) bremsstrahlung SPECT imaging is valuable for safety assessment and absorbed dose verification after \({}^{90}\)Y radioembolization in liver malignancies. However, SPECT imaging suffers from noise and limited spatial resolution due to the collimator response; the resulting reconstruction problem is hence ill-posed and challenging to solve.

Numerous reconstruction algorithms have been proposed for SPECT reconstruction, of which the most popular ones are model-based image reconstruction algorithms such as maximum likelihood expectation maximization (MLEM) [3] and ordered-subset EM (OSEM) [4]. These methods first construct a mathematical model for the SPECT imaging system, then maximize the (log-)likelihood for a Poisson noise model. Although MLEM and OSEM have achieved great success in clinical use, they have a trade-off between recovery and noise. To address that trade-off, researchers have proposed alternatives such as regularization-based (or maximum a posteriori in Bayesian setting) reconstruction methods [5, 6, 7]. For example, Panin _et al._[5] proposed total variation (TV) regularization for SPECT reconstruction. However, TV regularization may lead to "blocky" images and over-smoothing the edges. One way to overcome blurring edges is to incorporate anatomical boundary side information from CT images [8], but that method requires accurate organ segmentation. Chun _et al._[9] used non-local means (NLM) filters that exploit the self-similarity of patches in images for regularization, yet that method is computationally expensive and hence less practical. In general, choosing an appropriate regularizer can be challenging; moreover, these traditional regularized algorithms may lack generalizability to images that do not follow assumptions made by the prior.

With the recent success of deep learning (DL) and especially convolutional neural networks (CNN), DL methods have been reported to outperform conventional algorithms in many medical imaging applications such as in MRI [10, 11, 12], CT [13, 14] and PET reconstruction [15, 16, 17]. However, fewer DL approaches to SPECT reconstruction appear in the literature. Reference [18] proposed "SPECTnet" with a two-step training strategy that learns the transformation from projection space to image space as an alternative to the traditional OSEM algorithm. Reference [19] also proposed a DL method that
can directly reconstruct the activity image from the SPECT projection data, even with reduced view angles. Reference [20] trained a neural network that maps non-attenuation-corrected SPECT images to those corrected by CT images as a post-processing procedure to enhance the reconstructed image quality.

Though promising results were reported with these methods, most of them worked in 2D whereas 3D is used in practice [18, 19]. Furthermore, there has yet to be an investigation of end-to-end training of CNN regularizers that are embedded in unrolled SPECT iterative statistical algorithms such as CNN-regularized EM. End-to-end training is popular in machine learning and other medical imaging fields such as MRI image reconstruction [21], and is reported to meet data-driven regularization for inverse problems [22]. But for SPECT image reconstruction, end-to-end training is nontrivial to implement due to its complicated system matrix. Alternative training methods have been proposed, such as sequential training [23, 24, 25, 26] and gradient truncation [27]; these methods were shown to be effective, though they could yield sub-optimal reconstruction results due to approximations to the training loss gradient. Another approach is to construct a neural network that also models the SPECT system matrix, like in "SPECTnet" [18], but this approach lacks interpretability compared to algorithms like unrolled CNN-regularized EM, i.e., if one sets the regularization parameter to zero, then the latter becomes identical to the traditional EM.

As an end-to-end training approach has not yet been investigated for SPECT image reconstruction, this paper first describes a SPECT forward-backward projector written in the open-source and high performance Julia language that enables efficient auto-differentiation. Then we compare the end-to-end training approach with other non-end-to-end training methods.

The structure of this article is as follows. Section II describes the implementation of our Julia projector and discusses end-to-end training and other training methods for the unrolled EM algorithm. Section III compares the accuracy, speed and memory use of our Julia projector with Monte Carlo (MC) and a Matlab-based projector, and then compares reconstructed images with end-to-end training versus sequential training and gradient truncation on different datasets (XCAT and VP phantoms), using qualitative and quantitative evaluation metrics. Section IV and V conclude this paper and discuss future works.

_Notation:_ Bold upper/lower case letters (e.g., \(\mathbf{A}\), \(\mathbf{x}\), \(\mathbf{y}\), \(\mathbf{b}\)) denote matrices and column vectors, respectively. Italics (e.g., \(\mu,y,b\)) denote scalars. \(y_{i}\) and \(b_{i}\) denote the \(i\)th element in vector \(\mathbf{y}\) and \(\mathbf{b}\), respectively. \(\mathbb{R}^{N}\) and \(\mathbb{C}^{N}\) denote \(N\)-dimensional real/complex normed vector space, respectively. \((\cdot)^{*}\) denotes the complex conjugate and \((\cdot)^{\prime}\) denotes Hermitian transpose.

## II Methods

This section summarizes the Julia SPECT projector, a DL-based image reconstruction method as well as the dataset used in experiments and other experiment setups.

### _Implementation of Julia SPECT projector_

Our Julia implementation of SPECT projector is based on [28], modeling parallel-beam collimator geometries. Our projector also accounts for attenuation and depth-dependent collimator response. We did not model the scattering events like Compton scatter and coherent scatter of high energy gamma rays within the object. Fig. 1 illustrates the SPECT imaging system modeled in this paper.

For the forward projector, at each rotation angle, we first rotate the 3D image matrix \(\mathbf{x}\in\mathbb{R}^{n_{x}\times n_{y}\times n_{z}}\) according to the third dimension by its projection angle \(\theta_{l}\) (typically \(2\pi(l-1)/n_{\text{view}}\)); \(l\) denotes the view index, which ranges from 1 to \(n_{\text{view}}\) and \(n_{\text{view}}\) denotes the total number of projection views. We implemented and compared (results shown in Section III) both bilinear interpolation and 3-pass 1D linear interpolation [29] with zero padding boundary condition for image rotation. For attenuation correction, we first rotated the 3D attenuation map \(\mathbf{\mu}\in\mathbb{R}^{n_{x}\times n_{y}\times n_{z}}\) (obtained from transmission tomography) also by \(\theta_{l}\), yielding a rotated 3D array \(\tilde{\mu}(i,j,k;l)\), where \(i,j,k\) denotes the 3D voxel coordinate. Assuming \(n_{y}\) is the index corresponding to the closest plane of \(\mathbf{x}\) to the detector, then we model the accumulated attenuation factor \(\bar{\mu}\) for each view angle as

\[\bar{\mu}(i,j,k;l)=\mathrm{e}^{-\Delta_{y}\left(\frac{1}{2}\bar{\mu}(i,j,k;l) +\sum_{x=j+1}^{n_{y}}\bar{\mu}(i,s,k;l)\right)}, \tag{1}\]

where \(\Delta_{y}\) denotes the voxel size for the (first and) second coordinate. Next, for each \(y\) slice (an \((x,z)\) plane for a given \(j\) index) of the rotated and attenuated image, we convolve with the appropriate slice of the depth-dependent point spread function \(\mathbf{p}\in\mathbb{R}^{p_{k}\times p_{k}\times n_{y}\times n_{\text{view}}}\) using a 2D fast Fourier transform (FFT). Here we use replicate padding for both the

Fig. 1: SPECT imaging model for parallel-beam collimators, with attenuation and depth-dependent collimator point spread response.
\(i\) and \(k\) coordinates. The view-dependent PSF accommodates non-circular orbits. Finally, the forward projection operation simply sums the rotated, blurred and attenuated activity image \(\mathbf{x}\) along the second coordinate \(j\). Algorithm 1 summarizes the forward projector, where denotes a 2D convolution operation.

```
Input: 3D image \(\mathbf{x}\in\mathbb{R}^{n_{x}\times n_{y}\times n_{x}}\), 3D attenuation map \(\mathbf{\mu}\in\mathbb{R}^{n_{x}\times n_{y}\times n_{x}}\), 4D point spread function \(\mathbf{p}\in\mathbb{R}^{p_{x}\times p_{x}\times n_{y}\times n_{\rm view}}\), voxel size \(\Delta_{y}\). Initialize:\(\mathbf{v}\in\mathbb{R}^{n_{x}\times n_{x}\times n_{\rm view}}\) as all zeros. for\(l=1,...,n_{\rm view}\)do \(\tilde{\mathbf{x}}\leftarrow\) rotate \(\mathbf{x}\) by \(\theta_{l}\) \(\tilde{\mathbf{\mu}}\leftarrow\) rotate \(\mathbf{\mu}\) by \(\theta_{l}\) for\(j=1,...,n_{y}\)do \(\tilde{\mathbf{\mu}}\leftarrow\) calculate by (1) using \(\tilde{\mathbf{\mu}}\) \(\tilde{\mathbf{x}}(i,j,k)\leftarrow\tilde{\mu}(i,j,k)\) \(v(i,k,l)\leftarrow\tilde{x}(i,j,k)\)\(\triangleright\)\(p(i,k;j,l)\)  end for end Output: projection views \(\mathbf{v}\in\mathbb{R}^{n_{x}\times n_{x}\times n_{\rm view}}\)
```

**Algorithm 1**SPECT forward projector

All of these steps are linear, so hereafter, we use \(\mathbf{A}\) to denote the forward projector, though it is not stored explicitly as a matrix. As each step is linear, each step has an adjoint operation. So the backward projector \(\mathbf{A}^{\prime}\) is the adjoint of \(\mathbf{A}\) that satisfies

\[\langle\mathbf{A}\mathbf{x},\mathbf{y}\rangle=\langle\mathbf{x},\mathbf{A}^{\prime}\mathbf{y}\rangle, \quad\forall\mathbf{x},\mathbf{y}. \tag{2}\]

The exact adjoint of (discrete) image rotation is not simply a discrete rotation of the image by \(-\theta_{l}\). Instead, one should also consider the adjoint of linear interpolation. For the adjoint of convolution, we assume the point spread function is symmetric along coordinates \(i\) and \(k\) so that the adjoint convolution operator is just the forward convolution operator along with the adjoint of replicate padding. Algorithm 2 summarizes the SPECT backward projector.

```
Input: Array of 2D projection views \(\mathbf{v}\in\mathbb{R}^{n_{x}\times n_{x}\times n_{\rm view}}\), 3D attenuation map \(\mathbf{\mu}\in\mathbb{R}^{n_{x}\times n_{y}\times n_{x}}\), 4D point spread function \(\mathbf{p}\in\mathbb{R}^{p_{x}\times p_{x}\times n_{y}\times n_{\rm view}}\), voxel size \(\Delta_{y}\). Initialize:\(\mathbf{x}\in\mathbb{R}^{n_{x}\times n_{y}\times n_{x}}\) as all zeros. for\(l=1,...,n_{\rm view}\)do \(\tilde{\mathbf{\mu}}\leftarrow\) rotation \(\mathbf{\mu}\) by \(\theta_{l}\) for\(j=1,...,n_{y}\)do \(\tilde{\mathbf{\mu}}\leftarrow\) calculate by (1) using \(\tilde{\mathbf{\mu}}\) \(\tilde{\mathbf{v}}(i,k,l)\leftarrow\) adjoint of \(v(i,k,l)\otimes p(i,k,j,l)\) \(\tilde{x}(i,j,k)\leftarrow\tilde{v}(i,k,l)\cdot\tilde{\mu}(i,j,k;l)\)  end for\(\mathbf{x}+=\) adjoint rotate \(\tilde{\mathbf{x}}\) by \(\theta_{l}\)  end for Output:\(\mathbf{x}\in\mathbb{R}^{n_{x}\times n_{y}\times n_{x}}\)
```

**Algorithm 2**SPECT backward projector

To accelerate the for-loop process, we used multi-threading to enable projecting or backprojecting multiple angles at the same time. To reduce memory use, we pre-allocated necessary arrays and used fully in-place operations inside the for-loop in forward and backward projection. To further accelerate auto-differentiation, we customized the chain rule to use the linear operator \(\mathbf{A}\) or \(\mathbf{A}^{\prime}\) as the Jacobian when calling \(\mathbf{A}\mathbf{x}\) or \(\mathbf{A}^{\prime}\mathbf{y}\) during backpropagation. We implemented and tested our projector in Julia v1.6 ; we also implemented a GPU version in Julia (using CUDA.jl) that runs efficiently on a GPU by eliminating explicit scalar indexing. For completeness, we also provide a PyTorch version but without multi-threading support, in-place operations nor the exact adjoint of image rotation.

### _Unrolled CNN-regularized EM algorithm_

Model-based image reconstruction algorithms seek to estimate image \(\mathbf{x}\in\mathbb{R}^{N}\) from noisy measurements \(\mathbf{y}\in\mathbb{R}^{M}\) with imaging model \(\mathbf{A}\in\mathbb{R}^{M\times N}\). In SPECT reconstruction, the measurements \(\mathbf{y}\) are often modeled by

\[\mathbf{y}\sim\text{Poisson}(\mathbf{A}\mathbf{x}+\bar{\mathbf{r}}), \tag{3}\]

where \(\bar{\mathbf{r}}\in\mathbb{R}^{M}\) denotes the vector of means of background events such as scatters. Combining regularization with the Poisson negative log-likelihood yields the following optimization problem:

\[\hat{\mathbf{x}} =\operatorname*{arg\,min}_{\mathbf{x}\geq\mathbf{0}}f(\mathbf{x})+R(\mathbf{x}),\] \[f(\mathbf{x}) \triangleq\mathbf{1}^{\prime}(\mathbf{A}\mathbf{x}+\bar{\mathbf{r}})-\mathbf{y}^{ \prime}\log(\mathbf{A}\mathbf{x}+\bar{\mathbf{r}}), \tag{4}\]

where \(f(\mathbf{x})\) is the data fidelity term and \(R(\mathbf{x})\) denotes the regularizer. For deep learning regularizers, we follow [23] and formulate \(R(\mathbf{x})\) as

\[R(\mathbf{x})\triangleq\frac{\beta}{2}\|\mathbf{x}-\mathbf{g}_{\mathbf{\theta}}(\mathbf{x})\|_{2}^ {2}, \tag{5}\]

where \(\beta\) denotes the regularization parameter; \(\mathbf{g}_{\mathbf{\theta}}\) denotes a neural network with parameter \(\mathbf{\theta}\) that is trained to learn to enhance the image quality.

Based on (4), a natural reconstruction approach is to apply variable splitting with \(\mathbf{u}=\mathbf{g}_{\mathbf{\theta}}(\mathbf{x})\) and then alternatively update the images \(\mathbf{x}\) and \(\mathbf{u}\) as follows

\[\mathbf{u}_{k+1} =\mathbf{g}_{\mathbf{\theta}}(\mathbf{x}_{k}),\] \[\mathbf{x}_{k+1} =\operatorname*{arg\,min}_{\mathbf{x}\geq 0}f(\mathbf{x})+\frac{\beta}{2}\| \mathbf{x}-\mathbf{u}_{k+1}\|_{2}^{2}, \tag{6}\]

where subscript \(k\) denotes the iteration number. To minimize (6), we used the EM-surrogate from [30] as summarized in [23], leading to the following vector update:

\[\hat{\mathbf{x}}_{k} =\frac{1}{2\beta}\left(-\mathbf{d}(\beta)+\sqrt{\mathbf{d}(\beta)^{2}+4 \beta\mathbf{x}_{k}\odot\mathbf{e}(\mathbf{x}_{k})}\right), \tag{7}\] \[\mathbf{d}(\beta)\triangleq\mathbf{A}^{\prime}\mathbf{1}-\beta\mathbf{u}_{k}, \quad\mathbf{e}(\mathbf{x}_{k})\triangleq\mathbf{A}^{\prime}(\mathbf{y}\oslash(\mathbf{A}\mathbf{x}_{k}+ \bar{\mathbf{r}}))\,, \tag{8}\]

where \(\odot\) and \(\oslash\) denote element-wise multiplication and division, respectively. To compute \(\mathbf{x}_{k+1}\), one must substitue \(\hat{\mathbf{x}}_{k}\) back into \(\mathbf{e}(\cdot)\) in (8), and repeat. Hereafter, we refer to
(6) as one outer iteration and (7) as one inner EM iteration. Algorithm 3 summarizes the CNN-regularized EM algorithm.

```
Input: 3D projection measurements \(\mathbf{y}\), 3D background measurements \(\mathbf{\bar{r}}\),  system model \(\mathbf{A}\), initial guess \(\mathbf{x_{0}}\),  deep neural network \(\mathbf{g_{\theta}}\),  outer iterations \(K\). for\(k=0,...,K-1\)do \(\mathbf{u_{k+1}}=\mathbf{g_{\theta}}(\mathbf{x_{k}})\) \(\mathbf{x_{k+1}}\leftarrow\) repeat (7) until convergence tolerance or maximum # of inner iterations is reached  end for Output:\(\mathbf{x_{K}}\)
```

**Algorithm 3**SPECT CNN-regularized EM algorithm

To train \(\mathbf{g_{\theta}}\), the most direct way is to unroll Algorithm 3 and train end-to-end with an appropriate target; this supervised approach requires backpropagating through the SPECT system model, which is not trivial to implement with previous SPECT projection tools due to the memory issues. Non-end-to-end training methods, e.g., sequential training [23], first train \(\mathbf{u}_{k}\) by the target and then plug into (7) at each iteration. This method must use non-shared weights for the neural network per each iteration. Another method is gradient truncation [27] that ignores the gradient involving the system matrix \(\mathbf{A}\) and its adjoint \(\mathbf{A}^{\prime}\) during backpropagation. Both of these training methods, though reported to be effective, may be sub-optimal because they approximate the overall training loss gradients.

### _Phantom Dataset and Simulation Setup_

We used simulated XCAT phantoms [31] and virtual patient phantoms for experiment results presented in Section III. Each XCAT phantom was simulated to approximately follow the activity distributions observed when imaging patients after \({}^{177}\mathrm{Lu}\) DOTATATE therapy. We set the image size to \(128\times 128\times 80\) with voxel size \(4.8\times 4.8\times 4.8\mathrm{mm}^{3}\). Tumors of various shapes and sizes (5-100mL) were located in the liver as is typical for patients undergoing this therapy.

For virtual patient phantoms, we consider two radionuclides: \({}^{177}\mathrm{Lu}\) and \({}^{90}\mathrm{Y}\). For \({}^{177}\mathrm{Lu}\) phantoms, the true images were from PET/CT scans of patients who underwent diagnostic \({}^{68}\)Ga DOTATATE PET/CT imaging (Siemens Biograph mCT) to determine eligibility for \({}^{177}\mathrm{Lu}\) DOTATATE therapy. The \({}^{68}\)Ga DOTATATE distribution in patients is expected to be similar to \({}^{177}\mathrm{Lu}\) and hence can provide a reasonable approximation to the activity distribution of \({}^{177}\mathrm{Lu}\) in patients for DL training purposes but at higher resolution. The PET images had size \(200\times 200\times 577\) and voxel size \(4.073\times 4.073\times 2\) mm\({}^{3}\) and were obtained from our Siemens mCT (resolution is 5-6 mm FWHM [32]) and reconstructed using the standard clinic protocol: 3D OSEM with three iterations, 21 subsets, including resolution recovery, time-of-flight, and a 5mm (FWHM) Gaussian post-reconstruction filter. The density maps were also generated using the experimentally derived CT-to-density calibration curve.

For \({}^{90}\mathrm{Y}\) phantoms, the true activity images were reconstructed (using a previously implemented 3D OSEM reconstruction with CNN-based scatter estimation [33]) from \({}^{90}\mathrm{Y}\) SPECT/CT scans of patients who underwent \({}^{90}\mathrm{Y}\) microsphere radioembolization in our clinic.

In total, we simulated 4 XCAT phantoms, 8 \({}^{177}\mathrm{Lu}\) and 8 \({}^{90}\mathrm{Y}\) virtual patient phantoms. We repeated all of our experiments 3 times with different noise realizations. All image data have University of Michigan Institutional Review Board (IRB) approval for retrospective analysis. For all simulated phantoms, we selected the center slices covering the lung, liver and kidney corresponding to SPECT axial FOV (39cm).

Then we ran SIMIND Monte Carlo (MC) program [34] to generate the radial position of SPECT camera for 128 view angles. The SIMIND model parameters for \({}^{177}\mathrm{Lu}\) were based on \({}^{177}\mathrm{Lu}\) DOTATATE patient imaging in our clinic (Siemens Intevo with medium energy collimators, a 5/8" crystal, a 20% photopeak window at 208 keV, and two adjacent 10% scatter windows) [35]. For \({}^{90}\mathrm{Y}\), a high-energy collimator, 5/8" crystal, and a 105 to 195 keV acquisition energy window was modeled as in our clinical protocol for \({}^{90}\mathrm{Y}\) bremsstrahlung imaging. Next we approximated the point spread function for \({}^{177}\mathrm{Lu}\) and \({}^{90}\mathrm{Y}\) by simulating point source at 6 different distances (20, 50, 100, 150, 200, 250mm) and then fitting a 2D Gaussian distribution at each distance. The camera orbit was assumed to be non-circular (auto-contouring mode in clinical systems) with the minimum distance between the phantom surface and detector set at 1 cm.

## III Experiment Results

### _Comparison of projectors_

We used an XCAT phantom to evaluate the accuracy and memory-efficiency of our Julia projector.

#### Iii-A1 Accuracy

We first compared primary (no scatter events included) projection images and profiles generated by our Julia projector with those from MC simulation and the Matlab projector. For results of MC, we ran two SIMIND simulations for 1 billion histories using \({}^{177}\mathrm{Lu}\) and \({}^{90}\mathrm{Y}\) as radionuclide source, respectively. Each simulation took about 10 hours using a 3.2 GHz 16-Core Intel Xeon W CPU on MacOS. The Matlab projector was originally implemented and compiled in C99 and then wrapped by a Matlab MEX file as a part of the Michigan Image Reconstruction Toolbox (MIRT) [36]. The physics modeling of the Matlab projector was the same as our Julia projector except that it only implemented 3-pass 1D linear interpolation for image rotation. Unlike the memory-efficient Julia version, the Matlab version pre-rotates the patient attenuation map for all projection views. This strategy saves time during EM iterations for a single patient,
but uses considerable memory and scales poorly for DL training approaches involving multiple patient datasets.

Fig. 2 compared the primary projections generated by different methods without adding Poisson noise. Visualizations of image slices and line profiles illustrate that our Julia projector (with rotation based on 3-pass 1D interpolation) is almost identical to the Matlab projector, while both give a reasonably good approximation to the MC. Using MC as reference, the NRMSE of Julia1D/Matlab/Julia2D projectors were 7.9%/7.9%/7.6% for \({}^{177}\mathrm{Lu}\), respectively; while the NRMSE were 8.2%/8.2%/7.9% for \({}^{90}\mathrm{Y}\). We also compared the OSEM reconstructed images using Julia (2D) and Matlab projectors, where we did not observe notable difference, as shown in Fig. 3. The overall NRMSD between Matlab and Julia (2D) projector for the whole 3D OSEM reconstructed image ranged from 2.5% to 2.8% across 3 noise realizations.

#### Iii-A2 Speed and memory use

We compared the memory use and compute times between our Julia projector (with 2D bilinear interpolation) and the Matlab projector using different number of threads when projecting a \(128\times 128\times 80\) image. Fig. 4 shows that our Julia projector has comparable computing time for a single projection with 128 view angles using different number of CPU threads, while using only a very small fraction of memory (\(\sim\)5%) and pre-allocation time (\(\sim\)1%) compared to the Matlab projector.

#### Iii-A3 Adjoint of projector

We generated a set of random numbers to verify that the backprojector is an exact adjoint of the forward projector. Specifically, we generated the system matrix of size \((8\times 6\times 7)\times(8\times 8\times 6)\) using random (nonnegative) attenuation maps and random (symmetric) PSF. Fig. 5 compares the transpose of the forward projector to the backprojector. As shown in Fig. 5 (d), the Frobenius norm error of our backprojector agrees well with the regular transpose within an accuracy of \(10^{-6}\) across 100 different realizations, as expected for 32-bit floating point calculations. A more comprehensive comparison is available in the code tests at [https://github.com/JuliaImageRecon/SPECTrecon.jl](https://github.com/JuliaImageRecon/SPECTrecon.jl).

Fig. 2: Primary (scatter-free) projections generated by MC simulation, Matlab projector and our Julia projector with 3-pass 1D linear interpolation and 2D bilinear interpolation for image rotation, using \({}^{177}\)Lu and \({}^{90}\)Y radionuclides. Subfigure (i)-(l) show line profiles across tumors as shown in subfigure (a) and (e), respectively. MC projections were scaled to have the same total activities as the Matlab projector per field-of-view.
### _Comparison of CNN-regularized EM using different training methods_

This section compares end-to-end training with other training methods that have been used previously for SPECT image reconstruction, namely the gradient truncation and sequential training. The training targets were simulated activity maps on \({}^{177}\mathrm{Lu}\) XCAT phantoms and \({}^{177}\mathrm{Lu}\) & \({}^{90}\mathrm{Y}\) virtual patient phantoms. We implemented an unrolled CNN-regularized EM algorithm with 3 outer iterations, each of which had one inner iteration. Only 3 outer iterations were used (compared to previous works such as [27]) because we used the 16-iteration 4-subset OSEM reconstructed image as a warm start for all reconstruction algorithms. We set the regularization parameter (defined in (5)) as \(\beta=1\). The regularizer was a 3-layer 3D CNN, where each layer had a \(3\times 3\times 3\) convolutional filter followed by ReLU activation (except the last layer), and hence had 657 trainable parameters in total. We added the input image \(\mathbf{x}_{k}\) to the output of CNN following the common residual learning strategy [37]. End-to-end training and gradient truncation could also work with a shared weights CNN approach, but were not included here for fair comparison purpose, since the sequential training only works with non-shared weights CNN. All the neural networks were initialized with the same parameters (drawn from a Gaussian distribution) and trained on an Nvidia RTX 3090 GPU for 600 epochs by minimizing mean square error (loss) using AdamW optimizer [38] with a constant learning rate 0.002.

Besides line profiles for qualitative comparison, we also used mean activity error (MAE) and normalized root mean square error (NRMSE) as quantitative evaluation metrics, where MAE is defined as

\[\mathrm{MAE}\triangleq\left|1-\frac{\frac{1}{n_{p}}\sum_{j\in\mathrm{VOI}} \hat{\mathbf{x}}[j]}{\frac{1}{n_{p}}\sum_{j\in\mathrm{VOI}}\mathbf{x}_{\text{true}}[j] }\right|\times 100\%, \tag{9}\]

where \(n_{p}\) denotes number of voxels in the voxels of interest (VOI). \(\hat{\mathbf{x}}\) and \(\mathbf{x}_{\text{true}}\) denote the reconstructed image and the true activity map, respectively. The NRMSE is defined as

\[\mathrm{NRMSE}\triangleq\frac{\sqrt{\frac{1}{n_{p}}\sum_{j\in\mathrm{VOI}} \left(\hat{\mathbf{x}}[j]-\mathbf{x}_{\text{true}}[j]\right)^{2}}}{\sqrt{\frac{1}{n_{ p}}\left(\sum_{j\in\mathrm{VOI}}\mathbf{x}_{\text{true}}[j]\right)^{2}}}\times 100\%. \tag{10}\]

Fig. 4: Time and memory comparison between Matlab projector and our Julia projector for projecting 128 view angles of a \(128\times 128\times 80\) image. “time pre” denotes the time cost for pre-allocating necessary arrays before projection; “time proj” denotes the time cost for a single projection; “mem” denotes the memory usage. All methods were tested on MacOS with a 3.8 GHz 8-Core Intel Core i7 CPU.

Fig. 5: Accuracy of the backprojector. In subfigure (d), \(\mathbf{A}^{\prime}\) denotes regular transpose of \(\mathbf{A}\); \(\mathbf{A}_{b}\) denotes the backprojector.

Fig. 3: Comparison of one slice of the \(128\times 128\times 80\) OSEM reconstruction (16 iterations, 4 subsets) using Matlab and Julia (2D interpolation) projectors.
All activity images were scaled by a factor that normalized the whole activity to 1 MBq per field of view (FOV) before comparison. All quantitative results (Table I, Table II, Table III) were averaged across 3 different noise realizations.

#### Iv-B1 Loss function, computing time and memory use

We compared the training and validation loss using sequential training, gradient truncation and end-to-end training. We ran 1800 epochs for each method on \({}^{177}\mathrm{Lu}\) XCAT phantoms with the AdamW optimizer [38]. Fig. 6 shows that the end-to-end training achieved the lowest validation loss while it had comparable training loss with the gradient truncation (which became lower at around 1400 epochs). For visualization, we concatenated the first 600 epochs of each outer iteration for the sequential training method, as shown by the spikes in sequential training curve. We ran 600 epochs for each algorithm for subsequent experiments because the validation losses were pretty much settled at around 600 epochs.

We also compared the computing time of each training method. We found that for MLEM with 3 outer iterations and 1 inner iteration, where each outer iteration had a 3-layer convolutional neural network, sequential training took 48.6 seconds to complete a training epoch; while gradient truncation took 327.1 seconds and end-to-end training took 336.3 seconds. Under the same experiment settings, we found sequential training took less than 1GB of memory to backpropagate through one outer iteration; compared to approximately 6GB used in gradient truncation and end-to-end training that backpropagated through three outer iterations.

#### Iv-B2 Results on \({}^{177}\mathrm{Lu}\) XCAT phantoms

We evaluated the CNN-regularized EM algorithm with three training methods on 4 \({}^{177}\mathrm{Lu}\) XCAT phantoms we simulated. We generated the primary projections by calling forward operation of our Julia projector and then added uniform scatters with 10% of the primary counts before adding Poisson noise. Of the 4 phantoms, we used 2 for training, 1 for validation and 1 for testing.

Fig. 7 shows that the end-to-end training yielded incrementally better reconstruction of the tumor in the liver center over OSEM, sequential training and gradient truncation. Fig. 7 (g) also illustrates this improvement by the line profile across the tumor. For the tumor at the top-right corner of the liver, all methods had comparable performance; this can be attributed to the small tumor size (5mL) for which partial volume (PV) effects associated with SPECT resolution are higher; and hence its recovery is even more challenging.

Table I demonstrates that the CNN-regularized EM algorithm with all training methods (sequential training, gradient truncation and end-to-end training) consistently had lower reconstruction error than the OSEM method. Among all training methods, the proposed end-to-end training had lower MAE over nearly all lesions and organs than other training methods. The relative reduction in MAE by the end-to-end training was up to 32% (for lesion 3) compared to sequential training. End-to-end training also had lower NRMSE for most lesions and organs, and was otherwise comparable to other training methods. The relative improvement compared to sequential training was up to 29% (for lesion 3).

#### Iv-B3 Results on \({}^{177}\mathrm{Lu}\) VP phantoms

Next we present test results on 8 \({}^{177}\mathrm{Lu}\) virtual patient phantoms. Out of 8 \({}^{177}\mathrm{Lu}\) phantoms, we used 4 for training, 1 for validation and 3 for testing.

Fig. 8 shows that the improvement of all learning-based methods was limited compared to OSEM, which was also evident from line profiles. For example, in Fig. 8 (g), where the line profile was drawn on a small tumor. We found that OSEM yielded a fairly accurate estimate already, and

Fig. 6: Training and validation loss of three backpropagation methods.

Fig. 7: Qualitative comparison of different training methods and OSEM tested on \({}^{177}\mathrm{Lu}\) XCAT phantoms. Subfigure (a)-(c): true activity map, attenuation map and OSEM reconstruction (16 iterations and 4 subsets); (d)-(f): regularized EM using sequential training, gradient truncation, end-to-end training, respectively; (g) and (h): line profiles in (a).
we did not observe as much improvement as we had seen on \({}^{177}\mathrm{Lu}\) XCAT phantoms for end-to-end training or even learning-based methods. Table II also demonstrates this observation. The OSEM method had substantially lower MAE and NRMSE compared to the errors shown for \({}^{177}\mathrm{Lu}\) XCAT data (cf Table I). Moreover, the end-to-end training method had comparable accuracy with gradient truncation. For example, gradient truncation was the best on lesion, liver and lung in terms of MAE; end-to-end training had the lowest NRMSE on lesion, liver, lung, kidney and spleen. Perhaps this could be due to the loss function used for training, i.e., MSE loss was used in our experiments so that end-to-end training might yield lower NRMSE. A more comprehensive study would be needed to verify this conjecture.

#### Iv-B4 Results on \({}^{90}\mathrm{Y}\) VP phantoms

We also tested with 8 \({}^{90}\mathrm{Y}\) virtual patient phantoms. Of the 8 phantoms, we used 4 for training, 1 for validation and 3 for testing.

Fig. 9 compares the reconstruction quality between OSEM and CNN-regularized EM algorithm using sequential training, gradient truncation and end-to-end training. Visually, the end-to-end training reconstruction yields the closest estimate to the true activity. This is also evident through the line profiles (subfigure (m) and (n)) across the tumor and the liver.

Table III reports the mean activity error (MAE) and NRMSE for lesions and organs across all testing phantoms. Similar to the qualitative assessment (Fig. 9), the end-to-end training also produced lower errors consistently across all testing lesions and organs. For instance, compared to sequential training/gradient truncation, the end-to-end training relatively reduced MAE on average by 8.7%/7.2%, 18.5%/11.0% and 24.7%/16.1% for lesion, healthy liver and lung, respectively. The NRMSE was also relatively reduced by 6.1%/3.8%, 7.2%/4.1% and 6.1%/3.0% for lesion, healthy liver and lung, respectively. All learning-based methods consistently had lower errors than the OSEM method.

### _Results at intermediate iterations_

One potential problem associated with end-to-end training (and gradient truncation) is that the results at intermediate iterations could be unfavorable, because they are not directly trained by the targets [39]. Here, we examined the images at intermediate iterations and did not observe such problems as illustrated in Fig. 10, where images at each iteration gave a

\begin{table}
\begin{tabular}{|c|c|c|c|c|} \hline \multicolumn{5}{|c|}{MAE(\%)} \\ \hline Lesion/Organ & OSEM & Sequential & Truncation & End2end \\ \hline Lesion (6-152mL) & 11.1 \(\pm\) 2.5 & 9.4 \(\pm\) 3.2 & **6.7**\(\pm\) 2.4 & 7.3 \(\pm\) 2.8 \\ \hline Liver & 4.8 \(\pm\) 0.1 & 4.5 \(\pm\) 0.2 & **3.4**\(\pm\) 0.6 & 4.0 \(\pm\) 0.2 \\ \hline Healthy liver & 4.1 \(\pm\) 0.1 & 4.1 \(\pm\) 0.1 & **3.5**\(\pm\) 0.6 & 4.1 \(\pm\) 0.2 \\ \hline Lung & 3.4 \(\pm\) 0.1 & 3.0 \(\pm\) 0.2 & **2.4**\(\pm\) 0.7 & 3.0 \(\pm\) 0.5 \\ \hline Kidney & 5.2 \(\pm\) 0.3 & 4.3 \(\pm\) 0.1 & 2.6 \(\pm\) 0.1 & **2.3**\(\pm\) 0.2 \\ \hline Spleen & 0.8 \(\pm\) 0.2 & **0.6**\(\pm\) 0.1 & 1.3 \(\pm\) 0.6 & 1.2 \(\pm\) 0.4 \\ \hline \multicolumn{5}{|c|}{NRMSE(\%)} \\ \hline Lesion/Organ & OSEM & Sequential & Truncation & End2end \\ \hline Lesion (6-152mL) & 16.1 \(\pm\) 2.2 & 14.9 \(\pm\) 2.4 & 14.3 \(\pm\) 1.7 & **14.2**\(\pm\) 2.1 \\ \hline Liver & 15.9 \(\pm\) 0.2 & **15.3**\(\pm\) 0.1 & 15.5 \(\pm\) 0.6 & **15.3**\(\pm\) 0.1 \\ \hline Healthy liver & 16.8 \(\pm\) 0.1 & **16.6**\(\pm\) 0.1 & 17.3 \(\pm\) 0.5 & 17.1 \(\pm\) 0.3 \\ \hline Lung & 22.3 \(\pm\) 0.3 & 22.1 \(\pm\) 0.4 & 22.0 \(\pm\) 0.4 & **21.9**\(\pm\) 0.5 \\ \hline Kidney & 17.4 \(\pm\) 0.1 & 16.8 \(\pm\) 0.1 & 16.4 \(\pm\) 0.3 & **16.3**\(\pm\) 0.5 \\ \hline Spleen & 13.5 \(\pm\) 0.2 & 12.4 \(\pm\) 0.3 & **12.3**\(\pm\) 0.7 & **12.3**\(\pm\) 0.5 \\ \hline \end{tabular}
\end{table} TABLE II: The average(\(\pm\)standard deviation) MAE(%) and NRMSE(%) across 3 noise realizations of \({}^{177}\mathrm{Lu}\) VP phantoms.

Fig. 8: Qualitative comparison of different training methods and OSEM tested on \({}^{177}\mathrm{Lu}\) VP phantoms. Subfigure (g) and (h) correspond to line profiles marked in (a).

\begin{table}
\begin{tabular}{|c|c|c|c|c|} \hline \multicolumn{5}{|c|}{MAE(\%)} \\ \hline Lesion/Organ & OSEM & Sequential & Truncation & End2end \\ \hline Lesion (1 67mL) & 12.5 \(\pm\) 0.6 & 6.7 \(\pm\) 1.8 & 2.8 \(\pm\) 0.9 & **2.1**\(\pm\) 1.1 \\ \hline Lesion 2 (10mL) & 20.2 \(\pm\) 0.9 & 11.5 \(\pm\) 4.1 & 10.8 \(\pm\) 0.9 & **9.7**\(\pm\) 1.1 \\ \hline Lesion 3 (9mL) & 25.6 \(\pm\) 0.6 & 18.8 \(\pm\) 0.4 & 15.2 \(\pm\) 0.9 & **12.8**\(\pm\) 1.0 \\ \hline Lesion 4 (5mL) & 43.0 \(\pm\) 0.6 & 40.0 \(\pm\) 1.2 & 38.8 \(\pm\) 0.8 & **38.7**\(\pm\) 0.7 \\ \hline Liver & 6.4 \(\pm\) 0.7 & 6.2 \(\pm\) 1.5 & 4.6 \(\pm\) 1.1 & **3.7**\(\pm\) 1.2 \\ \hline Lung & 24.4 \(\pm\) 0.7 & 2.2 \(\pm\) 0.4 & **0.7**\(\pm\) 0.6 & 0.9 \(\pm\) 0.5 \\ \hline Spleen & 14.2 \(\pm\) 0.9 & 12.6 \(\pm\) 2.4 & **14.8**\(\pm\) 0.7 & 9.3 \(\pm\) 1.5 \\ \hline Kidney & 15.9 \(\pm\) 1.0 & 15.1 \(\pm\) 1.2 & 14.4 \(\pm\) 1.4 & **13.6**\(\pm\) 1.6 \\ \hline \multicolumn{5}{|c|}{NRMSE(\%)} \\ \hline Lesion/Organ & OSEM & OSEM & Sequential & Truncation & End2end \\ \hline Lesion 1 (67mL) & 27.3 \(\pm\) 0.3 & 21.7 \(\pm\) 1.3 & 18.9 \(\pm\) 0.6 & **18.3**\(\pm\) 0.6 \\ \hline Lesion 2 (10mL) & 26.8 \(\pm\) 0.6 & 19.2 \(\pm\) 2.2 & 16.4 \(\pm\) 0.4 & **16.3**\(\pm\) 0.8 \\ \hline Lesion 3 (9mL) & 28.4 \(\pm\) 0.4 & 22.8 \(\pm\) 0.8 & 18.3 \(\pm\) 0.7 & **16.3**\(\pm\) 0.7 \\ \hline Lesion 4 (5mL) & 43.5 \(\pm\) 0.5 & 41.1 \(\pm\) 1.3 & **40.0**\(\pm\) 0.7 & 40.2 \(\pm\) 0.6 \\ \hline Liver & 28.5 \(\pm\) 0.1 & 25.0 \(\pm\) 0.8 & **24.3**\(\pm\) 0.3 & 24.5 \(\pm\) 0.3 \\ \hline Lung & 32.1 \(\pm\) 0.1 & 31.2 \(\pm\) 1.1 & **29.5**\(\pm\) 0.3 & 30.4 \(\pm\) 0.4 \\ \hline Spleen & 25.7 \(\pm\) 0.3 & 22.8 \(\pm\) 1.1 & 20.4 \(\pm\) 0.4 & **19.9**\(\pm\) 0.6 \\ \hline Kidney & 40.8 \(\pm\) 0.3 & 39.7 \(\pm\) 0.4 & 39.7 \(\pm\) 0.2 & **39.2**\(\pm\) 0.3 \\ \hline \end{tabular}
\end{table} TABLE I: The average(\(\pm\)standard deviation) MAE(%) and NRMSE(%) across 3 noise realizations of \({}^{177}\mathrm{Lu}\) XCAT phantoms.
fairly accurate estimate to the true activity. Perhaps under the shallow-network setting (e.g., 3 layers used here, with only 3 outer iterations), the network for each iteration was less likely to overfit the training data. Another reason could be due to the non-shared weights setting so that the network could learn suitable weights for each iteration.

## IV Discussion

Training end-to-end CNN-based iterative algorithms for SPECT image reconstruction requires memory efficient forward-backward projectors so that backpropagation can be less computationally expensive. This work implemented a new SPECT projector using Julia that is an open-source, high performance and cross-platform language. With comparisons between Monte Carlo (MC) and a Matlab-based projector, we verified the accuracy, speed and memory-efficiency of our Julia projector. These favorable properties support efficient backpropagation when training end-to-end unrolled iterative reconstruction algorithms. Most modern DL algorithms process multiple data batches in parallel, so memory efficiency is of great importance for efficient training and testing neural networks. To that extent, our Julia projector is much more suitable than the Matlab-based projector.

Fig. 10: Visualization of intermediate iteration results of different training methods. Subfigure (d)-(f): sequential training; (g)-(i): gradient truncation; (j)-(l): end-to-end training.

Fig. 9: Qualitative comparison of different training methods and OSEM tested on \({}^{90}\)Y VP phantoms. Subfigure (a)-(f) and (g)-(l) show two slices from two testing phantoms. Subfigure (m) and (n) correspond to line profiles in (a) and (g), respectively.

\begin{table}
\begin{tabular}{|c|c|c|c|c|c|} \hline \multicolumn{6}{|c|}{MAE(\%)} \\ \hline Lesion/Organ & OSEM & Sequential & Truncation & End2end \\ \hline Lesion (3-356mL) & 32.5 \(\pm\) 1.3 & 25.3 \(\pm\) 1.3 & 24.9 \(\pm\) 1.0 & **23.1**\(\pm\) 1.8 \\ \hline Liver & 25.0 \(\pm\) 0.1 & 18.7 \(\pm\) 0.1 & 17.8 \(\pm\) 1.3 & **15.6**\(\pm\) 3.6 \\ \hline Healthy liver & 25.1 \(\pm\) 0.2 & 23.8 \(\pm\) 0.5 & 21.8 \(\pm\) 1.2 & **19.4**\(\pm\) 3.1 \\ \hline Lung & 88.4 \(\pm\) 2.1 & 64.9 \(\pm\) 1.6 & 58.3 \(\pm\) 6.6 & **48.9**\(\pm\) 8.4 \\ \hline \multicolumn{6}{|c|}{NRMSE(\%)} \\ \hline Lesion/Organ & OSEM & Sequential & Truncation & End2end \\ \hline Lesion (3-356mL) & 35.3 \(\pm\) 1.5 & 29.6 \(\pm\) 1.4 & 28.9 \(\pm\) 1.1 & **27.8**\(\pm\) 1.2 \\ \hline Liver & 29.9 \(\pm\) 0.4 & 22.7 \(\pm\) 0.1 & 22.1 \(\pm\) 0.9 & **21.2**\(\pm\) 1.5 \\ \hline Healthy liver & 31.6 \(\pm\) 0.4 & 27.9 \(\pm\) 0.3 & 27.0 \(\pm\) 0.9 & **25.9**\(\pm\) 2.0 \\ \hline Lung & 62.4 \(\pm\) 1.3 & 59.2 \(\pm\) 1.1 & 57.3 \(\pm\) 3.0 & **55.6**\(\pm\) 4.6 \\ \hline \end{tabular}
\end{table} TABLE III: The average(\(\pm\)standard deviation) MAE(%) and NRMSE(%) across 3 noise realizations of \({}^{90}\)Y VP phantoms.

Fig. 10: Visualization of intermediate iteration results of different training methods. Subfigure (d)-(f): sequential training; (g)-(i): gradient truncation; (j)-(l): end-to-end training.
We used the CNN-regularized EM algorithm as an example to test end-to-end training and other training methods on different datasets including \({}^{177}\mathrm{Lu}\) XCAT phantoms, \({}^{177}\mathrm{Lu}\) and \({}^{90}\mathrm{Y}\) virtual patient phantoms. Simulation results demonstrated that end-to-end training improved reconstruction quality on these datasets. For example, end-to-end training improved the MAE of lesion/liver in \({}^{90}\mathrm{Y}\) phantoms by 8.7%/16.6% and 7.2%/12.4% compared to sequential training and gradient truncation. This improvement could be attributed to the correct gradient was used in backpropagation. Although the end-to-end training yielded the lowest reconstruction error on both \({}^{177}\mathrm{Lu}\) XCAT phantoms and \({}^{90}\mathrm{Y}\) VP phantoms, the reconstruction errors on \({}^{177}\mathrm{Lu}\) VP phantoms were comparable with the gradient truncation. This could be due to the choice of loss functions and CNN architectures in the EM algorithm, which we will explore in the future. Also we noticed that the recovery of the nonuniform activity in VP phantoms was generally higher than activity for the XCAT phantom (MAE reported in Table I and Table II) because the assigned "true" activities at the boundaries of organs did not drop sharply, and instead, were blurred out. And therefore the OSEM algorithm was fairly competitive as reported in Table II; in \({}^{90}\mathrm{Y}\) VP results, the OSEM performed worse than learning-based methods, which could be attributed to the high downscatter associated with \({}^{90}\mathrm{Y}\) SPECT due to the continuous bremsstrahlung energy spectrum. We found all learning methods did not work very well for small tumors (e.g., 5mL), potentially due to the worse PV effect. Reducing PV effects in SPECT images has been studied extensively [40, 41]. Recently, Xie _et al._[42] trained a deep neural network to learn the mapping between PV-corrected and non-corrected images. Incorporating their network into our reconstruction model using transfer learning is an interesting future direction.

Although promising results were shown in previous sections, this work has several limitations. First, we did not test numerous hyperparameters and CNN architectures, nor with a wide variety of phantoms and patients for different radionuclides therapies. Secondly, our experiments used OSEM images as warm start to the CNN-regularized EM algorithm, where the OSEM itself was initialized with a uniform image. We did not investigate using other images such as uniform images as the start of the EM algorithm. Using a uniform image to initialize the network would likely require far more network iterations which would be very expensive computationally and therefore impractical. Additionally, this paper used fixed regularization parameter (\(\beta\) in (5)) rather than declaring \(\beta\) as a trainable parameter. We compared different methods for backpropagation, which requires using the same cost function (4) for a fair comparison. If one set \(\beta\) as a trainable parameter, then different methods could learn different \(\beta\) values, leading to different cost functions. However, the investigation of trainable \(\beta\) values is an interesting future work. Another limitation is that we did not investigate more advanced parallel computing methods such as distributed computing using multiple computers to further accelerate our Julia implementation of SPECT forward-backward projector. Such acceleration is feasible using existing Julia packages if needed. The compute times reported in Fig. 4 show that the method needs a few seconds per 128 projection views using 8 threads, which is already feasible for scientific investigation.

We also found there exists a trade-off between computational cost and reconstruction accuracy for different training methods. End-to-end training yielded reconstruction results with the lowest MAE and NRMSE because the correct gradient was used during backpropagation. Sequential training yielded worse results, but it was significantly faster and more memory efficient than the end-to-end training method. It is notably faster because it splits the whole training process and trains each of neural networks separately, and its backpropagation does not involve terms associated with the MLEM algorithm, so sequential training is actually equivalent to training that neural network alone without considering MLEM. Sequential training also used much less memory because the training was performed iteration by iteration, one network by one network, and hence the memory limitation did not depend on the number of unrolled iterations in the MLEM algorithm.

## V Conclusion

This paper presents a Julia implementation of backpropagatable SPECT forward-backward projector that is accurate, fast and memory-efficient compared to Monte Carlo (MC) and a previously developed analytical Matlab-based projector. Simulation results based on \({}^{177}\mathrm{Lu}\) XCAT phantoms, \({}^{90}\mathrm{Y}\) and \({}^{177}\mathrm{Lu}\) virtual patient (VP) phantoms demonstrate that: 1) End-to-end training yielded reconstruction images with the lowest MAE and NRMSE when tested on XCAT phantoms and \({}^{90}\mathrm{Y}\) VP phantoms, compared to other training methods (such as sequential training and gradient truncation) and OSEM. 2) For \({}^{177}\mathrm{Lu}\) VP phantoms, end-to-end training method yielded better results than sequential training and OSEM; but was rather comparable with gradient truncation. We also found there exists a trade-off between computational cost and reconstruction accuracy in different training methods (e.g., end-to-end training and sequential training). These results indicate that end-to-end training, which is feasible with our developed Julia projector, is worth investigating for SPECT reconstruction.

## Acknowledgement

All authors declare that they have no known conflicts of interest in terms of competing financial interests or personal relationships that could have an influence or are relevant to the work reported in this paper.
