Monte-Carlo-Based Scatter Correction for Quantitative SPECT Reconstruction Realization and Evaluation Rolf Bippus1 , Andreas Goedicke1 , Henrik Botterweck2 1 2 Philips Research Laboratories, Aachen Fachhochschule Lübeck, Bildgebende Verfahren in der Medizintechnik rolf.bippus@philips.com Abstract. Quantitative SPECT as well as simultaneous acquisition of multiple isotopes with SPECT in the clinical field, although clinically interesting, are still limited by reconstruction artifacts and computing power. As a considerable step in this direction, we have implemented an efficient reconstructor with variance reduced Monte-Carlo-simulation in the forward and/or backward projection of an OS-EM iteration. Apart from a quantitatively accurate scatter estimation the integrated MC simulation allows us to include all effects that are relevant for the multiple isotope task. The reconstruction problem is explicitly formulated as a combined maximum-likelihood estimation for all unknowns in one step and implemented efficiently on the basis of the OS-EM algorithm using the Effective Scatter Source approach. The algorithm has been evaluated on simulated as well as measured phantom data. Two isotopes (Tc99m, Tl201) were tested on the voxelized NCAT phantom and evaluated quantitatively. In addition we performed phantom measurements to evaluate the method on a Philips CardioMDT M with a VantageT M line source for attenuation correction. A clear advantage of the proposed approach is its robustness and generalizability. It is currently being evaluated in clinical applications like simultaneous dual isotope cardiac imaging. 1 Introduction Absolute quantification as well as simultaneous acquisition of multiple isotopes are receiving increasing interest in clinical applications (e.g. [1], [2]). However patient scatter, down-scatter plus collimator effects such as penetration and Pbfluorescence result in data contamination. Basically one can distiguish between energy and/or spatial distribution based scatter estimation and correction from reconstruction based methods [3]. The OS-EM (Ordered Subset Expectation Maximization) approach we follow with our implementation falls into the second category. Scatter is estimated based on the current activity estimate in the forward projection step. Using the approximative scheme of effective scatter source (ESS) Estimation, originally published by Frey et al. [4], scatter estimation is effectively separated from other effects in the algorithm. The original approach uses precomputed scatter kernels to Quantitative SPECT Reconstruction 371 obtain the ESS via convolution. In contrast we apply variance reduced MonteCarlo simulation to estimate the ESS. An attenuation map is used as a physical model of the patient. Using one of the reconstruction based scatter estimation methods one can differentiate sequential and simultaneous reconstruction as depicted in Fig. 1. We follow the scheme of simultaneous reconstruction (Fig. 1(b)) where the cross contamination estimation is integrated into the forward projector and all estimates isotope activities are forward projected into all energy windows and updated in each sub-iteration of the OS-EM algorithm. 2 Methods The full update equations for the ordered subset EM algorithm for all isotope activities (typicall 1 or 2) fǫ,i (isotope i in voxel ǫ) on all available data are given by equations (1) and (2). The forward projection step computes (1) the contributions of isotopes i to all projection pixels p in each energy-window e, based on the recent activity estimate (k) ǫ,i is the normalized system operator. In the backprojection step fǫ,i , where Hp,e (2), the correction factors are obtained from the ratio of the observed and the estimated projection values. X (k) X X (k) X ǫ,i ǫ,i ν (k) (p, e) = νi (p, e) = , fǫ,i Hp,e Hp,e =1 (1) i (k+1) fǫ,i 2.1 (k) = fǫ,i i ǫ 1 X X ǫ,i ν(p, e) Hp,e (k) nǫ,i p e ν (p, e) p,e , nǫ,i = X ǫ,i Hp,e (2) p,e Approximations / Implementation Using brute force counting statistics with MC-simulation to obtain the ν (k) (p, e) in equation (1) is not feasible. Therefore we applied the principle of convolution based forced detection. Simulating the photon tracks, at each interaction of the photons with the patient the effective cross section is calculated for the photon being redirected towards the detector for a set of intermediary energy intervals. These are accumulated in 3D and are treated as secondary emission distributions, the ESS. These are then projected onto all camera positions and energy windows, taking (c) (d) Fig. 1. Sequential (a) vs. simultaneous (b) dual isotope reconstruction scheme. 372 Bippus, Goedicke & Botterweck into account attenuation and collimator response functions (CRF). The CRFs, simulated beforehand and stored, are taking into account all major physical degradation effects. The principle is illustrated in Fig. 2. According to the Dual-Matrix aproach, the backprojector is approximated ǫ ) as well as the relative weights by the reversed PSFs and the attenuation (Bp,e ωi (p, e) of the isotope i contributing to pixel (p, e) in the forward projection, thus ignoring the spatial distribution of the scatter. (k+1) fǫ,i 2.2 (k) = fǫ,i 1 X X ǫ (k) ν(p, e) B ω (p, e) (k) nǫ,i p e p,e i ν (p, e) (3) Activity Estimation Quantitative results are obtained by defining volumes of interest (VOI) and calculating the mean activity within the VOI. We evaluated quantitative activities after 30 to 50 full iterations of the OS-EM algorithm. Corresponding images are not suited for visual inspection or VOI delineation due to high noise. Thus images for visual inspection or delineation use 4 full iterations of the OS-EM algorithm. For the simulated data we defined the VOIs directly from true activity images. 3 3.1 Results NCAT Phantom Results The voxelized NCAT phantom was simulated such that the activities represent a patient dose of 25 mCi Tc99m and 3.5 mCi of Tl201 . In the Tl201 distribution a heart lesion was included, which was not present in the Tc99m image. The lesion activity was chosen 15 % of the normal myocardial Tl201 activity. Fig. 2. Principle of down-scatter estimation in the forward projection step. Quantitative SPECT Reconstruction 373 For quantitative comparison of true and reconstructed activity, we measured absolute activity in the 17 segments of a standard 17 segment polar plot representation. Heart re-orientation and segmentation were performed on the true activity image and kept constant for the reconstructed image. As before the values represent absolute activity estimates averaged over each segment of the heart. 3.2 Phantom Measurement Results Phantom measurements were performed using a Data Spectrum Cooperation Anthropomorphic Torso PhantomT M with a fillable heart insert. The phantom was filled according to an estimated biodistribution resulting from 90 MBq Tl201 and 650 MBq Tc99m administered to an adult. The heart contained 2 lesions of 5.6 ml (lesion 1) and 12 ml (lesion 2), each filled with Tc99m only, which gives no Tl201 and approx. 65% Tc99m (due to the walls of the lesions) in the lesions. Resulting images and polar plots are shown in Fig. 4. 4 Conclusion Quantitative MC reconstruction for SPECT and multiple isotopes shows to be feasible. The more exact modeling of scatter within the object and collimator results in reduced image artifacts. The computational effort is reasonable, especially compared to other state of the art scatter corrections like ESSE. An advantage of the proposed approach is its robustness and generalizability: Fast adaptation to new isotope combinations is possible. Next steps will include the evaluation of its potential in clinical application of dual isotope cardiac imaging. (a) truth (b) 30 iterations Fig. 3. NCAT, T c99m -T l201 SDI, T l201 mean activity in kBq/ml in each of 17 segments of standard polar plot. 374 Bippus, Goedicke & Botterweck Fig. 4. Phantom measurements on the CardioMD. Cardiac SDI filled with Tl201 and Tc99m and two lesions filled with Tc99m only. (a) (b) References 1. Berman D, Kang X, Tamarappoo B, et al. stress thallium-201/rest technetium- 99m sequential dual isotope high-speed myocardial perfusion imaging. JACC Cardiovasc Imaging. 2009;2(3):273–82. 2. Sgouros G. Dosimetry of internal emitters. J Nucl Med. 2005;46(1):273–82. 3. King J M A ad Glick, Pretorius PH, Wells RG, et al. Attenuation, Scatter and Spatial Resolution Compensation in SPECT. In: Wernick M, Aarsvold JS, editors. Emission Tomography. The Fundamentals of PET and SPECT. Acedemic Press; 2004. 4. Frey EC, Tsui BMW. A new method for modelling the spatially variant, objectdependent scatter response function in SPECT. In: Proc IEEE Nucl Sci Symp; 1996. p. 1082–6.