Speckle Noise Reduction in Medical Ultrasound Images Using Modelling of Shearlet Coefficients as a Nakagami Prior

The diagnosis of UltraSound (US) medical images is affected due to the presence of speckle noise. This noise degrades the diagnostic quality of US images by reducing small details and edges present in the image. This paper presents a novel method based on shearlet coefficients modeling of logtransformed US images. Noise-free log-transformed coefficients are modeled as Nakagami distribution and speckle noise coefficients are modeled as Gaussian distribution. Method of Log Cumulants (MoLC) and Method of Moments (MoM) are used for parameter estimation of Nakagami distribution and noise free shearlet coefficients respectively. Then noise free shearlet coefficients are obtained using Maximum a Posteriori (MaP) estimation of noisy coefficients. The experimental results were presented by performing various experiments on synthetic and real US images. Subjective and objective quality assessment of the proposed method is presented and is compared with six other existing methods. The effectiveness of the proposed method over other methods can be seen from the obtained results.


Introduction
UltraSound (US) images are used for clinical diagnosis of patients suffering from diseases related to internal organs of the human body.The real-time operation, non-ionization and non-invasiveness properties of this imaging modality make it popular among other imaging modalities, such as X-ray, CT-scan and MRI.US images are generated by the recording of received echo signal from internal body parts when US waves are transmitted inside the human body.In its generation process, US images are corrupted by a noise, which arises due to interference between transmitted and received signal.This noise is of multiplicative nature and appears in granular pattern in the formed US image popularly known as "speckle noise".This noise degrades the visible quality of US images and hence makes diagnosis difficult.Fine details and edges of US images are lost due to the presence of speckle [1].Hence, speckle noise reduction is an important aspect and may also be the preprocessing step of various image processing algorithms.Speckle denoising techniques are divided into two broad categories: • compounding techniques, • post-processing techniques.
Compounding techniques are quite old and are used during US image acquisition process whereas postprocessing techniques are applied after generation of US images.Post-processing methods are broadly classified as: • spatial domain methods, • transform domain methods.
In spatial domain techniques, arithmetic and logical operations are performed on the image pixels directly while in transform domain methods operations are performed on transformed coefficients obtained after operation of an appropriate transform on the image [2].A variety of spatial domain techniques is available in c 2018 ADVANCES IN ELECTRICAL AND ELECTRONIC ENGINEERING the literature for speckle noise reduction in US images [3], [4], [5], [6], [7], [8], [9], [10], [11], [12] and [13].Earlier methods include various spatial domain filters, such as mean [3], median [4], adaptive weighted median [5] and Wiener filter [6].Although, these filters show good denoising performance, but failed to preserve useful information and are also unable to differentiate the boundaries between gray levels where the difference between gray levels is not sufficient.Some popular spatial domain techniques, based on computation of local statistics of noisy image are Lee [7], Frost [8], Kuan filter [9], PMAD [10], Speckle Reducing Anisotropic Diffusion (SRAD) [11], Detail Preserving Anisotropic Diffusion (DPAD) [12] filters.The major issue associated with these filters is that the performance of these filters is dependent upon window size and number of iterations.Currently, hybrid spatial domain filters are also used for speckle noise reduction because of the availability of fast processors.A hybrid spatial domain filter is presented in [13], wherein a combination of SRAD, DPAD and Bilateral filter is used for speckle noise reduction in US images.Transform domain techniques are yet powerful for speckle denoising.These techniques are broadly classified as thresholding techniques and Bayesian estimation shrinkage techniques.In thresholding techniques, a proper threshold is chosen for retaining noise free coefficients from noisy coefficients.The challenges associated with thresholding based techniques are: • proper selection of thresholding function, • matching of distribution of signal and noise at various scales.
In Bayesian Shrinkage (BS) techniques, a prior distribution of transformed coefficients is to be assumed and original noise-free signal is estimated from noisy signal in Bayesian environment.Also, an accurate estimation of signal and noise characteristics is required to be achieved for both speckle suppression and edge preservation.Thus, the performance of these filters is mainly dependent on the correct selection of prior for modelling of transformed coefficients.A large number of BS techniques based on wavelets and its extensions are available in literature.
A Minimum Mean Square Error (MMSE) estimator based on symmetric alpha-stable as a prior for modelling of log-transformed US image is presented in [14].Heavy-tailed behavior has been shown when wavelet coefficients are statistically modelled [15].This is because of the non-Gaussian statistics exhibited by wavelet coefficients in different subbands.A Maximum Likelihood Estimation (MLE) approach based on modelling of wavelet coefficients of speckle noise as Rayleigh distribution is presented in [16].A similar approach but with log-transformed image was presented by Gupta et al. [17] wherein wavelet coefficients corresponding to log-transformed speckle noise are modelled as Rayleigh distribution.This method utilizes Maximum a Posteriori (MAP) criterion for estimating wavelet coefficients corresponding to signal from noisy data.Bhuiyan et al. proposed a MAP estimator [18] using Normal Inverse Gaussian (NIG) as a prior for modelling of wavelet coefficients corresponding to signal.This method presents a simple and fast approach for estimating the parameters of prior distribution by local statistics applied to the neighbourhood of wavelet coefficients.A method used to capture heavytailed marginal distribution of wavelet coefficients and their interscale dependencies based on modelling of 2D-GARCH Generalised Gaussian (2D-GARCH GG) is presented in [19].The approach is more flexible and uses MAP criterion for the estimation of noise-free coefficients.A Levy Shrink algorithm based on Dual-Tree Complex Wavelet Transform (DTCWT) is presented in [20] for speckle noise reduction in US images.In this technique, heavy-tailed Levy distribution is used as a prior model for modelling of wavelet coefficients corresponding to signal and noise.Parameters of levy distribution are estimated using fractional moments and noise-free signal is obtained using MAP estimation.A combination of MAP estimation and Laplace mixture prior is developed by Lu et al. [21].Other methods based on Cauchy prior model and bivariate Cauchy-Rayleigh mixture prior model are presented in [22] and [23].
Although, wavelet transform based techniques are very popular from last two decades for speckle noise reduction.However, the wavelets are useful for handling only point discontinuity.US images consist of lines and curves, so other popular transforms, such as Ridglet, Curvelet, Contourlet, Ripplet and Shearlet Transform were used for speckle noise reduction.All these transforms are extensions of Wavelets.Modelling of contourlet coefficients as a NIG prior and logtransformed contourlet coefficient as a Cauchy prior is presented in [24] and [25], respectively.A method based on MMSE and MAP estimators in the nonsubsampled contourlet domain is presented in [26].Also, a technique for speckle noise reduction based on the combination of diffusion filtering and MAP estimation of noiseless coefficients in Curvelet domain is presented in [27].
In this paper, a new transform domain method for speckle noise reduction in US images is presented.Discrete Shearlet Transform (DST) is used for the purpose of signal decomposition.Noise-free log-transformed shearlet coefficients are modelled as Nakagami distribution and speckle noise coefficients are modelled as Gaussian distribution.The major contribution of the presented work lies in the estimation of parameters of Nakagami distribution using Method of Log Cumulants (MoLC) from noisy observations.Thereafter, a method based on Bayesian framework popularly known as MAP estimation is employed to obtain speckle-free coefficients.Section 2. , presents the speckle noise model in medical US images.Section 3. , presents the brief description of DST which is used to obtain the transformed coefficients.In Sec.
4. , parameters of Nakagami distribution are estimated using MoLC and noise-free shearlet coefficients using absolute central moments.MAP estimation used for estimating noisefree shearlet coefficients is presented in Sec.
5. In Sec. 6. , results are obtained for the proposed method and are compared with other state of the art methods.Finally, conclusion of the proposed work is presented in Sec. 7.

Speckle Noise Model
Speckle noise in US images can be represented as Eq. ( 1) given below [1].
where, x is original US image which is corrupted by speckle noise component n to yield noisy US image y.
The additive noise can be obtained from multiplicative speckle noise by taking log transformation of noisy image.Suppose, f , g and n are log values of x, y and n, respectively.It can be represented as given in Eq. ( 2).
Then N-level Discrete Shearlet Transform (DST) is applied on the log-transformed noisy image to obtain the low pass and high pass subbands.For any arbitrary subband, the DST of the log-transformed image can be represented as given in Eq. (3).
S x , S y and S n are the DST of the log-transformed version of noiseless, noisy and noise component, respectively.All detail subbands are despeckled while approximation band remain untouched.The logtransformed speckle noise can be modelled as zeromean Gaussian distribution.

Discrete Shearlet Transform
Wavelet transform does not possess the directionality property due to the association of only scaling and translation parameters.The idea behind Shearlets is to provide directional information along with retaining useful information related to wavelets.Shearlets show an efficient representation of the images containing edges as compared to wavelets which are only efficient in dealing with point discontinuity.Discrete Shearlet Transform (DST) can be realized by combining Laplacian Pyramid (LP) and Shear Filter (SF) or Directional Filter (DF) [28].Figure 1 shows the twolevel decomposition of DST.Equation ( 4) presents the multidimensional representation of Shearlets.A parabolic scaling matrix A α,2 j is used to change the resolution of generating functions.The S k given in Eq. ( 5) is the shearing matrix used to change the orientation of parabolic scaling matrix.The value of α ∈ (0, 2). and The Discrete Shearlet Transform (DST) can be expressed as given in Eq. (6).
for, j = 0, 1..., J − 1 where Ψ d j,k , is a support function, p j is Fourier coefficient of 2-D filter and M c is the sampling matrix, chosen according to [29].

Modelling of Noise-Free Shearlet Coefficients
Noise-free shearlet coefficients of US image can be modelled as Nakagami distribution due to its flexibility to describe the amplitude distribution of noise.The probability density function of Nakagami distribution [30] is given in Eq. ( 8).
In Eq. ( 8), L is the shape parameter (L > 0), λ is the scale parameter (λ > 0) and Γ(.) is gamma function.Nakagami distribution is used for modelling right-skewed, positive datasets (x ≥ 0).From statistical properties of Nakagami distribution, MoLC based speckle reduction method can be appropriately derived on MAP estimation framework.The parameters of Nakagami distribution are estimated using MoLC and that of noise-free shearlet coefficients are estimated using absolute central moments.Then shearlet transformed coefficients corresponding to clean image are estimated using Bayesian MAP estimation.

Parameter Estimation of Nakagami Distribution
MoLC method is used for parameter estimation of Nakagami distribution from the observed noisy data.This method is widely used for parameter estimation because of its low variance and high computational efficiency.Absolute values of shearlet coefficients are used for parameter estimation as MoLC method works only with a random variable having positive values ( + ).MoLC method is based on Mellin transform defined as given in Eq. ( 9) below.
where φ x (s) is the first characteristic function of second kind and f x (x) is the density function of Nakagami distribution.The First characteristic function of second kind for Nakagami distribution is obtained using Eq. ( 9) and is given in Eq. ( 10) below.
The second characteristic function of second kind is obtained by taking the natural logarithm of the first characteristic function of second kind given in Eq. (11).
The n th order cumulants or log cumulants can be easily derived by taking the n th order derivative of Eq. ( 13) and can be represented as given in Eq. ( 14).
The first two log cumulants can be obtained by putting n = 1 and n = 2 in Eq. ( 14), respectively, and are given by Eq. ( 15) and Eq. ( 16) respectively.
Here, Ψ(n, L) is the Polygamma function.These first two log cumulants can also be estimated empirically from observed noisy data S y = S y,i |i ∈ (0, N 1 ) and are given in Eq. ( 17) and Eq. ( 18), respectively.Here, N 1 denotes the number of wavelet coefficients in i th subband.
Parameters L and λ corresponding to noisy shearlet coefficients can be estimated by putting log cumulants values in Eq. ( 15) and Eq. ( 16).

Parameter Estimation for Noise-Free Shearlet Coefficients
For parameter estimation of noise-free shearlet coefficients, it is assumed that the distribution model of noisy shearlet coefficients is also Nakagami distribution.Speckle noise is modelled as zero-mean Gaussian noise.Moments of Nakagami distribution are obtained by putting s = n + 1; in first characteristic function of second kind (Eq.( 11)), i.e.
resulting in (Eq.( 20)): Then 2 nd and 4 th order moments corresponding to noisy shearlet coefficients are obtained by putting n = 2 and n = 4 in Eq. ( 20), yield Eq. ( 21) and Eq.(22). and As noise is zero-mean Gaussian so absolute central moments for the noise-free shearlet coefficients can be calculated from noisy shearlet coefficients using the expressions given in Eq. ( 23) and Eq. ( 24): and where σ N is the noise standard deviation which can be estimated from the shearlet coefficients in different detail bands using robust Median Absolute Deviation (MAD) [31] as given in Eq. ( 25) below: 2 nd and 4 th order moments of noise-free shearlet coefficients can be estimated from Eq. ( 23), Eq. ( 24) and Eq.(25).As the distribution of noise-free shearlet coefficients is Nakagami distribution, so the ratio of their moments can also be computed employing Eq. ( 21) and Eq. ( 22) using S x = S y , yields Solving Eq. ( 26), parameter L of noise-free shearlet coefficients can be estimated.Then, parameter λ can be estimated in Eq. ( 27) using Eq. ( 21), i.e.

Goodness of Fit for Nakagami Distribution with Shearlet Coefficients
Figure 3 shows the histogram and estimated pdfs of tenth detailed band shearlet coefficients for synthetic kidney US image.It should be noted that three levels of 2-D DST are used for the purpose of image decomposition resulting in 32 detail bands and 1 approximation band.It can be seen that the data corresponding to shearlet decomposed detail band has been found to fit well with Nakagami distribution as compared to Gamma, Inverse Gaussian and Rayleigh Distributions.Also, for objective evaluation of data fitting with various distributions Kolmogorov-Smirnov (KS) statistic [32] is used.The KS distance obtained for Nakagami, Gamma, Inverse Gaussian and Rayleigh Distributions are presented in Tab. 1 at various decomposition levels.Table 1 also proves that the best-fitted distribution is Nakagami distribution as it gives the lowest value of KS distance at all decomposition levels as compared to the KS distances obtained from other distributions.

Proposed Method
The flow chart of the proposed method is presented in Fig. 4 Log-transformed noisy image is used as input image.2-D DST is applied on the input image to obtain detail bands and approximation band.Then modelling of all detail band shearlet coefficients is done assuming the distribution of shearlet coefficients as a Nakagami prior.Parameters of Nakagami distribution are estimated using MoLC method and then for noiseless coefficients using absolute central moment method.MAP estimation is then applied on the noisy shearlet coefficients to obtain the noise-free shearlet coefficients.Then 2-D inverse DST followed by exponential operation is applied on denoised shearlet coefficients to obtain the final denoised image.A Bayesian MAP estimator is designed to estimate the noise-free shearlet coefficients Sx . Also, Applying Bayes theorem on Eq. ( 30), we have Using Eq. ( 31) in Eq. ( 34), we get: where, f Sx (S x ) is the priori distribution of noise-free coefficients, assumed here Nakagami distribution and f Sn (.) is the probability density function of noise which is assumed to be Gaussian.Now substituting the value of f Sx (S x ) (after taking logarithm) in Eq. ( 35) to yield Eq. ( 28) To maximize function F (S x ) given in Eq. ( 36), derivative is the simplest way But Eq. ( 37) is not linear for S x .Therefore, another equation proposed by Hyvarinen [33] is taken.
Substituting the value of derivative given in Eq. (37) into Eq.(38) where S x = S y , resulting Eq. ( 29).Thus, estimated clean shearlet coefficients can be obtained using Eq. ( 29).'T ' is a tuning parameter used in the MAP expression as the noise variance of highfrequency subbands is variable.

Experimentation and Result Analysis
Experiments were performed on synthetic and real US images.Objective and subjective performance analysis of our proposed method is carried out against six existing methods viz.Wiener [6], Speckle Reducing Anisotropic Diffusion (SRAD) [11], Wavelet [34], Tab. 2: Comparison of various image quality metric for synthetic US image of kidney for noise variance σ 2 n = 0.1, 0.2, 0.3.BayesShrink [35], Soft thresholding [36] and SWT-Cauchy [37].For objective assessment of synthetic images, parameters Peak Signal to Noise Ratio (PSNR) [34], Signal to Noise Ratio (SNR) [34], Structure Similarity Index Measures (SSIM) [38], Figure of Merit (FOM) [39] and Correlation Coefficient (CC) [40] were used.Due to the non-availability of noise-free real US image as reference image, a commonly used parameter Mean to Variance Ratio (MVR) [41] is used for objective quality assessment of real US images.Subjective quality assessment of both synthetic and real US images was carried out using visual quality inspection.

Experiment on Synthetic US Image
The synthetic image of a kidney (469×522) was created using 'simulator Field II' software [42].Then speckle noise of variance σ 2 n = 0.1, 0.2 and 0.3 is introduced in synthetic US image using MATLAB.
Table 2 shows the different parameter values obtained at noise variance (σ 2 n ) = 0.1, 0.2 and 0.3.Results for Wiener filter are taken from the model available in Matlab for a window size of 5×5.The chosen parameters for SRAD filter for experiment are 100 iterations each for σ 2 n = 0.1, 0.2 for time step size (∆t) = 0.7, 1 and 1.2 respectively.For Wavelet filter, two levels (J = 2) of db4 wavelet decomposition is used and the reconstructed image is obtained after removal of all detail bands.In BayesShrink method denoising is performed using db4 wavelet with three level (J = 3) of decomposition.Soft thresholding method uses three levels (J = 3) of db1 wavelet decomposition for its operation.For the proposed method, results are taken on different values of tuning parameter T varying from T = 0.2 to T = 1 in a step size of 0.2 for all noise variances.The above-chosen filter parameter values are found to be optimal after performing several experiments on synthetic US images.The experiments were conducted for all the above-mentioned filters.It is also evident from Tab. 2 that the proposed method outperforms other methods in terms of parameter values obtained for a range of noise variance.It can be seen clearly that the noise is eliminated and the edges are also retained in the obtained denoised image using the proposed method.From subjective analysis of denoised images, it is evident that the quality of image obtained using Wiener filter is the worst.Almost, similar quality of images are obtained using SRAD and Wavelet filters but are better than the one obtained from Wiener filter.Also, similar performance is shown by both BayesShrink and soft thresholding methods and also the obtained denoised images are better than the one obtained from Wiener, SRAD and Wavelet filters in terms of fine detail preservation.The quality of image obtained using SWT-Cauchy method is the best among all existing methods but it is less pleasing as compared to the image obtained from the proposed method.

Experiments on Real US Images
A set of real US images was collected from an online US image database [43].Experiments were performed on these images and subjective quality assessment was carried out.For objective quality assessment, a parameter known as MVR is used.MVR is the ratio of mean to variance of a selected region in an image.A different value of MVR is obtained for a differently selected region in each image.A higher value of MVR parameter is desired in case of a better-denoised image.8 presents the MVR plot obtained for Abdomen image.For all the methods presented in Fig. 7, MVR values are calculated in all highlighted regions.It is evident from Fig. 8 that a higher value of MVR is obtained for the proposed method as compared to other methods in all four regions.Similarly, a total number of 120 different regions were chosen from a set of 50 real US images for MVR calculation.Table 3 shows the average MVR values along with their standard deviations for 120 measurements.MVR values are calculated in all these regions for all the denoising methods under consideration.

Conclusion
In this paper, a new speckle reduction method based on shearlet coefficient modelling of US image is presented.Detailed band shearlet coefficients are modelled as a Nakagami prior.The denoised shearlet coefficients are obtained using MAP estimation in each detail subband.Subjective and objective evaluation is done for synthetic and real US images.From objective evaluation of the proposed method, it is evident that the higher values of all assessment parameters are obtained for both synthetic and real US images.The effectiveness of the proposed method can also be seen from the subjective evaluation of all denoised images obtained from various methods.Hence, the proposed method outperforms other existing methods and is better in view of US image denoising and fine detail preservation.
Engineering College, Ghaziabad for providing the necessary research facilities.

Figure 2 (
Figure 2(a) and Fig. 2(b) show the tiling of the spatial frequency plane induced by shearlets and frequency support of shearlets [29].

Figure 5 (
Figure 5(a), Fig. 5(b), Fig. 5(c), Fig. 5(d), Fig. 5(e), Fig. 5(f), Fig. 5(g) and Fig. 5(h) present synthetic kidney US images for subjective quality evaluation for existing methods and the proposed method.It can be seen that the visual quality of the image obtained using the proposed method is the best among all other denoised images obtained from other existing methods.It can be seen clearly that the noise is eliminated and

Figure 7
Figure7presents a real US image of the Abdomen.Rectangular boxes of four colours viz.blue, red, green and yellow show the various regions denoted as Region 1 to Region 4. These regions are selected for calculation of MVR values which are represented as MVR1, MVR2, MVR 3 and MVR4.Figure8presents the MVR plot obtained for Abdomen image.For all the methods presented in Fig.7, MVR values are calculated in all highlighted regions.It is evident from Fig.8that a higher value of MVR is obtained for the proposed method as compared to other methods in all four regions.

Figure
Figure7presents a real US image of the Abdomen.Rectangular boxes of four colours viz.blue, red, green and yellow show the various regions denoted as Region 1 to Region 4. These regions are selected for calculation of MVR values which are represented as MVR1, MVR2, MVR 3 and MVR4.Figure8presents the MVR plot obtained for Abdomen image.For all the methods presented in Fig.7, MVR values are calculated in all highlighted regions.It is evident from Fig.8that a higher value of MVR is obtained for the proposed method as compared to other methods in all four regions.

Fig. 6 :
Fig. 6: Real head US reconstructed images obtained for various denoising methods.

Tab. 3 :
Comparison of MVR values for different methods.METHOD Average MVR Values with Std.Deviation Noisy Image 11.57 ± 5.89 Wiener [6]