A LIST-MODE OSEM-BASED ATTENUATION AND SCATTER COMPENSATION METHOD FOR SPECT.

Reliable attenuation and scatter compensation (ASC) is a pre-requisite for quantification and beneficial for visual interpretation tasks in SPECT. In this paper, we develop a reconstruction method that uses the entire SPECT emission data, i.e. data in both the photopeak and scatter windows, acquired in list-mode format and including the energy attribute of the detected photon, to perform ASC. We implemented a GPU-based version of this method using an ordered subsets expectation maximization (OSEM) algorithm. The method was objectively evaluated using realistic simulation studies on the task of estimating uptake in the striatal regions of the brain in a 2-D dopamine transporter (DaT)-scan SPECT study. We observed that inclusion of data from the scatter window and using list-mode data yielded improved quantification compared to using data only from the photopeak window or using binned data. These results motivate further development of list-mode-based ASC methods that include scatter-window data for SPECT.

INTRODUCTION

Single-photon emission computed tomography (SPECT) has an important role in the diagnosis and therapy of several diseases such as Parkinson’s disease, coronary artery disease, and many cancers. A major image-degrading process in SPECT is the scatter and resultant attenuation of photons as they traverse through the tissue before they reach the detector. Reliable attenuation and scatter compensation (ASC) is a pre-requisite for quantification tasks, such as quantifying biomarkers from SPECT images [1] or performing SPECT-based dosimetry [2, 3]. Also, ASC has been observed to be beneficial for visual interpretation tasks [4]. Thus, there is an important need for reliable ASC methods.

In this paper, we focus on developing ASC methods when an attenuation map is available from a transmission scan, typically a CT scan. Several such methods have been developed [5, 6, 7, 8]. However, typically existing methods do not use the precise value of the energy attribute of the detected photon, as is available when the data are stored in list-mode (LM) format. Using this precise value provides an avenue to improve the ASC. This was a major motivation for ASC approaches based on extensive spectral analysis and modeling [9, 10]. Prior work in PET imaging has shown that using LM data and incorporating energy information led to improved ASC [11, 12]. More recently, we observed that SPECT LM emission data containing the energy attribute contains information to jointly estimate the activity and attenuation distributions [13]. Also, several studies have shown that LM data can yield improved reconstruction and quantification compared to binned data in SPECT imaging [14, 15, 16, 17, 18, 19]. These investigations motivate the use of SPECT LM data containing the energy attribute for ASC.

Existing SPECT reconstruction methods also typically use only the photo-peak (PP) window data for estimating the activity. In this context, recent studies have shown that addition of data from scatter window can provide more information to estimate the activity uptake compared to using data from only the PP window [13, 20]. Using data from PP and scatter windows for activity estimation also increases the effective sensitivity of SPECT systems, that otherwise, is typically very low (< 100 counts/million emitted counts). Based on these scientific premise, we hypothesize that processing the entire SPECT data from PP and scatter windows in LM format containing the energy attribute can provide improved quantification compared to using binned data or using only PP window data.

To investigate this hypothesis, we developed a method to perform ASC using the SPECT emission data acquired in LM format and containing the energy attribute. The method is inspired by the LM reconstruction approach for PET imaging proposed by Parra et al. [21] and extends upon theory originally briefly proposed by Jha et al. for jointly estimating the activity and attenuation distribution from SPECT LM data [22]. We developed this theory specifically for estimating the activity distribution with known attenuation and derived a maximum-likelihood expectation-maximization (MLEM) algorithm for this task. An ordered-subsets version of this algorithm was then developed, and implemented on GPU-based hardware for computational efficiency. The method was objectively evaluated using simulation studies in the context of a quantitative 2D dopamine transporter (DaT)-scan SPECT study.

METHODS

THEORY

Consider a preset-time scintillation-detector-based SPECT system imaging an activity distribution denoted by the vector . The system acquires and stores data in LM format over a fixed acquisition time, T. Let J denote the number of detected events. Note that the proposed technique is also applicable to a preset-count system. Denote attributes collected for the jth LM event by the attribute vector Â. This vector contains attributes such as the position of interaction with the scintillator, energy deposited in the scintillator, time of interaction, and the angular orientation of the detector that interacted with the photon. Denote the full LM dataset as a set of attribute vector . Since the detected LM events are independent, the LM data likelihood is given by Our approach to develop the reconstruction technique is to estimate that maximizes the probability of . While obtaining the expression for Pr(J|) is easy since J is Poisson distributed, deriving an expression for pr(Âj|) is complicated. To address this issue, we use the fact that each detected photon traverses through a specific discrete sub-unit of space after being emitted. We refer to this sub-unit as a path [13, 22]. For example, in Fig.1, P1, P2, and P3 denote 3 such paths. Note that we do not know in advance of the path that a photon takes. To address this issue while deriving our reconstruction method, we define a latent variable as follows: Defining this latent variable enables developing an expectation-maximization (EM) technique to perform the reconstruction. Further, it enables deriving the expression for . More specifically, we can expand the term in (1) in terms of a mixture model as The number of detected events J is Poisson distributed with mean βT where β is the mean rate of detected photons. Using this fact and starting from Eq. 1, we can write the log-likelihood of the acquired LM data, denoted by as The term represents the radiation that is transmitted through the path . The expression for this term is given by [13] where denotes the activity of the starting voxel of the path, . The term is independent of object activity, and models the sensitivity of the path to the detector surface. This term models the attenuation and scatter of photons in the tissue as well as the transmission of photons through the collimator. Using these expressions, we can derive the log-likelihood of the observed LM data and the latent variables, denoted by , as Having defined the log-likelihood, we put forth the LM MLEM (LM-MLEM) technique. In the expectation (E) step of iteration (t+1), we take the expectation of (5) conditioned on observed data using the previous estimate of object activity distribution, (t) The result is equivalent to replacing in (5) with its expected value conditioned on the observed data. In the maximization (M) step, we maximize this conditional expectation of the log likelihood. This yields the following iterative update equation: where f and denote the activity in the q voxel and the paths that originate from the q voxel, respectively.

Fig. 1:

A schematic of SPECT system demosnstrating the definition of a path.

THE ORDERED-SUBSET LM-MLEM (OS-LM-MLEM)

The proposed method is computationally expensive. To reduce the compute time, we developed an ordered-subsets (OS) version of the technique, and implemented this version on parallelized computing hardware. We first describe the developed OS version.

The technique of OS is widely used to achieve faster convergence of EM-based reconstruction methods [23]. To develop the OS version of the developed LM-MLEM technique, similar to conventional OSEM-based algorithms, we divided the LM data into subsets based on the detector angle of each LM event. In each sub-iteration, an estimate of the activity uptake is reconstructed using all the events in a subset. Denote the subset by the index i. Also, let and denote the set of paths that reach the detector and the set of events that are detected at angles that are elements of the subset i, respectively. Then, starting from (6), the iterative update corresponding to the i’th subset and (t+1)’th iteration is derived to be where N denote the number of global iterations and N is the number of sub-iterations or equivalently the number of subsets in (7) can be computed as follows: where is a 2-D vector denoting the iteration state given by

At any iteration, the number of calculation scales linearly as number of non-zero voxels, N. The reconstruction time was reduced by mimimizing N in the initial estimate of the activity map. This was done by generating an initial crude estimate of phantom boundary using OSEM reconstruction from binned sinogram. Further, the computational complexity increases exponentially as the order of scatter. To reduce the number of paths, we considered only up to first-order scatter events and all the scatter was assumed to occur in plane.

To further increase the computational speed up, the method was parallelized and implemented on NVIDIA GPU hardware. In each iteration of the OSEM, the evaluation of the denominator of (8) and the numerator of (7) for a fixed voxel q were performed in parallel using reduction algorithm.

OBJECTIVE EVALUATION OF PROPOSED TECHNIQUE

The proposed method was evaluated in the context of a quantitative 2-D DaT-scan SPECT study, where the task was to estimate the mean activity uptake in the caudate and putamen regions of the reconstructed SPECT image. There is much interest in exploring whether uptakes in these regions can help with improved diagnosis of Parkinson’s disease. A 2D clinical SPECT scanner with a geometry similar to the Optima 640 parallel-hole collimator and imaging uptake of Ioflupane (I-131) tracer within the brain was simulated [13]. Patient anatomy and physiology were modeled using the Zubal digital brain phantom [24]. The LM data acquisition was modeled using a Monte Carlo-based code, where for each photon, we collected the position of interaction, energy of the photon, and the angular orientation of the detector. To simulate a low-dose setting, we collected around one-third of the number of photons typically acquired clinically.

The proposed OS-LM-MLEM reconstruction method was used to estimate the activity map. The experiments were repeated for multiple noise realizations. The normalized root-mean-square error (RMSE) of the estimated activity uptakes were computed.

We evaluated the effect of including data from the scatter window on estimating the activity map. For this purpose, we considered three configurations, namely, using data only from PP window, using data only from photons with energy higher than 120 keV, and using data from entire energy spectrum.

We also evaluated the effect of binning the energy and position attributes of the LM data on quantification performance. The position attribute was binned into 64 bins and the energy attribute into 2 and 3 bins in different experiments.

RESULTS

In Fig. 2, the normalized RMSE plots for the activity uptake in the caudate and putamen regions are shown as a function of iteration for the different configurations. Figs. 2a and 2b show that as we increased the range of energies considered, the RMSE reduced. For example, including data from the entire energy spectrum resulted in approximately 5% decrease in the RMSE for activity uptake in both caudate and putamen compared to using only data from PP window. Figs. 2c and 2d show that using LM data yielded lower RMSE for both caudate and putamen regions. For example, when using LM data, the RMSE of activity uptake decreased by 6% in caudate and 5% in putamen compared to using two energy bins. Representative reconstructed images with different configurations are shown in Fig. 3. Visually also, the results with LM data and using the entire energy spectrum appear to have improved quality. Overall, these results demonstrate that data processed in LM format and encompassing all the emission and scattered photons yielded superior performance.

Fig. 2:

Normalized RMSE plot as a function of iteration number for (a,b) different configurations of energy window and (c,d) different number of enrgy bins. The left and right side plots are for caudate and putamen, respectively.

Fig. 3:

(a) True activity and reconstructed image using OS-LM-MLEM (b) with LM data and the entire energy spectrum; (c) with the entire energy spectrum but binned data with 2 energy bins; (d) with LM data only from PP window. The reconstructed images (b,c,d) are represented on the same scale.

A pure LM-MLEM based technique was also developed and compared to the OS-LM-MLEM technique. We found that the reconstruction results were similar and as expected, the OSEM technique with four subsets yielded close to four-times computational speedup.

CONCLUSIONS AND FUTURE WORK

In this manuscript, we proposed an OS-LM-MLEM-based reconstruction method that uses both the scattered and photopeak data acquired in LM format to perform ASC in SPECT. Results from realistic simulation studies conducted in the context of measuring regional activity uptakes in a 2-D DaT-scan SPECT study demonstrated that inclusion of data in the scatter window yielded improved quantification compared to using only PP data. Further, processing data in LM format yielded improved quantification compared to binning the energy and position attributes. A challenge with the proposed method is the large computational time. To address this, improved optimization strategies will need to be developed. Additionally, further development of the method for 3D imaging, and application of the method to myocardial perfusion SPECT as well as other clinical SPECT imaging applications are important research frontiers. Overall, these results motivate further development of LM-based ASC methods that include data from the scatter window.