Comparison of noise estimation methods used in denoising ^99mTc-sestamibi parathyroid images using wavelet transform

Abstract

The objective of this study was to compare the performance of variance, median absolute deviation, and the square of median absolute deviation methods of noise estimation in denoising of ^99mTc-sestamibi parathyroid images using wavelet transform. Sixty-eight ^99mTc-sestamibi parathyroid images including 33 images acquired at zoom 1.0 and 35 acquired at zoom 2.0 were denoised using the wavethresh package in R. The image decomposition and reconstruction method discrete wavelet transform, wavelet filter db4, shrinkage method hard, and thresholding policy universal were used. The noise estimation in the process was made using var, mad and madmad functions, which use variance, mean absolute deviation, and the square of mean absolute deviation, respectively. The quality of denoised images was assessed both qualitatively and quantitatively. A nonparametric two-sample Kolmogorov–Smirnov test was applied to find whether the difference in image quality produced by these three noise estimation methods was significant at 95% confidence. Noise estimation using madmad function produced the best quality denoised image. Further, the quality of the denoised image using madmad function was significantly better than the quality of the denoised image obtained with var or mad function (P = 1). The estimation of noise using madmad functions in wavelet transforms provides the best-denoised image for both zoom 1.0 and zoom 2.0 ^99mTc-sestamibi parathyroid images.

Keywords

^99mTc-sestamibi parathyroid images - noise estimation methods - wavelet transform

Introduction

Nuclear medicine images are noisy. Several approaches have been proposed to filter the noise from nuclear medicine images.[1],[2],[3],[4],[5] Filtering in spatial domain and frequency domain are frequently used in nuclear medicine, which involves a trade-off between reducing the spatial resolution in the image and noise reduction. The application of filter also causes image blurring or edge removal along with the removal of noise from the image.[6] Although image filtering in the frequency domain is routinely used in nuclear medicine, recently image denoising in the wavelet domain is becoming popular.[7],[8],[9],[10],[11]

The wavelet transform of the image (signal plus noise) is very sparse-the image gets concentrated in a few wavelet coefficients (WCs) but the noise remains spread out. It is easy to separate the signal from noise by keeping large coefficients (which corresponds to true image) and delete the small ones (which correspond to noise) and then apply inverse wavelet transform to get back the image. However, this requires some a-priory information about the noise-level in the image.

To estimate noise in the image various statistical functions such as var (variance), mad (median absolute deviation), and madmad (square of median absolute deviation) functions have been proposed.[12] Variance (var) gives the idea about how “spread out” the data is. The variance represents the average squared deviation from the mean and is calculated using the formula:

The variance is not robust to outliers, that is, a few values that are separate from the main body of the data can increase the value of the statistic by an arbitrarily large amount. Median absolute deviation (mad) is a robust measure of the variability of a univariate sample of quantitative data. To compute the mad, we first compute the median, and then for each data-point we compute the distance between its value and the median. The mad is defined as the median of these distances. Thus mad is not sensitive to the presence of outliers. The square of median absolute deviation (madmad function) is equal to the square of the median absolute deviation (mad) function.

In this study we have estimated noise-level using var, mad, and madmad functions and compared the performance of these function in denoising ^99mTc-sestamibi parathyroid image in wavelet domain.

Methods

A total of 68 images acquired between 18 July 2016 and 09 August 2017 were selected, of which, 33 images were acquired with zoom 1.0 and 35 images were acquired with zoom 2.0. The images were acquired on Siemens Symbia T6 (Siemens Medical Solutions, Illinois, USA). Symbia T6 is a dual head gamma camera with single photon emission computed tomography (SPECT) computed tomography. Planar images (included in this study) of the neck and mediastinum were the images acquired at 2 h post intravenous administration of ^99mTc-sestamibi (500MBq-700MBq) with low-energy high-resolution collimator. However, the routine protocol used in the department for ^99mTc-sestambi images is the acquisition of images at 15 min and 2 h postadministration of ^99mTc-sestamibi. The digital data were acquired for 700 Kilo counts in a 256 × 256 matrix, with zoom 1.0 and zoom 2.0.

The procedure used for denoising ^99mTc-sestamibi image:

1. Compute the two-dimensional (2D) discrete wavelet transform (DWT) of the image: The 2D DWT algorithm is essentially the application of many 1D filters. First, the columns are passed through the smoothing (H) and bandpass (G) filters and the rows of each of these resultant images are again passed through each of G and H, this results in four images. Three of them, genioglossus, geniohyoid, and hyoglossus correspond to the highest resolution WCs. The HH image is a smoothed version of the original and can be further attacked in exactly the same. After each attack, the dimension of the images is halved due to subsampling.[13] The 8-levels wavelet decomposition of a 256 × 256 ^99mTc-sestamibi parathyroid image using db4 wavelet filter is shown in [Figure 1a]

2. Estimate the noise level from WCs using var, mad and madmad functions: The noise level was estimated from 3rd level to 8th level of decomposition of details of WCs. The estimated value of noise was used for thresholding WCs. The universal policy that computes a threshold based on Donoho and Johnstone's “universal thresholds.” The threshold is sqrt (2*log (n)) *noise, where noise is equal to the noise estimated from the above-mentioned three functions, and n is the number of data points (or number of WCs)

3. Threshold the WCs (hard thresholding rule) using universal threshold policy: The universal policy uses the formula: Sqrt (2*log (n)) *noise, where noise is equal to the noise estimated from the above-mentioned three functions, and n is the number of data points (or number of WCs), to compute the threshold. Hard threhsolding rule compares a WC with the threshold, if WC is larger in absolute magnitude it is left alone, if WC is smaller it is set to zero

4. Apply 2D inverse DWT on the modified WCs to get the image in the spatial domain.

Figure 1 (a) Wavelet coefficient, (b) threshold coefficient using var, (c) threshold coefficient using mad, (d) threshold coefficient using madmad

The above mentioned procedures were performed using the wavethresh package.[12]

Performance comparison

The following five methods were used to assess the quality of denoised images: (1) nonreference image quality evaluator called BRISQUE score.[14] A smaller score indicates better perceptual quality. The BRISQUE score is usually in the range (0, 100), (2) quantifying the difference between the input and denoised image (SumofSquareDiff: Sum of square difference between the input image and its corresponding denoised image);[15] (3) visual inspection of denoised image and its input image placed side by side; (4) visual inspection of residual image and histogram of residual image: The residual image is defined as the difference between the original (always slightly noisy) image and its denoised version. If a denoising method performs well, the residual image must look like a noise even with nonnoisy images and should contain as little structures as possible.[16],[17] and (5) visual inspection of the image as a result of the local Pearson correlation coefficient test between the denoised image and its residual image.[17] A 7 × 7 sliding windows was used to compute the local correlation.

Statistical analysis

For the purpose of increase in perceptual quality for a given test image Xi and its corresponding enhanced image Yi, one expects that

The following hypothesis was evaluated using a non-parametric two-sample Kolmogorov-Smirnov (KS) test as:[18]

Software package

For reading, writing, and plotting images we have used EBImage[19] and Imager package.[20] Base R functions were used for statistical analysis.[21] MATLAB R2019b was used for calculating BRISQUE score. The choice of these software packages was based on the availability of well debugged built-in functions for the accomplishing the required tasks in this study.

Results

The human eye is the only one able to decide if the quality of the image has been improved by the denoising method. [Figure 2] and [Figure 3] show the results of an experiment on ^99mTc-sestamibi zoom 1.0 and zoom 2.0 images, respectively. Visually the quality of denoised image reconstructed using madmad method of estimation of noise is superior. The accurate reconstruction of edges, texture, and details and absence of artifacts can be seen in case of madmad function in comparison to the other two methods of noise estimation. [Figure 2] and [Figure 3] can be zoomed in to verify these findings. It may be re-iterated that all other variables are same except the method of estimation of noise in the images so that one can compare the visual quality of the denoised images, the non-presence of artifacts, and the correct reconstruction of edges, texture, and details.

Figure 2 Denoising experience with zoom 1.0 images. (a) noisy image zoom 1.0, Brisque score = 38.19, (b) denoised image reconstructed using var, Brisque score = 82.80, (c) denoised image reconstructed using mad, Brisque score = 84.16, (d) denoised image reconstructed using madmad, Brisque score = 43.54. The removed or distorted details must be compared with the corresponding residual images shown in Figure 5

Figure 3 Denoising experience with zoom 2.0 images. (a) noisy image zoom 2.0, Brisque score = 48.61, (b) denoised image reconstructed using var, Brisque score = 89.21, (c) denoised image reconstructed using mad, Brisque score = 93.50, and (d) denoised image reconstructed using madmad, Brisque score = 69.10. The removed or distorted details must be compared with the corresponding residual images shown in Figure 4

The visual inspection of the residual image helps to understand the performance and limitations of the denoising algorithm. If the denoising algorithm removes details or texture, then the residual image has details and structure (i.e., quantitatively large noise). In such case, the blurred or degraded structures of the denoised images coincide with the noticeable structures of its residual image [compare denoised image [Figure 3b] obtained with var function with residual image [Figure 4b] and also compare denoised image [Figure 3c] obtained with mad function with residual image [Figure 4c], where [Figure 4a] is the input image]. The comparison of denoised image [obtained with madmad function [Figure 3d] and residual image [Figure 4d] shows that best-denoised image was obtained when madmad function was used to estimate the noise because only with this function the residual image looks like Gaussian white noise. The comparison of denoised image [Figure 2] with its corresponding residual image [Figure 5] also shows that the performance of noise estimation using madmad function resulted in the best denoised image. We can see in [Figure 5] that the residual images corresponding to denoised images using var and mad function contain details or texture while the residual image corresponding to denoised image with madmad function does not have any noticeable geometrical structures and most of the regions appear as white noise except removal of details corresponding to the salivary gland. The histogram of the residual image corresponding to madmad function appears to follow normal distribution [Figure 6i] however [Figure 6c] and [Figure 6f] do not follow normal distribution.

Figure 4 Denoised image shown in Figure 2. (a) input image Figure 2a, (b) residual image for Figure 2b, (c) residual image for Figure 2c, (d) residual image for Figure 2d

Figure 5 Input noisy image zoom 2.0 and residual images for denoised image shown in Figure 3. (a) Input image Figure 3a, (b) residual image for Figure 3b, (c) residual image for Figure 3c, (d) residual image for Figure 3d

Figure 6 Input image was same for creating both residual and correlation image. (a) Residual image var, (b) Pearson correlation of denoised and residual image Figure 6a, (c) Histogram of the residual image Figure 6a, (d) Residual image mad, (e) Pearson correlation of denoised and residual image Figure 6d, (f) Histogram of the residual image Figure 6d, (g) Residual image madmad, (h) Pearson correlation of denoised and residual image Figure 6g, (i) Histogram of the residual image Figure 6g

[Figure 6b], [Figure 6e] and [Figure 6h] show the result of a correlation coefficient test between denoised image and residual image reconstructed using the var, [Figure 6a] mad [Figure 6d] and madmad [Figure 6g] function. Here “white” represents the rejection of independent hypothesis and “black” represents acceptance of independent hypothesis at α = 0.05. From the result of the Pearson correlation coefficients test also, madmad function was found to be the best as the residual image in the case of madmad function were found to be independent of the denoised image [Figure 6h].

The residual was quantified by finding the sum of squares of the difference between the input image and denoised image (named as “SumofSquareDiff”). [Figure 7] displays the boxplot of “SumofSquareDiff” for the residual corresponding to var, mad and madmad functions for zoom 1.0 and zoom 2.0 images. The mean of the “SumofSquareDiff” for madmad function was found to be significantly less than that of both var and mad function (at P = 1, KS test) for zoom 1.0 and zoom 2.0 images. The smallest median value of “SumofSquareDiff” for madmad function also supports the best performance of madmad function in estimating noise for denoising ^99mTc-sestamibi parathyroid images using wavelet transform.

Figure 7 (a) Boxplot of Brisque score for zoom 2.0 images, (b) Boxplot of Brisque score for zoom 1.0 images, (c) Boxplot of SumofSquareDiff for zoom 2.0 images and (d) Boxplot of SumofSquareDiff for zoom 1.0 images

Based on BRISQUE score none of the denoised image series were perceptually better than the input image [Figure 7]. The minimum distortions were introduced in denoised image, in case of madmad function for both zoom 1.0 (P = 5.952861e-09, at α = 0.05 and zoom 2.0 images (P = 9.000494e-10, at α = 0.05.

Discussion

We have compared the performance of three methods for estimation of noise on denoising ^99mTc-sestamibi parathyroid images using wavelet transform. The performance comparison was based on the quality of the denoised image each method produced. Based on the five different methods of image quality assessment, the denoised image produced using madmad function was found to have the best image quality having smoothness with a pleasant visual appearance. The use of var and mad function resulted in over-smooth images with relatively large reconstruction error compared to madmad function.

The over-smooth image is the consequence of universal thresholding policy producing large threshold. The threshold is sqrt (2*log(n)) *noise, where noise is equal to the noise estimated from the above mentioned three functions, and n is the number of data points (or the number of wavelet coefficients). The comparatively large value estimated by var and mad function resulted in large thresholds.

In zoom 2.0 category images, mad function estimated the maximum noise and threshold [Figure 8a] and [Figure 8b], and hence yielded overly smoothed images with the largest reconstruction error and artifacts [Figure 3c]. A large threshold value (for mad function) makes more number of coefficients as zero, which leads to smooth signal and destroys details that may have caused blur and artifacts [Figure 3c]. A small threshold [for madmad function, [Figure 8a] and [Figure 8b] had made very small number of coefficient equal to zero [Figure 1d] for madmad function, compared to [Figure 1b] for var, and [Figure 1c] for mad function which resulted in a relatively less smooth image and less reconstruction error and artifacts [Figure 3d] in comparison to denoised image yielded by var function [Figure 3b] and mad functions [Figure 3c]. The same explanation holds good for denoised zoom 1.0 images in which var function estimated the maximum value of noise [Figure 8c] and threshold [Figure 8d] and hence smoothest image, whereas madmad function estimated minimum value of noise and threshold and hence the least smooth image.

Figure 8 (a) Boxplot of noise for zoom 2.0 images, (b) Boxplot of estimated threshold for zoom 2.0 images, (c) Boxplot of noise for zoom 1.0 image, and (d) Boxplot of estimated threshold for zoom 1.0 image

Our result is similar to the results reported.[12] Their recommendation is to estimate the noise using madmad function as this function provides a better estimate of noise. The mad function is the method of noise estimation from the finest resolution sub-band.[22] They have found mad function as a robust method for noise estimation. The noise estimation is different from Donoho study because we have estimated noise from 3rd level to 7th level details sub-band, applied the estimated threshold to the 3rd level to 7th level details sub-band, and also that we have compared the performance of the var, mad, madmad function for denoising ^99mTc-setsamibii parathyroid images.

The medical image denoising using wavelet transform have been first attempted by Waver et al.[23] They performed denoising of MRI images, and their results were encouraging except that the method eliminates small structures that were confused with noise. In nuclear medicine imaging, Ogawa et al.[9] performed denoising of Scintigraphic images, the method of decomposition and reconstruction were translational invariant wavelet transform. They proposed the algorithm to denoise scintigraphic images and demonstrated the effectiveness of their method by denoising brain phantom, staircase phantom, clinical ^99mTc-MDP Bone Scan, and67Ga-Scan images.

Nawres et al.[8] performed denoising of scintigraphic images on planar bone and heart images acquired over different duration of acquisition, consequently increasing count levels. A modified version of the Bayesian threshold was applied for threshold value estimations. Wavelet-based denoising has also been performed in nuclear medicine in the fields of SPECT and PET images. Afef et al.[10] used a 2D preprocessing Daubechies wavelet transformed for removal of Poison noise in the acquired projections and reconstructed the data with Ordered subset expectation maximization (OSEM) algorithm for a tomographic bone SPECT image reconstruction. Bal et al.[11] proposed a PET denoising technique on simulated phantom and ¹⁸F-FDG clinical images using wavelet and curvelet domains. The performance of the denoising method was compared with VisuShrink, BayesShrink, NeighShrink, and ModineighShrink. The algorithm efficiently denoised the Phantom images and the clinical images of Gaussian noise, Poisson noise, and Mixed Gaussian-Poisson noise.

Methods to reduce Poisson noise have also been proposed by many researchers,[24],[25] and these methods sometimes work well. Various image denoising techniques have been used to remove noise from scintigraphic images. These include linear filters[26] and order statistic filters such as a median filter in the spatial domain,[27] and Butterworth filter[28] and Wiener filter in the frequency domain,[29] deep convolution neural network,[30] median modified Wiener filter technique,[31] a blind-deconvolution framework after a noise-reduction algorithm based on a nonlocal mean[32] etc.

The significance of this study is that it clearly compares and suggests that madmad should be used for best results when denoising parathyroid Scintigraphic images, in comparison to var and mad.

There are many variables that affect the success of denoising in the wavelet transform domain. Exploration of all for different types of Scintigraphic images will be an enormous work. We have explored only the effect of method of noise estimation in this study.

Conclusions

The noise estimation using madmad function provides best denoised image for both zoom 1.0 and zoom 2.0 ^99mTc-sestamibi parathyroid images.

Conflict of Interest

There are no conflicts of interest.

Financial support and sponsorship

Nil.

• References

• 1 Pizer SM, Todd-Pokropek AE. Improvement of scintigrams by computer processing. Semin Nucl Med 1978;8:125-46.
• 2 Cinotti L, Meignan M, Usdin JP, Vasile N, Castaigne A. Diagnostic value of image processing in myocardial scintigraphy. J Nucl Med 1983;24:768-74.
• 3 Hannequin P, Mas J. Statistical and heuristic image noise extraction (SHINE): A new method for processing Poisson noise in scintigraphic images. Phys Med Biol 2002;47:4329-44.
• 4 Wesolowski CA, Yahil A, Puetter RC, Babyn PS, Gilday DL, Khan MZ. Improved lesion detection from spatially adaptive, minimally complex, Pixon® reconstruction of planar scintigraphic images. Comput Med Imaging Graph 2005;29:65-81.
• 5 Takalo R, Hytti H, Ihalainen H. Adaptive autoregressive model for reduction of poisson noise in scintigraphic images. J Nucl Med Technol 2011;39:19-26.
• 6 Gonzalez RC, Woods RE. Intensity transformations and spatial filtering. In: Digital Image Processing. 3rd ed. Delhi: Indian Edition Published by Dorling Kindersley India Pvt. Ltd., licensees of Pearson Education in South Asia; 2009. p. 126-220.
• 7 Nawres K, Kamel H, Noureddine E. Image denoising using wavelets: A powerful tool to overcome some limitations in nuclear imaging. In: Vol. 1. 2006 2nd International Conference on Information & Communication Technologies; 2006. p. 1114-8.
• 8 Ogawa K, Sakata M, Li Y. Adaptive noise reduction of scintigrams with a wavelet transform. Int J Biomed Imaging 2012;2012:130482.
• 9 Afef H, Bechir L, Dorra BS. A based daubechies wavelet transform approach to enhance the OSEM bone SPECT image reconstruction. CEIT-2017. Proceed Eng Technol PET 2018;33:94-100.
• 10 Bal A, Banerjee M, Sharma P, Maitra M. An efficient wavelet and curvelet-based PET image denoising technique. Med Biol Eng Comput 2019;57:2567-98.
• 11 Ali R, Yunfeng P, Amin RU. A novel bayesian patch-based approach for image denoising. IEEE Access 2020;8:38985-94.
• 12 Guy Nason. Wavethresh: Wavelets Statistics and Transforms. Package Version 4.6.8; 2016. Available from: https://CRAN.R-project.org/package=wavethresh. [Last accessed on 2020 May 11].
• 13 Mallat SG. A theory for multiresolution signal decomposition: The wavelet representation. IEEE Trans Pattern Analysis Machine Intelligence 1989;11:674-93.
• 14 Mittal A, Moorthy AK, Bovik AC. No-reference image quality assessment in the spatial domain. IEEE Trans Image Process 2012;21:4695-708.
• 15 Available from: https://earlglynn.github.io/RNotes/package/waveslim/modwt-2d.html. [Last accessed on 2020 May 11].
• 16 Buades A, Coll B, Morel JM. A review of image denoising algorithms, with a new one. Multiscale Modeling Simulation 2005;4:490-530.
• 17 Brunet D, Vrscay ER, Wang Z. The use of residuals in image denoising. In: Kamel M, Campilho A, editors. Image Analysis and Recognition. ICIAR. Lecture Notes in Computer Science. Vol. 5627. Berlin: Springer, Heidelberg; 2009. p.
• 18 Sheskin DJ. Handbook of Parametric and Nonparametric Statistical Procedures. 4th ed. London, U.K.: Chapman & Hall; 2007.
• 19 Pau G, Fuchs F, Sklyar O, Boutros M, Huber W. EBImage-an R package for image processing with applications to cellular phenotypes. Bioinformatics 2010;26:979-81.
• 20 Barthelme S. Imager: Image Processing Library Based on 'CImg'. R Package Version 0.41.1; 2018.
• 21 Core Team R. A language and environment for statistical computing. Foundation for Statistical Computing, Vienna, Austria. URL; 2018. Available from: https://www.R-project.org/. [Last accessed on 2020 May 11].
• 22 Donoho DL, Johnstone IM. Adapting to unknown smoothness via wavelet shrinkage. J Am Stat Assoc 1995;90:1200-24.
• 23 Weaver JB, Xu YS, Healy DM Jr., Cromwell LD. Filtering noise from images with wavelet transforms. Magn Reson Med 1991;21:288-95.
• 24 Wang L, Lu J, Li Y, Yahagi T, Okamoto T. Noise Reduction Using Wavelet with Application to Medical X-Ray Image. In2005 IEEE International Conference on Industrial Technology; 2005. p. 33-8.
• 25 Wang L, Lu J, Li Y, Yahagi T, Okamoto T. Noise removal for medical X-ray images in wavelet domain. Elect Eng Japan 2008;163:37-46, 22-27.
• 26 Gonzalez RC, Elwoods R. Digital Image Processing. Reading, MA: Addison-Wesley; 1993. p. 14-23.
• 27 Pandey AK, Sharma A, Sharma PD, Saroha K, Parida GK, Bal CS, et al. Denoising of iodine-131 images using a median filter. Nucl Med Commun 2019;40:308-16.
• 28 Pandey AK, Pant GS, Malhotra A. Standardization of SPECT filter parameters. IJNM 2004;19:30-5.
• 29 Pandey AK, Santosh Y, Sharma PD, Yadav D, Bal C, Kumar R. Restoration of I-131 whole body image using a Wiener filter. Nucl Med Commun 2020;41:426-35.
• 30 Minarik D, Enqvist O, Trägårdh E. Denoising of scintillation camera images using a deep convolutional neural network: A monte carlo simulation approach. J Nucl Med 2020;61:298-303.
• 31 Park CR, Kang SH, Lee Y. Median modified Wiener filter for improving the image quality of gamma camera images. Nucl Eng Technol 2020. Volume 52 Issue 10 / Pages.2328-2333.
• 32 Kim K, Lee MH, Lee Y. Investigation of a blind-deconvolution framework after noise reduction using a gamma camera in nuclear medicine imaging. Nucl Eng Technol 2020. [In press].

Publication History

Received: 12 May 2020

Accepted: 25 May 2020

Article published online:
30 March 2022