Research Article  Open Access
Yan Yan, Gengsheng L. Zeng, "Scatter and Blurring Compensation in Inhomogeneous Media Using a Postprocessing Method", International Journal of Biomedical Imaging, vol. 2008, Article ID 806705, 11 pages, 2008. https://doi.org/10.1155/2008/806705
Scatter and Blurring Compensation in Inhomogeneous Media Using a Postprocessing Method
Abstract
An efficient postprocessing method to compensate for the scattering and blurring effects in inhomogeneous medium in SPECT is proposed. A twodimensional point spread function (2DPSF) was estimated in the image domain to model the combination of these two physical effects. This 2DPSF in the inhomogeneous medium is fitted with an asymmetric Gaussian function based on Monte Carlo simulation results. An efficient further blurring and deconvolution method was used to restore images from the spatially variant 2DPSF kernel. The compensation is performed using a computersimulated NCAT phantom and a flanged Jaszczak experimental phantom. The preliminary results demonstrate an improvement in image quality and quantity accuracy with increased image contrast (25% increase compared to uncompensated image) and decreased error (40% decrease compared to uncompensated image). This method also offers an alternative to compensate for scatter and blurring in a more time efficient manner compared to the popular iterative methods. The execution time for this efficient postprocessing method is only a few minutes, which is within the clinically acceptable range.
1. Introduction
Single photon emission computed tomography (SPECT) images are degraded by attenuation, collimator and detector blurring, and photon scatter. Several studies have shown that compensations for these degradations can improve the quantitative accuracy and clinical lesion detectability [1–6]. The goal of this study is to develop a new method that can compensate for the scatter and blurring effects and improve the quantitative and qualitative accuracy of clinically realistic SPECT images. Currently, the stateoftheart compensation method is to model the scatter and blurring effects in the projector/backprojector pair of an iterative reconstruction algorithm [7–20]. The main problem with this iterative compensation method is its heavy computational burden. Also, preprocessing procedures have been investigated to compensate for these physical degradations. The blurring is compensated in a preprocessing procedure such as using the frequencydistance principle [21–24]. The scatter is corrected using energydistributionbased methods [25–27].
Previously, we have proposed a postprocessing method to compensate for the scattering in homogeneous media [28]. In this paper, we extend the postprocessing method to compensate for the scatter and blurring in inhomogeneous scattering media. We first reconstruct a raw image using an efficient analytical or iterative algorithm that corrects for attenuation only. We then model the scatter and blurring using a spatially variant twodimensional point spread function (2DPSF) in the inhomogeneous scattering media and parameterize the 2DPSF based on Monte Carlo simulations. Finally, we use an efficient further blurring and deconvolution method to restore the image.
2. Method
2.1. Monte Carlo Simulations
Monte Carlo simulations have been widely used in different areas of medical physics with the advantage of powerful computing systems [29, 30]. These Monte Carlo modeling techniques are ideal for SPECT because of the stochastic nature of radiation emission, transport, and detection processes. However, they require very long computational times. In this paper, we used the Monte Carlo simulation package SIMSET [31] to generate SPECT data with scatter contamination and detector response. The Monte Carlo data was used as a standard for scatter and blurring modeling. In the simulation, the collimator was modeled as a parallel hole collimator with a thickness of 2 cm and hole diameter of 0.14 cm. The detection energy window was centered at 140 keV with a width of 10%. The radius of rotation was 20 cm. For each phantom study, two sets of projection data were simulated. The first dataset was primary photons representing an ideal data acquisition, in which scattered photons were perfectly rejected. The second dataset contained scattered photons only. In each simulation, one billion photon histories were generated to yield lownoise projection data.
2.2. Computer Simulation Phantom
An NCAT phantom [32] was used in computer simulations. The attenuation map and activity distribution are shown in Figure 1. The intensity ratio of the activity in the myocardium versus background tissues was 5:1. The 40 cm 40 cm object region was digitized onto a 129 129 array with a pixel size of 0.31 cm (the array size of 129 129 was chosen to allow the placement of the point source in the center of the object). The object is centered on the SPECT camera's rotation axis. The projection data was collected with 300 view angles over a full 360^{°}.
(a)
(b)
2.3. Experimental Phantom
A flanged Jaszczak hotrod/coldsphere phantom was scanned for one hour using a Philips IRIX SPECT system. The phantom was filled with water and 21.6 mCi of Tc99m. Three lowenergy highresolution parallelhole collimators were used during data acquisition. The rotation radius of the collimators was 24 cm. The data were collected with 180 view angles over a full 360^{°}. The image was reconstructed in a 128 128 array with an image pixel size of 0.28 cm.
2.4. The Postprocessing Method
Our postprocessing method consists of three steps: (1) an efficient analytical or iterative algorithm that corrects for attenuation only is used to reconstruct a raw image; (2) a spatially variant twodimensional point spread function (2DPSF) in the inhomogeneous scattering medium is estimated; (3) an efficient, noniterative method is developed to restore the image.
2.4.1. Reconstruction Algorithm
We started with the raw SPECT projection data. They were contaminated by attenuation, scattering, and collimator blurring. Instead of trying to subtract the estimated scattered data from the projections, we can directly reconstruct the raw image from these projection data using either an analytical reconstruction algorithm [33–36] or an iterative MLEM reconstruction algorithm [37–39]. Here, we used the iterative algorithm as follows: where represents one pixel in the image space, is the measured SPECT emission data, and is the known coefficient that represents the contribution of image pixel i to projection bin j with the attenuation map . The summation over k is the projector, and the summation over j is the backprojector. This algorithm reconstructs a raw image with attenuation compensation (AC). However, the scattering and collimator blurring are not corrected. Our postprocessing method was applied to this raw image .
2.4.2. The TwoDimensional Point Spread Function (2dPsf)
The raw image can be modeled as a blurred version of the original image f: where the blurring kernel is what we call a 2DPSF in the image domain, and is a small positive number (i.e., in our study), representing half the size of the kernel h. The discrete version of this relation can be written as
For each image pixel , h is a 2D blurring matrix with the size of . The true image f can be solved if the kernel h is known. However, this 2DPSF h contains the effects of collimator blurring and scattering and it is normally hard to obtain. Furthermore, the 2DPSF is spatially variant, which means that it changes for every image pixel .
This 2DPSF models the scattering effect and collimator blurring in the 2D image domain instead of in the conventional 1D projection domain. It is also different from the “effective scatter source image” as proposed by Frey and Tsui [11]. Both being in the image domain, the effective scatter source image is different for each projection view and when a projection is applied to this effective image, the estimated scattered projection at this view is obtained; our proposed 2DPSF is a kernel that relates the true image and the raw reconstructed image. Also, we can obtain any projected blurring kernel by performing an attenuated projection operator on the 2DPSF.
We model the 2DPSF h in (3) as a Gaussian function with five variables (a short explanation of the reason why we are able to use the simple Gaussian function is discussed in the appendix): the magnitude of the Gaussian , full width at halfmaximum in the longaxis direction , full width at halfmaximum in the shortaxis direction , and the center of the Gaussian: where , and are the functions of the point source position .
Here, we demonstrate a similarity comparison of a measured 2DPSF and its corresponding Gaussian fitting function. Figure 2(a) is the 2DPSF calculated from a point source at a distance of 7.5 cm to the center of the rotation using Monte Carlo simulation. Figure 2(b) is the Gaussian function with parameters fitted to the 2DPSF. The comparison indicates that the twodimensional Gaussian distribution is a good fit for the 2DPSF.
(a)
(b)
2.4.3. Parameterization of 2dPsf
In order to observe the variations of the 2DPSF in different locations of the object, we perform Monte Carlo simulations for point source at eight different locations. A uniform cylinder phantom with elliptical crosssections is used. The raw reconstructed image of each point source is related to the 2DPSF at the same location. The locations of the point sources are displayed in Figure 3.
In our previous study of the 2DPSF [40], we discovered that in homogeneous scattering media, the 2DPSF is rotationally symmetric with respect to the rotation center, which means that the 2DPSF with a constant radial distance has the same shape for all angles but is rotated by a certain angle. Therefore, the 2DPSF is estimated only on the positive xaxis (Figure 3). Also, because of the localized character (i.e., small width) of the 2DPSF, it is convenient to have this assumption in the inhomogeneous case except for the variations from the different attenuation coefficients. For rotationally symmetric point sources in inhomogeneous media, if the attenuation coefficients are the same for two different point locations, we assume that the 2DPSFs are the same; if the attenuation coefficients are different, we estimate the 2DPSF according to different attenuation factors. To estimate the variations of the 2DPSF with local attenuation coefficients, three sets of Monte Carlo simulations have been performed. In each of these three simulations, we use uniform attenuation maps with the same shapes but with different coefficients: c representing the bone, c representing the tissue, and c representing the lung. All the other configurations are the same for these three simulations. In Figures 4–6, we show the variations of the amplitude , FWHM on the short axis, and FWHM on the long axis as a function of d, the distance from the point source to the rotation center. The variations of the amplitude are shown in Figure 4. It decreases with d and also varies for different attenuators. The values of are larger in highly attenuated objects: largest in the bone and smallest in the lung. This distribution agrees with the scatter probability derived from the KleinNishima formula [41]. We fit as a function of both d and the attenuation distribution : in which where is dimensionless and represents the relative magnitude increase due to scattering. The secondorder polynomial of the attenuation factors is chosen to get reasonably wellfitting results.
The full width at halfmaximum (FWHM) of the Gaussian function is determined by the combined effects of collimator blurring and scatter blurring. Because the amplitude of the reconstructed image from the scattered data is small compared to the reconstructed image from the primary photons (see the appendix), the FWHM of the h is mainly determined by the FWHM of the reconstructed image from the primary photons. Therefore, the blurring in the 2DPSF is most affected by collimator blurring and is independent of the attenuators. This is verified by the simulation results as shown in Figures 5 and 6. As the source moves from the center of the object toward the edge, the value of decreases (Figure 5). Little change is observed among three different objects. Figure 6 shows a plot of , the full width at halfmaximum on the long axis of the Gaussian function. It is observed that the value of barely changes as the source moves from the center toward the edge. Also, note that the largest difference is a little over one pixel. Similar variations are proposed by Zeng and Huang [42]. These two parameters are fitted based on a distancedependent model [43], which will not change for different attenuators: where r is the radius of the collimator hole, l is the thickness of the collimator, d is the distance from the point source to the center of the object , and D represents the distance from the center of rotation to the detector in cm, which equals to 20 cm in our simulations.
With (4)–(7), we can calculate the 2DPSF for any point source inside the object. These empirical formulae eliminate the need for extensive Monte Carlo simulations for each point source. Also, the estimation of the 2DPSF using these formulae is independent of the raw image of the simulation phantom. It is derivable from the attenuation map and configurations of the collimator and can be calculated before reconstruction.
2.4.4. Image Restoration
As the 2DPSF is a spatially variant, a normal deconvolution algorithm cannot be used. In order to avoid the long computational time of using iterative restoration algorithms, we used a furtherblurring and deconvolution method to restore the image [42]. This method converted the raw image with a spatially variant point spread function into a further blurred image with a spatially invariant point spread function, and then used an efficient technique (e.g., a frequency domain filtering) for deblurring.
The further blurring was implemented by using a rotational convolution. Let the raw reconstructed image be . We rotated the image counterclockwise by a small angle about the axis of the detector rotation obtaining and rotated clockwise by obtaining . When necessary, we rotated the image counterclockwise by obtaining and rotated clockwise by obtaining and so on. A weighted sum of these rotated images gives a further blurred image : where the weighting factors form a convolution kernel. The sum of is normalized to 1 to assure the consistency of the image intensity. The weighting factors are chosen empirically, so that the 2DPSF is spatially invariant. is the amplitude of 2DPSF discussed in Section 3, which is used to normalize the amplitude of the raw image with different radial distances d.
Now, this further blurred image has an approximately spatially invariant point spread function and this is nothing but the 2DPSF at the center of the object. Then, we perform an efficient inverse filtering (e.g., the Wiener Filtering) on this image to obtain the restored image f in (1).
2.5. Assessment of Restored Images
Several measurements were performed in the computer simulation results to evaluate the improvement of the image quality using the proposed compensation method. (a)Sumsquared error (SSE) was used to measure the average discrepancy of the restored image with respect to the original image. It is defined as the averaged sum of the squared pixel difference as follows: where and represent the restored value and the true value for pixel i, respectively, and N is the total number of pixels calculated.(b)Contrast (CR) between myocardium and background is defined as where FG represents the average pixel value in the myocardium, and BG represents the average pixel value in the background.(c)Noise was measured as the standard deviation of pixel counts in the uniform background, normalized by the mean activity of that region.
3. Results
3.1. Computer Simulations
The Monte Carlo simulation data for the NCAT phantom was used to reconstruct the raw image using the MLEM iterative algorithm. The attenuation correction was performed in the MLEM reconstruction using a blurred attenuation map. The map was blurred with a Gaussian function to match the resolution of the emission data. The parameters of this Gaussian function are empirically chosen to obtain the optimal reconstructed image with least crosstalk artifacts. Figure 7(b) demonstrates the raw reconstructed image without using the blurred attenuation map, in which there were crosstalk effects due to the unmatched resolution between the attenuation map and the emission data. Figure 7(c) shows the raw image reconstructed using the blurred attenuation map. The 2DPSF was then estimated using a smeared attenuation map and collimator parameters. This smeared attenuation map used in the 2DPSF estimation is different from the blurred map used in the reconstruction. This smearing is to integrate the influence of neighborhood pixels into the point source because the 2DPSF represents the total scattered photons that originate from a point source and interact with its neighborhood pixels. For restoration, the further blurred image was obtained by rotational convolution: where represents the raw image, and represents the image rotated by an angle . The weighting factors were determined empirically to get image so that it has a spatially invariant point spread function. Wiener filtering was then performed to restore the further blurred image Figure 7(d) shows the restored image. It is observed that there exists a dishing effect in the liver and boundary area. This is caused by the over filtering in the inverse filtering step. The tradeoff between over filtering and the effectiveness of the inverse filtering is a limitation of our method and needs further investigation in the future. A horizontal profile through the center of the images is shown in Figure 8.
(a)
(b)
(c)
(d)
The contrast, SSE, and noise were calculated for all images to illustrate the improvement of image quality in the restored image (Table 1). The raw image here is reconstructed with a blurred attenuation map. It is observed that after compensation, the quantitative accuracy and contrast were improved, and noise was controlled from being elevated.

3.2. Phantom Experiment
The MLEM iterative algorithm with 50 iterations was used for raw image reconstruction and attenuation correction. The 2DPSF was estimated based on the waterfilled uniform attenuator and the lowenergy highresolution collimator. For image restoration, the further blurred image was obtained by rotational convolution as follows: The restored image was obtained after application of Wiener filtering. Figure 9(c) shows the restored image. It is observed that the restored image is less noisy than the raw image, as shown in Figure 9(c) compared to Figure 9(b). This may be due to the fact that the Wiener filter is a bandpass filter, and the highfrequency noise is suppressed.
(a)
(b)
(c)
4. Discussions
The goal of this postprocessing method is to develop a time efficient compensation method and overcome the heavy computational burden in the iterative reconstructionbased method. There are several issues to mention in this section.
4.1. Accuracy and Generality of The 2dPsf Estimation
The estimation of the 2DPSF is the main challenge of the proposed method. The 2DPSF models the scattering and collimator blurring in the image domain instead of the conventional projection domain. We used a Gaussian function to approximately model the 2DPSF. The validity of using the Gaussian function was discussed. We derived empirical formulae for the parameters of the Gaussian function from the Monte Carlo simulations. As the 2DPSF is objectdependent, we need to pay attenuation to the generality of the estimations. As discussed in this study, the parameter of the Gaussian function depends on the local attenuation coefficient, and the FWHMs stay the same for different attenuation distributions and only depend on collimator configurations. More Monte Carlo simulations for objects with various sizes and shapes are desired to further determine a more accurate and general 2DPSF model.
4.2. “further Blurring and Deconvolution” Restoration
This method efficiently restores images with spatially variant point spread functions. The advantage of this approach is its fast implementation compared with the conventional iterative algorithms. However, this is an approximate method, and the rotation angle and the weighting factors in (8) are currently determined empirically. As the goal of the postprocessing method is to cut down the computational time for the compensation, an efficient restoration method like this one is desired. Other efficient restoration methods with the capability of deblurring spatially variant point spread functions can also be used.
4.3. Computational Time
The computational time for this postprocessing compensation method was reduced compared to the iterative reconstructionbased method. The current time for fast implementation of the reconstructionbased method for a image array is in the range of thirty minutes [1, 17]. In this proposed postprocessing method, the computer time for getting a twodimensional attenuationcorrected raw image was in the order of seconds for fast reconstruction algorithms [34, 35, 39]. We then precalculated the 2DPSF and stored it in the computer memory. In the last step, all we needed was a few more seconds for image restoration (i.e., rotational convolution and Wiener filtering). Therefore, the total computer time for this postprocessing compensation method was only few minutes and is acceptable for clinical applications.
5. Conclusions
We have presented an efficient postprocessing method to compensate for scattering and blurring effects in inhomogeneous media. The major challenge of the method is to accurately estimate the 2DPSF in the image domain. Empirical formulae are proposed to model the 2DPSF variations with the various locations within nonuniformly attenuated objects. From the clinical aspect, the implementation of our method is faster (within several minutes) than the iterative reconstructionbased compensation method. One limitation of this study is that it is developed in two dimensions and does not consider scattered photons from outofplane sources. Our future work includes modeling the scattering with a 3DPSF.
Appendix
Monte Carlo simulations are performed to generate two sets of projection data: one is the primary scatterfree data, denoted by p, and the other is the scattered data only, denoted by s. Both datasets are contaminated by blurring effect. We use and to represent the raw reconstructed images from primary and scattered photons, respectively. The sum of these two images (denoted by ) is the raw image with both scatter and blurring contaminations. As defined in (1), the amplitude of the 2DPSF is determined by the total volume change in the image with repect to Figure 10 shows the profiles of the , and of a point source. Although the ratio of total volume from the image over the total volume from the image is around 2:1, the amplitude of the is much larger than that of the . Therefore, the shape of the 2DPSF is determined by the shape of point source image . Also in Figures 11(a) and 11(b), we demonstrate the comparison of the 2DPSF (related to the ) with the fitted Gaussian function in both linear and logarithmic scales. It is observed that the Gaussian function fits the overall shape of the 2DPSF very well except in the tails of the 2DPSF. The discrepancy in the tails is magnified in Figure 11(b) in logarithmic scale. Compared with the amplitude and total area of , the discrepancy in the tails from the Gaussian fittings is very small and neglectable. Therefore, the Gaussian function is a good fit for the 2DPSF.
(a)
(b)
(a)
(b)
Another point worth mentioning is the asymmetry shape along the long axis direction of the 2DPSF. Previous study [9] discovered the asymmetric shape of the projectiondomain scatter response function and modeled it with two different Gaussians on either side of the point source. In our approach, it is not necessary to model the scatter with two Gaussians. One reason for that is our 2DPSF is a reconstructed image from all the asymmetric projections. The asymmetry in the projections is balanced out and not significant observed in the 2DPSF. Figures 12(a) and 12(b) show the plots of the 2DPSF on the long axis for point sources with three different radial distances. The asymmetry can be barely observed even in the logarithmic scale. Furthermore, as the amplitude of is much smaller than that of as observed in Figure 10, the asymmetric distribution of the scatter image contributes little to the image , which represents the overall effects of both blurring and scatter. Therefore, it is reasonable to model the scatter and blurring with one Gaussian function in (4).
(a)
(b)
References
 F. J. Beekman, H. W. A. M de Jong, and S. van Geloven, “Efficient fully 3D iterative SPECT reconstruction with Monte Carlobased scatter compensation,” IEEE Transactions on Medical Imaging, vol. 21, no. 8, pp. 867–877, 2002. View at: Publisher Site  Google Scholar
 S. Sankaran, E. C. Frey, K. L. Gilland, and B. M. W. Tsui, “Optimum compensation method and filter cutoff frequency in myocardial SPECT: a human observer study,” Journal of Nuclear Medicine, vol. 43, no. 3, pp. 432–438, 2002. View at: Google Scholar
 M. V. Narayanan, M. A. King, P. H. Pretorius et al., “Humanobserver receiveroperatingcharacteristic evaluation of attenuation, scatter, and resolution compensation strategies for ${}^{99\text{m}}\text{TC}$ myocardial perfusion imaging,” Journal of Nuclear Medicine, vol. 44, no. 11, pp. 1725–1734, 2003. View at: Google Scholar
 M. V. Narayanan, P. H. Pretorius, S. T. Dahlberg et al., “Evaluation of scatter compensation strategies and their impact on human detection performance Tc99 m myocardial perfusion imaging,” IEEE Transactions on Nuclear Science, vol. 50, no. 5, part 2, pp. 1522–1527, 2003. View at: Publisher Site  Google Scholar
 X. He, E. C. Frey, J. M. Links, K. L. Gilland, W. P. Segars, and B. M. W. Tsui, “A mathematical observer study for the evaluation and optimization of compensation methods for myocardial SPECT using a phantom population that realistically models patient variability,” IEEE Transactions on Nuclear Science, vol. 51, no. 1, part 1, pp. 218–224, 2004. View at: Publisher Site  Google Scholar
 H. Zaidi and K. F. Koral, “Scatter modelling and compensation in emission tomography,” European Journal of Nuclear Medicine and Molecular Imaging, vol. 31, no. 5, pp. 761–782, 2004. View at: Publisher Site  Google Scholar
 C. E. Floyd, Jr., R. J. Jaszczak, S. H. Manglos, and R. E. Coleman, “Compensation for collimator divergence in SPECT using inverse Monte Carlo reconstruction,” IEEE Transactions on Nuclear Science, vol. 35, no. 1, part 12, pp. 784–787, 1988. View at: Publisher Site  Google Scholar
 B. M. W. Tsui, H.B. Hu, D. R. Gilland, and G. T. Gullberg, “Implementation of simultaneous attenuation and detector response correction in SPECT,” IEEE Transactions on Nuclear Science, vol. 35, no. 1, part 12, pp. 778–783, 1988. View at: Publisher Site  Google Scholar
 E. C. Frey and B. M. W. Tsui, “A practical method for incorporating scatter in a projectorbackprojector for accurate scatter compensation in SPECT,” IEEE Transactions on Nuclear Science, vol. 40, no. 4, part 12, pp. 1107–1116, 1993. View at: Publisher Site  Google Scholar
 E. C. Frey and B. M. W. Tsui, “Modeling the scatter response function in inhomogeneous scattering media for SPECT,” IEEE Transactions on Nuclear Science, vol. 41, no. 4, part 12, pp. 1585–1593, 1994. View at: Publisher Site  Google Scholar
 E. C. Frey and B. M. W. Tsui, “A new method for modeling the spatiallyvariant, objectdependent scatter response function in SPECT,” in Proceedings of IEEE Nuclear Science Symposium, vol. 2, pp. 1082–1086, Anaheim, Calif, USA, November 1996. View at: Publisher Site  Google Scholar
 B. C. Penney, M. A. King, and K. Knesaurek, “A projector, backprojector pair which accounts for the twodimensional depth and distance dependent blurring in SPECT,” IEEE Transactions on Nuclear Science, vol. 37, no. 2, part 1, pp. 681–686, 1990. View at: Publisher Site  Google Scholar
 F. J. Beekman, E. G. J. Eijkman, M. A. Viergever, G. F. Borm, and E. T. P. Slijpen, “Object shape dependent PSF model for SPECT imaging,” IEEE Transactions on Nuclear Science, vol. 40, no. 1, pp. 31–39, 1993. View at: Publisher Site  Google Scholar
 G. L. Zeng, G. T. Gullberg, B. M. W. Tsui, and J. A. Terry, “Threedimensional iterative reconstruction algorithms with attenuation and geometric point response correction,” IEEE Transactions on Nuclear Science, vol. 38, no. 2, part 12, pp. 693–702, 1991. View at: Publisher Site  Google Scholar
 A. Welch, G. T. Gullberg, P. E. Christian, F. L. Datz, and H. T. Morgan, “A transmissionmapbased scatter correction technique for SPECT in inhomogeneous media,” Medical Physics, vol. 22, no. 10, pp. 1627–1635, 1995. View at: Publisher Site  Google Scholar
 A. Welch and G. T. Gullberg, “Implementation of a modelbased nonuniform scatter correction scheme for SPECT,” IEEE Transactions on Medical Imaging, vol. 16, no. 6, pp. 717–726, 1997. View at: Publisher Site  Google Scholar
 D. J. Kadrmas, E. C. Frey, S. S. Karimi, and B. M. W. Tsui, “Fast implementations of reconstructionbased scatter compensation in fully 3D SPECT image reconstruction,” Physics in Medicine and Biology, vol. 43, no. 4, pp. 857–873, 1998. View at: Google Scholar
 R. G. Wells, A. Celler, and R. Harrop, “Analytical calculation of photon distributions in SPECT projections,” IEEE Transactions on Nuclear Science, vol. 45, no. 6, part 3, pp. 3202–3214, 1998. View at: Publisher Site  Google Scholar
 C. Bai, G. L. Zeng, and G. T. Gullberg, “A slicebyslice blurring model and kernel evaluation using the Klein Nishina formula for 3D scatter compensation in parallel and converging beam SPECT,” Physics in Medicine and Biology, vol. 45, no. 5, pp. 1275–1307, 2000. View at: Google Scholar
 H. W. A. M. De Jong, E. T. P. Slijpen, and F. J. Beekman, “Acceleration of Monte Carlo SPECT simulation using convolutionbased forced detection,” IEEE Transactions on Nuclear Science, vol. 48, no. 1, part 1, pp. 58–64, 2001. View at: Publisher Site  Google Scholar
 R. M. Lewitt, P. R. Edholm, and W. Xia, “Fourier method for correction of depthdependent collimator blurring,” in Medical Imaging III: Image Processing, vol. 1092 of Proceedings of SPIE, pp. 232–243, Newport Beach, Calif, USA, JanuaryFebruary 1989. View at: Google Scholar
 W. G. Hawkins, N.C. Yang, and P. K. Leichner, “Validation of the circular harmonic transform (CHT) algorithm for quantitative SPECT,” Journal of Nuclear Medicine, vol. 32, no. 1, pp. 141–150, 1991. View at: Google Scholar
 S. J. Glick, B. C. Penney, M. A. King, and C. L. Byrne, “Noniterative compensation for the distancedependent detector response and photon attenuation in SPECT imaging,” IEEE Transactions on Medical Imaging, vol. 13, no. 2, pp. 363–374, 1994. View at: Publisher Site  Google Scholar
 W. Xia, R. M. Lewitt, and P. R. Edholm, “Fourier correction for spatially variant collimator blurring in SPECT,” IEEE Transactions on Medical Imaging, vol. 14, no. 1, pp. 100–115, 1995. View at: Publisher Site  Google Scholar
 R. J. Jaszczak, K. L. Greer, C. E. Floyd, Jr., C. C. Harris, and R. E. Coleman, “Improved SPECT quantification using compensation for scattered photons,” Journal of Nuclear Medicine, vol. 25, no. 8, pp. 893–900, 1984. View at: Google Scholar
 M. A. King, G. J. Hademenos, and S. J. Glick, “A dualphotopeak window method for scatter correction,” Journal of Nuclear Medicine, vol. 33, no. 4, pp. 605–612, 1992. View at: Google Scholar
 K. Ogawa, T. Ichihara, and A. Kubo, “Accurate scatter correction in single photon emission CT,” Annals Nuclear Medicine and Sciences, vol. 7, no. 3, pp. 145–150, 1994. View at: Google Scholar
 Y. Yan and G. L. Zeng, “A postprocessing method for scatter compensation in SPECT,” in Proceedings of the International Meeting on Fully ThreeDimensional Image Reconstruction in Radiology and Nuclear Medicine, pp. 186–189, Lindau, Germany, July 2007. View at: Google Scholar
 H. Zaidi, “Relevance of accurate Monte Carlo modeling in nuclear medical imaging,” Medical Physics, vol. 26, no. 4, pp. 574–608, 1999. View at: Publisher Site  Google Scholar
 I. Buvat and I. Castiglioni, “Monte Carlo simulations in SPET and PET,” Quarterly Journal of Nuclear Medicine, vol. 46, no. 1, pp. 48–61, 2002. View at: Google Scholar
 R. L. Harrison, S. D. Vannoy, D. R. Haynor, S. B. Gillispie, M. S. Kaplan, and T. K. Lewellen, “Preliminary experience with the photon history generator module of a publicdomain simulation system for emission tomography,” in Proceedings of the IEEE Nuclear Science Symposium and Medical Imaging Conference (NSS/MIC '93), pp. 1154–1158, San Francisco, Calif, USA, OctoberNovember 1993. View at: Google Scholar
 W. P. Segars, Development of a new dynamic NURBSbased cardiactorso (NCAT) phantom, Ph.D. dissertation, The University of North Carolina, Chapel Hill, NC, USA, 2001.
 F. Natterer, “An inversion of the attenuated Radon transform,” Inverse Problems, vol. 17, no. 1, pp. 113–119, 2001. View at: Google Scholar
 R. G. Novikov, “An inversion formula for the attenuated Xray transformation,” Arkiv for Matematik, vol. 40, no. 1, pp. 145–167, 2002. View at: Publisher Site  Google Scholar  MathSciNet
 H. Zaidi and K. F. Koral, “Scatter modeling and compensation in emission tomography,” European Journal of Nuclear Medicine and Molecular Imaging, vol. 31, no. 5, pp. 761–782, 2004. View at: Publisher Site  Google Scholar  MathSciNet
 Q. Huang, G. L. Zeng, J. You, and G. T. Gullberg, “An FDKlike conebeam SPECT reconstruction algorithm for nonuniform attenuated projections acquired using a circular trajectory,” Physics in Medicine and Biology, vol. 50, no. 10, pp. 2329–2339, 2005. View at: Google Scholar
 L. A. Shepp and Y. Vardi, “Maximum likelihood reconstruction for emission tomography,” IEEE Transactions on Medical Imaging, vol. 1, no. 2, pp. 113–122, 1982. View at: Publisher Site  Google Scholar
 K. Lange and R. Carson, “EM reconstruction algorithms for emission and transmission tomography,” Journal of Computer Assisted Tomography, vol. 8, no. 2, pp. 306–316, 1984. View at: Google Scholar
 G. T. Gullberg, R. H. Huesman, J. A. Malko, N. J. Pelc, and T. F. Budinger, “An attenuated projectorbackprojector for iterative SPECT reconstruction,” Physics in Medicine and Biology, vol. 30, no. 8, pp. 799–816, 1985. View at: Publisher Site  Google Scholar
 Y. Yan and G. L. Zeng, “A postprocessing method for scatter and collimator blurring compensation using spatially variant point spread function,” in Proceedings of IEEE Nuclear Science Symposium and Medical Imaging Conference (NSS/MIC '06), San Diego, Calif, USA, OctoberNovember 2006. View at: Google Scholar
 O. Klein and T. Nishina, “Über die Streuung von Strahlung durch freie Elektronen nach der neuen relativistischen Quantendynamik von Dirac,” Zeitschrift für Physik A Hadrons and Nuclei, vol. 52, no. 1112, pp. 853–868, 1929. View at: Publisher Site  Google Scholar
 G. L. Zeng and Q. Huang, “Compensation for collimator blurring using rotational and axial convolution,” in Proceedings of the International Meeting on Fully ThreeDimensional Image Reconstruction in Radiology and Nuclear Medicine, pp. 329–332, Lindau, Germany, July 2007. View at: Google Scholar
 H. W. A. M. de Jong, E. T. P. Slijpen, and F. J. Beekman, “Acceleration of Monte Carlo SPECT simulation using convolutionbased forced detection,” IEEE Transactions on Nuclear Science, vol. 48, no. 1, part 1, pp. 58–64, 2001. View at: Publisher Site  Google Scholar
Copyright
Copyright © 2008 Yan Yan and Gengsheng L. Zeng. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.