Evaluation of six scatter correction methods based on spectral analysis in (99m)Tc SPECT imaging using SIMIND Monte Carlo simulation.

Compton-scattered photons included within the photopeak pulse-height 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 cold-sphere hot-background phantom. For evaluation of different scatter correction methods, three assessment criteria including image contrast, signal-to-noise ratio (SNR) and relative noise of the background (RNB) are considered. Except for the dual-photopeak window (DPW) method, the image contrast of the five cold spheres is improved in the range of 2.7-26%. 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 9-10% 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 Compton-scattered 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.[456789101112] Some of them are based on a spectral analysis that is performed by setting additional energy windows in 99mTc 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 dual-energy window (DEW) method, photopeak energy distribution analysis (PEDA) method, dual-photopeak window (DPW) method, cannel ratio method (CRM), and three energy window (TEW) methods (trapezoidal and triangular approximation). A cold-sphere hot-background 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, signal-to-noise 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 126-154 keV energy window centered on 140 keV is used as photopeak energy window (Wpk ).

Dual-energy window (DEW) method

This method[45] relies on the assumption that the spatial distribution of the counts acquired in a secondary window placed in Compton region of the energy spectrum (WC = 92-125 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:

where Spk is the scatter projection obtained from the photopeak window, TC 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, USpk , can be deduced by:

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:

where S1 and S2 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:

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 cut-off 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.

Dual-photopeak 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 (wlw ) to the upper subwindow (wuw ), R, as shown in the following equation:

Then, the number of scatter photon events in the photopeak window can be estimated as:

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:

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:

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 upper-edge energies of the photopeak window divided by the width of the window, respectively, that is:

where Tnw1 and Tnw2 denote the total counts detected in the left and right narrow energy windows, respectively. Two narrow windows with the equal widths (wnwl = wnw2 = 2 keV), centered on E1 = 126 keV and E2 = 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 E1 = 126 keV. In this method, it is assumed:

There isn’t any scattered photon in the narrow window centered on the upper-edge energy of the photopeak window (E2 = 154 keV), and therefore the photons detected in this window are only unscattered photons:

The photopeak window is symmetric around the energy E0 (140 keV). Therefore,

From these assumptions, it can be concluded that

Therefore, the number of scattered photons in the photopeak window can be estimated using following equation:

The width of the lower and upper narrow energy windows used for this method consist of wnw1 = 6 keV and wnw2 = 8 keV.

Simulation

For evaluation of the different scatter correction methods, a cold-sphere hot-background phantom containing six water-filled spheres (labeled as S1-S6 in Figure 1) with diameters 3.2, 2.6, 2, 1.6, 1.2, and 1 cm placed in a 99mTc 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 half-height 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

(a) An illustration of the cross-section of the cold-sphere hot-background phantom used in this study. (b) A projection simulated by SIMIND Monte Carlo program. (c) A slice (slice 64) from reconstructed image by MATLAB software

In this study, the simulated SPECT system has been equipped with a low-energy high-resolution (LEHR) collimator. One hundred and twenty-eight 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 back-projection 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:

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:

Signal-to-noise 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:

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 pixel-by-pixel 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, Spk (i,j), and the pixel counts of the Compton projection, TC (i,j), has been shown in Figure 2.

Figure 2

A plot of the photopeak scattered counts, Spk(i,j), against the total counts acquired through Compton window, TC(i,j), for a typical projection image. Each point corresponds to a pixel count. The solid line corresponds to the linear fitting curve

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 multiple-scattered 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 reconstructed image resulting from the photopeak scatter component (left) and the estimated scatter component (right). These images indicate that the spatial distribution of the scattered photons detected in the photopeak window cannot be appropriately approximated by the spatial distribution of the photons in Compton window

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

Comparison of the image contrast resulting from six scatter correction method for five cold spheres labeled as S1-S5

Table 2

Comparison of the SNR resulting from six scatter correction method for the five cold spheres

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 cut-off energy between two subwindows is 133.5 keV. Therefore, the subwindows w1 = 126-133 keV and w2 = 134-154 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 S1 (i,j) to S2 (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 (w1 ) 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 mean value, standard deviation, minimum and maximum values of the ratio S1(i,j) to S2(i,j) for seven projections into the region of interest (ROI; 33:96, 30:99)

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 scatter-free 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

Plot of SUR(i,j) against R(i,j) for a typical projection. Each point corresponds to a pixel (i,j) in the projection image. The solid line corresponds to the nonlinear regression fitting SUR(i,j) = A[R(i,j)]B + C

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

Illustration of the pixel values of the calculated SUR (*) and the true SUR (O) for (a) the region of interest (ROI) near to the edge and (b) the ROI near to the center of the phantom in a typical projection

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

The reconstructed image of the photopeak scatter component estimated by the DPW method that indicates well the structure information of the cold-sphere hot-background phantom

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

The mean value, standard deviation, minimum and maximum values of the parameters H and G for a typical projection in three different regions of interest (ROIs). It should be noted that the maximum value of H is reduced significantly for the smaller ROIs

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

The spectrum of the scattered photons in the photopeak window along with the area estimated using the trapezoidal and triangular approximations. The trapezoidal approximation is an overestimation for the scattered photons in the photopeak window, whereas the triangular approximation is an underestimation

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 99mTc 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 cut-off 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 cut-off 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%), and

The 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 R-square 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

Plots of the true unscatter counts against the estimated unscatter counts of the photopeak window. Each point corresponds to a pixel (i,j) and solid line corresponds to the linear fitting between true US(i,j) and estimated US(i,j). The R-square is statistic parameter that measures how successful the fit is in explaining the variation of the data. A value closer to 1 indicates a better fit

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.