Application of Noise Invalidation
Denoising in MRI
Pegah Elahi 2012-09-18
Application of Noise Invalidation
Denoising in MRI
Pegah Elahi 2012-09-18
LiTH-IMT/MASTER-EX12/019SE
Supervisor: Soosan Beheshti, Associate Professor,
Department of Electrical and Computer Engineering, Ryerson University Examiner: Hans Knutsson, Professor,
Abstract
Magnetic Resonance Imaging (MRI) is a common medical imaging tool that have been used in clinical industry for diagnostic and research purposes. These images are sub-ject to noises while capturing the data that can eect the image quality and diagnos-tics.Therefore, improving the quality of the generated images from both resolution and signal to noise ratio (SNR) perspective is critical. Wavelet based denoising technique is one of the common tools to remove the noise in the MRI images. The noise is eliminated from the detailed coe cients of the signal in the wavelet domain. This can be done by applying thresholding methods. The main task here is to nd an optimal threshold and keep all the coe cients larger than this threshold as the noiseless ones. Noise Invalidation Denoising technique is a method in which the optimal threshold is found by comparing the noisy signal to a noise signature (function of noise statistics). The original NIDe approach is developed for one dimensional signals with additive Gaussian noise. In this work, the existing NIDe approach has been generalized for applications in MRI images with dierent noise distribution. The developed algorithm was tested on simulated data from the Brainweb database and compared with the well-known Non Local Mean l-tering method for MRI. The results indicated better detailed structural preserving for the NIDe approach on the magnitude data while the signal to noise ratio is compatible. The algorithm shows an important advantageous which is less computational complexity than the NLM method. On the other hand, the Unbiased NLM technique is combined with the proposed technique, it can yield the same structural similarity while the signal to noise ratio is improved.
KEYWORDs: Magnetic Resonance Imaging, Noise Invalidation Denoising, Unbiased Non Local Mean ltering, Wavelet Transform Function
Acknowledgements
I would like to take an oppurtunity to express my appriciation for all the people who helped me through this work. At srt, I am grateful to my supervisor Prof. Soosan Beheshti at Ryerson University not only for giving me the opportunity to work under her supervison but also for her support and guidance throughout this thesis work.
I would like to acknowlege Prof. Hans Knutsson for accepting being my examiner and evaluating my work which means a lot based on his knowldge and experience in medical imaging. Also, many thanks to Prof. Göran Salerud director of studies in Biomical engineering department for his patience and support during my thesis work.
Last but not least, I would like to express my deepest gratitude to my parents for their unconditional love and support not only during this project but in every steps of my life. I could not achieve my goals without them standing by my side and their encouragement.
Contents
1 Introduction 1
2 Background and literature review 3
2.1 MRI imaging . . . 3
2.1.1 creating a MRI image . . . 3
2.2 MRI signal distribution . . . 5
2.2.1 Signal distribution of the raw data . . . 6
2.2.2 Signal distribution of the magnitude data . . . 6
2.2.3 Signal distribution of the squared magnitude data . . . 7
2.3 Noise estimation . . . 7
2.3.1 Estimating the noise variance from the background . . . 7
2.3.2 Estimating the noise variance based on local statistics . . . 8
2.4 Denoising methods . . . 9
2.4.1 Gaussian ltering . . . 9
2.4.2 Bilateral ltering . . . 9
2.4.3 Anisotropic Diusion Filtering . . . 10
2.4.4 Wavelet denoising . . . 10
2.4.5 Non Local Mean (NLM) lters . . . 15
3 Method and Simulation Results 17 3.1 Method and material . . . 17
3.1.1 Noise Invalidation Denoising for Rician distribution . . . 17
3.1.2 Noise Invalidation Denoising for chi-square distribution . . . 18
3.1.3 Noise variance estimation . . . 19
3.1.4 Evaluation criteria . . . 19
3.1.5 Material . . . 20
3.2 Simulation Results . . . 21
3.2.1 Results for NIDe for Rician distributed data . . . 21
3.2.2 Results for NIDe for chi-squared distributed data . . . 31
3.2.3 Results for combining NIDe with Richian distribution with UNLM 39 4 Discussion and future works 41 4.1 Comparison based on evaluation criteria . . . 41
4.2 Limitations . . . 53
4.3 Conclusion and future trends . . . 53
List of Figures
2.1 Spin and its related angular momentum [15] . . . 4
2.2 Generated current signal because of applying RF wave [2] . . . 5
2.3 Local second order distribution for a: normal image and b: MRI image. The dashed line shows the distribution related to the noisy images[20] . . 8
2.4 Eect of sorting the absolute values: samples of unsorted noise, sorted noise values and changed axes for top to bottom respectively[21] . . . 14
2.5 Finding threshold with NIDe. Solid line is the noisy data and dashed lines represent the noise only boundaries[21] . . . 15
3.1 Probability density function plots for noisy image, data in the background and signal area top to bottom respectively . . . 22
3.2 ANS function on noise only data with Rayleigh distribution . . . 23
3.3 Noise con dence area for the Rician distributed Noise signature. The blue line shows the mean while the two red lines shows the upper and lower boundaries . . . 24
3.4 Finding the optimal threshold for the NIDe with Rician distribution method. The blue line shows the sorted noisy data and how it exists the con dence area . . . 24
3.5 Denoising result of NIDe for Rician distribution with 3% noise . . . 25
3.6 Denoising result of NIDe for Rician distribution with 5% noise . . . 26
3.7 Denoising result of NIDe for Rician distribution with 7% noise . . . 27
3.8 NIDe approach based on Rayleigh distribution results on T2 weighted image with 3% noise level . . . 28
3.9 NIDe approach based on Rayleigh distribution results on T2 weighted image with 5% noise level . . . 29
3.10 NIDe approach based on Rayleigh distribution results on T2 weighted image with 7% noise level . . . 30
3.11 PDF plot of the squared data . . . 31
3.12 Eect of sorting on the noisy data with chi-squared distribution . . . 32
3.13 Noise con dence area for the chi-squared distribution noise. The blue line shows the mean while the red lines are the limits of this area. . . 33
3.14 Finding the threshold in the NIDe approach based on chi-squared distri-bution. The blue line shows the sorted noisy data. . . 33 3.15 Final Result of the denoising with NIDe on the squared data with 3% noise 34 3.16 Final Result of the denoising with NIDe on the squared data with 5% noise 35 3.17 Final Result of the denoising with NIDe on the squared data with 7% noise 36
3.18 Final result of NIDe algorithm based on chi-squared distribution on the T2 weighted MRI with 3% noise. . . 37 3.19 Final result of NIDe algorithm based on chi-squared distribution on the
T2 weighted MRI with 5% noise. . . 38 3.20 Final result of NIDe algorithm based on chi-squared distribution on the
T2 weighted MRI with 7% noise. . . 39 3.21 T1-weighted MRI image with 5% noise level ltered with NIDe-UNLM . 40 4.1 Comparing gures of each methods for T1 weighted with 3% noise . . . . 44 4.2 Comparing gures of each methods for T1 weighted with 5% noise . . . . 45 4.3 Comparing gures of each methods for T1 weighted with 7% noise . . . . 46 4.4 Comparing gures of each methods for T2 weighted with 3% noise . . . . 48 4.5 Comparing gures of each methods for T2 weighted with 5% noise . . . . 49 4.6 Comparing gures of each methods for T2 weighted with 7% noise . . . . 50 4.7 Final results of NIDe-Rician combined with UNLM . . . 52
List of Tables
3.1 Threshold values T1/T2 weighted MRI images at dierent noise level for Rician distribution . . . 30 3.2 Threshold values T1/T2 weighted MRI images at dierent noise level for
chi squared distribution . . . 39 4.1 Comparison table for noise level of 3% on T1 weighted MRI image . . . . 41 4.2 Comparison table for noise level of 5% on T1 weighted MRI image . . . . 41 4.3 Comparison table for noise level of 7% on T1 weighted MRI image . . . . 42 4.4 Speed dierence . . . 42 4.5 Comparison table for noise level of 3% on T2 weighted MRI image . . . . 47 4.6 Comparison table for noise level of 5% on T2 weighted MRI image . . . . 47 4.7 Comparison table for noise level of 7% on T2 weighted MRI image . . . . 47 4.8 The improved statistics for NIDe-Rician combined with UNLM . . . 51
1 Introduction
Magnetic Resonance Imaging (MRI) is one of the eective medical equipments that has been proven to be less harmful for patients compared to other medical modalities. This method is more accurate for imaging soft tissues such as brain and muscles. It is known that better resolution in MRI images results in less SNR (signal to noise ratio) and vice versa [2]. Therefore, applying a proper denoising technique provides a better opportunity for accessing to an acceptable image quality and higher SNR.
MRI image is generated by using a reconstruction method based on the sampling process (e.g. using the inverse Fourier transform of the raw data in k-space in case of regular Cartesian grid sampling). However, the resulted data is still complex values due to the magnetic eld or hardware de ciencies. Therefore, the nal image is calculated as the magnitude of the complex data[2] .
The noise in MRI signal is a white Guassian noise added to each part of the com-plex data. The magnitude image which is computed as the square root of two squared Gaussian variables added together has Rician distribution. One can demonstrate that the distribution becomes Rayleigh in high SNR areas and Gaussian in low SNR (back-ground) regions of the image [5]. So, the developed ltering methods which are applied on the magnitude image should be compatible with this type of signal. The signal depen-dent nature of the Rician distribution results in low contrast images, so a new method in which the ltering method is applied on the squared image is introduced. The squared signal has a non-central chi-square distribution and a signal independent noise bias. The noise bias can be easily removed after applying the lters [13].
There are two general categories for the denoising methods: acquisition-based and post acquisition ltering methods. The former methods involve increasing the scan times, per-forming several similar measurements and averaging them, increasing the voxel dimension or using more advanced hardware. However, using this type of noise reduction increases both the acquisition time and the expenses. Therefore, the second group can be con-sidered as an e cient and aordable technique for MRI image denoising [10]. Several ltering methods have been developed in past decades. The most classical lters have been used are Gaussian lters. However, these lters result in blurring edges. Hence, more edge preserving lters such as bilateral lter which functions as a weighted aver-aging lter have been developed. This lter has better results compared to linear lters but still removes some detailed information [23]
The other type of lter that is applied for MRI images is Anisotropic Diusion lter which keeps edges by averaging pixels in the orthogonal dimension of the local gradient. However, it still eliminates small featuressmall features [23].
Recently, a new method has been proposed by Buades called Non-Local Mean ltering (NLM) which uses a weighted averaging method based on the similarities between
dif-ferent neighborhoods in the image. This method takes advantage of the fact that there are many repeated structures in most natural images and removes the noise by nding and averaging these structures [13].
One of the most common methods that have widely been used for MRI imaging is wavelet based ltering method. In this procedure, approximation and detailed subbands of the image are generated in the wavelet domain and then the detailed subbands which contain the high frequency data including small features and mostly noise are modi ed by thresholding method. Here, the most important factor in success of the denoising method is to nd the optimal threshold by which the noise can be removed while the detailed structures and specially the edges are preserved. The most common methods for nding the optimal threshold are VisuShrink, SureShrink, BaysShrink and NeighShrink. Although, they still remove some detailed information [5].
The aim of this project is to develop new method of noise removal. The main concen-tration is on wavelet optimal thresholding for denoising MRI images. For this purpose, a new method called Noise Invalidation Denoising (NIDE) is investigated. This method uses a noise signature extracted from noise order statistics for denoising the image. It does not require making any assumptions about the noise free signal and can be applied on noisy data in any type of basis transformation (orthogonal or non-orthogonal) [21]. This method is originally proposed for additive Gaussian noise.
Therfore, in this thesis the NIDe approach is furthur developed to be compatible with the MRI images. For this purpose, the original method is modi ed to be applicable on the non-central chi squared and Rician distribution of the magnitude and squared magnitude data.
2 Background and literature review
This chapter introduces some basic information about MRI imaging and the available denoising methods for this type of the image. In order to describe the denoising methods, it is necessary to provide some theoretical information about the MRI signal distribution in presence of noise. All these are explained in more details as following.
2.1 MRI imaging
Mangentic Resonance imaging (MRI) is one of the common methods have been used from 1946 in wide range of applications. Some of these areas are molecular investigations, chemical information and imaging of living tissues and organs. For the imaging purposes, MRI imaging technique provides images from internal parts of the body from received signals. This application of MRI is so important for the pathological purposes [7].
Although, there are other biomedical imaging techniques available but MRI imaging has some properties that makes it more suitable method. By this method, multi dimen-sional images can be produced. The main dierence in using MRI comparing to other available techniques here is its capability to produce two dimensional, volumetric three or four dimensional images in any direction. These images can also provide temporal information without any speci c mechanical adjustments [7].
The technique creates images from measured signals coming directly from the object without any needs to inject contrast agents. It does not have the problems regarding the radiation as it works in the radio frequency range. Finally, it can creates dierent images focusing on one speci c property of a same object by just changing some intrinsic parameters of the system. The problem with this technique is related to application of magnetic elds in the encoding process. Therefore, there would be some limitation for using this technique for example in treating patients with metal implants in their bodies [7, 2].
As the main focus of this work is on improving produced MRI images, presenting a general information about the process of the MRI images can be useful. This is described in 2.1.1 section.
2.1.1 creating a MRI image
The MRI imaging is based on the behavior of Nuclei consisting odd number of protons or neutrons when an object placed in magnetic eld. These particles have a spin angular momentum. They are called spins and the basics of the nuclear magnetic resonance (NMR). A spin and it's angular momentum is illustrated in gure 2.1 [15].
Figure 2.1: Spin and its related angular momentum [15]
The process of creating a MRI image based on spins can be categorized in four steps. Each of these steps are explained in more details as following.
Placing the object in a strong magnetic eld
At rst the object will be placed in a magnetic eld B0. This causes the spins align
themselves to the magnetic eld by precessing around that direction. This alignment can be parallel or anti-parallel to the magnetic eld based on the energy level of the spin which is low or high respectively. The dierence between these two energy levels is related to the magnetic eld and creates a bulk magnetization vector in the direction of the magnetic eld B0 [2, 7].
Applying RF waves to the object
In the next step an RF signal is applied to the object. This causes some of the spins absorb the energy and become exited. Therefore, they go from the low energy state to the high energy state. This transformation can be seen in bulk magnetization vector in this way that it ips than to the x-y plane while precessing. The precessing frequency equals to Larmor frequency (w = γB0) in which γ is the gyromagnetic ratio. This phenomena
creates a current signal in the receiver [15, 2]. This process is illustrated in gure 2.2. Spatial encoding the signal
In order to determine the spatial location of the signal, several gradient elds are applied. This changes the magnetic eld locally. The local magnetization now can be presented by the equation 2.1 [2].
m (x, t) = m0(x) exp −i w0+ γgxTn t (2.1)
where g (x)is the applied gradient eld. The received signal is the sum over the object which equals to sampling in the frequency domain or k-space.
Figure 2.2: Generated current signal because of applying RF wave [2] The Received signal can be formulated as equation 2.2 to 2.3 [2].
F {m (x)} = ˆ ∞ −∞ m (x) exp −i2πxTu dx (2.2) where u = (kx, ky, kz)and kx(t) = γ 2π ˆ 2π 0 gx(τ ) dτ ky(t) = γ 2π ˆ 2π 0 gy(τ ) dτ kz(t) = γ 2π ˆ 2π 0 gz(τ ) dτ (2.3) Image reconstruction
At the nal step, the image is reconstructed based on the applied sampling method. The simplest reconstruction technique is an inverse FFT which corresponds to the sampling on Cartesian grid. The obtained data from the reconstruction is still complex valued. This can be explained due to de ciencies of the magnetic elds or used equipments. Therefore, the magnitude or phase of the complex data is used commonly for representing the nal image [2].
2.2 MRI signal distribution
As it was mentioned in last section, the recorded data before reconstruction is complex valued. This data is also aected by an additive complex Gaussian noise with zero mean and variance of σ2. The inverse Fourier transform does not change the characteristics
of the noise due to the transform linearity and orthogonality. A MRI signal aected by noise can be formulated as the equation 2.4 [5, 7] .
y = (p + nre) + i (q + nim) (2.4)
where p ,q, nreand nimare the image and noise real and imaginary parts.
As the denoising methods are performed mostly on the raw data, magnitude data or squared magnitude data, the characteristics of the MRI signal distribution is explored based on this classi cation.
2.2.1 Signal distribution of the raw data
As it was mentioned, the inverse Fourier transform does not change the properties of the complex data. Therefore, the data after reconstruction still has a format the same as the equation 2.4. Therefore, the probability density function (PDF) of the signal can be presented by the equation 2.5 which describes a Gaussian PDF[7].
p (wr, wi|A, ϕ, σ ) = 1 2πσ2exp − (wr− A cos ϕ)2 σ2 exp − (wi− A sin ϕ)2 σ2 (2.5)
where wrand wiare the real and imaginary parts of the data with amplitude A and
phase value of ϕ.
2.2.2 Signal distribution of the magnitude data
The MRI image is most commonly the magnitude of the raw data and can be formulated with equation 2.6 [5]. Using the magnitude data makes the data unaected by some factors such as phase variations because of RF inhomogeneity, delay and the sampling method [7].
|y| = q
(p + nre)2+ (q + nim)2 (2.6)
The PDF of the magnitude data is not anymore Gaussian because of the non-linearity of the square root function. It can be derived by nding the joint PDF of the real and imaginary data and using the polar coordinates[7]. Therefore, The PDF can be represented as equation 2.7 [5]. p (y |A, σ ) = y σ2exp − y2+ A2 2σ2 I0 Ay σ2 (2.7) where A is the noiseless signal magnitude and I0is the zeroth order modi ed Bessel
function [5].
One important attribute of the Rician distribution is its signal dependency. It means that the distribution is changed to Rayleigh distribution when the SNR is low. This can be related to the areas of the image like background. The PDF of a Rayleigh distribution can be represented by equation 2.8 [5, 7].
p (y |σ ) = y σ2 exp
− y2
On the other hand, when SNR is high i.e. in image areas the signal distribution changes to Gaussian [5, 7].
2.2.3 Signal distribution of the squared magnitude data
The magnitude data can be considered as the sum of two Gaussian distributed data each squared. Therefore, its distribution is non central chi square and can be shown as equation 2.9 [7]. p (y) = 1 2σ2 y µ2 N −12 exp −y + µ 2 2σ2 IN −1 √yµ σ2 (2.9) It should be noted that the estimation of noise in the magnitude image is signal depen-dent, so it is di cult to remove the bias from the data. But for the squared magnitude data this bias becomes additive. Therefore, it is easy to remove the bias from the squared magnitude data. The bias equals to 2σ2which is shown in equation 2.10 [5].
Ey2 = E p2 + E q2 + 2σ2 (2.10)
2.3 Noise estimation
Dierent methods has been proposed to estimate the noise variance from the MRI mag-nitude data. They either use the background part of the image or some image statistics to nd the noise variance. These methods are explained here brie y.
2.3.1 Estimating the noise variance from the background
The second moment of the Rayleigh distribution (the signal distribution of the back-ground) equals to E y2 = 2σ2. Therefore, the noise variance can be obtained from the
equation 2.11 [20]. σ2= 1 2N N X i=1 Mi2 (2.11)
Also the noise can be estimated from the rst moment of the Rayleigh distribution. This can be seen in equation 2.12 [20].
σ = r 2 π 1 N N X i=1 Mi (2.12)
However, estimating the noise variance based on background areas has some drawbacks. In these methods, the main assumption is that the signal in the background area is always zero and the background area is not selected automatically. These factors aects the result and make it inaccurate [20].
2.3.2 Estimating the noise variance based on local statistics
The rst approach for estimating the noise based on local statistics is using the second order moment of the whole image (not only the background). It is proved that for a normal image the shape of the second-order moment distribution is not changed when it gets noisy by Rician noise. In MRI images the second order moment distribution has an unique property which is a maximum in the origin. This maximum appears in the image due to the presence of background in the image. Again, when the image is noisy a shift in the maximum can be seen. This principle is illustrated in gure 2.3 [20].
Figure 2.3: Local second order distribution for a: normal image and b: MRI image. The dashed line shows the distribution related to the noisy images[20]
Therefore, the position of the maximum can be used to estimate the noise variance. Equation 2.13 shows this relationship [20].
σ2= 1
2mode (µ2ij) (2.13)
where µ2ij = E
n
Iij2ois a selected neighborhood.
The same reasoning can also be applied for the local mean distribution of the MRI image. Again, due to the background area a maximum can be seen in origin. This max-imum will be shifted when the image is noisy. Therefore, the variance can be estimated by nding the maximum of the local mean. This is presented in equation 2.14 [20].
σ = r
2
πmode (µ1ij) (2.14)
2.4 Denoising methods
Dierent denoising methods have been used to improve the quality and increasing the SNR of MRI images. The SNR of an MR image can be represented by the equation 2.15 [23].
SN R = C.f (Ob) .g (Im) (2.15)
Where C is a constant depending on the physical properties of the system like magnetic ux, f (Ob) is related on the dimensions of the target, and g (Im) is based on selected imaging factors. The last function can be measured knowing readout frequency ω0, voxel
dimension Vh and the imaging time T (g (Im) = ω0Vh
√
T). It is obvious that having better resolution leads to less SNR. Therefore, applying a proper denoising technique provides the opportunity of accessing to an image with acceptable quality along with higher SNR [23].
There are two general categories for the denoising methods: acquisition-based and post acquisition ltering methods. The former methods involve increasing the scan times, per-forming several similar measurements and averaging them, increasing the voxel dimension or using more advanced hardware. However, using this type of noise reduction methods increases both the acquisition time and the expenses. Therefore, the second group can be considered as an e cient and aordable technique for MRI image denoising [10].
One important dierence in post acquisition-based denoising methods is the data they are applied to (raw, magnitude, squared magnitude). Some common methods of MRI denoising are described here.
2.4.1 Gaussian ltering
This is one of the most basic methods for denoising MRI images. However, it removes sharp edges and ne details as it is a smoothing lter [23].
2.4.2 Bilateral ltering
This method can be described as weighted averaging method. Two weight factors are used in this method. One of them is based on the spatial distance dierence while the other is based on the intensity dierence. An output of such a lter can be described as equation 2.16 [23, 22]. h (x) = P f () c (, x) .S (f () , f (x)) P c (, x) .S (f () , f (x)) (2.16) where c (, x)is the spatial distance measuring factor, f (x)is the original image and S (f () , f (x))is the intensity dierence measuring factor. The similarity measuring fac-tors can be described by equations 2.17 and 2.18 in which both of them are Gaussian [23, 22].
c (, x) = exp1/2 k − f (x)k σd 2 (2.17) S (, x) = exp1/2 kf− f (x)k σr 2 (2.18) where σdand σr are spatial and photometric spread respectively.
The bilateral ltering performs better than linear ltering methods in case of smoothing and keeping the edges. However, it still removes some details [23].
2.4.3 Anisotropic Diusion Filtering
This method is based on Linear Minimum Mean Square Error (LMMSE) estimator that is combined with a Partial Dierential Equation (PDE). The LMMSE is de ned as equation 2.19 [14, 16, 24].
ˆ
A2(x) =DM (x)2E− 2σ2+ K (x)M2(x) −DM (x)2E (2.19)
where M (x)is the magnitude data, A (x)is the noiseless signal value and K (x)is a factor de ned as equation 2.20 [14].
K (x) = 1 − 4σ2 D M (x)2 E − 2σ D M (x)4 E −DM (x)2 E2 (2.20)
where h.idescribes the expected value.
The PDE function is also can be described as equation 2.21 [14]. (
u (x, 0) = g
∂u(x,t)
∂t = div (c∇u (x, t))
(2.21) where c (x, t) = 1 − K (x, t)is the diusion function.
Although this lter preserve edges by averaging orthogonal to the local gradient, but it still removes detailed information. It also modi es the image statistics because of its edge improving characteristic [13].
2.4.4 Wavelet denoising
The wavelet transform function is used to analyze a signal in its dierent scales i.e. detailed and approximation structures. In contrast to the Fourier transform that is used to present the spectral content of a signal, wavelet transform can be used to evaluate the local properties of a signal. Furthermore, the wavelet transform can be used to detect edges[18].
To perform a wavelet transform, a speci c wavelet is selected. The wavelet is brought to a speci c scale and also shifted. Then the correlation of the analyzed signal and the
wavelet is explored. The result determines the properties of the signal in that scale. The whole process can also be de ned as a digital lter bank consisting low and high pass lters. These lters are applied on the low pass results until the desired level. The output of the nal low pass lter presents the approximation of the analyzed signal while the output of the high pass lter can be considered as the detailed part (the high frequency part) of the signal [18].
The wavelets are generated by dilation and translation of a general wavelet function called mother wavelet ψ (x). This is shown by equation 2.22 [18] .
ψa,b(x) = 1 √ aψ x − b a (2.22) The simplest form of the wavelet transform is the continuous wavelet transform. How-ever, this method is redundant and shift invariant. They are used often for signal char-acterization. Therefore, the other class of wavelet transforms i.e. Discrete wavelet trans-form (DWT) are introduced. This class of wavetrans-forms can be considered as simple sam-pling of the Continuous one. A common DWT wavelet coe cient and approximation coe cient basis can be represented in equation 2.23 and 2.24 respectively [18].
ψj,k(x) = 1 √ 2jψ x − 2jk 2j (j,k)∈Z2 (2.23) ϕj,k(x) = 1 √ 2jϕ x − 2jk 2j kZ (2.24) A signal can be decomposed with this basis according equation 2.25 [18].
f (x) = ∞ X j=−∞ ∞ X k=−∞ ωj,kψj,k(x) + ∞ X k=−∞ sJ,kϕJ,k(x) (2.25)
where the wavelet coe cients ωj,k and approximation coe cients sj,kcan be obtained
as equation 2.26 [18]. ωi,j = D f, gψi,j E sj,k= hf,ϕgi,ji (2.26)
thee.function represent the dual basis function.
The DWT method is not shift invariant. Therefore, it is not suitable for some appli-cations such as pattern identi cation and makes some problems for denoising due to not enough redundancy. Therefore, Non-decimated Discrete Wavelet Transform (nDWT) can be used to ful ll the requirements for denoising and pattern recognition [18].
In nDWT, the number of wavelet coe cients is not decreased in each scale. They are achieved by sampling the CWT at all integer locations at each translation. The wavelet coe cient basis for this function are de ned as 2j/2ψ 2−j(x − k)
and the signal is achieved from the equation 2.27 [18].
f (x) = ∞ X j=−∞ ∞ X k=−∞ D f, ˜ψj,k E ψj,k(x) (2.27)
The basis here are not anymore linearly independent and they form a frame. This method is more redundant and also shift invariant. Therefore it could be a good option for image denoising [18].
Denoising in Wavelet domain
Dierent method have been used to denoise an image in wavelet domain. They either based on estimation techniques or thresholding [11, 19, 4, 1]. In the thresholding method, the coe cients that are less than a given threshold are eliminated. There are two general types of thresholding techniques: Soft and Hard Thresholding. The hard thresholding that is not optimal for discontinuities and remains unchanged for observations is de ned as equation 2.28 [6].
ηH(ω, t) =
(
ω |ω| ≥ t
0 |ω| < t (2.28)
On the other hand, the soft thresholding approach in which the data shrinks in each try and is a better option for discontinuities can be be de ned as equation 2.29 [6].
ηS(ω, t) = w − t ω ≥ t 0 |ω| < t ω + t ω ≤ −t (2.29) The main task in this part is to nd an optimal threshold. Dierent methods have been used. In the next section one of the methods that is used to nd an optimal thresholding and is a basis for the purpose of this project is introduced.
Noise Invalidation denoising Noise invalidation Denoising can be considered as a method which works in the wavelet domain. As it was mentioned before, several methods are available for noise ltering in wavelet domain. They are designed to estimate a threshold value and then use a thersholding method to remove the noise coe cients from an image. However, the main logic of these methods are based on some characteristics of the noise-free signal and they are often application based. On the contrary, NIDe approach for denoising relies on the noise properties. It de nes a con dence noise area and consider the data outside of this area as a noise free signal. It should be noted that the original NIDe method performs on signals contaminated with additive Gaussian noise[21]. The basic principles of this method can be described as following.
A signature function can be de ned for any variables v and z as g (z, v)with nite mean GE(z) and variance Gvar(z) over the random noise vector V . Therefore, the signature
for N samples of the random noise and its mean and variance can be de ned as equations 2.30 to 2.32 respectively[21]. g z, vN = 1 N N X i=1 g (z, vi) (2.30)
E g z, vN = GE(z) (2.31)
var g z, vN = 1
NGvar(z) (2.32)
The noise signature function that is used for this method is based on Absolute Noise Sorting (ANS) de ned as equation 2.33 [21].
g (z, vi) =
(
1 |vi| ≤ z
0 |vi| > z (2.33)
The mean and variance of this signature can be formulated as equations 2.34and 2.35 respectively [21].
E g z, VN = F (z) (2.34)
var g z, VN = 1
N F (z) (1 − F (z)) (2.35)
where F (z) = 2φ z
σ − 1 is the cumulative distribution function of absolute value of
the additive noise.
It can be realized from the proposed signature function, that for each z , g z, vN = m N
in which m is the number of samples with absolute values less than z. On the other hand, the mthvalue in a sorted array is the largest v
ithat is less than z can be formulated as
m = N g z, vN
[21]. The eect of ANS function can be illustrated in gure 2.4.
As it is obvious in the gure 2.4, the sorted data place in such a denser area (middle image) than the original data (top image). Also, changing the axes shows that values are around mean NF (z)with variance F (z) (1 − F (z)). Therefore it is possible to de ne an area in which the data can be considered as noise. In other words, the probability of data being placed in a limited area around the mean value should be high (equation 2.36) [21].
P r {LN(z) ≤ g (z, vN) ≤ UN(z)} (2.36)
in which the up and down boundaries can be calculated based on the equations 2.38and 2.37respectively[21]. LN(z) = F (z) − λ r 1 NF (z) (1 − F (z)) (2.37) UN(z) = F (z) − λ r 1 NF (z) (1 − F (z)) (2.38)
According to the equations 2.37 and 2.38, the boundaries are functions of mean and standard deviation of the sorted noise only data. λis the controlling factor to increase the probability and at the same time to prevent having very loose boundaries [21].
Figure 2.4: Eect of sorting the absolute values: samples of unsorted noise, sorted noise values and changed axes for top to bottom respectively[21]
Considering a noisy signal with coe cients θi = vi+ ¯θi with ¯θias the noiseless data,
the mean and variance of samples of a random process noisy data can be represented by the equations 2.39 and 2.40 [21].
E g z, ΘN = 1 N N X i=1 H z, ¯θi (2.39) var g z, ΘN = 1 N2 N X i=1 H z, ¯θi 1 − H z, ¯θi (2.40) where H z, ¯θi
is the expected value of the signature function of noisy data.
As the noise only data that was described before, sorting the noisy data leads to a dense area of data. Comparing this area with the noise only con dence area will determine a point at which the the noisy data leave the noise con dence area. This concept can be formulated as equation 2.41 [21].
T = max
z ∀z ≤ x : LN(x) g x, θ N
UN(x) (2.41)
The explained behavior is also illustrated in gure 3.2.
Figure 2.5: Finding threshold with NIDe. Solid line is the noisy data and dashed lines represent the noise only boundaries[21]
2.4.5 Non Local Mean (NLM) lters
The NLM lters is a weighted averaging lter based on the intensity similarity of the pixels. As it was mentioned in section 2.4.2, Bilateral ltering behaves also like this. The factor that can discriminate between these two methods or other similar methods is the comparing region. In contrast to bilateral ltering that nd the similarity of a pixel with its neighborhood, the NLM lter is based on region comparison. Therefore, the pattern redundancy is not local anymore and the far pixels are not ignored [13, 3].
For a NLM lter on a noisy image Y , the output can be described as equation 2.42. N LM (Y (p)) = X
∀q∈Np
w (p, q) Y (q) (2.42)
Here, the weights are calculated according to the similarity of the squared window neighborhoods of two pixels p and q with Radius of Rsim. The weights w (p, q)can be
driven based on the equations 2.43 and 2.44 [13]. w (p, q) = 1
Z (p)exp − d (p, q)
Z (p) =X
∀q
exp −d (p, q)
h2 (2.44)
where h is a controlling factor and d which is the Euclidean distance of pixels in two neighborhoods can be described as equation 2.45 [13].
d(p, q) = Gσ kY (N p) − Y (Nq)k2Rsim (2.45)
Here, Gσis a Gaussian function with zero mean and variance of σ2.
The other important parameter in this method is Rsearchwhich is de ned as the radius
of a window to compare the pixels on that neighborhood (Np) instead of the whole image
for a speci c pixel. This will decrease the complexity of the method [13].
As it was mentioned before, the noise bias in the magnitude image is hard to remove. But it can be removed easily from the squared magnitude image. Based on this assump-tion, a version of NLM ltering called Unbiased Non Local Mean ltering (UNLM) is developed. It is described based on the equation 2.46 [13].
U N LM (Y ) = q
3 Method and Simulation Results
Based on the knowledge about MRI images, their signal distribution and also denoising methods that were explained in the previous chapter, It is possible to introduce the generalized NIDe method for MRI denoising. The new method improves the signal to noise ratio while keeping the image quality as good as possible. Therefore, in this chapter at rst the proposed method for denoising MRI images is introduced. Then after introducing the criteria for evaluating the new method, used material are described. Finally the results of running implemented algorithm on two dierent type of MRI images (T1/T2 weighted) are investigated.
3.1 Method and material
The NIDe approach that is described in chapter two is originally implemented for one dimensional signals aected with additive Gaussian noise. Although, some of its interest-ing characteristics such as usinterest-ing noise statistics instead of the noise free signal properties implies that this method can also be used on two dimensional MRI signals with dierent noise distribution after applying some modi cations. These modi cations are applied to adjust the original method with MRI dierent signal and noise distribution. The exam-ined MRI signal can either be an image signal (amplitude signal) or squared image with Rician or non-central chi-squared distribution respectively. Therefore, two approaches based on NIDe method are developed for MRI signals which are described in the follow-ing.
It should be noted that the NIDe approach for MRI images works in a wavelet domain. Therefore, it can be considered as one of the wavelet domain methods used for MRI denoising. In this method, a non-decimated wavelet transform is used. The wavelet family used here is Haar and the image signal is decomposed to three level. Then the NIDe approach will be applied on the detailed coe cients to nd the optimal threshold and soft thresholding is the method used here to eliminate the unwanted coe cients.
3.1.1 Noise Invalidation Denoising for Rician distribution
As it was mentioned in chapter 2, MRI signals are most often displayed as magnitude signal. This signal has a Rician distribution which transforms to a Rayleigh distribution in background areas and a Gaussian distribution in the image area with high SNR. Therefore, it is possible to lter the image with the NIDe approach if the noise signature has a dierent distribution. This distribution should be compatible with the signal properties of the MRI image.
Accordingly, the ANS function in this case is de ned for a random process VN with
independent identically distributed members. These members have a Rayleigh distribu-tion. Therefore the mean and variance of the ANS function g z, VN
are functions of a Rayleigh distribution CDF (instead of Gaussian distribution).
In this case the CDF in equations 2.34 and 2.35 can be formulated as equation 3.1[7]. F (z) = 1 − exp − z 2 2σ2 (3.1) Like the original NIDe approach with additive Gaussian noise, the sorting of absolute values of the noise only data with Rayleigh distribution causes the data gather in a condense area around mean value. Therefore, it is possible to de ne a noise con dence area around the mean.. This change of CDF also aects the up and down limits (equations 2.37 and 2.38) of the noise con dence region.
The new threshold can be described as the point at which noisy signal (MRI signal with noise) will have a dierent distribution than Rayleigh. In the case of the MRI image signal, this new distribution is Gaussian. Therefore by eliminating the coe cients under the obtained threshold, the signal will be noise free.
Using NIDe approach for Rician distribution, the noise bias in the image is still signal dependent and can not easily be removed. However, the bias can be easily removed from the squared signal. Therefore, the ltered image signal is transformed to the time domain and squared. Now the bias can be removed easily from the squared image. This method is based on the method proposed in [13] to create an unbiased NLM ltering method. After removing the noise bias, the square root of the signal is taken and the image can be displayed. This process can be formulated based on the following formula.
ˆ I =
q
N IDe (I)2− 2σ2 (3.2)
3.1.2 Noise Invalidation Denoising for chi-square distribution
It is also possible to apply the the NIDe approach to the squared image. As it was explained in chapter 2, the signal distribution in a squared MRI image is not Rician anymore. It changes to the non-central chi-squared distribution as it can be considered as sum of the two squared Gaussian distribution with a nite mean and variance.
To apply the NIDe method it is necessary to have a noise distribution with chi-square distribution. Therefore, the comulitative distribution function can be de ned as equation 3.3[12]. F (z) = 1 − Qk 2 √ λ,√z (3.3) where QM(a, b) is the Marcum Q-function, λis the non-centrality parameter and k is
the degree of freedom.
Again this change of CDF will aect the mean and variance of the ANS function in equations 2.34 and 2.35. Therefore, the sorted absolute value of noise coe cients are populated around the new mean which is a function of the chi-squared CDF. The lower
and upper limits are also can be calculated as the equations 2.37 and 2.38 respectively with the dierence that the new CDF of the sorted absolute values should substituted in the equations.
The algorithm nd the optimal threshold for de ning noiseless signal by comparing the noisy squared data and the noise con dence area. It is the point at which the sorted absolute valued signal leave the noise con dence area.
The noise bias in this case is no longer signal dependent and as it was mentioned before, it can be removed by easily subtracting it from the squared data coe cients directly. Finally, the root square of the data is the image can be displayed.
3.1.3 Noise variance estimation
The noise variance estimation method that is used in the implemented algorithm is based on the Local second order moment approach. This approach is described in detail in chapter 2. However, the size of the window for the neighborhood in which mean is calculated is seven.
3.1.4 Evaluation criteria
Dierent metrics have been used for the purpose of evaluating. These criteria can be categorized in two main group of quantitative and qualitative metrics. Among the rst group Root Mean Squared Error (RMSE), SSIM and PSNR and from the second category contrast value are common measures that have been used in many research[5]. Therefore, these metrics are also used in this research to compare and validate the proposed method of denoising a MRI image.
Root Mean Squared Error
RMSE calculates the dierence between original and estimated data to measure how the applied method was useful. The less the result is, the method is more valid. RMSE can be formulated as equation 3.4 [5]. RM SE = v u u t P i P j h I (i, j) − ˆI (i, j)i2 M × N (3.4)
where the I (i, j)is the noiseless data, I (i, j)ˆ is the estimated data and M and N represent the dimension of the image [5].
Contrast
It is a qualitative metric that indicates the dierence in brightness between two regions. The more the contrast is, the method can keep more details. Therefore, it is important to show that the denoising method keeps the details and also remove noise. To measure the contrast, a region of interest (ROI) should be selected. The contrast can be formulated as equation 3.5 [5].
C = Smax− Smin
Smax+ Smin (3.5)
Smaxand Sminrepresent the maximum and minimum intensity of the pixels in the
selected area[5].
Structural similarity index (SSIM)
To be able to investigate the similarity in the visual sense, structural similarity index is a good criteria to be explored. As RMSE is not an optimal measure for measuring if the denoising method keeps the detailed structures, the SSIM index is measured in this work. It can be calculated to show the similarity between the denoised (Y) and noiseless (X) image in a neighborhood with radius 11 based on the equation 3.6 [5, 25].
SSIMN(x, y) = (2µxµy+ c1) (2σxy+ c2) µ2 x+ µ2y+ c1 σ2 x+ σ2y+ c2 (3.6)
where µis the mean, σis the standard deviation and c1 and c2 are some constant related to the maximum intensity in the image. The nal SSIM is calculated by achieving the mean of all the SSIMs for the N locals. The mean and standard deviation can be calculated based on the equations 3.7 to 3.9[5, 25].
µx= N X i=1 ωixi (3.7) σx = N X i=1 ωi(xi− µx)2 !1/2 (3.8) σxy = N X i=1 ωi(xi− µx) (yi− µy) (3.9)
The value of the SSIM is between -1 to 1 and the more is this value the more the similarity is [25, 5].
3.1.5 Material
The proposed method is implemented in Matlab. The data for the simulation are taken from the Brainweb database [9]. The data consist of T1 weighted and T2 weighted MRI image of normal brain with resolution 217 × 181. Moreover, data are contaminated by adding Rician noise to the raw data taken from the brainweb site. The noise is added according to the equation 3.10 [5].
ˆ I =
q
(I + n1)2+ n22 (3.10)
where n1and n2are two random Gaussian distributed variables with N 0, σ2
and I is the original data [5].
The level of noise varies between 3%, 5% and 7% on each of the T1 and T2 weighted MRI images to show the eect of dierent noise level on dierent types of MRI images.
3.2 Simulation Results
Two NIDe methods i.e. NIDe with Rician and NIDe with chi-squared distribution are implemented and runned on the existing MRI data. In the next step, The Nide-Rician approach is combined with UNLM method to investigate the eect of combing two meth-ods. The results are described separately for each method as following.
3.2.1 Results for NIDe for Rician distributed data
The NIDe algorithm at rst executed on each MRI data with dierent noise level. As it was described before, the magnitude MRI data has a Rician distribution which is signal dependent i.e. it is Rayleigh in low SNR and Gaussian in high SNR areas respectively. This fact is illustrated in gure 3.1.
Figure 3.1: Probability density function plots for noisy image, data in the background and signal area top to bottom respectively
After estimating the noise variance and transforming to the wavelet domain, noise signature is produced. Figure 3.2 shows how ANS function causes most of the noise coe cients are gathered around mean.
Figure 3.2: ANS function on noise only data with Rayleigh distribution
After calculating the up and down boundaries, a noise con dence area is determined. The value of the landa which determine how wide is this con dence area is 5 < λ < 9. However, the best estimate for landa value in the this experiment is set to eight. The noise con dence area is illustrated in gure 3.3.
Figure 3.3: Noise con dence area for the Rician distributed Noise signature. The blue line shows the mean while the two red lines shows the upper and lower boundaries Finally the optimal threshold is obtained by comparing the sorted noisy image data with the con dence noise area. This process is illustrated in gure 3.4.
Figure 3.4: Finding the optimal threshold for the NIDe with Rician distribution method. The blue line shows the sorted noisy data and how it exists the con dence area
The other important point here is the noise bias. As it was mentioned before, the noise bias is removed by squaring the denoised image and subtracting it from this data. Then
again the nal image is obtained by taking the square root of the data.
The nal results of NIDe denoising on a T1 weighted MRI image with 3%, 5% and 7% noise, is illustrated in gures 3.5to 3.7.
Figure 3.7: Denoising result of NIDe for Rician distribution with 7% noise To illustrate the eect of the proposed denoising method on dierent types of MRI image with dierent noise contribution, The algorithm was performed on T2 weighted with dierent levels on noise too. Figures 3.8 to 3.10 show the nal result of NIDe method on T2 weighted image with 3%, 5% and 7% noise level.
Figure 3.8: NIDe approach based on Rayleigh distribution results on T2 weighted image with 3% noise level
Figure 3.9: NIDe approach based on Rayleigh distribution results on T2 weighted image with 5% noise level
Figure 3.10: NIDe approach based on Rayleigh distribution results on T2 weighted image with 7% noise level
The value of the threshold for T1/ T2weighted images with 3-7% noise level are listed in Table 3.1.
Noise level 3% 5% 7%
T1 weighted MRI image 39.9719 94.7414 147.9547 T2 weighted MRI image 223.0496 420.2845 434.4847
Table 3.1: Threshold values T1/T2 weighted MRI images at dierent noise level for Ri-cian distribution
Other statistics regarding to other levels of noise are mentioned in next chapter for the purpose of comparing the proposed method with other existing approaches. The comparison is made in terms of evaluation criteria explained in section 3.1.4.
3.2.2 Results for NIDe for chi-squared distributed data
In the second type of NIDe i.e. with the chi-squared distribution, on the squared data of the image. As it was mentioned before, the data now have non-central chi squared distribution. This is illustrated in gure3.11.
Figure 3.11: PDF plot of the squared data
To apply the NIDe method, again the noise variance estimated and wavelet transform with non decimated wavelet transform approach is applied on the data. Then the noise bias equal to 2σ2is subtracted from the approximation coe cients. The bias removal can
be done in this step because in the case of squared data, the noise bias is not anymore signal dependent. It is additive and can be easily removed by subtraction.
Then, the noise signature for chi-squared distribution noise is generated based on the ANS method. Figure 3.12shows how ANS function behave on the noise only data. By applying this function, the major part of the noise only data gather around the mean.
Figure 3.12: Eect of sorting on the noisy data with chi-squared distribution The noise con dence area for noise with chi-squared distribution with regard to its up and down boundaries is determined. The value of landa which de nes the largeness of noise con dence are in this approach is 3 < λ < 5. However, the best t for the experiment is 3.5. The result noise con dence area based on the mentioned values is illustrated in gure 3.13.
Figure 3.13: Noise con dence area for the chi-squared distribution noise. The blue line shows the mean while the red lines are the limits of this area.
Finally, the optimal threshold is found by comparing the sorted noisy squared data with the noise con dence area. When the sorted noisy data leaves the con dence area, the amplitude of data on that point is considered as the threshold for ltering the image. The process of nding this threshold is illustrated in gure 3.14.
Figure 3.14: Finding the threshold in the NIDe approach based on chi-squared distribu-tion. The blue line shows the sorted noisy data.
After soft thresholding the detailed coe cients of the squared data with the optimal threshold found in the previous step, data will be transformed to the spatial area. This is done by inverse non-decimated wavelet transform. Finally, squared root of the data are used for the purpose of displaying the image.
The nal ltered image of a T1-wieghted MRI image with 3%, 5% and 7% noise level is illustrated in gure 3.15 to 3.17.
Figure 3.17: Final Result of the denoising with NIDe on the squared data with 7% noise The algorithm is also performed on T2-wieghted MRI image with dierent noise levels. gure 3.18 to 3.20 illustrate the results of ltering on a T2 weighted MRI image with 3%, 5% and 7% noise levels respectively.
Figure 3.18: Final result of NIDe algorithm based on chi-squared distribution on the T2 weighted MRI with 3% noise.
Figure 3.19: Final result of NIDe algorithm based on chi-squared distribution on the T2 weighted MRI with 5% noise.
Figure 3.20: Final result of NIDe algorithm based on chi-squared distribution on the T2 weighted MRI with 7% noise.
The threshold values for T1/ T2weighted images with 3-7% noise level found by NIDe approach based on chi-squared distribution are listed in Table 3.2.
Noise level 3% 5% 7%
T1-weighted MRI image 1.5436e+04 8.3615e+04 2.4690e+05 T2-weighted MRI image 0.9303e+05 1.6124e+06 1.2995e+06
Table 3.2: Threshold values T1/T2 weighted MRI images at dierent noise level for chi squared distribution
Other statistics regarding the evaluation of this method is explained in the next chapter comparing with other available methods of MRI denoising.
3.2.3 Results for combining NIDe with Richian distribution with UNLM
In this part, the Nide method for image data with rician distribution is combined with UNLM. The combination is in a way that rst the NIDe approach applied and the detailed
coe cients are soft thresholded with the achieved threshold. Then to take advantage of the weighted averaging eect of the UNLM, it is applied on the thresholded detailed coe cients.
The result for T1 image 5% level is illustrated in gure 3.21.
4 Discussion and future works
In this chapter, the proposed method is compared with a well-known method of NLM and UNLM for MRI denoising. Then the advantageous and limitations of the method is explaind. Finally, some ideas that would be good to be considered for the future works on MRI image denoising are presented.
4.1 Comparison based on evaluation criteria
For the purpose of evaluating the proposed method, the NIDe approach (for both chi-squared and Rician distribution) is compared with Non Local Mean and Unbiased Non Local Mean technique. The NLM tool box on Mathwork is used for the basic NLM [17]. However, for the UNLM some modi cations are applied.
As it was mentioned in 3.1.4 segment, four evaluation criteria are used. Tables 4.1 to 4.3 show the result of each criteria for four dierent denoising methods i.e. NIDe-Rician and chi-squared, UNLM and NLM on T1 wieghted MRI image with 3%, 5% and 7% noise levels respectively.
Comparison RMSE PSNR Contrast SSIM
NIDe-Rician 21.56 34.84 0.53 0.91
NIDe-chi square 25.12 33.51 0.53 0.88
UNLM 19.91 35.53 0.51 0.90
NLM 25.56 33.36 0.52 0.85
Table 4.1: Comparison table for noise level of 3% on T1 weighted MRI image
Comparison RMSE PSNR Contrast SSIM
NIDe-Rician 29.72 32.06 0.55 0.87
NIDe-chi square 32.05 31.40 0.55 0.86
UNLM 28.45 32.41 0.53 0.86
NLM 40.03 29.47 0.52 0.79
Comparison RMSE PSNR Contrast SSIM
NIDe-Rician 37.65 30.00 0.50 0.84
NIDe-chi square 39.29 29.63 0.50 0.82
UNLM 37.10 30.13 0.54 0.81
NLM 55.19 26.68 0.49 0.73
Table 4.3: Comparison table for noise level of 7% on T1 weighted MRI image The statistics show that the UNLM method has the better results in terms of RMSE and PSNR for low levels of noise . However the structural similarity is always better for NIDe approach. This comparision proves that the visual perception in terms of detailed structures is better in NIDe-Rician method. Although the images are more blurred which is not a case with UNLM method. The blurring eect could be due to the wavelet transformation. It was mentioned before that wavelet trnasformation acts as a lter bank consisting low and high pass lters. The low pass lter can results in blurring in sharp transitions such as edges. However, in UNLM the wieghts are functions of Guassian distribution. Therefore, they prevent blurring eects. This eect of blurring can be seen in the diernece image in gure 3.6 in the chapter 3. It is obvious that the NIDe approch can not preserve edges as good as UNLM but it keeps the detailed information whereas the UNLM can not preserve them.
Furthuremore, the results from the NIDe method on chi-squared data is worse than NIDe-Rician and UNLM method. The reseaon can be due to the high degree of similarity between chi-squared noise signature and non-central chi squared distribution of the image. This results in smaller optimal threshold. Therefore, the noise can not be removed properly with this method. This problem is more obvious with the lower levels of noise that the similarity between two distributions is more.
The NIDe Rician has better results in all aspects comparing to the basic NLM method. However, the NIDe chi gets better results than NLM by increasing the noise level.
One obvious advantageous of running NIDe approach comparing to the NLM and UNLM is the algorithm speed. Although the NLM method does not do the similarity check for each pixel on the whole image (just in a searching neighborhood), it is still slow. It has to investigate for all similar patterns for all the pixels in the image. Therefore, the speed is decreasing noticably for higher resolution images. This dierence in calculation speed can be represented in table 4.4 for two methods on NIDe and UNLM in three dierent sizes.
Image size original Size half size image quarter sized image
UNLM 18.02s 4.37s 1.03s
NIDe approach 2.8s 0.83s 0.37s
Table 4.4: Speed dierence
This advantageous can also be explained in a dierent manner. The complex complex-ity can be represented by Time computational complexcomplex-ity (O(N)). It is proved that the
time computational complexity of the UNLM method is O (N × B × S) where N is the size of the image, B is the size of each neigborhood and S is the size of the search neigh-borhood [8] while for the NIDe approach the computational complexity is only depends on the length of the data (size of the image) and can be represented by the O (N).
The results for each level of noise in a T1-weighted MRI image is illustrated in gures 4.1 to 4.3 .
The same conclusion can be obtained by invastigating the T2-wieghted MRI images. The statistics are presented in Tables 4.5 to 4.7while gures 4.4 to 4.6 illustrated the resulted images for all four methods with 3%, 5% and 7% noise levels.
Comparison RMSE PSNR Contrast SSIM
NIDe-Rician 110.57 32.44 0.71 0.93
NIDe-chi square 113.94 32.17 0.70 0.91
UNLM 108.72 32.58 0.72 0.92
NLM 116.4008 31.99 0.71 0.91
Table 4.5: Comparison table for noise level of 3% on T2 weighted MRI image
Comparison RMSE PSNR Contrast SSIM
NIDe-Rician 167.64 28.82 0.72 0.86
NIDe-chi square 170.11 28.69 0.71 0.86
UNLM 143.96 30.14 0.67 0.85
NLM 189.93 27.74 0.69 0.83
Table 4.6: Comparison table for noise level of 5% on T2 weighted MRI image
Comparison RMSE PSNR Contrast SSIM
NIDe-Rician 197.64 27.39 0.72 0.82
NIDe-chi square 201.38 27.23 0.70 0.81
UNLM 178.26 22.84 0.70 0.80
NLM 216.90 26.58 0.69 0.81
As some of the noise is not removed by the NIDe approach, the UNLM method is ap-plied in wavelet domain after thresholding the detailed coe cients with optimal thresh-old. The statistics resulted from this method for a T1 wieghted MRI image with 5% noise level is illustrated in table.
It is obvious that the PSNR and RSME is improved while the criteria regarding the detailed structures are still the same. These statistics are listed in table 4.8.
Comparison RMSE PSNR Contrast SSIM
NIDe-Rician 29.72 32.06 0.55 0.87
UNLM 28.45 32.41 0.53 0.86
NIDe-Rician with UNLM 28.48 32.47 0.54 0.87 Table 4.8: The improved statistics for NIDe-Rician combined with UNLM However the image is still blurred. This can be illustarted in gure 4.7.
It should be noted that applying UNLM method on the NIDe approach decreases the speed of the algorithm.
4.2 Limitations
As it was mentioned in the 4.1 section, the following list can be considered as some limitations of proposed method.
• Blurring eect: This eect can be due to wavelet transform function.
• Lower PSNR: It is specially obvious for the NIDe with chi-squared distribution due to the high level of similarity between noise signature and MRI signal distribution. However, the latter drawback can be improved by using UNLM method or some similar methods after the soft thresholding. To decrease the blurring eect, using some other wavelet transform functions chould be useful.
4.3 Conclusion and future trends
MRI images are one of the useful and most common modalities have been used for decades in the eld of medicine. Its special characteristics such as its harmlessness, its capability to generate good visualization of some internal tissues and also functional behaviors emphasizes its critical role. However, the noise and artifacts in the nal images can aect the diagnosis. Therefore, developing signal processing methods compared to hardware based techniques or methods such as averaging can be cost eective.
NIDe method developed for MRI images tries to provide alternative solution for this problem. There are two methods based on NIDe proposed. In the rst method (NIDe Rician) the NIDe approach which is originally based on additive Gaussian noise is further developed to be consistent with the Rician noise. The second method (NIDe-chi squared) is compatible with the chi-squared distributed noise for the squared magnitude MRI data. The results of applying NIDe approach on simulated data shows that the NIDe-Rician method gains better results than the other method (NIDe-chi squared). While the results have less PSNR and are more blurred comparing the well-known UNLM method, it keeps more detailed structures. Combining the NIDe and UNLM improved the PSNR while the blurring eect still exists. The important advantageous method over exstings method is its less computational complexity.
Therefore, investigating some other wavelet transform functions can be useful to de-crease the blurring eect. On the other hand, the proposed method is only performed on the simulated data. Hence, implementing the method on some real clinical data can be a good idea to prove the eectiveness of NIDe approach.
Bibliography
[1] Ignace Lemahieu Aleksandra Pizurica, Wilfried Philips and Marc Acheroy. A ver-satile wavelet domain noise ltration technique for medical imaging. IEEE Trans-actions on Medical Imaging, 22(3):323331, 2003.
[2] Hans Knutsson Anders Eklund, Mats Andersson. Medical image analysis, mri, fmri, image registration, image segmentation. Technical report, Division of Medical In-formatics, Department of Biomedical Engineering, Linköping University, 2010. [3] Jean-Michel Morel Antoni Buades. A non-local algorithm for image denoising. In
Computer Vision and Pattern Recognition, 2005.
[4] Paul Bao and Lei Zhang. Noise reduction for magnetic resonance images via adap-tive multiscale products thresholding. IEEE Transactions on Medical Imaging, 22(9):10891099, 2003.
[5] Jyotinder S. Sahambi C. Shyam Anand. Wavelet domain non-linear ltering for mri denoising. Magnetic Resonance Imaging, 28:842861, March 2010.
[6] Iain M. Johnstone David L. Donoho. Threshold selection for wavelet shrinkage of noisy data. In Engineering in Medicine and Biology Society, 1994. Engineering Advances: New Opportunities for Biomedical Engineers. Proceedings of the 16th Annual International Conference of the IEEE, 1994.
[7] A.J. den Dekker and J. Sijbers. Advanced Image Processing in Magnetic Resonance Imaging. CRC Press, Taylor & Francis Group, 2005.
[8] Michel Kocher Dimitri Van De Ville. Sure-based non-local means. IEEE Signal Processing Letters, 16(11):973976, 2009.
[9] V. Kollokian J.G. Sled N.J. Kabani C.J. Holmes A.C. Evans D.L. Collins, A.P. Zij-denbos. Design and construction of a realistic digital brain phantom. IEEE Trans-actions on Medical Imaging, 17(3):463468, June 1998.
[10] Ron Kikinis. Guido Gerig, Olaf Kubler and Ferenc A. Jolesz. Nonlinear anisotropic ltering of mri data. IEEE Transactions on Medical Imaging, 11(2):221232, June 1992.
[11] Saeed Gazor Hossein Rabbani, Reza Nezafat. Wavelet-domain medical image de-noising using bivariate laplacian mixture model. IEEE Transactions on Biomedical Engineering, 56(12):28262837, 2009.
[12] http://en.wikipedia.org/wiki/Noncentralchi squareddistribution.
[13] Juan J. Lull Gracian Garci-Marti Luis Marti-Bonmati Montserrat Robles Jose V. Manjo n, Jose Carbonell-Caballero. Mri denoising using non-local means. Medical Image Analysis, 12:514523, February 2008.
[14] Karl Krissian and Santiago Aja-Fernandez. Noise-driven anisotropic diusion l-tering of mri. IEEE Transactions on Image Processing, 18(10):22652274, October 2009.
[15] Tianhu Lei. Statistics of Medical Imaging. Chapman and Hall/CRC, 2011.
[16] P. Perona and J. Malik. Scale-space and edge detection using anisotropic diusion. IEEE Transactions on Pattern Analysis and Machine Intelligence, 12(7):629639, 1990.
[17] Gabriel Peyre. Toolbox non-local means. Mathwork Website, January 2007. A toolbox for the non-local means algorithm.
[18] Aleksandra Pizurica. Image denoising using wavelets and spatial context modeling. Master's thesis, UNIVERSITEIT GENT, 2002.
[19] Jimmy Latt-Sara Brockstedt Ronnie Wirestam, Adnan Bibic and Freddy Stahlberg. Denoising of complex mri data by wavelet-domain ltering: Application to high-b-value diusion-weighted imaging. Magnetic Resonance in Medicine, 56:11141120, 2006.
[20] Carlos Alberola-Lopez Santiago Aja-Fernandez and Carl-Fredrik Westin. Noise and signal estimation in magnitude mri and rician distributed images: A lmmse ap-proach. IEEE Transactions on Ima, 17(8):13831398, AUGUST 2008.
[21] Xiao-Ping Zhang Soosan Beheshti, Masoud Hashemi and Nima Nikvand. Noise invalidation denoising. IEEE Transactions on Signal Processing, 58(12):60076016, December 2010.
[22] Manduchi R. Tomasi, C. Bilateral ltering for gray and color images. In Computer Vision, 839- 846.
[23] Andrew P. Bradley Kerry McMahon Dominic Kennedy Yaniv Gal, Andrew J. H. Mehnert and Stuart Crozier. Denoising of dynamic contrast-enhanced mr images using dynamic nonlocal means. IEEE Transactions on Medical Imaging, 29(2):302 310, February 2010.
[24] S.T. Yongjian Yu, Acton. Speckle reducing anisotropic diusion. IEEE Transactions on Image Processing, 11(11):1260 1270, 2002.
[25] Hamid R. Sheikh Eero P. Simoncelli Zhou Wang, Alan C. Bovik. Image quality assessment: From error visibility to structural similarity. IEEE Transactions on Image, 13(4):114, April 2004.
Nomenclature
ANS Absolute Noise Sorting DWT Discrete wavelet transform FFT Fast Fourier Transform
LMMSE Linear Minimum Square Error MRI Magnetic Resonance Imaging
nDWT Non-decimated Discrete Wavelet Transform NLM Non Local Mean
NLM Non Local Mean
NMR nuclear magnetic resonance PDE Partial Dierential Equation PDF probability density function RMSE Root Mean Squared Error ROI region of interest
SNR Signal to Noise Ratio SSIM Structural similarity index
