A generalization of the maximum likelihood expectation maximization (MLEM) method: Masked-MLEM

Objective. In our previous work on image reconstruction for single-layer collimatorless scintigraphy, we developed the min–min weighted robust least squares (WRLS) optimization algorithm to address the challenge of reconstructing images when both the system matrix and the projection data are uncertain. Whereas the WRLS algorithm has been successful in two-dimensional (2D) reconstruction, expanding it to three-dimensional (3D) reconstruction is difficult since the WRLS optimization problem is neither smooth nor strongly-convex. To overcome these difficulties and achieve robust image reconstruction in the presence of system uncertainties and projection noise, we propose a generalized iterative method based on the maximum likelihood expectation maximization (MLEM) algorithm, hereinafter referred to as the Masked-MLEM algorithm. Approach. In the Masked-MLEM algorithm, only selected subsets (‘masks’) from the system matrix and the projection contribute to the image update to satisfy the constraints imposed by the system uncertainties. We validate the Masked-MLEM algorithm and compare it to the standard MLEM algorithm using experimental data obtained from both collimated and uncollimated imaging instruments, including parallel-hole collimated SPECT, 2D collimatorless scintigraphy, and 3D collimatorless tomography. Additionally, we conduct comprehensive Monte Carlo simulations for 3D collimatorless tomography to further validate the effectiveness of the Masked-MLEM algorithm in handling different levels of system uncertainties. Main results. The Masked-MLEM and standard MLEM reconstructions are similar in cases with negligible system uncertainties, whereas the Masked-MLEM algorithm outperforms the standard MLEM algorithm when the system matrix is an approximation. Importantly, the Masked-MLEM algorithm ensures reliable image reconstruction across varying levels of system uncertainties. Significance. With a good choice of system uncertainty and without requiring accurate knowledge of the actual system matrix, the Masked-MLEM algorithm yields more robust image reconstruction than the standard MLEM algorithm, effectively reducing the likelihood of erroneously reconstructing higher activities in regions without radioactive sources.


Introduction
Generating high-quality images is an essential procedure for obtaining structural or functional information in many imaging scenarios, such as computed tomography (CT), single photon emission computed tomography (SPECT), and positron emission tomography (PET) (Defrise & Gullberg 2006).
In emission and transmission tomography, image reconstruction is performed using either analytical or iterative methods (Hsieh et al. 2013, Zeng 2001).Iterative reconstruction has several advantages over analytical reconstruction.These include better image quality, flexibility for reconstruction using incomplete data or irregular sampling, and the potential to reduce patient dose, as imaging physics and projection noise can be integrated into the iterative process (Hsieh et al. 2013, Zeng 2001, Beister et al. 2012, Yan et al. 2013).A popular iterative reconstruction algorithm aimed at finding the most likely image given the observed data is maximum likelihood expectation maximization (MLEM) (Shepp & Vardi 1982, Lange et al. 1984), or its variant, ordered subsets expectation maximization (OSEM) (Hudson & Larkin 1994).MLEM-based algorithms are particularly suitable for imaging problems with Poisson noise characteristics and enforce the non-negativity of reconstructed images, which are desirable properties in PET or SPECT reconstructions (Lange et al. 1984, Zeng 2001).
Mathematically, image reconstruction addresses a linear ill-posed inverse problem to recover an unknown image from incomplete and noisy measurements (Hsieh et al. 2013, Yan et al. 2013).Typically, this process is structured as a convex optimization problem using Tikhonov regularization, where the data fidelity term is a least-squares (LS) optimization problem and the additive regularization term imposes some prior knowledge on the image to be reconstructed (Tikhonov 1963, Panin et al. 1998, Yan et al. 2013).With a growing research interest in image priors, various regularization terms including total variation (TV), nonlocal means (NLM), and deep-learning-based priors, have been successfully integrated into iterative algorithms whereas assuming a known and accurate system response in the data fidelity term (Panin et al. 1998, Sawatzky et al. 2008, Chen et al. 2016, Gong et al. 2018).In reality, the system matrix used for image reconstruction suffers from uncertainties due to complicated physical effects (Liu et al. 2012).Previous studies has investigated the impact of noise in the system matrix and shown that these uncertainties in the system matrix can propagate to the reconstructed images (Qi & Huesman 2004, Rafecas et al. 2004).Therefore, in addition to incorporating projection statistical models and image priors, incorporating constraints on system uncertainties is expected to further improve the fidelity of the reconstructed images and yield robust reconstructions.
Due to the huge dimensionality and ill-posedness of the system matrix, considering the system uncertainties explicitly during the reconstruction process is difficult and computationally expensive.Because of this, only a few studies have investigated system uncertainties in image reconstruction (Qi & Huesman 2004, Rafecas et al. 2004, Liu et al. 2012, Kucharczak et al. 2018, Zheng et al. 2020, Lunz et al. 2021).Instead of directly addressing the system uncertainties, previous studies have adopted alternative approaches to account for system uncertainties indirectly (Rico et al. 2009, Rico & Strauss 2010, Kucharczak et al. 2018, Lunz et al. 2021).By transforming the system uncertainties to the projection confident intervals, the non-additive interval-based EM algorithm has been developed to reconstruct image confidence intervals based on the projection confident intervals, providing an estimate of the uncertainties in the reconstructed images (Rico et al. 2009, Rico & Strauss 2010, Kucharczak et al. 2018).An alternative approach to transforming system uncertainties into projection uncertainties is by training a neural network to learn a correction for the projection data, which can effectively provide a correction for the system matrix (Lunz et al. 2021).Deep-learningbased approaches take advantage of the learning capabilities of neural networks to model the complex system uncertainties and the projection uncertainties, but they require a large amount of training data to effectively learn and generalize.
In our previous work, we developed a min-min weighted robust least squares (WRLS) algorithm to directly address the uncertainties present in both the projection data and the system matrix (Zheng et al. 2020).We demonstrated the success of the WRLS algorithm in two-dimensional (2D) collimatorless scintigraphy using the CVXPY toolkit for iterative reconstruction (Zheng et al. 2020, Diamond & Boyd 2016).However, the WRLS optimization problem is neither smooth nor strongly-convex (Boyd et al. 2004).When applying some traditional gradient-based algorithms to solve the WRLS optimization problem, the non-smoothness will introduce difficulties in differentiating the objective function, and the lack of strong convexity will lead to multiple local optima, making it difficult to find the global optimum and resulting in slow convergence.As image dimensionality increases, the optimization landscape becomes more complex, and finding the global optimum is more challenging.For example, in three-dimensional (3D) tomographic reconstruction, the CVXPY toolkit is computationally memory-intensive and time-consuming to handle such a large-scale WRLS optimization problem even on a dedicated server.Without an efficient iterative solver, the applications of the WRLS algorithm will be limited.
In this study, we propose the Masked-MLEM algorithm as a generalized iterative MLEM approach to achieve robust image reconstruction in the presence of system uncertainties and projection noise.We use the Masked-MLEM algorithm to reconstruct images from experimental data in three types of imaging scenarios: parallel-hole collimated SPECT, 2D collimatorless scintigraphy, and 3D collimatorless tomography.The images reconstructed using the Masked-MLEM algorithm are compared with those using the standard MLEM algorithm to examine the robustness of the Masked-MLEM algorithm.

Masked-MLEM
The Masked-MLEM algorithm is derived based on the WRLS algorithm, where "Masked" refers to selecting subsets in both the system matrix and the projection in order to satisfy the constraints on the system uncertainties, and using the selected subsets for the image update.In our previous work, the WRLS algorithm took into account the uncertainties in both the projection data y and the system matrix A by assuming a box uncertainty set such that element-wisely A lb ⪯ A ⪯ A ub , where A lb and A ub are the lower and upper bounds of A respectively (Zheng et al. 2020).The box uncertainty set of A indicates the constraints such that element-wisely A lb x ⪯ y ⪯ A ub x, where x is the image to be reconstructed.By relaxing the strict equality constraint Ax = y to these additional soft constraints based on the system uncertainties, the feasible reconstructed image space can be expanded and the reconstruction process has more flexibility to explore a larger solution space when the exact solution for Ax = y may not be feasible.
We first reformulate the WRLS optimization problem in (Zheng et al. 2020) as where Σ is a diagonal matrix representing an additive Gaussian noise in y, and µ is a slack variable introduced to make the objective function convex (Boyd et al. 2004).Instead of applying any optimization techniques like the CVXPY toolkit to directly solve the WRLS optimization problem, we decompose it into smaller subproblems and address the challenges posed by non-smoothness and lack of strong convexity within these subproblems.Looking at Eq. 1, if we can find some x * satisfying A u x * − y u ⪯ 0, then the optimal solutions for Eq. 1 will be µ * = 0 and x * , i.e. x * is the final reconstructed image.Otherwise, when some entries in the vector A u x − y u are positive, the objective function in Eq. 1 is to minimize the ℓ ∞ norm of the positive entries in A u x − y u (i.e.∀i = 1, . . ., m, max{(a lb ) T i x − y i , y i − (a ub ) T i x, 0}, see (Zheng et al. 2020)).Since the ℓ ∞ norm optimization is neither smooth nor strongly-convex, solving Eq. 1 in large-scale image reconstruction problems is computationally intractable and prone to convergence issues.
To overcome these challenges, one optimization approach is to relax the ℓ ∞ norm with a smooth convex surrogate, such as the ℓ 2 norm (i.e.∀i = 1, . . ., m, max{(a lb ) T i x − y i , 0} 2 + max{y i − (a ub ) T i x, 0} 2 ).This relaxation extends the objective function in Eq. 1 to a larger and smooth function.After relaxation, instead of optimizing the original objective function in Eq. 1, we optimize the upper bound of Eq. 1 given by min where M = diag{A u x − y u ≻ 0} is a "mask" applied to A u and y u to select the positive subsets in A u x − y u and nullify the non-positive subsets (see proof below).Given that some elements in M A u and M y u may be negative, which mathematically makes the constraints in Eq. 1 convex but is physically trivial in Eq. 2, we element-wisely take the absolute values of M A u and M y u as M |A u | and M |y u | to make the selected subsets in A u and y u non-negative.Proof.
Eq. 1 See (Zheng et al. 2020) The form of Eq. 2 is a generalized data fidelity term that takes into account the system uncertainty A u and the additive Gaussian noise in y.It can be solved efficiently by a number of iterative algorithms used for LS optimization, where the modification is to apply a mask to select the positive subsets in A u x − y u for the image update at each iteration, as if the actual system matrix is M |A u | and the projection is M |y u | for LS optimization.Consequently, the optimal solution to Eq. 2 will be different from Eq. 1.We have tested that the difference between the reconstructed images using Eq. 1 and Eq. 2 can be negligible considering the intrinsic image resolution and projection noise.
Inspired by Eq. 2 where we assume the actual system matrix is M |A u | and the projection is M |y u |, when we assume Poisson noise in y, we derive the Masked-MLEM algorithm based on the MLEM algorithm as shown in Algorithm 1 (see Appendix A for derivation details).The Masked-MLEM algorithm is a generalization of the standard MLEM algorithm, i.e.MLEM reduces to Masked-MLEM (Dasgupta et al. 2008).The validity of this reduction is straightforward: if we know how to solve Masked-MLEM, by setting A lb = A = A ub and running Masked-MLEM, we will get the solution of MLEM.
Algorithm 1 Masked-MLEM (N out , N in ) 1: Initialize x 0 to be a positive vector.2: for s = 1, . . ., N out do ▷ We find the optimal solution of Eq. 2. 5: x k ← x s 9: for k = 1, . . ., N in do 10: ▷ | • | is the element-wise absolute value.12: x s ← x k 13: return x s In Algorithm 1, for each iteration in the outer loop N out , we select a mask M such that only the selected subsets in A u and y u will contribute to the update of the image x s .To speed up reconstruction, we can use the same M to update the image for N in iterations, rather than updating M whenever the image is updated, since updating M may be computationally expensive for huge matrices and the selected masks at subsequent iterations may be similar.We have checked that the Masked-MLEM algorithm converges when using appropriate A u , N out , and N in by monitoring the loss curve during the reconstruction process.The reconstruction loss is defined as , where x k+1 and x k represent the reconstructed images at iteration k + 1 and k respectively, and serves as the metric for convergence throughout the iterations.For the convergence of the Masked-MLEM algorithm, N in should be determined empirically through experimentation and should take a small value (usually less than 10 depending on reconstruction problems), otherwise the mask M will be outdated and not suitable for future image updates.
The Masked-MLEM algorithm is also compatible with ordered subsets (OS), i.e.Masked-OSEM.If we update the image using OS from the mask M , it does not guarantee convergence since M varies with the image update.For the purpose of convergence, we use the OS to calculate masks, followed by the Masked-MLEM algorithm.The detailed Masked-OSEM algorithm is shown in Algorithm 2.

System Uncertainties
To apply the Masked-MLEM algorithm, we design the system matrix upper and lower bounds for both collimated and uncollimated imaging modalities.For the parallel-hole collimated SPECT, we use the standard ray-tracing algorithm to generate the system matrix A S , which does not account for unwanted photon transmission through the collimator (Huesman et al. 1977, Siddon 1985, Gullberg et al. 1989).Because of the Algorithm 2 Masked-OSEM (N out , N in , N subset ) 1: Initialize x 0 to be a positive vector.2: for s = 1, . . ., N out do 3: ▷ We find the optimal solution of Eq. 2.
13: return x s 14: ▷ Perform the EM step.17: x s ← x k 19: return x s imperfect collimator performance of SPECT systems in reality, A S is too idealistic and we assume that the point spread function (PSF) of the system response is a Gaussian function rather than an ideal δ function.We blur each PSF in A S with a 5 by 5 Gaussian kernel and add the Gaussian smear to A S as the penetration component, resulting in a new system matrix A G .Thus the system matrix lower and upper bounds for the parallel-hole collimated SPECT are calculated by where α is a hyperparameter adjusting the magnitude of the Gaussian smear in A S .
For collimatorless imaging, we follow the same approach as in (Zheng et al. 2020) to calculate the system matrix lower and upper bounds, and add more hyperparameters to make the bounds more flexible (see Fig. 1).Without any collimation, the system matrix A for collimatorless gamma cameras is theoretically calculated based on solid angles.The element a ij in A, which represents the probability that photons from the image voxel j reach the detector pixel i, is calculated by where p is the detector pixel size, R ij and r ij are the distance and the normal distance between the image voxel j and the detector pixel i.However, A is a weak approximation to the real system matrix for collimatorless imaging, as it only takes into account the solid angle.The detailed procedures to obtain the system matrix lower and upper bounds based on A are shown in Algorithm 3, where ϵ, η, θ, and ζ are hyperparameters used to tune A lb and A ub .The matrices Â, A max , A ′ , Ā, and W in Algorithm 3 are illustrated in Fig. 1 and follow the same definitions in (Zheng et al. 2020).
The system matrix A for collimatorless imaging may not be accurate and we assume the actual system matrix is in the neighborhood of A (shown as the blue cloud).
The system matrix Â corresponds to the ideal scenario where a perfect parallel-hole collimator is used.The system matrix A max defines the upper limit that the current system can achieve, i.e. it represents the scenario where every image voxel has an equally maximum probability of reaching every detector pixel.The system matrix A ′ determines the characteristics of the PSF and the range Ā ± W is desired to be close to the actual system matrix.
The main idea of Algorithm 3 is to change the characteristics of the PSF by assuming that the largest response in A is accurate, while considering that the uncertainties of the remaining a ij s are linearly proportional to their respective magnitudes.The hyperparameter ϵ determines the proportion of the largest response in A that is assumed to be accurate and retained.The hyperparameter η determines the location of A ′ and the hyperparameter θ determines the location of Ā in Fig. 1.Both η and θ contribute to shaping the characteristics of the PSF.Additionally, the hyperparameter ζ controls the size of the uncertainty set.These hyperparameters jointly define a system uncertainty set and can be individually adjusted.The value of ϵ is determinate when we have knowledge of the ideal system response Â using a perfect parallel-hole collimator.If A is already close to the actual system matrix, a Algorithm 3 System uncertainty (ϵ, η, θ, ζ) smaller θ can be set, and vice versa.Similarly, if the uncertainty set is expected to be small, a smaller ζ can be selected accordingly.Among these hyperparameters, η is the most sensitive and typically needs to be set close to zero.Finding the optimal value for η often involves trying out various smaller values through experimentation.We have tested in our previous study that using these hyperparameters to adjust the system uncertainties works for collimatorless image reconstruction (Zheng et al. 2020).
For reconstruction in collimatorless imaging, the system matrix for the standard MLEM algorithm is A, and the system matrix uncertainty A u for the Masked-MLEM algorithm is obtained by manually adjusting ϵ, η, θ, and ζ in Algorithm 3.Both A and A u set zero responses for the elements corresponding to the dead detector pixels or strips in collimatorless imaging.

Phantom Studies
We obtained data from three imaging instruments utilizing both collimated and uncollimated imaging modalities to validate the Masked-MLEM algorithm: 1) Parallelhole collimated SPECT; 2) 2D collimatorless scintigraphy; and 3) 3D collimatorless tomography.Due to limited access to these imaging instruments, we used three unique phantoms filled with both different radionuclides and activities in different imaging modalities.The focus of this study is to compare the standard MLEM and the Masked-MLEM approaches, not the imaging instruments themselves.
2.3.1.SPECT For SPECT imaging, a NEMA SPECT triple line source phantom (NEMA SPECT Triple Line Source Phantom n.d.), filled with about 5 mCi/mL 99m Tc pertechnetate by following a standard filling procedure, was scanned by the GE Infinia Hawkeye 4 SPECT/CT for 30 s per projection (GoldSeal Infinia Hawkeye 4 n.d.).The line sources are positioned in a right triangle configuration with one at the center and the other two at a distance of 75 mm from the center.Each line source has a diameter of 1 mm and a useful height of 184 mm.The SPECT scanner uses a NaI detector consisting of 128 by 128 detector pixels and a parallel-hole collimator (LE-HR).To obtain the projection data, 10% energy window was used for the 140 keV photopeak of 99m Tc and no scatter correction was applied.We used the projection data from 60 projection angles separated by 6 • for image reconstruction.The system matrix for the standard MLEM algorithm is A S , and the system matrix uncertainty A u for the Masked-MLEM algorithm is obtained by adjusting α in Eq. 3. Due to the huge dimensionality of A u , we used 6 subsets for the OSEM and the Masked-OSEM reconstructions.We ran both the OSEM algorithm and Masked-OSEM algorithm for 50 iterations until convergence (the reconstruction loss was less than 0.004) with N out = 50 and N in = 1 for the Masked-OSEM reconstruction.

Collimatorless Tomography
For 3D collimatorless tomography, a fillable micro hollow sphere phantom filled with 0.53 kBq/µL of 225 Ac aqueous solution in its three different-sized spheres, was rotated at 8 angles separated by 45 • and scanned by a 3D position-sensitive HPGe double-sided strip detector (DSSD) for 30 min per projection (see Fig. 3) (Frame, Barnowski, Gunter, Mihailescu & Vetter 2022, Frame, Bobba, Gunter, Mihailescu, Bidkar, Flavell & Vetter 2022).The detector is part of a Compton camera and consists of 37 by 37 orthogonal stripes with an active pixel size of 2 mm × 2 mm × 15 mm (Frame, Barnowski, Gunter, Mihailescu & Vetter 2022, Frame, Bobba, Gunter, Mihailescu, Bidkar, Flavell & Vetter 2022).The cylinder phantom has a diameter of 40 mm and a height of 82 mm, and was filled with water.The micro hollow spheres, from smallest to largest, have diameters of 4.3 mm, 6.2 mm, and 7.8 mm; therefore, the total activities of 225 Ac were approximately 0.6 mCi, 1.8 mCi, and 3.6 mCi, respectively.The center of the micro hollow sphere phantom to the surface of the detector was about 52.5 mm.The detailed experiment setup for measuring the micro hollow sphere phantom is shown in (Frame, Bobba, Gunter, Mihailescu, Bidkar, Flavell & Vetter 2022).Given that the 3D position-sensitive detector could record the depth of interactions, we summed the counts of the single interaction events as the projection data, whose interaction depths were within 3 mm from the surface of the first detector of Compton cameras, and energies were in the 18% energy window of the characteristic X-ray peak centering at 87 keV, the 1.1% energy window of the 218 keV photopeak of 221 Fr, and the 0.6% energy window of the 440 keV photopeak of 213 Bi.With the constraint on the interaction depth, the system matrix A is applicable for projections coming from a mixture of energies.We ran both the MLEM algorithm and Masked-MLEM algorithm for 500 iterations until convergence (the reconstruction loss was less than 0.002) with N out = 500 and N in = 1 for the Masked-MLEM reconstruction.
For all three experiment studies, we perform the 3D image reconstruction and apply a 3 by 3 by 3 median filter to the 3D reconstructed images for denoising.In some cases, we collapse the 3D reconstructed images to the corresponding 2D images for comparison between the standard MLEM algorithm and the Masked-MLEM algorithm.For example, in collimatorless scintigraphy, we sum the 3D reconstructed image along the normal direction of the detector surface to obtain a 2D coronal image of the mouse phantom, since the depth reconstruction of the mouse phantom using only one projection is inaccurate.In this study, we do not consider image regularizations and compare the reconstructed images using the standard MLEM algorithm and the Masked-MLEM algorithm.To evaluate the quality of reconstructions, we normalize both the reconstructed images and the ground truth images by their respective maximum values, and utilize three figures of merit (FOMs) that are not influenced by scaling factors: normalized root mean squared error (NRMSE), peak signal-to-noise ratio (PSNR), and structural similarity (SSIM).NRMSE and PSNR quantitatively measure the accuracy between the reconstructed images and the ground truth images, while SSIM is used to assess the similarity between them.

Results
In In collimatorless scintigraphy, the projection data coming from the characteristic X-rays and the 218 keV photopeak of 221 Fr are shown in the first column of Fig. 5.We use ϵ = 0.04, η = 0.02 for X-rays projection and η = 0.04 for 221 Fr projection, θ = 0.2, and ζ = 0.5 in Algorithm 3 to obtain A u for the Masked-MLEM reconstruction.We reconstruct the 3D images of the mouse phantom using the standard MLEM algorithm and the Masked-MLEM algorithm based on the single-angle projection data from X-rays and 221 Fr, respectively.The collapsed 2D coronal images overlaid on a coronal slice of In collimatorless tomography, the projection data of the micro hollow sphere phantom filled with 225 Ac at 8 angles are shown in Fig. 6.We use ϵ = 0.04, η = 0.5, θ = 0.2, and ζ = 0.5 in Algorithm 3 to obtain A u for the Masked-MLEM reconstruction.The representative transverse, coronal and sagittal slices from the center of the reconstructed spheres are shown in Fig. 7.The FOMs of all reconstructed images in Fig. 4, 5, and 7 are shown in Table 1.

Discussion
In this study, the Masked-MLEM algorithm is developed and validated on both collimated and uncollimated imaging instruments and compared with the standard MLEM algorithm.Due to limited access, the use of different phantoms, radionuclides, and total activities for imaging affects the reconstruction performance of the three imaging instruments, which is one limitation of this study.However, it will not be a major concern of this study as we are comparing the reconstruction results across different algorithms, rather than across different imaging instruments.
In SPECT imaging, there is almost no difference between the reconstruction of the triple line source phantom using the standard OSEM algorithm and the Masked-OSEM algorithm (see Fig. 4 and Table 1).Similar results were obtained in previous studies which selected the worst-case system matrix in its uncertainty set for robust reconstruction (Liu et al. 2012).The reasons may have the following two aspects.First, the system matrix A S for the simple-structured triple line source phantom is already the optimal or suboptimal system matrix in the uncertainty set.Second, our definition of the system matrix upper bound, which adds an adjustable Gaussian smear to A S as  the penetration component, will only make the reconstructed images more restricted, but there may be other better designs of the system uncertainties which require further exploration.
In contrast to SPECT imaging, the Masked-MLEM algorithm outperforms the standard MLEM algorithm in collimatorless imaging as shown in Fig. 5 and 7, and in Table 1.Image resolution is a limiting factor in collimatorless imaging because we do not have any collimation and only rely on the solid angle model, which sacrifices the image resolution for higher sensitivity.Therefore, the three key factors for collimatorless imaging are: 1) Proximity, i.e. how close the phantom is to the detector, which determines the resolution of the reconstructed image and to what extent the system matrix is ill-conditioned; 2) Sensitivity, i.e. how noisy the projection is, which determines the signal-to-noise ratio (SNR) of the reconstruction; 3) Uniformity of source distribution, i.e. whether the source is uniformly distributed throughout the phantom or dominates in certain locations within the phantom, which determines the reconstruction accuracy.
The collimatorless scintigraphy experiment is an example showing good proximity, good sensitivity, and a uniform source distribution.In collimatorless scintigraphy, the Masked-MLEM algorithm successfully reconstructs the positions of the brain, liver, and tumor of the mouse phantom with fewer artifacts and better FOMs compared to the standard MLEM algorithm (see Fig. 5 and Table 1).Compared to the Masked-MLEM reconstructions, the standard MLEM reconstructions distribute more activities to adjacent organs that should not be radioactive, for example, the kidneys of the mouse phantom in Fig. 5. Given that the mouse phantom was placed on the surface of the detector with good proximity and the source organs were distant from each other with similar total activities of 225 Ac (i.e.uniform source distribution), the source organs can be easily distinguished from the projection data in Fig. 5, which makes the image reconstruction easier and improves the reconstruction accuracy.Whether using the MLEM algorithm or the Masked-MLEM algorithm, the images reconstructed from the 221 Fr projection are noisier than those from the X-rays projection in Fig. 5, since the projection data of 221 Fr have fewer counts and thus worse statistics than those of Xrays in Fig. 5. Therefore, it is important to obtain sufficient counts for collimatorless projection to improve the image reconstruction quality.
Good proximity, good sensitivity, and a uniform source distribution are ideal conditions for collimatorless reconstruction; however, in practice, such conditions are not always met.The collimatorless tomography experiment is an example of this.In collimatorless tomography, although the Masked-MLEM algorithm outperforms the standard MLEM algorithm in Fig. 7, neither of them obtains reconstructions with high resolution in Fig. 7, or good FOMs in Table 1.This can be explained by the poor proximity of the micro hollow sphere phantom, and the non-uniform activity distribution among the spheres.In regard to the latter, the total activities in the two largest spheres were much higher than those in the smallest one.Due to this non-uniform distribution and the poor image resolution of collimatorless tomography, only the largest two spheres appear in the reconstruction (see Fig. 7).Furthermore, due to the dead stripes in the detector (see Fig. 6), useful projection information is missing, degrading the image quality.We also notice a shift in position towards the boundaries of the imaging region of interest (ROI) in collimatorless tomography (see Fig. 7), indicating that some collimation should be required to provide more accurate position information for reconstruction.Though reduced reconstruction quality, we experimentally demonstrate the feasibility of performing proximity reconstruction with a Compton camera system, as discussed in (Caravaca et al. 2022).How to improve the image reconstruction quality in collimatorless imaging will be investigated in our future studies.
In all experiment studies, the Masked-MLEM algorithm consistently produces robust image reconstructions by effectively minimizing the likelihood of reconstructing higher activities in regions without radioactive sources (see Fig. 4, 5, and 7).Furthermore, the Masked-MLEM algorithm offers the capability to reconstruct images without accurate knowledge of the actual system matrix.Although it is preferable to have an accurate system matrix, in many cases, we are only able to obtain an approximation to the actual system matrix due to complicated factors such as scattering, attenuation, and nonuniform detector responses.Additionally, computing an accurate system matrix can often be time-consuming or impractical.The power of the Masked-MLEM algorithm lies in leveraging A lb and A ub based on the current system matrix A to avoid the need for time-consuming computations of an optimal system matrix.This allows for efficient image reconstruction without the burdensome computational overhead, regardless of the availability or accuracy of the system matrix.In cases where an accurate system matrix is available, the Masked-MLEM reconstruction is similar to the standard MLEM reconstruction, as observed in SPECT imaging.
However, in situations where the system matrix is an approximation, such as in collimatorless imaging, the Masked-MLEM reconstruction outperforms the standard MLEM reconstruction, which employs the straightforward system matrix A that considers only solid angles.
The success of the Masked-MLEM algorithm depends heavily on the choice of the system uncertainty A u .Although the Masked-MLEM algorithm will choose the best system matrix in its uncertainty set, it does not mean that we can set the system uncertainty set as large as possible.On the other hand, the optimal choice of A u should be as restrictive as possible, including all possible good system matrices and excluding all impossible system matrices.We have tested that η is the most sensitive hyperparameter for tuning A u in collimatorless imaging, and all other hyperparameters can be kept the same for different reconstruction problems.In this study, A u is largely defined from a mathematical point of view, and it will be better to understand some prior knowledge of the system uncertainties to define A u according to the physical effects in real applications.The box uncertainty of the system matrix provides the flexibility to element-wisely design the system uncertainties and the power to include any model for system uncertainties, which has the potential to be widely used in other inverse optimization problems involving system uncertainties.Although loading the lower and upper bounds of the system matrix occupies more storage memory, the computation time of each iteration in the Masked-MLEM algorithm is comparable to that of the standard MLEM algorithm, which is acceptable.

Conclusion
In this study, we develop the Masked-MLEM algorithm as a generalization of the standard MLEM algorithm for robust image reconstruction.The Masked-MLEM algorithm incorporates a box uncertainty set in the system matrix and accounts for Poisson noise in the projection data.We validate the Masked-MLEM algorithm using data from both collimated and uncollimated imaging instruments by carefully designing the corresponding system matrix uncertainties.In SPECT imaging, the Masked-MLEM and standard MLEM reconstructions are similar, and in collimatorless imaging, the Masked-MLEM algorithm outperforms the standard MLEM algorithm.In both collimated and uncollimated imaging, the Masked-MLEM algorithm enhances the robustness of image reconstruction, yielding more reliable results in the presence of system uncertainties, and effectively reducing the likelihood of reconstructing higher activities in regions without radioactive sources.Future studies will focus on designing better system uncertainties for SPECT imaging and proposing new imaging modalities for further improving the reconstructed images in collimatorless imaging.
For 2D collimatorless scintigraphy, a fillable mouse phantom with its brain, liver, and tumor filled with 261 nCi, 261 nCi, and 239 nCi of 225 Ac aqueous solution, respectively, was put on the surface of a CZT detector and measured for 20 min (Fillable Mouse/rat phantom 2022).The experiment setup is shown in Fig.2.The dimensions of the mouse phantom are 97.86 mm in length, 27.48 mm in width, and 22.70 mm in height.The detector is composed of 5 stacked CZT crystals, with each crystal consisting of 16 by 16 pixels.This results in a total of 80 by 16 detector pixels.Each pixel has a size of 1.6 mm × 1.6 mm × 5.0 mm.The bottom of the mouse phantom to the surface of the detector was about 2 mm.The mouse phantom was scanned on only one projection and we used a 20% energy window for the characteristic X-ray peak centering at 90 keV and a 6% energy window for the 218 keV photopeak of 221 Fr, one daughter of 225 Ac, to obtain the projection data for both MLEM reconstruction and Masked-MLEM reconstruction.We ran both the MLEM algorithm and Masked-MLEM algorithm for 200 iterations until convergence (the reconstruction loss was less than 0.002) with N out = 200 and N in = 1 for the Masked-MLEM reconstruction.

Figure 2 .
Figure 2. Experiment setup of the mouse phantom for collimatorless scintigraphy.

Figure 4 .
Figure 4. Reconstruction of the triple line source phantom in SPECT imaging.The top row shows the transverse images of the triple line source phantom obtained by longitudinally summing the 3D images reconstructed by the standard OSEM algorithm and the Masked-OSEM algorithm.The ground truth of the triple line source phantom is reconstructed from the CT image.The bottom row shows the reconstruction profiles of the triple line source phantom in the top row.

Figure 5 .
Figure 5. Reconstruction of the mouse phantom in collimatorless scintigraphy with its brain, liver, and tumor filled with 261 nCi, 261 nCi, and 239 nCi of 225 Ac.The first column shows the projection data coming from the characteristic X-rays and the 218 keV photopeak of 221 Fr where the zero response in certain positions was due to the dead detector pixels.It is noteworthy that the dead detector pixels appeared at the same positions across all five CZT crystals, and were not the result of a multistage acquisition.The second column shows the normalized concentrations of the source organs as the ground truth.The last two columns show the 2D coronal images overlaid on the CT image and reconstructed using the standard MLEM algorithm and the Masked-MLEM algorithm, respectively.

Figure 6 .
Figure 6.Projection of the micro hollow sphere phantom filled with 225 Ac at 8 angles.The projection counts came from the single interaction events near the detector surface (within 3 mm) with energies emitted from the characteristic X-rays, the 218 keV photopeak of 221 Fr, and the 440 keV photopeak of 213 Bi.The zero response was caused by the dead detector strips.

Figure 7 .
Figure 7. Ground truth and reconstruction of the micro hollow sphere phantom filled with 0.6 mCi, 1.8 mCi, and 3.6 mCi of 225 Ac from the smallest to the largest micro hollow spheres.From top to bottom, it shows the transverse, coronal and sagittal slices of the reconstructed spheres.

Table 1 .
FOMs of reconstructed images in Fig.4, 5, and 7 using the standard MLEM algorithm and the Masked-MLEM algorithm.