Abstract
An efficient method for tomographic imaging in nuclear medicine is PET. Higher sensitivity, higher spatial resolution, and more accurate quantification are advantages of PET, in comparison to SPECT. However, a high noise level in the images limits the diagnostic utility of PET. Noise removal in nuclear medicine is traditionally based on the Fourier decomposition of the images. This method is based on frequency components, irrespective of the spatial location of the noise or signal. The wavelet transform presents a solution by providing information on frequency contents while retaining spatial information, alleviating the shortcoming of Fourier transformation. Thus, wavelet transformation has been extensively used for noise reduction, edge detection, and compression. Methods: In this research, SimSET software was used for simulation of PET images of the nonuniform rational Bspline–based cardiactorso phantom. The images were acquired using 250 million counts in 128 × 128 matrices. For a reference image, we acquired an image with high counts (6 billion). Then, we reconstructed these images using our own software developed in a commercially available program. After image reconstruction, a 250millioncount image (noisy image or test image) and a reference image were normalized, and then root mean square error was used to compare the images. Next, we wrote and applied denoising programs. These programs were based on using 54 different wavelets and 4 methods. Denoised images were compared with the reference image using root mean square error. Results: Our results indicate stationary wavelet transformation and global thresholding are more efficient at noise reduction than are other methods that we investigated. Conclusion: Wavelet transformation is a useful method for denoising simulated PET images. Noise reduction using this transform and loss of highfrequency information are simultaneous with each other. It seems we should attend to mutual agreement between noise reduction and visual quality of the image.
PET is an efficient technique to determine 3dimensional distributions of radiotracers in a patient's body. The technique is used to map the biologic function and metabolic changes of the organs under investigation. PET has a good sensitivity and specificity in diagnosis and differentiation of malignant from benign tumors (1,2).
Although PET has made crucial progress, it nevertheless bears the main weakness of nuclear medicine: poor count density in images. No matter which organ is imaged, noise is always present in the nuclear medicine images and always causes error in quantification. Signaltonoise ratio, although considerably higher in PET than in SPECT, is yet much lower than in other tomography techniques such as CT and MRI. The inherent noise of the PET images considerably increases if improvement techniques such as scatter or randoms correction are applied (3).
Conventionally, finite impulse response filters are used to improve the signaltonoise ratio of nuclear medicine data (3). These filters are mainly of the lowpass and resolution recovery (adaptive) types. Butterworth and gaussian filters are common lowpass filters, and Metz and Wiener filters are the main resolution recovery filters (4–6).
Noise removal in nuclear medicine is traditionally based on Fourier decomposition of images. However, this method has 2 major drawbacks. The Fourier transform decomposes the data into a series of sine and cosine functions, representing the frequency components of the data. Because the sinusoidal functions are periodic and of infinite length, the object domain information (spatial information, in the case of an image) is ignored. Therefore, noise removal based on the Fourier transform affects exclusively the frequency components, irrespective of the spatial location of the noise or signal. The result is uniform distortion of the image, regardless of the local signaltonoise ratio. Moreover, this type of decomposition requires the frequency of the data to be uniform over the entire image, which is not always the case in nuclear medicine data (7).
An alternative method is the wavelet transform, which successfully overcomes these shortcomings. The wavelet transform provides a timefrequency representation of the signals. The wavelet was developed as an alternative to the shorttime Fourier transform because high frequencies are better resolved in time and low frequencies are better resolved in the frequency domain (8,9).
During the last decade, wavelet transformation has gradually been replacing the conventional Fourier transform in many applications. Wavelet transformation has several advantages over the conventional Fourier transform that are potentially advantageous in nuclear medicine. The main advantage is flexibility in the selection of base functions for decomposing the data. Although the Fourier transform uses just the sine and cosine functions, wavelet transforms may use an infinite set of possible base functions. Thus, wavelet analysis provides access to information that may be obscured by Fourier analysis. Another advantage of the wavelet is localization of the base functions in space (8,9). Although the sine and cosine functions have infinite length, the wavelets are quite short, creating an opportunity for analyzing signals that contain discontinuities and sharp spikes—a common situation in nuclear medicine images. Moreover, wavelet transformation also brings the possibility of extracting edge information, which provides essential visual cues in the clinical interpretation of images.
The wavelet has the ability to approximate an image with just a few coefficients independent of the original image resolution and, thus, makes possible the comparison of images of different resolutions.
These excellent features make the wavelet transform an exceptionally powerful tool to detect and encode important elements of an image in terms of wavelet coefficients. Adaptive thresholding of the coefficients corresponding to undesired components removes the noise from the image much more efficiently than through the conventional Fourier method. Wavelet denoising is now an accepted and widely used method of noise reduction (8,10,11) but has not yet received considerable attention in nuclear medicine. Regardless of some reports in favor of wavelet denoising (12), it is nevertheless on the waiting list of nuclear medicine.
MATERIALS AND METHODS
The nonuniform rational Bspline–based cardiactorso phantom was used to generate a torso of a typical human as the virtual object to be imaged (13,14). Adjustments of activity distribution in the phantom were based on ^{18}FFDG uptake in the organs of a healthy human. Adjustments of attenuation coefficients for the tissues in the phantom were based on the Zubal phantom (15,16) attenuation coefficients. The phantom was constructed in a 256 × 256 × 256 matrix corresponding to 1.6 mm^{3} voxels. To fully resemble a human torso, normal cardiac and respiratory motion was also considered in the creation of the phantom.
The SimSET PET simulator, version 2.6.2.6, was used to simulate a Discovery LS PET scanner (GE Healthcare). The SimSET package uses Monte Carlo techniques to model physical processes (17). Validation of the PET Monte Carlo simulator has been done already (18).
We simulated Monte Carlo 2dimensional PET using ^{18}FFDG as the radiopharmaceutical. The energy range according to 511 keV ± 30% was adjusted on 350–650 keV (the standard window) (19).
The imaging system was adjusted to fully cover the cardiac and liver regions, as these organs are the main subjects of the PET imaging. The total number of slices after reconstruction was 20. Intentionally, the simple backprojection technique was used to reconstruct the images in crude format without any manipulation.
The simulation first generated a history of 6 billion photons (6,000 million photons) to create the noisefree data to be used as a reference image. Then, the simulation was repeated generating 250, 500, 1,000, 1,500, 2,000, 2,500, and 3,000 million photons as test images of different signaltonoise ratios.
Images were simulated with the phantom representing a person under normal conditions: the lengths of the beating heart cycle and of the respiratory cycle were adjusted to 1 and 5 s, respectively. Maximum diaphragm motion in normal breathing was specified as 2 cm. The phantom length was 40 cm (−20 to +20), and we adjusted scan range to 14 cm (−2 to +12) to display the entire heart and liver. The number of phantom slices was specified as 256, and the number of image slices after reconstruction was specified as 20. We acquired 2 images: one with 6 billion counts as a reference image and another with 250 million counts as a noisy image.
MATLAB software (The MathWorks) was used for programming. In the first step, we reconstructed simulated PET images using our programs. Before denoising, reference and noisy images were compared using the line profile and the root mean squared error (RMSE). Using the line profile, we displayed pixel values in the specific row of the image for both the reference and the noisy images. This line profile can show noise in the image approximately. We plot pixel indices along the xaxis and pixel values along the yaxis. RMSE can indicate the difference between the images. The range of RMSE is from zero to infinity. The best value, zero, is achieved if the 2 images are quite similar. We calculated the RMSE value between the reference and noisy images and after denoising between the reference and denoised images.
Images were denoised by 4 different methods of wavelet denoising: singlelevel discrete wavelet, which transforms in 2 dimensions (DWT); singlelevel discrete stationary wavelet, which transforms in 2 dimensions (SWT); global thresholding, which uses a positive real number for uniform threshold; and leveldependent thresholding, which thresholds according to the decomposition level in 2 methods—hard thresholding and soft thresholding. In all methods, the noisy image was decomposed to the approximation, horizontal detail, vertical detail, and diagonal detail. Then, the image was reconstructed using the techniques explained in the “Results” section.
Fiftyfour wavelets (Haar [1 wavelet], Daubechies [9 wavelets], Symlets [8 wavelets], Coiflets [5 wavelets], BiorSplines [15 wavelets], ReverseBior [15 wavelets], and DMeyer [1 wavelet]) were used in each of the denoising methods. The optimum wavelet in each method was determined in terms of minimum RMSE between the denoised test images and the reference image. All calculations were performed using software developed in MATLAB, version 7.1.
The test images were decomposed using wavelet transformation into the approximation, horizontal detail, vertical detail, and diagonal detail. Then images were reconstructed again using different combinations of the approximation and details or using different thresholding methods.
We compared the noisy and denoised images globally (RMSE) and locally using profiles. These are the methods conventionally used in nuclear medicine. There are methods that are more complicated but are not suitable for nuclear medicine. Identical profiles were drawn on the images. In simulated images, the reproducibility is high. Subtraction of 2 profiles may be helpful but is not necessary.
RESULTS
In the first method (Table 1; Fig. 1), DWT was used for denoising, and then we reconstructed images using 4 procedures: only with approximation, with approximation and horizontal detail, with approximation and vertical detail, and with approximation and diagonal detail. We would lose important information in the reconstructed image if approximation were eliminated. To show the change, we calculated relative RMSE (RMSE between the reference and denoised images divided by the RMSE between the reference and noisy images). Under the best conditions, the relative RMSE is about 0.31 when we used solely approximation. Adding vertical or horizontal detail to the approximation gave a relative RMSE of about 0.58, showing no significant differences between the choices.
Figure 1 showed that image quality improved after denoising, but as the line profiles indicate, noise reduction was more significant when only approximation was used as the reconstruction procedure. The best wavelets in this method of denoising are Daubechies.
Then, the image was reconstructed using other procedures: with approximation, horizontal detail, and vertical detail; with approximation, horizontal detail, and diagonal detail; and with approximation, vertical detail, and diagonal detail. Afterward, we compared these reconstructed images with the image that was reconstructed using only approximation. The result was similar to that of the previous procedure.
In the second method (Table 2; Fig. 2), SWT was used for denoising. We did decomposition to 3 levels and we used only approximation for reconstruction, because the RMSE value increased when we added horizontal, vertical, or diagonal detail to the approximation. The minimum RMSE values were observed when only approximation was used for reconstruction. At levels 2 and 3 of decomposition, relative RMSE values were about 0.09 and 0.12, respectively, and the difference between them was not statistically significant.
Figure 2 indicates that at levels 2 and 3 of decomposition, noise reduction is more considerable than at level 1. Level 2 shows the smallest relative RMSE and is most appropriate. At level 1, the best wavelet in this method of denoising is Haar.
In the third method (Table 3; Fig. 3), a global threshold was used for denoising. In this method, a positive real number is used for uniform threshold. The same threshold is applied at all levels for all subimages. We compared the results at 3 levels of decomposition. Our results show that noise reduction is more significant at levels 2 and 3 of decomposition than at level 1. However, their difference is not statistically significant. The best wavelets in this method of denoising are Daubechies.
In the fourth method (Table 4; Fig. 4), a leveldependent threshold was used for denoising. This method uses different thresholds for each transformation level. The threshold is obtained using a wavelet coefficients selection rule based on the strategy of Birge et al. (20). We compared the results of hard and soft thresholding at 1 level of decomposition. Our results showed that leveldependent soft thresholding more efficiently reduces noise than does hard thresholding.
Moreover, four methods of leveldependent soft thresholding at 1 level of decomposition were compared. (These 4 methods are different from the 4 methods of denoising.) In these methods, the threshold values and numbers of coefficients were determined on the basis of the αparameter (21). Threshold values and numbers of coefficients to be kept for denoising correspond to the αparameter. The MATLAB default value for α equals 3. In other procedures, penalized threshold values were used. In the penalhi procedure, a high value of α is used (2.5 ≤ α < 10). In the penalme and penallo procedures, respectively, a medium value (1.5 < α < 2.5) and a low value (1 < α < 2) of α are used. The difference between the default procedure and other procedures was statistically significant, but the results of the penalhi, penalme, and penallo procedures were similar. The best wavelet in both methods of denoising (hard and soft) is Haar.
DISCUSSION
In this study, we compared 4 different methods for wavelet denoising in PET images. All these methods were effective for denoising, but their results were different.
The results of DWT indicated that this method improves the RMSE value. However, the best result was obtained by using only approximation. In nuclear medicine images, highfrequency noise exists. When we eliminate details that contain highfrequency information, noise reduces significantly. Adding one of the details to the approximation degrades RMSE, and this degradation is more pronounced for horizontal and vertical details.
Classic DWT is not shiftinvariant: the DWT of a translated version of a signal is not the same as the DWT of the original signal (11). SWTs, or undecimated algorithms, apply the lowpass and highpass filters without any decimation. The SWT is similar to the DWT except the signal is never subsampled and instead the filters are upsampled at each level of decomposition. The SWT is an inherently redundant scheme, as each set of coefficients contains the same number of samples as the input; therefore, for a decomposition of n levels there is a redundancy of 2n. Therefore, SWT is a timeconsuming method and needs more space in the computer for processing but is more accurate than DWT because it is a shiftinvariant method.
DWT looks much the same as the undecimated DWT (or SWT) except for downsampling by 2 and upsampling by 2. The downsampling by 2 is often referred to as decimation by 2. With the downsampling, the detail coefficients and approximation coefficients are each about half the length of the original signal.
We have to be careful using the DWT instead of the SWT for 2 main reasons: First, downsampling by 2 in the DWT can produce aliasing (throwing away half the samples can lead to false signals). Second, this transform is not shiftinvariant (sometimes called timeinvariant) (22,23).
In SWT, we studied the results of a reconstruction procedure that uses only approximation because when we add one detail to approximation for reconstruction, RMSE increases and image quality degrades. The degradation can be due to intense distortion, which happens after decomposition. In this method, the best result for RMSE is provided by the procedure that uses only approximation and the third level of decomposition. However, the visual quality of the image degrades because at the third level, in comparison to previous levels, more information is omitted.
It seems that in the reconstruction procedure using only approximation, the SWT is more effective than the DWT for noise reduction. At the first level of decomposition using SWT, relative RMSE is about 0.196 and is comparable to the best state of DWT.
In global thresholding, the best result for RMSE is obtained by thresholding at levels 2 and 3 of decomposition. However, the visual quality of the image is better at level 2. At level 3, noise is reduced but important edge information is also lost from the image.
Our results indicate that the global thresholding method is more effective than the leveldependent thresholding method of denoising. Using leveldependent thresholding, the relative RMSE is about 0.31 at best and is comparable to the relative RMSE in global thresholding at the third level. Therefore, uniform thresholding at all levels of decomposition is more efficient in denoising.
Of the methods that we examined on simulated PET images, the best denoising results were generated by SWT using only approximation at all 3 levels and global thresholding at levels 2 and 3. However, the reduction in noise using wavelet transformation was accompanied by a loss of highfrequency information. We should attend to mutual agreement between noise reduction and visual quality of the image. Therefore, SWT at level 1 has better results than the other methods when approximation is used for reconstruction.
Applying the methods to real patient data does not add much to the paper. Clinical application of the results will be useful but difficult and is not possible for us at this time.
CONCLUSION
Wavelet transformation is a useful method to reduce noise in PET images and therefore enhance the images. The characteristic of timefrequency localization of the wavelet helps us to shift and scale the signals and keep important information for analysis after transformation. Our results show that the resolution and contrast of images are almost not influenced by the wavelet denoising methods but that much noise is removed and image quality is improved.
Footnotes

COPYRIGHT © 2009 by the Society of Nuclear Medicine, Inc.