ORIGINAL ARTICLE Year : 2013  Volume : 38  Issue : 4  Page : 189197 Evaluation of six scatter correction methods based on spectral analysis in ^{99m} Tc SPECT imaging using SIMIND Monte Carlo simulation Mahsa Noori Asl^{1}, Alireza Sadremomtaz^{1}, Ahmad BitarafanRajabi^{2}, ^{1} Department of Physics, Faculty of Sciences, University of Guilan, Rasht, Iran ^{2} Department of Nuclear Medicine, Rajaei Cardiovascular, Medical and Research Center, Iran University of Medical Sciences, Tehran, Iran Correspondence Address: Comptonscattered photons included within the photopeak pulseheight window result in the degradation of SPECT images both qualitatively and quantitatively. The purpose of this study is to evaluate and compare six scatter correction methods based on setting the energy windows in ^{99m} Tc spectrum. SIMIND Monte Carlo simulation is used to generate the projection images from a coldsphere hotbackground phantom. For evaluation of different scatter correction methods, three assessment criteria including image contrast, signaltonoise ratio (SNR) and relative noise of the background (RNB) are considered. Except for the dualphotopeak window (DPW) method, the image contrast of the five cold spheres is improved in the range of 2.726%. Among methods considered, two methods show a nonuniform correction performance. The RNB for all of the scatter correction methods is ranged from minimum 0.03 for DPW method to maximum 0.0727 for the three energy window (TEW) method using trapezoidal approximation. The TEW method using triangular approximation because of ease of implementation, good improvement of the image contrast and the SNR for the five cold spheres, and the low noise level is proposed as most appropriate correction method.
Introduction One of the major problems in single photon emission computed tomography (SPECT) is the inclusion of scattered photons within the photopeak energy window used for acquisition of the projection images. Scattering can occur as Compton scatter and coherent scatter. In the energy range of interest in nuclear medicine, the probability of occurrence for Compton scatter is much more than that for coherent scatter. In Compton scatter, not only the photon direction is changed but also photon energy dependent on the scatter angle is reduced. [1] Because of the poor energy resolution of the NaI(Tl) scintillation crystal used in the imaging system (about 910% full width at half maximum (FWHM) at 140 keV), it is unavoidable to detect some scattered photons in the photopeak window. [2] Because of the changed direction of Comptonscattered photons, they carry inaccurate spatial information, and therefore the detection of them leads to reduce image contrast and degrade image quality. Thereby, it seems that scatter correction is necessary to improve the image contrast and enhancement of the diagnostic accuracy. [3] Many methods for scatter correction of SPECT images have been proposed by several investigators. [4],[5],[6],[7],[8],[9],[10],[11],[12] Some of them are based on a spectral analysis that is performed by setting additional energy windows in 99m Tc spectrum, and others rely on a spatial analysis that is implemented by convolution and deconvolution procedures. In this study, six scatter correction methods from the first category are individually evaluated. These methods consist of dualenergy window (DEW) method, photopeak energy distribution analysis (PEDA) method, dualphotopeak window (DPW) method, cannel ratio method (CRM), and three energy window (TEW) methods (trapezoidal and triangular approximation). A coldsphere hotbackground phantom including six spheres with different diameters is simulated to generate the projection images for the different energy windows and evaluate each of the scatter correction methods based on the change resulted in image contrast, signaltonoise ratio (SNR), and relative noise of the background (RNB) for each of the cold spheres. Materials and Methods Scatter correction methods The underlying assumptions for each of the scatter correction methods are described below. In all of these methods, a 126154 keV energy window centered on 140 keV is used as photopeak energy window (W pk ). Dualenergy window (DEW) method This method [4],[5] relies on the assumption that the spatial distribution of the counts acquired in a secondary window placed in Compton region of the energy spectrum (W C = 92125 keV) is same with the spatial distribution of the scatter counts included in the photopeak window, and that they only quantitatively differ in a factor k, that is: [INLINE:1] where S pk is the scatter projection obtained from the photopeak window, T C is the total projection from the Compton window, and (i,j) denotes a given pixel position in the projection image. Then, the estimated projection of unscatter photons in the photopeak window, US pk , can be deduced by: [INLINE:2] It is common to use a k valve of 0.5 proposed originally by Jaszczak et al., [4] but since the k value depends on the phantom and imaging situation used, we individually calculate this value here. Photopeak energy distribution analysis (PEDA) method In this method, [6] it is assumed that the photopeak energy window can be divided into two subwindows such that for each of the projection image pixels, the number of scatter events included within these two subwindows are approximately equal, that is: [INLINE:3] where S 1 and S 2 denote the number of scatter photon events in the lower and upper subwindows, respectively. Therefore, the scatter correction can be performed by subtracting the total projections of the lower subwindow from that of the upper subwindow, that is: [INLINE:4] It is obvious that this subtraction results in removing some unscatter events as well as the scatter events. The essential step in this method is to determine the cutoff energy between subwindows 1 and 2 so that for each pixel, the ratio of the scatter counts detected in two subwindows is as close as possible to unit. Dualphotopeak energy (DPW) method In this correction method, [7] the photopeak window is divided into two equal width subwindows and is assumed that for any image pixel, there is a regression relation between the ratio of the scatter counts to the unscatter counts in the photopeak window, SUR, and the ratio of the total counts detected in the lower subwindow (w lw ) to the upper subwindow (w uw ), R, as shown in the following equation: [INLINE:5] Then, the number of scatter photon events in the photopeak window can be estimated as: [INLINE:6] Cannel ratio method (CRM) The same two subwindows used in the DPW method are used again in this correction method. [8] This method assumes that both the ratio of unscatter counts and the ratio of scatter counts detected in two subwindows are constant: [INLINE:7] In this way, the parameters G(i,j) and H(i,j) are calculated for each pixel in the projection images. Then, the number of unscatter photons in the photopeak window can be estimated by following equation: [INLINE:8] where G and H are the mean values calculated for the two parameters G(i,j) and H(i,j). In this method, the scatter correction is performed directly and without need to estimate the photopeak scatter component. Three energy window (TEW) method using trapezoidal approximation In this method [9] (referred to as TEW1), it is assumed that the photopeak scatter spectrum can be estimated by the area of a trapezoid that its left and right heights are equal to the total number of photons acquired in two narrow energy windows centered on lower and upperedge energies of the photopeak window divided by the width of the window, respectively, that is: [INLINE:9] where T nw1 and T nw2 denote the total counts detected in the left and right narrow energy windows, respectively. Two narrow windows with the equal widths (w nwl = w nw2 = 2 keV), centered on E 1 = 126 keV and E 2 = 154 keV, respectively, are used for this correction method. Three energy window (TEW) method using triangular approximation In this method [10] (referred to as TEW2), similar to trapezoidal approximation, two narrow energy windows located on both sides of the photopeak window are used for the scatter correction. However, instead the trapezoidal area, the spectrum of the scattered counts of the photopeak window is estimated by the area of a right triangle that its height is equal to the estimated scatter counts in the left narrow energy window centered on energy E 1 = 126 keV. In this method, it is assumed: There isn't any scattered photon in the narrow window centered on the upperedge energy of the photopeak window (E 2 = 154 keV), and therefore the photons detected in this window are only unscattered photons: [INLINE:10] Simulation For evaluation of the different scatter correction methods, a coldsphere hotbackground phantom containing six waterfilled spheres (labeled as S1S6 in [Figure 1]) with diameters 3.2, 2.6, 2, 1.6, 1.2, and 1 cm placed in a 99m Tc uniformly filled cylindrical phantom (diameter 20 cm and height 22 cm) is simulated using SIMIND Monte Carlo simulation program. The centers of these spheres are located at the plane passing from the halfheight of the cylindrical phantom and the radial distance of them from the axis of the cylinder is equal to 5.4 cm [see [Figure 1].{Figure 1} In this study, the simulated SPECT system has been equipped with a lowenergy highresolution (LEHR) collimator. One hundred and twentyeight projections (matrix size 128 × 128 and pixel size 0.3 cm) are acquired at equidistant angles in a 360° camera rotation with a radius of rotation 20 cm that is same for all of the simulations in this study. Using the SIMIND Monte Carlo simulation, the projections related to the different energy windows required for the six scatter correction methods are generated. For any energy window sitting, the SIMIND simulation is able to create the scatter projection images containing only the scattered photons, as well as the total (scatter + unscatter) projection images. The unscatter projections can be obtained by subtracting the scatter projections from the total projections. The scatter correction methods are programmed in MATLAB environment. The corrected projections are reconstructed by filtered backprojection method using Hann filter. Finally, the reconstructed images from the different scatter correction methods are evaluated using assessment criteria described below. Assessment criteria Before the use of the assessment criteria, the regions of interest (ROIs) used for the calculation of these criteria are must defined. In this study, the ROI for each cold sphere is defined as all of the pixels that have located entirely inside and far from the edges of the circle of the sphere. Based on this definition, the ROIs for cold spheres with diameters 3.2, 2.6, 2, 1.6, 1.3, and 1 cm consist of 52, 30, 16, 12, 6, and 2 pixels, respectively. For background, the ROI consists of 256 pixels ranging from row 57 to 72 and column 57 to 72. Image contrast For each cold sphere, the image contrast is defined as deference between the mean of counts in the ROI of the sphere, th N c , and the mean of counts in the ROI of the background, th N b , divided by the latter: [INLINE:11] Relative noise of the background (RNB) The RNB is defined as the ratio of the standard deviation (SD) in the background's ROI, δb , to the mean of counts in this ROI, that is: [INLINE:12] Signaltonoise ratio (SNR) According to two definitions given above, for each cold sphere, the SNR is defined as the ratio of the image contrast to the RNB, that is: [INLINE:13] The program related to each assessment criterion is also written in the MATLAB environment, and the evaluation and comparison between the different correction methods is performed based on these criteria. Results In this section, we describe the results obtained from each of the scatter correction methods individually. DEW method The essential step in the DEW method is the calculation of the k value. This value can be calculated by pixelbypixel dividing the photopeak scatter projections by the corresponding projections acquired through the Compton window. A typical projection's ROI (33:96, 30:99), covering the entire area of the cylindrical phantom in the projection, is used for the calculation of the k value. Since the k values calculated for the different pixels in a projection image are somewhat difference, therefore, it is necessary to average the k values calculated for the different pixels and projections. The k values obtained in this study have been ranged from the minimum value 0.1615 to the maximum value 0.7444 with a mean value 0.4252 and a SD of 0.0515. Also, a linear fitting between the pixel counts of the photopeak scatter projection, S pk (i,j), and the pixel counts of the Compton projection, T C (i,j), has been shown in [Figure 2].{Figure 2} The results of the simulation indicate that 98.6% of the photons detected in the Compton window are scattered photons and 44.4% of them undergo multiple scattering; whereas, only 13.1% of scatter photons detected in the photopeak window are multiplescattered photons. Also, the reconstructed images from the true scatter component in the photopeak window and the scatter component estimated by the DEW method indicate that there are some differences between these two images [see [Figure 3]. The estimated scatter component image is approximately uniform so that it is not possible to detect the location of the spheres; whereas, the true scatter component image is noisy and the location of the spheres is somewhat detectable. Therefore, the subtraction of estimated scatter projections from corresponding photopeak projections results in some errors in the corrected images.{Figure 3} The results obtained show that the DEW method improves the image contrast of cold sphere 1, 2, 3, 4, and 5 about 22.6, 21, 17, 10.5, and 6%, respectively [see [Table 1]. The RNB for the uncorrected and the corrected images is 0.0312 and 0.0454, respectively that show the scatter correction by the DEW method results in the increase of the image noise. Also, the SNR, for cold sphere 1, 2, and 3 is slightly lower and for cold sphere 4 and 5 slightly higher than the uncorrected images [see [Table 2].{Table 1}{Table 2} PEDA method Based on the description offered for this method, two subwindows with the approximately equal scatter counts for any pixel must be determined within the photopeak window. The consideration of the different subwindows indicates that the optimal cutoff energy between two subwindows is 133.5 keV. Therefore, the subwindows w 1 = 126133 keV and w 2 = 134154 keV with the number of scattered photons 436 and 444, respectively are chosen as most appropriate for this purpose. The mean value, SD, minimum and maximum values of the ratio S 1 (i,j) to S 2 (i,j) for seven projections in the projection's ROI (33:96, 30:99) are given in [Table 3]. 47.7% of detected photons in the lower subwindow (w 1 ) are the unscattered photons, including ~11.3% of the photopeak unscattered photons. Therefore, the subtraction in the PEDA method removes some unscattered photons that results in the significant reduction of the SNR in the corrected image.{Table 3} The PEDA method improves the image contrast of cold sphere 1, 2, 3, 4, and 5 about 26, 19, 12, 8.5, and 8.6%, respectively. The RNB for the corrected images is 0.0543. The SNR, for all of the cold spheres except for sphere 5 is lower than the uncorrected images. DPW method Using the photopeak scatter and scatterfree projections generated by SIMIND Monte Carlo simulation, it is possible to calculate the ratio of the photopeak scatter counts to the photopeak unscatter counts (scatter fraction) for each pixel, SUR(i,j). The plot of the SUR(i,j) against R(i,j) for a typical projection has been shown in [Figure 4]. The values of parameters A, B, and C resulting from the nonlinear regression fitting for the projection's ROI are −0.276, −1.429, and 0.4595, respectively.{Figure 4} As shown in [Figure 5], for a typical projection, the majority of the pixel values of the calculated SUR in the regions near to the edges of the phantom are higher than that of the true SUR (overestimation); whereas in the central regions of the phantom, the state is vice versa (underestimation). It is also obvious that the amount of variation among the pixel values of the calculated SUR (SD ~0.02) is lesser than that of the true SUR (SD ~0.05). Therefore, it seems that the nonlinear regression fitting [Eq. 5] cannot estimate the scatter fraction of the photopeak window correctly.{Figure 5} From the data of [Table 1], it is obvious that the improvement in the image contrast of the cold spheres is negligible so that it is only 4% for the largest sphere. As shown in [Figure 6], the structure of the reconstructed image from the photopeak scatter component estimated by the DPW method is similar to the structure of the reconstructed image from the total counts of the photopeak window, and therefore the subtraction in the DPW method cannot reduce the number of counts within the sphere's ROI significantly. Thus, it can be understood why the DPW method does not improve the image contrast greatly. On the other hand, the SNR for all of the cold spheres has been increased [see [Table 2] that shows the DPW method results in reduction of variation among pixel values and thereby reduction of SD in the background's ROI. For this reason, the RNB for the corrected image (RNB ~0.03) is slightly lower than that for the uncorrected image.{Figure 6} CRM The first step for use of the CRM is the calculation of the parameters G and H. The mean value, SD, minimum and maximum values calculated for these two parameters in three different ROIs have been given in [Table 4]. The data of this table show that the SD of the parameter H is relatively big showing the big variations in the H values from a pixel to another pixel that can lead to significant errors in the corrected images using the mean H value.{Table 4} As the previous scatter correction methods, the mean values of G and H calculated for first ROI (33:96, 30:99) are used for the scatter correction. From [Table 1], it can be understood that the use of CRM leads to improve the image contrast of cold sphere 1, 2, 3, 4, and 5, about 24.5, 13.9, 10, 11, and 4.9%, respectively. The improvement in the image contrast of sphere 4 is somewhat larger than that of sphere 3 that shows the nonuniform performance of this scatter correction method. The RNB for the corrected images is 0.0538 that causes the SNR, for all of the cold spheres except for sphere 5, become lower than the uncorrected images. TEW methods Since the SIMIND can generate the projections containing only the scattered photons in the photopeak window, it is possible to study the true scatter spectrum of the photopeak window and thereby estimate the area of this spectrum using some approximate methods. The spectrum of the photopeak true scattered photons along with the estimated spectra using two approximations, trapezoidal and triangular approximations, is shown in [Figure 7]. The number of scattered photons, corresponding to the area under the true total scatter spectrum of the photopeak window is equal to 937.{Figure 7} TEW1. The number of scatter photons estimated using the trapezoidal area defined by this approximation is 1,323.7 representing an overestimation about 41.27%. This overestimation results in the high image noise (RNB = 0.0727), and thereby the significant reduction of the SNR for all of the cold spheres in the corrected images as compared to that for the uncorrected images. The improvement in the image contrast resulting from this method for the sphere 1, 2, 3, 4, and 5 is about 24.6, 17, 4, 14, and 2.7% respectively, showing a nonuniform correction similar to the CRM. TEW2. Based on the narrow windows used in this approximation, 4.6% of the photons detected at 154 keV are the scattered photons and the relative difference between the number of unscattered photons at 124 keV and 154keV is 10.67%. The number of scatter photons estimated using this triangular area is 838.25 representing an underestimation about 10.54%. The RNB for the corrected images is 0.0412. The improvement in the image contrast resulting from this method for the sphere 1, 2, 3, 4, and 5 is about 24, 18, 14, 10, and 4.6% respectively. The SNR, for all of the cold spheres in the corrected images is higher than the uncorrected images. Conclusion In this study, the six scatter correction methods based on setting two or three energy windows in 99m Tc spectrum were evaluated using the three criteria. Among these correction methods, two methods DEW and CRM involve the calculation of the mean values for one (k) or two parameters (G and H). Since each pixel in the projection image has an individual value for these parameters, that is, lower or higher than the calculated mean value, therefore the use of these mean values results in undercorrected or overcorrected images. This is important especially for parameter H with a relatively big SD. On the other hand, the PEDA and DPW methods need to determine an optimal cutoff energy between the two subwindows and the values A, B, and C for the regression relation between R(i,j) and SUR(i,j), respectively. According to the results obtained, it is not possible to determine the cutoff energy in the PEDA method so that the number of scattered photons for all of the pixels in the scatter projections of the two subwindows is closely equal. Also, the regression relationship obtained for the DPW method cannot approximate the relation between R(i,j) and SUR(i,j) correctly. Finally, the TEW methods do not need to calculate any additional parameter, and therefore the correction can be implemented directly on the pixels of projections. But these methods also leads to over or under estimate the scatter counts in the photopeak window. Thereby, it can be concluded that none of the six scatter correction methods investigated can estimate the number of scattered photons in the photopeak window closely. The comparison of the image contrast improvement resulting from the six scatter correction methods [see [Table 1] shows that: Highest image contrasts are resulting from: The PEDA method for cold sphere 1 (an increase of 26%) and cold sphere 5 (an increase of 8.6%),The DEW method for cold sphere 2 (an increase of 21%) and cold sphere 3 (an increase of 17%), andThe TEW1 method for cold sphere 4 (an increase of 14%). For the DPW method, contrast improvement of the five cold spheres is lowest.The CRM and TEW1 method perform a nonuniform correction (the contrast improvement of cold sphere 3 is lower than that of cold sphere 4), this is clear specially for latter. The comparison of the SNR resulting from the six scatter correction methods [see [Table 2] shows that: For the two methods TEW1 and DPW, the SNR of the five cold spheres is lowest and highest respectively.The SNR of the five cold spheres resulting from the DEW and TEW2 is very close to the SNR for uncorrected images. Plots of the pixel counts of the true unscatter projection against the pixel counts of the projection corrected by different scatter correction methods have been shown in [Figure 8]. From the Rsquare of these plots, it is obvious that the DEW and TEW1 methods indicate the best and worst match between the estimated US(i,j) and the true US(i,j).{Figure 8} From the results obtained in this study, among the six correction methods evaluated, the TEW method using triangular approximation (TEW2), because of the ease of implementation and uniform performance, the good improvement of the image contrast, the low RNB, the SNR higher than the uncorrected image, and the good match between the true unscatter counts and the corrected counts, is suggested as the most favorite scatter correction method. However, if a presupposed k value (for example originally suggested k value of 0.5) is used for the DEW method, this method will be also an appropriate choice for scatter correction. References


