# Improving the Performance of the Prony Method Using a Wavelet Domain Filter for MRI Denoising

Authors: Rodney Jaramillo; Marianela Lentini; Marco Paluszny
Journal: Computational and Mathematical Methods in Medicine
 (2014)

---

## Abstract

The Prony methods are used for exponential fitting. We use a variant of the Prony method for abnormal brain tissue detection in sequences ofT2 weighted magnetic resonance images. Here, MR images are considered to be affected only by Rician noise, and a new wavelet domain bilateral filtering process is implemented to reduce the noise in the images. This filter is a modification of Kazubek’s algorithm and we use synthetic images to show the ability of the new procedure to suppress noise and compare its performance with respect to the original filter, using quantitative and qualitative criteria. The tissue classification process is illustrated using a real sequence of T2 MR images, and the filter is applied to each image before using the variant of the Prony method.


---

## Body

## 1. Introduction

The main goal of this paper is to detect abnormal tissues inT2 weighted magnetic resonance image sequences. The identification of the distribution of relaxation times in T2 weighted magnetic resonance brain images has been used to classify tissues and as a tool for the segmentation of tumors [1–3]. This is possible because the tissue corresponding to a given pixel in a sequence of T2 weighted MRI images is characterized by a relaxation rate and therefore allows for differentiation of the tissues that conform to the region determined by a set of pixels. In this paper, we assume that the intensity at a pixel, in a sequence of T2 weighted MRI brain images, can be approximated by
(1)yi=b+∑j=1k‍Cje-iλjΔt,i=1,…,n,
where Δt is the echo time interval. Here, yi represents the intensity of gray in the pixel, n is the number of images in the sequence, and k is the number of tissues considered in the pixel, n≥2k+1. The exponents λj are called the nonlinear parameters and correspond to the relaxation rates associated with the different tissues present in the images, and the coefficients b,C1,…,Ck are called the linear parameters and are related to the background noise of the image and the proportions of each one of the k tissues in the pixel, respectively. The background noise is the noise which comes from all possible sources, some of which are unknown, and usually corresponds to patient movement and the image acquisition process. Martín-Landrove et al. [3, 4] proposed the model and used a variant of Prony’s method to find the linear and nonlinear parameters. We have used two methods to compute the parameters in the exponential model (1): the variant of the Prony method already mentioned [3–6] and the variable projection method of Golub and Pereyra [7–9]. Prony-type methods are widely used for estimating a signal frequency component; a discussion of such methods can be found in [10–13]. On the other hand, the variable projection method has been used in many applications in fields like electrical engineering, medical image processing, chemical sciences, environmental sciences, magnetic resonance applications, and so forth. A comprehensive bibliographic reference can be found in [8], (144 items).In [5] Paluszny et al., comparing the solution obtained by the variant of the Prony method and the solution given by the variable projection method, it was noted that both methods successfully converge in the case of low level noisy data, and none produces a satisfactory solution for high noise levels; the variable projection method proved to be more robust in the presence of higher noise levels but required ten times more computational time to get results comparable to those of the variant of the Prony method.The exponential fit given by (1) is an ill-conditioned problem, and this means that relatively small changes in the data may result in large variations in the linear and nonlinear parameters; therefore, a filter is required to modify the data before using any of the considered methods. In particular, [6] provides an analysis of the sources of instability for the variant of the Prony method.In the literature, there are a number of different filters to enhance MRI. We chose to modify the filter designed by Kazubek [14], to get a wavelet domain filter for removing noise on each image, before applying the Prony method, and trying to keep the computation time low.


## 2. Wavelet Domain Filters for MRI Denoising

### 2.1. The Two-Dimensional Wavelet Transform

Let us consider an imagef, and let f(x,y) be the gray intensity corresponding to the pixel given by coordinates (x,y). The discrete wavelet transform (DWT) is a decomposition of f in two functions:
(2)f(x,y)≈Fs(x,y)+Fw(x,y),
where Fs is a representation of f in terms of shifted and dilated versions of a low-pass scaling function ϕ(x,y) and Fw is a representation of f in terms of shifted and dilated versions of a prototype bandpass wavelet function ψ(x,y). The coefficients which determine the components Fs and Fw are called scaling coefficients and wavelet coefficients, respectively. Once the functions ϕ(x,y) and ψ(x,y) have been chosen, the components Fs and Fw are uniquely determined by the coefficients; therefore, it is said that the discrete wavelet transform consists of the two sets of coefficients. The scaling coefficients allow for a low resolution reconstruction of the image and the wavelet coefficients deal with the reconstruction of the details.There is a large amount of bibliographic information about the theoretical as well as computational aspects of the DWT; we will give, here, a brief description to point out the relevant issues regarding the design of a filter, in the wavelet domain, for MR images.


### 2.2. Wavelet Domain Filters for MRI

When we consider a noisy image and its DWT, it is supposed that the noise in the image is determined by the noise associated with the scaling and wavelet coefficients. A wavelet domain filter corresponds to two procedures, one for the noise reduction in the scaling coefficients and the other for dealing with the noise associated with the wavelet coefficients. The filtered image is the result of applying the inverse wavelet transform to the modified coefficients. The differences between various proposals for wavelet domain filters lie in the algorithms to be applied to each one of the two sets of coefficients; see, for example, [14–22].Although the noise in MRI comes from various sources, we consider that the perturbation in magnitude follows a Rice distribution, [15, 23]. Kazubek considers that most of the energy used in the image generation corresponds to the scaling coefficients and hence these contribute the most to the perturbation of the image in the presence of noise; therefore, Kazubek proposes to apply a filter to reduce the noise in the scaling coefficients under the hypothesis that the noise is Rician.According to some authors, numerical experiments show that noise in the wavelet coefficients can be approximated by a Gaussian distribution; for this reason we made a satisfactory noise reduction in these coefficients using a Wiener filter (consult [14, 15, 21]) in addition to the Kazubek filter in the scaling coefficients.The arithmetic mean of a random variable which follows a Rician distribution is a biased estimator of the value without noise. The next section shows a brief description of the Rician noise and a formula for reducing the bias between the value without noise and the arithmetic mean. We will use this formula later in the filter, on the scaling coefficients.


## 2.1. The Two-Dimensional Wavelet Transform

Let us consider an imagef, and let f(x,y) be the gray intensity corresponding to the pixel given by coordinates (x,y). The discrete wavelet transform (DWT) is a decomposition of f in two functions:
(2)f(x,y)≈Fs(x,y)+Fw(x,y),
where Fs is a representation of f in terms of shifted and dilated versions of a low-pass scaling function ϕ(x,y) and Fw is a representation of f in terms of shifted and dilated versions of a prototype bandpass wavelet function ψ(x,y). The coefficients which determine the components Fs and Fw are called scaling coefficients and wavelet coefficients, respectively. Once the functions ϕ(x,y) and ψ(x,y) have been chosen, the components Fs and Fw are uniquely determined by the coefficients; therefore, it is said that the discrete wavelet transform consists of the two sets of coefficients. The scaling coefficients allow for a low resolution reconstruction of the image and the wavelet coefficients deal with the reconstruction of the details.There is a large amount of bibliographic information about the theoretical as well as computational aspects of the DWT; we will give, here, a brief description to point out the relevant issues regarding the design of a filter, in the wavelet domain, for MR images.


## 2.2. Wavelet Domain Filters for MRI

When we consider a noisy image and its DWT, it is supposed that the noise in the image is determined by the noise associated with the scaling and wavelet coefficients. A wavelet domain filter corresponds to two procedures, one for the noise reduction in the scaling coefficients and the other for dealing with the noise associated with the wavelet coefficients. The filtered image is the result of applying the inverse wavelet transform to the modified coefficients. The differences between various proposals for wavelet domain filters lie in the algorithms to be applied to each one of the two sets of coefficients; see, for example, [14–22].Although the noise in MRI comes from various sources, we consider that the perturbation in magnitude follows a Rice distribution, [15, 23]. Kazubek considers that most of the energy used in the image generation corresponds to the scaling coefficients and hence these contribute the most to the perturbation of the image in the presence of noise; therefore, Kazubek proposes to apply a filter to reduce the noise in the scaling coefficients under the hypothesis that the noise is Rician.According to some authors, numerical experiments show that noise in the wavelet coefficients can be approximated by a Gaussian distribution; for this reason we made a satisfactory noise reduction in these coefficients using a Wiener filter (consult [14, 15, 21]) in addition to the Kazubek filter in the scaling coefficients.The arithmetic mean of a random variable which follows a Rician distribution is a biased estimator of the value without noise. The next section shows a brief description of the Rician noise and a formula for reducing the bias between the value without noise and the arithmetic mean. We will use this formula later in the filter, on the scaling coefficients.


## 3. Removing the Bias from Rice Distributed Data

In a magnetic resonance image, the gray intensity corresponding to a given pixel is the magnitude of a complex valued function whose real and imaginary components are distributed as Gaussians with the same standard deviation. The different approaches for MRI noise removal can be divided into three classes; some methods act on the real and imaginary components, most of the methods deal with the square of the magnitude, and a third kind of methods acts on the magnitude of the complex image. We are interested in the latter class of methods because our data are in magnitude.In [14], Kazubek implements a wavelet domain filter which acts on the magnitude of the image and whose relevant issues are the use of the Haar system of orthogonal functions in the wavelet transform and a formula to apply on the scaling coefficients that is used to suppress the bias discussed in Section 2.Next, we present an analysis of Kazubek's formula for bias removal and propose a modification of his algorithm by applying the bilateral filter to the scaling coefficients [24].LetA1 and A2 be the mean values of two random variables μ1 and μ2, respectively, which are the real and imaginary components of a complex image. If μ1 and μ2 are distributed as Gaussians with standard deviation σ, then the probability density function for the magnitude, μ=μ12+μ22, is a Rice distribution given by
(3)pμ(μ∣A)=μσ2e-((μ2+A2)/2σ2)I0(Aμσ2),
where I0 is the modified Bessel function of the first kind and zeroth order and A=A12+A22 is the underlying noise-free signal amplitude. Let m be the mean value of μ; then,
(4)mσ=v(x)∶=π2[(1+x22)I0(x24)+x22I1(x24)]e-x2/4,
where x=A/σ and I1 is the modified Bessel function of the first kind and first order [15].Letz=m/σ; we can get a numerical approximation of the inverse x=v-1(z) for 0≤x≤50.We have thatx∈(0,50] implies -∞<20log⁡10x<20log⁡1050≈33.98; then, the inverse function v-1 corresponds to a signal whose noise measurement runs from -∞ to 34 decibels, (dB).Bessel functionsI0 and I1 satisfy I0′=I1; then,
(5)v′(x)=π2e-x2/4[3x2I1(x24)+x2I0(x24)+x34I1′(x24)π2e-x2/4w-x34I0(x24)]
and using the power series expansion of functions I0, I1, and I1′, we get
(6)v′(x)=π2e-x2/4[∑k=0∞‍x4k+12(64k)(k!)2π2e-x2/4e+∑k=0∞‍x4k+316(k+1)(64k)(k!)2].
(Details of the proof can be found in the Appendix.) Therefore, v(x) is a monotonically increasing function in the interval [0,+∞) and v′(0)=0.Using polynomial approximations of the functions(7)x24I0(x24)e-x2/4,x24I1(x24)e-x2/4
for x>15 (consult Abramowitz and Stegun [25]), we get
(8)lim⁡x→∞v(x)x=L,
where
(9)L=0.999999+Er
and |Er|≤5.139(10-7). Then, for large x, v(x) behaves as the function
(10)H(x)=L2x2+v(0)2.On the other hand, we have(11)H(0)=v(0),H′(0)=v′(0)=0,lim⁡x→∞H(x)x=L,H′(x)>0,∀x>0,
and a straightforward computation shows that
(12)H-1(z)=z2L2-v(0)2L2.This function is a particular case of functionaz2+b with a>0 and b<0.The idea to get a numerical approximation ofv-1(z) consists in considering a function Φ(z) such that function
(13)F(z)=az2+b+Φ(z)
and its derivative satisfy similar asymptotic features of v-1(z) and [v-1]′(z), respectively. Let c<0 and let d<0; if we consider Φ(z)=cedz, then
(14)F′(z)>0,lim⁡z→∞F(z)z=a<∞.H
1
(
z
)
=
a
z
2 is a monotonically increasing function in the interval [0,∞) with H1(0)=0. Function H2(z)=-b-cedz takes positive values and is a monotonically decreasing function in the same interval. Then, there is a unique positive root, z0, of
(15)az2=-(b+cedz)
and if z≥z0, then
(16)az2+b+cedz≥az02+b+cedz≥az02+b+cedz0=0.It follows that[z0,∞) is the domain of
(17)F(z)=az2+b+cedz,F(z0)=0,lim⁡z→z0+F′(z)=∞.To get the parametersa, b, c, and d, we perform a least square fit to the discretized curve F(z) for xj=j(0.1) and zj=v(xj), j=1,…,500. This problem was solved numerically using the MATLAB functionlsqcurvefit. With MATLAB R2008a and taking the initial values [a0,b0,c0,d0]=[1,-1,-2,-1], we get the optimal solution: a=1.0000108, b=-1.0122372, c=-2.7102422, and d=-1.2598921. Figure 1(a) shows the function F(z) as a numerical approximation of V-1(z) and Figure 1(b) shows the quality if the solution given by the MATLAB functionlsqcurvefit. In MRI, as in acoustic, electricity, and telecommunications, the signal to noise ratio (SNR) is frequently used to measure, in decibels, the noise impact on a signal. Let r and t be two gray level images, a reference image r(x,y) and a noisy image t(x,y), of size nxny. The SNR is computed as
(18)SNR=10log⁡10(∑1nx∑1ny[r(x,y)]2∑1nx∑1ny[r(x,y)-t(x,y)]2).(a) shows the functionF as a numerical approximation of v-1. (b) is a plot of the function xj-F(zj).
(a)
(b)Nowak [15] suggests that a SNR greater than 10 dB corresponds to a noisy image, whose Rice distribution can be approximated by a Gaussian distribution. We will consider a method to suppress the bias in our MR images using the inverse function discussed in this section when SNR is less than 34 dB and a simple wavelet domain filter for the other cases.


## 4. Implementing a New Wavelet Domain Bilateral Filter for MRI

### 4.1. Formula for the Scaling Coefficients

The componentFs(x,y) in (2) is a linear combination of functions that belong to a specific family of orthogonal functions; each scaling coefficient depends on three parameters: L, the level of decomposition in the discrete wavelet transform, and the pair (k1,k2) which characterizes the support of the corresponding orthogonal function. Let f(x,y) be the function in (2); if we apply the discrete wavelet transform using the Haar system of orthogonal functions, then the scaling coefficients are given by
(19)cL,k1,k2=12L∑i=02L-1‍∑j=02L-1‍f(k12L+i,k22L+j).Details can be found in [26].


### 4.2. Wiener Filter for the Wavelet Coefficients

As the scaling coefficients, the wavelet coefficients depend on three parametersdl,h1,h2, 1≤l≤L. For simplicity, let us write d=dl,h1,h2 and suppose that the noise in the wavelet coefficients may be approximated by a Gaussian distribution [14]. Then, we attenuate the contribution of this coefficient by
(20)d^=αdd,0≤αd≤1.According to Nowak, we assume that the noise wavelet coefficient is an unbiased estimator of the value of the wavelet coefficient in the noise-free case and denote its mean byδ=E[d]. The filter weight αd is chosen by minimizing E[(δ-αdd)2]. In fact, if σd2 is the variance of the wavelet coefficient d, then
(21)E[(δ-αdd)2]=δ2-2αdδ2+αd2E[d2]=δ2-2αdδ2+αd2(σd2+δ2);
therefore αd is defined by
(22)αd=δ2δ2+αd2.Letσ*2 be an estimate of σd2; then,
(23)αd≈α*=E[d2]-σ*2E[d2].Usually the parameterσ*2 is approximated by τσ2, where τ≥1 is a new parameter to be chosen for each image and σ is as defined in Section 2. In previous works [15, 18, 21], the value τ=2 has proven to be convenient. Hence, the filter weight in (20) is defined by
(24)αd=(E[d2]-2σ2E[d2])+={E[d2]-2σ2E[d2],ifE[d2]>2σ2,0,otherwise.


### 4.3. Bilateral Filtering

The bilateral filter is a nonlinear filter for images, proposed by Tomasi and Manduchi [24], whereby contours are successfully recovered from noisy images. Let x∈I be a pixel in image I; the bilateral filter takes a sum with weights on the pixels in a local neighborhood of x, Nx. These weights depend on the spatial distance and the intensity in each pixel in the neighborhood. The response of the bilateral filter in the pixel is given by
(25)I~(x)=1C∑y∈Nx‍WS(x,y)WR(x,y)I(y),
where
(26)WS(x,y)=exp⁡-∥x-y∥22ρd2,WR(x,y)=exp⁡-|I(x)-I(y)|22ρr2,C=∑y∈Nx‍WS(x,y)WR(x,y).The performance of the bilateral filter depends on the choice ofρd, ρr and the size of the neighborhood Nx. Let x be a pixel located close to an edge which separates two regions of the image. When the pixel y is in the same region as x, WR(x,y) is close to one; on the other hand, when y is in the other region, then WR(x,y) is close to zero if the intensities I(x) and I(y) differ greatly and then the filter tends to preserve the pixel intensity due to the effect of the WR(x,y) components.In our proposal the bilateral filter is applied to the two-dimensional array consisting of the scaling coefficients of the wavelet transformI(k)=I(k1,k2)=cL,k1,k2. We follow Anand and Sahambi [20] and choose the neighborhood Nx to be a 15×15 matrix centered at x and the parameters ρd=5 and ρr=1.5σ, where σ is the standard deviation given in Section 3.


### 4.4. Algorithm for MRI Denoising

(1)
Compute an approximation of the variance,σ2, choosing a rectangular q1×q2 neighborhood Nη, in the background of the image [15]:
(27)σ2≈σ~2=12q1q2∑u∈Nη‍[I(u)]2.
(2)
Perform a levelL=3 wavelet decomposition using the Haar system of orthogonal functions.
(3)
Perform a bias correction of the level 3 scale coefficients using the inverse function considered in Section3.
From (19) we see that
(28)12LcL,k1,k2
is an approximation to the expected valuem given in (4) so that
(29)zN=1σ~2LcL,k1,k2
is an approximation of thez value defined in Section 2. Using (4) and (17), we get
(30)xσ~≈F(zN).
Now we compute the new scaling coefficientsc¯L,k1,k2 by
(31)c¯L,k1,k2=σ~2LF(zN).
(4)
Apply the bilateral filter to the two-dimensional arrayJ(k1,k2), given by the new scaling coefficients: J(k1,k2)=c¯L,k1,k2.
(5)
Obtain a provisional denoised image,IN, computing the inverse wavelet transform using the new scaling coefficients and the unchanged wavelet coefficients.
(6)
Perform a levelL=4 wavelet decomposition of IN using the Daubechie system with four vanishing moments.
(7)
Perform a correction of the wavelet coefficients given in (6) using (20) and (24).
(8)
Get a denoised image by computing the inverse wavelet transform on the set of coefficients given by the scaling coefficients in (6) and the new wavelet coefficients in (7).


## 4.1. Formula for the Scaling Coefficients

The componentFs(x,y) in (2) is a linear combination of functions that belong to a specific family of orthogonal functions; each scaling coefficient depends on three parameters: L, the level of decomposition in the discrete wavelet transform, and the pair (k1,k2) which characterizes the support of the corresponding orthogonal function. Let f(x,y) be the function in (2); if we apply the discrete wavelet transform using the Haar system of orthogonal functions, then the scaling coefficients are given by
(19)cL,k1,k2=12L∑i=02L-1‍∑j=02L-1‍f(k12L+i,k22L+j).Details can be found in [26].


## 4.2. Wiener Filter for the Wavelet Coefficients

As the scaling coefficients, the wavelet coefficients depend on three parametersdl,h1,h2, 1≤l≤L. For simplicity, let us write d=dl,h1,h2 and suppose that the noise in the wavelet coefficients may be approximated by a Gaussian distribution [14]. Then, we attenuate the contribution of this coefficient by
(20)d^=αdd,0≤αd≤1.According to Nowak, we assume that the noise wavelet coefficient is an unbiased estimator of the value of the wavelet coefficient in the noise-free case and denote its mean byδ=E[d]. The filter weight αd is chosen by minimizing E[(δ-αdd)2]. In fact, if σd2 is the variance of the wavelet coefficient d, then
(21)E[(δ-αdd)2]=δ2-2αdδ2+αd2E[d2]=δ2-2αdδ2+αd2(σd2+δ2);
therefore αd is defined by
(22)αd=δ2δ2+αd2.Letσ*2 be an estimate of σd2; then,
(23)αd≈α*=E[d2]-σ*2E[d2].Usually the parameterσ*2 is approximated by τσ2, where τ≥1 is a new parameter to be chosen for each image and σ is as defined in Section 2. In previous works [15, 18, 21], the value τ=2 has proven to be convenient. Hence, the filter weight in (20) is defined by
(24)αd=(E[d2]-2σ2E[d2])+={E[d2]-2σ2E[d2],ifE[d2]>2σ2,0,otherwise.


## 4.3. Bilateral Filtering

The bilateral filter is a nonlinear filter for images, proposed by Tomasi and Manduchi [24], whereby contours are successfully recovered from noisy images. Let x∈I be a pixel in image I; the bilateral filter takes a sum with weights on the pixels in a local neighborhood of x, Nx. These weights depend on the spatial distance and the intensity in each pixel in the neighborhood. The response of the bilateral filter in the pixel is given by
(25)I~(x)=1C∑y∈Nx‍WS(x,y)WR(x,y)I(y),
where
(26)WS(x,y)=exp⁡-∥x-y∥22ρd2,WR(x,y)=exp⁡-|I(x)-I(y)|22ρr2,C=∑y∈Nx‍WS(x,y)WR(x,y).The performance of the bilateral filter depends on the choice ofρd, ρr and the size of the neighborhood Nx. Let x be a pixel located close to an edge which separates two regions of the image. When the pixel y is in the same region as x, WR(x,y) is close to one; on the other hand, when y is in the other region, then WR(x,y) is close to zero if the intensities I(x) and I(y) differ greatly and then the filter tends to preserve the pixel intensity due to the effect of the WR(x,y) components.In our proposal the bilateral filter is applied to the two-dimensional array consisting of the scaling coefficients of the wavelet transformI(k)=I(k1,k2)=cL,k1,k2. We follow Anand and Sahambi [20] and choose the neighborhood Nx to be a 15×15 matrix centered at x and the parameters ρd=5 and ρr=1.5σ, where σ is the standard deviation given in Section 3.


## 4.4. Algorithm for MRI Denoising

(1)
Compute an approximation of the variance,σ2, choosing a rectangular q1×q2 neighborhood Nη, in the background of the image [15]:
(27)σ2≈σ~2=12q1q2∑u∈Nη‍[I(u)]2.
(2)
Perform a levelL=3 wavelet decomposition using the Haar system of orthogonal functions.
(3)
Perform a bias correction of the level 3 scale coefficients using the inverse function considered in Section3.
From (19) we see that
(28)12LcL,k1,k2
is an approximation to the expected valuem given in (4) so that
(29)zN=1σ~2LcL,k1,k2
is an approximation of thez value defined in Section 2. Using (4) and (17), we get
(30)xσ~≈F(zN).
Now we compute the new scaling coefficientsc¯L,k1,k2 by
(31)c¯L,k1,k2=σ~2LF(zN).
(4)
Apply the bilateral filter to the two-dimensional arrayJ(k1,k2), given by the new scaling coefficients: J(k1,k2)=c¯L,k1,k2.
(5)
Obtain a provisional denoised image,IN, computing the inverse wavelet transform using the new scaling coefficients and the unchanged wavelet coefficients.
(6)
Perform a levelL=4 wavelet decomposition of IN using the Daubechie system with four vanishing moments.
(7)
Perform a correction of the wavelet coefficients given in (6) using (20) and (24).
(8)
Get a denoised image by computing the inverse wavelet transform on the set of coefficients given by the scaling coefficients in (6) and the new wavelet coefficients in (7).


## 5. Validation of the MRI Filter Applied to Simulated Images

We compare the performances of the modified and the original Kazubek's filter using a synthetic image generated with MATLAB, see Figure2, and introduce the noisy images by adding Rician noise to the synthetic image. The Rician noise was generated as
(32)Je(m,n)=(J(m,n)+e1)2+e22,
where J(m,n) is the true signal and e1 and e2 are random numbers from a Gaussian distribution with mean zero and standard deviation σ; five levels of σ were used σ=[1,2,5,8,12], and the gray values in the image are between 0 and 88. To compare the performances of the two filters, we compute five measurements that appear more frequently in the literature: signal to noise ratio (SNR), peak signal to noise ratio (PSNR), root mean square error (RMSE), mean absolute error (MAE), and the structural similarity index (SSIM). Given two gray level images, a reference image r(x,y) and noise image t(x,y), of size nxny, the quantities PSNR, RMSE, and MAE are given by
(33)PSNR=10log⁡10(nxnymax⁡⁡[r(x,y)]2∑1nx∑1ny[r(x,y)-t(x,y)]2),RMSE=1nxny∑1nx‍∑1ny‍[r(x,y)-t(x,y)]2,MAE=1nxny∑1nx‍∑1ny‍|r(x,y)-t(x,y)|.Figure 2
Synthetic image, generated with MATLAB.The structural similarity index in a regionΩ is estimated as
(34)SSIMΩ(r,t)=(2μtμr+c1)(2θr,t+c2)(μr2+μt2+c1)(θr2+θt2+c2),
where μr and μt are the means over Ω of r and t, respectively, θr and θt are the corresponding variances, θr,t is the covariance of r and t, and the constants c1 and c2 are given by c1=0.01*255=2.55 and c2=0.03*255=7.65; see [27]. If the images are divided in s regions, Ω1,…,Ωs, the SSIM for the two images r(x,y) and t(x,y) is defined by
(35)SSIM=1s∑j=1s‍SSIMΩj(r,t).The resultant SSIM index is a value between −1 and 1, and value 1 is only reachable in the case of two identical data sets. Quantities SNR, PSNR, RMSE, and MAE are currently used in computer and medical sciences to compare two images; however, a good quantitative performance could be not consistent, in some cases, with visual perception. SSIM has become a more convenient way to compare two images, when PSNR and RMSE show indistinguishable results for the perception of the human eye.Tables1, 2, 3, 4, and 5 show the errors between the noisy images and the filtered images using the quantities SNR, PSNR, RMSE, MAE and SSIM, respectively. The first row in each table corresponds to the noisy image and various σ values, the second row corresponds to the filtered image using Kazubek's filter, and the third row corresponds to the filtered image using the modified filter. For each σ level, our modified Kazubek's algorithm shows an improved performance with respect to the original algorithm using all the criteria to measure the error. The performance of both filters is nearly indistinguishable in the cases σ=1,2 and we get better filtered images as the noise level increases to σ=5, 8, and 12, see Figure 3.Table 1
SNR between the original and the denoised images corresponding to the two filters andσ values.
σ
=
1
σ
=
2
σ
=
5
σ
=
8
σ
=
12
Original errors
29.77391
23.7840
16.0170
12.0372
8.9663
Kazubek's algorithm
33.7311
27.9346
20.4575
17.0542
14.1045
Modification to Kazubek's algorithm
34.0501
28.3961
21.4209
18.0140
15.3410Table 2
PSNR between the original and the denoised images corresponding to the two filters andσ values.
σ
=
1
σ
=
2
σ
=
5
σ
=
8
σ
=
12
Original errors
37.2277
31.5732
24.4287
21.4391
18.5411
Kazubek's algorithm
41.1923
35.7037
28.6968
25.4961
22.6839
Modification to Kazubek's algorithm
41.5063
36.1647
29.7469
26.9269
24.5610Table 3
RMSE between the original and the denoised images corresponding to the two filtersσ values.
σ
=
1
σ
=
2
σ
=
5
σ
=
8
σ
=
12
Original errors
1.2495
2.4935
6.1831
10.0105
14.9444
Kazubek's algorithm
0.7914
1.5401
3.6277
5.3316
7.4354
Modification to Kazubek's algorithm
0.7631
1.4618
3.2590
4.8087
6.5236Table 4
MAE between the original and the denoised images corresponding to the two filters andσ values.
σ
=
1
σ
=
2
σ
=
5
σ
=
8
σ
=
12
Original errors
1.0522
2.1056
5.2262
8.4441
12.6008
Kazubek's algorithm
0.5491
1.0710
2.4354
3.5928
5.1168
Modification to Kazubek's algorithm
0.5247
1.0195
2.2510
3.3490
4.6331Table 5
SSIM between the original and the denoised images corresponding to the two filters andσ values.
σ
=
1
σ
=
2
σ
=
5
σ
=
8
σ
=
12
Original errors
0.9094
0.7705
0.5543
0.4482
0.3640
Kazubek's algorithm
0.9864
0.9567
0.8482
0.7489
0.6532
Modification to Kazubek's algorithm
0.9878
0.9614
0.8636
0.7669
0.6775(a) The noisy image corresponding toσ=2. (b) The filtered image using the modification of Kazubek's algorithm.
(a)
(b)


## 6. Numerical Results

To see the effect introduced by the filter on a real magnetic resonance image, as a part of a tissue segmentation process, we consider the set of images referenced in [3]. In this case, we have a set of eight echoes (magnetic resonance images of the same axial section) from a sequence of T2-MRI, with an echo time of 44 milliseconds. According to [2, 3] we can see a part of the brain affected by tumoral tissues (ROI1) and another region showing no visible impairment (ROI2); regions ROI1 and ROI2 are called the region of lesion and region of control, respectively. Figure 4 shows the echo corresponding to n=6, from the set of images. For each pixel in a region the variant of Prony method is applied, to find the linear and nonlinear parameters defined in (1) with n=8 and ΔT=0.044; this process is performed for k=1,2, and 3, and we choose the ones minimizing the residual. In Table 6, we show the computation times of the classification process, using our implementation of the variant of the Prony method.Table 6
Computation times in seconds corresponding to the tissue classification process.
Method
ROI1
ROI2
Prony without filtering
9.39
9.29
Prony with filtering
11.40
11.29(a) ROI1 and (b) ROI2.
(a)
(b)For each region of interest, the frequency diagram is computed as follows: we divided the interval[0,30] into 100 bins of equal length; ςj=[ςj,ςj+1], 1≤j≤100, since the brain tissues relaxation rates belong to this interval, as we can see in [2]. For the pixels in a specific region and a selected subinterval ςj, we add the linear parameters whose corresponding nonlinear parameters belong to ςj, and this process is done for 1≤j≤100, to get a diagram of frequencies. Each frequency diagram is normalized to obtain a probability density function.In Figure5(a), we see the solution obtained in [3], which was calculated by the variant of the Prony method. Figure 5(b) shows the results using our implementation of the variant of the Prony method, in red the probability density corresponding to the control region, and in blue the corresponding to the region of lesion. The differences in the shape of these graphics are caused because our implementation of the method is not the same as that reported in [3]. Figure 5(c) shows the result of applying the MRI filter designed in this paper, on each one of the eight images, before running the variant of the Prony method. The filtering process produces a curve with less dispersion for the region ROI2 and a softer edge in the curve corresponding to region ROI1.Probability densities corresponding to ROI1 and ROI2, in blue and red, respectively. (a) Graphic reported in [3]. (b) Solution obtained by Prony’s method without a filtering process. (c) Solution obtained when we apply the wavelet domain filter to the images before using the Prony method.
(a)
(b)
(c)


## 7. Conclusion

This paper takes as its starting point the wavelet domain filter for MRI developed by Kazubek. We improved the original filter introducing the bilateral filter in the wavelet domain. The resulting filter exhibits a better performance and requires only a small increase in computational time of the whole tissue classification process. The new filter is applied to real brain MRI in a process of brain tissue classification proposed by Paluszny et al. [3], getting different results to those previously reported, which lead to an improvement in tissue identification with T2 relaxation techniques. Our results are validated with the main quality criteria for image denoising.



---

*Source: 810680-2014-04-14.xml*