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:
-
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]
-
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)
-
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
-
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 18F-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.
