Iterative Wavelet thresholding for suppression of speckle noise from magnetic resonance images.

Sudha, S. ; Suresh, G.R. ; Sukanesh, R. 等

Introduction

Introductory section furnishes brief idea about different available denoising schemes. Magnetic resonance imaging is a widely used medical imaging procedure because it is economical, comparatively safe, transferable, and adaptable. Though, one of its main shortcomings is the poor quality of images, which are affected by speckle noise. The existence of speckle is unattractive since it disgrace image quality and it affects the tasks of individual interpretation and diagnosis. Accordingly, speckle filtering is a central pre-processing step for feature extraction, analysis, and recognition from medical imagery measurements. Previously a number of schemes have been proposed for speckle mitigation.

An appropriate method for speckle reduction is one which enhances the signal to noise ratio while conserving the edges and lines in the image. Filtering techniques are used as preface action before segmentation and classification. On the whole speckle reduction can be divided roughly into two categories. The first one recovers the image by summing more than a few observations of the same object which suppose that no change or motion of the object happened during the reception of observations. D. Hillery et al [13] adopted filtering in the spectral domain, but the classical Wiener filter is not adequate while it is designed primarily for additive noise suppression Javier Portilla [3].To address the multiplicative nature of speckle noise, Jain developed a homomorphic approach, which by obtaining the logarithm of the image, translates the multiplicative noise into additive noise, and consequently applies the Wiener.

Adaptive filter takes a moving filter window and estimates the statistical information of all pixels' grey value, such as the local mean and the local variance. The central pixel's output value is dependent on the statistical information. Adaptive filters adapt themselves to the local texture information surrounding a central pixel in order to calculate a new pixel value. Adaptive filters implemented by Shi, Z. et.al [12],Li. C et.al [16] generally incorporate the Kuan filter Lee filter Frost filter Gamma MAP filters. These filters made obvious their superiority measured up to low pass filters, since they have taken into account the local statistical properties of the image. Adaptive filters present much better than low-pass smoothing filters, in preservation of the image sharpness and details while suppressing the speckle noise. X. Zong et.al [9]

Recently many challenges have been made to reduce the speckle noise using wavelet transform as a multi-resolution image processing tool. Speckle noise is a high-frequency component of the image and appears in wavelet coefficients. One widespread method exploited for speckle reduction is wavelet shrinkage. When multiplicative contamination is concerned; multiscale methods engage a preprocessing step consisting of a logarithmic transform to separate the noise from the original image. Then different wavelet shrinkage approaches are employed, which are based on Donoho's novel work. The well-known technique of wavelet shrinkage has been systematized by Donoho [11]. Achim et.al [4], Abrishami Moghaddam. H.et.al [1], suggested the speckle reduction through wavelet transform based on Bayesian approach by means of the statistical models of both noise and signal. Gagnon and Jouan [10] performed a comparative study between a complex wavelet coefficient shrinkage filter and several standard speckle filters that are largely used by imaging scientists, and show that the wavelet-based approach is deploy among the best for speckle removal. Fodor. I et.al [2], MaartenJanse [5] proposed wavelet thresholding method, in proportion to this larger wavelet coefficients of the image signify the signal and smaller ones signify the noise. The threshold (T) is computed based on statistical properties of the noisy data using different shrinkage rules.Thitimajshima.Pet.al [8] proposed wavelet thresholding function based on a shrinkage function such as hard-thresholding or soft-thresholding applies this threshold (T) to modify the wavelet coefficients. The main difficulty with this method is to find out the optimal threshold value.

In our work, we recommend a novel thresholding algorithm for denoising speckle in MRI with wavelets. We favor our approach by Bayes Shrinkage function proposed by Grace Chang.S. et.al [6], [7] .The statistical analysis process is exactly the same for all data sets. Carrying out the statistical test in the wavelet domain necessitates an inverse wavelet transform subsequently, which spreads out the activation in the final statistical map. The data sets used in this study require simple statistical test.

The paper is organized as follows: Section I, depicts the mathematical model for speckle noise and available speckle filters. Section II, briefly highlights the main features of wavelets and the wavelet decomposition. In section III an image adaptive threshold imposed on the wavelet coefficient is calculated to identify the significant structures. Section IV reviews how to parameterize the compactly supported threshold. Denoising procedure is explained in section V. Experimental results are given in Section VI in comparison with some existing wavelet thresholding schemes. Finally, Section VII concludes the paper.

Speckle Noise in MRI Images

MRI is the widely used medical imaging modality for diagnosing, visualizing and evaluating number of conditions of a patient. Most imaging modalities use injectable contrast, or dyes, for certain procedures. MRI contrast works by altering the local magnetic field in the tissue being examined. Normal and abnormal tissue will respond differently to this slight alteration, giving us differing signals. These varied signals are transferred to the images, allowing us to visualize many different types of tissue abnormalities and disease processes better than we could without the contrast. The main source of noise in MRI images is the Rician noise, and speckle noise. This work aims to suppress speckle in MRI.

Speckle noise affects all coherent imaging systems including medical MRI. Within each resolution cell a number of elementary scatterers reflect the incident wave towards the sensor. The backscattered coherent waves with different phases undergo a constructive or a destructive interference in a random manner. The acquired image is thus corrupted by a random granular pattern, called speckle that delays the interpretation of the image content.

A speckled image V = {[v.sub.1], [v.sub.2], [v.sub.3]......[v.sub.n]}is commonly modeled as [v.sub.1] = [f.sub.1][??]

where f = [f.sub.1], [f.sub.2], [f.sub.3],....[f.sub.n]} is a noise-free ideal image, and [??] = {[[??].sub.1] [[??].sub.2],.....[[??].sub.n]} is a unit mean random field..

Speckle filters

Some of the best known standard despeckling filters are Lee, Frost and Kuan. These filters use the second-order sample statistics within a minimum mean squared error estimation approach. Enhanced Lee and the enhanced Frost filters merge the filtering with a preliminary classification step: First the image pixels are consigned into one of the three classes: homogeneous, weakly textured or highly heterogeneous. Homogeneous image segments are simply averaged, while the highly heterogeneous ones are kept untouched; and only the remaining image segments are adaptively filtered. Another common despeckling approach is the homomorphic Wiener filter where the speckle noisy image is first subjected to a logarithmic transform and then filtered with an adaptive filter for additive Gaussian noise. Studies that compare different speckle filters in the image domain and in the wavelet domain usually show that wavelet domain filters are better able to preserve image details.

In the medical literature, speckle noise is referred as "texture", and may possibly contain useful diagnostic information. The desired grade of speckle smoothing preferably depends on the specialist's knowledge and on the application like enhancement for visual inspection or preprocessing for automatic segmentation. For automatic segmentation, sustaining the sharpness of the boundaries between different image regions is usually preferred while smooth out the speckled texture. For visual interpretation, smoothing the texture may be less desirable.

In most natural images counting medical images, there in general exists a context models like Markov random fields, or simply on adapting certain filter parameters based on measurements from a local window around each pixel. Most of the speckle filters assume speckle which is modeled as a multiplicative noise and perform logarithmic operation transforms speckle into additive white Gaussian noise. Because of different reasons such speckle model appears to be adolescent in the case of medical images. Hence the noise differs noticeably from the frequently assumed multiplicative model in the displayed medical MRI images.

Physicians generally have a preference of the original noisy images more willingly than the smoothed versions because the filters even if they are more sophisticated can destroy some relevant image details. In many cases noise suppression appreciably enhances the visibility of some image features and it certainly facilitates automatic image processing tasks such as segmentation. Thus it is essential to develop noise filters which can secure the conservation of those features that are of interest to the physician. The wavelet transform has recently entered the field of image denoising and it has firmly recognized its stand as a dominant denoising tool.

Wavelet Domain Noise Filtering

The wavelet transform is an extensive tool for processing an image [14], [15]. Contrasting the Fourier sinusoids, which offer a sharp frequency characterization of a given signal but are unable to categorize transient events, wavelets realize stability between localization in space or time and localization in the frequency domain. This balance is essential to multiresolution, which lets the analysis to deal with image features at any scale. As the discrete wavelet transform (DWT) corresponds to basis decomposition, it provides a non redundant and unique representation of the signal. More specifically, methods based on multiscale decompositions consist of three main steps: First, the raw data are decomposed by means of the wavelet transform, then the empirical wavelet coefficients are shrunk through a thresholding mechanism, and finally, the denoised signal is synthesized from the processed wavelet coefficients through the inverse wavelet transform. The discrete wavelet transform translates the image content into an approximation sub band and a set of detail sub band at different orientations and resolution scales. Typically, the band-pass content at each scale is divided into three orientation sub band characterized by horizontal, vertical and diagonal directions. The approximation sub band consists of the so-called scaling coefficients and the detail sub band is composed of the wavelet coefficients.

Several properties of the wavelet transform, which make this representation attractive for denoising, are

* Multiresolution--image details of different sizes are analyzed at the appropriate resolution scales

* Sparsity--the majority of the wavelet coefficients are small in magnitude.

* Edge detection--large wavelet coefficients coincide with image edges.

* Edge clustering--the edge coefficients within each sub band tend to form spatially connected clusters.

During a two level of decomposition of an image using a scalar wavelet, the two-dimensional data is replaced with four blocks. These blocks correspond to the sub bands that represent either low pass filtering or high pass filtering in each direction. The procedure for wavelet decomposition consists of consecutive operations on rows and columns of the two-dimensional data. The wavelet transform first performs one step of the transform on all rows. This process yields a matrix where the left side contains down sampled lowpass coefficients of each row, and the right side contains the high pass coefficients. Next, one step of decomposition is applied to all columns; this results in four types of coefficients.

(1) HH represents the diagonal features of the image and is formed by high pass filtering in both directions.

(2) HL represents the horizontal features of the image and is formed by low pass filtering the rows and then high pass filtering the columns.

(3) LH represents the vertical features of the image and is formed by high pass filtering the rows and then low pass filtering the columns.

(4) LL represents the coefficients that will be further decomposed in the next step. It is formed by low pass both filtering the rows and the columns.

Parameter Computation for Threshold (T)

Generally a small threshold value will leave behind all the noisy coefficients and subsequently the consequential denoised signal may still be over excited. A large threshold value alternatively, makes more number of coefficients as zero which directs to smooth signal and destroys details and the resultant image may cause blur and artifacts. And so optimum threshold value should be found out, which is adaptive to different sub band characteristics. Here, we describe an efficient method for fixing the threshold value for denoising by analyzing the statistical parameters of the wavelet coefficients. The innovative aspects of the present work consist of the estimating appropriate threshold.

Following:

In Bayes Shrink it is determined that the threshold for each sub band assuming a Generalized Gaussians distribution (GGD). The GGD is given by

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (1)

The parameter, [sigma]x is the standard deviation and [beta] is the shape parameter. Assuming such a distribution for the wavelet coefficients, we estimate [sigma]x and [beta] for each sub band and try to find the threshold T which minimizes the Bayesian Risk, i.e., the expected value of the mean square error.

[tau](T) = E[(X - X).sup.2] = [E.sub.X] [E.sub.X/Y] [(X - X).sup.2] (2)

The optimal threshold T is then given by

T* ([sigma]x, [beta]) arg Min (T) (3)

This is a function of the parameters [sigma]x and [beta] since there is no closed form solution for T*, numerical calculation is used to find its value. It is observed that the threshold value set by

T = [beta] [[sigma].sup.2] / [[sigma].sub.X] (4)

The normalized threshold, [T.sub.B]/[sigma] is inversely proportional to the standard deviation of X ([sigma]), and proportional to the noise standard deviation ([sigma]x). When [sigma]/[sigma]x [greater than or equal to] 1 the signal is much stronger than the noise, T is chosen to be small in order to preserve most of the signal and remove some of the noise; when [sigma]/[sigma]x [greater than or equal to] 1, the noise dominates and the normalized threshold is chosen to be large to remove the noise which has weighed down the signal. Thus, this threshold choice adapts to both the signal and the noise characteristics as reflected in the parameters [sigma] and [sigma]x.The parameters, [sigma]x and [beta], need to be estimated to compute T ([sigma]x).

This section focuses on the estimation of the GGD parameters, [[sigma].sub.X] and [beta] which in turn yields a data-driven estimate of T([[sigma].sub.X])that is adaptive to different sub band characteristics. The noise variance [[sigma].sup.2] needs to be estimated first. It may be possible to measure [[sigma].sup.2] based on information other than the corrupted image and it is estimated from the sub band [HH.sub.1] by the robust median estimator,

[[sigma].sup.2] = [[median[sub.|m,n|]/0.6745].sup.2] (5)

Where [[sigma].sup.2] is the noise variance, [[sigma].sub.X] is the signal standard deviation. [[sigma].sub.X] is the signal standard deviation needs to be estimated.

From the observation model Y = X + [epsilon], with X and [epsilon] independent o each other we have

[[sigma].sup.2.sub.Y] = [[sigma].sup.2.sub.X] + [[sigma].sup.2] (6)

Where [[sigma].sup.2.sub.Y] is the variance of Y. It can be found by

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

From this [[sigma].sub.X] can be derived as

[[sigma].sub.X] = [square root of max ([[sigma].sup.2.sub.Y] - [[sigma].sup.2], 0)]. (7)

The parameter [beta] is weighted variance, which involves neighboring coefficients of the wavelet decomposition for the estimation of the local variance. It is based on the estimation of the local weighted variance [[sigma].sub.w][[m, n].sup.2] of each wavelet coefficient [Y.sup.(l,o)][m, n] at level l and orientation o[member of] {V,H,D} using a window N, which covers a 5 x 5 neighborhood of [Y.sup.(l,o)][m, n] . The weighted variance of a coefficient [Y.sup.(l,o)] [m, n] with respect to the window 5 x 5 and a corresponding set of weights w = {[w.sup.(l,o).sub.j,k] | j, k [member of] N} is defined by:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (8)

The value weighted variance([beta]) of a given wavelet coefficient is determined by the weight in a local window containing neighboring coefficients as well as corresponding parent coefficients in the next higher decomposition level of the multiresolution decomposition, where the weights depend on the decomposition level and the orientation of given sub band.

The suggestion behind the recommend method is to iterate the context-based thresholding process on the denoised wavelet representation. Better visual quality substantially can be acquired by varying the context weights in each iteration step appropriately.

Image Denoising Algorithm

A general procedure is: (i) calculate the discrete wavelet transform; (ii) remove noise from the wavelet coefficients and (iii) reconstruct a denoised signal or image by applying the inverse wavelet transform. The noise-free component of a given wavelet coefficient is typically estimated by wavelet shrinkage , the idea of which is to heavily suppress those coefficients that represent noise and to retain the coefficients that are more likely to represent the actual signal or image discontinuities. This section depicts the image-denoising algorithm, which achieves near optimal soft thresholding in the wavelet domain for recovering original signal from the noisy one. The algorithm is very simple to implement and computationally more efficient. The wavelet transform employs Daubechies' least asymmetric compactly supported wavelet with eight vanishing moments with four scales of orthogonal decomposition. It has the following steps.

1. Perform the DWT of the noisy image up to 2 levels (L=2) to obtain seven sub bands, which are named as L[L.sub.1], HH1, L[H.sub.1], H[L.sub.1], H[H.sub.2], L[.sub.H2], H[L.sub.2] and L[L.sub.2].

2. Compute the threshold value T for each sub band, except the L[L.sub.2] band using equation

(i) Obtain the noise variance using the equation (5)

(ii) Calculate the signal standard deviation [[sigma].sub.X]by the equation (7)

(iii) Find out the parameter [beta] from equation (8)

3. Threshold the all sub band coefficients using Soft Thresholding given in equation (4) by substituting the threshold value obtained from the step 8.

4. Perform the inverse DWT to reconstruct the denoised image.

Experimental Results and Discussion

The performance of the wavelet thresholding method that has been proposed in this paper is investigated with simulations. White Gaussian noise, speckle noise with variance [sigma]2 = 0.03, [[sigma].sup.2] = 0.04, [[sigma].sup.2] = 0.05 is added to a 256x256 MRI image and Barbara and denoising by soft thresholding. For the same images denoising is carried out with the hard threshlding Bayes thresholding and proposed thresholding. The simulation results can be estimate impartially and instinctively. For objective evaluation, the signal to noise ratio (SNR) of each denoised image has been calculated using Signal to Noise Ratio (SNR) , which is defined as

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Where S, S represent the original and denoised images, respectively.

Table I presents the quantitative results of different noise suppression techniques for various noises. It can be observed that the proposed method has the best results from the signal to noise ratio and edge preservation points of view. We have also made comparisons with the Wiener filter, the best linear filtering possible. The version used is the adaptive filter, wiener2, in the MATLAB image processing toolbox, and they are considerably worse than the nonlinear thresholding methods, especially when [[sigma].sup.2] is large. The image quality is also not as good as those of the thresholding methods.

[FIGURE 3 OMITTED]

Conclusion

A speckle noise reduction method in wavelet domain using an Iterative thresholding based on Bayes approach is commenced. Major offerings of this work are the optimal selection of the wavelet threshold parameter. We have introduced a relatively simple context-based model for adaptive threshold selection within a wavelet thresholding framework. Estimations of local weighted variance with appropriately chosen weights are used to adapt the threshold. The proposed thresholding technique outperforms the soft thresholding, Bayes thresholding and the filters like wiener. The proposed method removes noise significantly and remains within 3% of the Bayes thresholding. Moreover the computational time is 3% faster for proposed method. However, by visual inspection it is evident that the denoised image, while removing a substantial amount of noise, suffers practically no degradation in sharpness and details. By iterating the thresholding operation, more accurate reconstruction can be realized. Experimental results show that our proposed method yields significantly improved visual quality as well as lower mean squared error compared to the other techniques in the denoising literature.

References

[1] Abrishami Moghaddam. H, Valadan Zouj, M.J. Dehghani. M. (2004)." Bayesian-based speckle reduction using wavelet transforms," Accepted in the Conference on applications of Digital Image Processing in 49th SPIE annual meeting.

[2] Denver. Fodor I. K, Kamath. C. (2003). "Denoising Through Wavelet Shrinkage: An Empirical Study, Journal of Electronic Imaging", 12, pp.151-160.

[3] Javier Portilla, Vasily Strela, Martin J.Wainwright and Eero P. Simoncelli. (2002) "Adaptive Wiener Denoising using a Gaussian Scale Mixture Model in the wavelet Domain", Proceedings of the 8th International Conference of Image Processing Thessaloniki, Greece.

[4] Achim Bezerianos A. and Tsakalides.P (2001). "Novel Bayesian Multiscale for Speckle Removal in Medical Ultrasound Images", IEEE Transactions. Medical Imaging Journal, 20[8], pp. 772-783.

[5] Maarten Janse.(2001)- "Noise Reduction by Wavelet Thresholding", Volume 161, Springer Verlag, United States of America, I edition

[6] Grace Chang. S., Bin Yu and M. Vattereli. (2000)--"Adaptive Wavelet Thresholding for Image denoising and Compression", IEEE Transaction, Image Processing, vol. 9, pp. 1532-1546.

[7] Grace Chang. S., Bin Yu and M. Vattereli.(2000)-"Spatially Adaptive Wavelet Thresholding with Context Modeling for Imaged noising", IEEE Transaction Image Processing, volume 9, pp. 1522-1530.

[8] Thitimajshima.P, Rangsanseri.Y, and Rakprathanporn.P(1998),"A Simple SAR Speckle Reduction by Wavelet Thresholding", Proceedings of the 19th Asian Conference on remote Sensing ACRS98, pp. P-14-1--P-14-5.

[9] X. Zong, A. F. Laine and E. A. Geiser, (1998) "Speckle reduction and contrast enhancement of echocardiograms via multiscale nonlinear processing", IEEE Transactions on. Medical Imaging, vol. 17, pp. 532-540.

[10] L. Gagnon and A. Jouan(1997)"Speckle filtering of SAR images a comparative study between complex wavelet based and standard filters, "Proc. SPIE,vol.3169, pp. 80-91.

[11] D. L. Donoho, (1995), "Denoising by soft-thresholding," IEEE Trans. Inform. Theory, vol. 41, pp. 613-627.

[12] Shi. Z and Fung, K. B(1994)" A Comparison of Digital Speckle Filters", Proceedings of IGRASS'94, August, Pasadena.

[13] D.Hillery et al. (1991)," Iterative wiener filters for images restoration", IEEE Transaction on SP, 39, pp. 1892-1899.

[14] I. Daubechies,(1990), "The wavelet transform, time-frequency localization and signal Analysis," IEEE Trans. Inform. Theory, vol. 36, no. 5, pp. 961-1005,

[15] S. Mallat (1989), "A theory for multiresolution signal decomposition and wavelet decomposition", IEEE Trans. Pattern Anal. Mach. Intell., vol. 11, no.7, pp.674-693,.

[16] Li. C, (1987)," Two Adaptive Filters for Speckle Reduction in SAR Imagery by Using .The variance Ratio", International Journal of Remote Sensing, 9 (4), pp. 641-653.

S. Sudha *, G.R. Suresh * and R. Sukanesh ($)

* : Faculty, Sona College of Technology, Salem--6360005, TamilNadu, India.

($) : Professor, Thiagarajar College of Engineering, Madurai--15, TamilNadu, India.

E-mail: sudha76s@yahoo.com, sureshgr@rediffmail.com

Table 1: SNR vs. noise variance [[sigma].sup.2] for three
standard test images for two different types of noises
obtained by Iterative wavelet thresholding (proposed)
method and existing methods.

            Peak Signal To Noise Ratio in dB (PSNR)

Images      [[sigma].sup.2]   Noise      Soft           Bayes
                              type       thresholding   thresholding

Cameraman   0.03              Gaussian   29.1324        30.4326
MRI                                      30.3216        30.352
Barbara                                  31.1621        32.1621
Cameraman   0.03              Speckle    33.2386        34.1572
MRI                                      42.4231        43.5621
Barbara                                  34.0234        36.0267
Cameraman   0.04              Gaussian   28.5286        29.5889
MRI                                      28.1794        29.1382
Barbara                                  35.1256        35.3255
Cameraman   0.04              Speckle    34.3127        33.3811
MRI                                      40.6524        41.6516
Barbara                                  33.5281        35.2831
Cameraman   0.05              Gaussian   27.0024        28.8521
MRI                                      28.1532        28.384
Barbara                                  29.4216        30.437
Cameraman   0.05              Speckle    30.7361        32.7841
MRI                                      40.2143        42.1171

Images      [[sigma].sup.2]   Noise      Weiner         Iterative
                              type                      wavelet
                                                        thresholding
                                                        (Proposed
                                                        method)

Cameraman   0.03              Gaussian   30.5786        31.4481
MRI                                      31.2547        32.4045
Barbara                                  32.7521        32.3842
Cameraman   0.03              Speckle    35.1374        36.8341
MRI                                      42.5834        44.5581
Barbara                                  36.1826        37.2692
Cameraman   0.04              Gaussian   29.8395        30.3761
MRI                                      29.2458        29.2684
Barbara                                  34.8234        36.1277
Cameraman   0.04              Speckle    34.2379        35.546
MRI                                      41.4392        42.9089
Barbara                                  35.9631        36.1864
Cameraman   0.05              Gaussian   28.3214        29.4641
MRI                                      28.5248        28.4418
Barbara                                  29.8761        30.6227
Cameraman   0.05              Speckle    32.1293        33.2592
MRI                                      41.2746        42.4282

Figure :1: 5*5 window with variable weight for
calculating weighted variance

[W.sub.8]   [W.sub.6]   [W.sub.4]   [W.sub.6]   [W.sub.8]
[W.sub.7]   [W.sub.3]   [W.sub.1]   [W.sub.3]   [W.sub.7]
[W.sub.5]   [W.sub.2]   [W.sub.0]   [W.sub.2]   [W.sub.5]
[W.sub.7]   [W.sub.3]   [W.sub.1]   [W.sub.3]   [W.sub.7]
[W.sub.8]   [W.sub.6]   [W.sub.4]   [W.sub.6]   [W.sub.8]