A Wavelet Multiscale De-Noising Algorithm Based on Radon Transform

Xueling Zhu1,*, Xiaofeng Yang2,*, Qinwu Zhou3, Liya Wang4, Fulai Yuan5 and Zhengzhong Bian3 1The School of Humanities and Social Sciences, National University of Defense Technology, Changsha, 2Department of Radiation Oncology, Emory University, Atlanta, GA 3Department of Biomedical Engineering, School of Life Science and Technology, Xi'an Jiaotong University, Xi'an, Shaanxi, 4Department of Radiology and Imaging Sciences, Emory University, Atlanta, GA, 5Department of Stomatology, Xiangya Hospital, Central South University, Changsha, 1,3,5China 2,4USA


Introduction
The rapid development of medical imaging technology and the introduction of new imaging modalities, call for new image processing methods which include specialized noise filtering, enhancement, classification and segmentation techniques. Denoising is a particularly delicate and difficult task in medical images. A tradeoff between noise reduction and the preservation of actual image features has t o b e m a d e i n a w a y t h a t e n h a n c e s t h e diagnostically relevant image content [1].
Single scale filtering methods, like linear (Wiener) or nonlinear (median) are commonly used to eliminate noise. They do not take into account different noise levels and count distributions. Linear filters can reduce the noise variance increasing the signal-to-noise ratio (SNR). However, they smooth images to the point of degrading the image contrast and detail. Nonlinear filters also reduce the noise variance, while preserving the edges to some extent, but image contrast and detail degradation still occurs. In order to avoid such problem many edge preserving filters, like the Anisotropic Diffusion (AD) Filter [2,3] have been used. Such filters respect edges by averaging pixels in the orthogonal direction of the local gradient. More recently, the trilateral filter [4], an evolution of the bilateral filter [5], has been proposed to take into account local structure in addition to intensity and geometric features.
On the other hand, multiscale filtering methods have also been proposed for image filtering. Unlike single scale filtering methods, multiscale expansions offer the possibility of separating features of interest and noise components into distinct sub-band coefficients.
In this paper we present a medical image filtering method based on the Radon and wavelet transforms. We perform Radon transform for input images to get sinograms. Then we apply 1D non-orthogonal wavelets transform along s in sinograms, and use threshold-based methods to filter it. Dissimilar to the traditional threshold schemes that apply the same threshold to the wavelet coefficients at every scale, the proposed method can get robust and adaptive noise threshold at every scale. We can take advantage of the interscale dependency information between wavelet scales to evaluate original noise variance. Finally we apply the inverse Radon transform algorithm to reconstruct the original images. Our method has been validated using MR images. The detailed steps of our method and its evaluative results are reported in the following sections.

Radon transform
The Radon transform of a 2D function is defined as: where s is the perpendicular distance of a line from the origin and  is the angle formed by the distance vector. Fourier slice theorem states that for a 2D function (,) f xy , the 1D Fourier transforms of the Radon transform along s , are the 1D radial samples of the 2D Fourier transform of (,) f xy at the corresponding angles [14].

Noise robustness
Since the Radon transform is a line integral of the image, the Radon transform of noise is constant for all of the points and directions and is equal to the mean value of the noise.
www.intechopen.com Therefore, this means zero-mean noise has no effect on the Radon transform of the image.
Here we suppose the image is a square with side length 2 ma  (0 a  ), where a is the half diagonal length of the square area in terms of pixels. We add up intensity values of the pixels to calculate the Radon transform as shown in Fig. 1. Suppose (,) f xy is a 2D discrete image whose intensity values are random variables with mean  and variance 2  From equation (7) we can get 22

The results for noise reducing
From above, we can see that the Radon transform is very robust in reducing the effect of additive noise. Here we add 3% and 10% Gaussian noise with a zero mean to the original phantom and at the same time we perform the Radon transform to get corresponding sinograms as shown in Fig. 2. Fig. 3(a) shows the profiles through the original image and the noised image (the profile position is shown in the top row of Fig. 2), Fig. 3(b) shows the profiles through their corresponding sinograms (the profile position is shown in the bottom row of Fig. 2). In Fig. 3(a) it is easy to see that the image with 10% noise has already been contaminated severely, and we cannot see any image details in the noised image. In Fig. 3 shows the transform has greatly decreased the noise in the original image. Since the Radon transform is a linear transform and is invertible, it does not lose any texture information.
And rotation of the input image corresponds to the translation of the Radon transform along  , which is more tractable.

Wavelet theory
Wavelets are mathematical functions that decompose data into different frequency components. Then each component is performed with a resolution matched to its scale [15]. The wavelet analysis decomposes a signal into a hierarchy of scales ranging from the coarsest scale to the finest one. Hence, wavelet transforms which provide representation of an image at various resolutions are a better tool for feature extraction from images. Wavelet coefficients of a signal are the projections of the signal onto the multiresolution subspaces. Wavelets are functions generated from basis function, called the mother wavelet, by dilations and translations in time (frequency) domain. If the mother wavelet is denoted by can be represented as where A and B are two arbitrary real numbers. The variables A and B represent the parameters for dilations and translations respectively in the time axis.

Dyadic wavelet transform
The discrete wavelet transform is an implementation of the wavelet transform using a discrete set of the wavelet scales and translation which obey some defined rules. For practical computations, it is necessary to discretize the wavelet transform. The scale parameter A is discretized on a logarithmic grid. The translation parameter B is then discretized with respect to the scale parameter, i.e. sampling is done on the dyadic sampling grid. A dyadic wavelet transform is a semi-discrete wavelet transform, makes scale factor binary discrete, but the displacement factor to change continuously. The discretized scale and translation parameters are given by  are the scaling functions and wavelet functions respectively.
As shown in Fig. 4, for any coarse scale 2 j , a discrete signal sequence Sf Wf  is named as the discrete wavelet transform of the original signal f . In Fig. 4, ( ) jj HG is taken for the 2 j th scale expansion of 11 () HG , that is, inserted 2 1 j  zeros among the filter coefficients.
www.intechopen.com Reconstruction is the reverse process of decomposition, so we only need to exchange the decomposition filters for the reconstruction filters. The decomposition and reconstruction of () f t is completed by the Mallat tower algorithm [16]. Fast wavelet reconstruction (16) and inverse fast wavelet reconstruction (17) are given below ,1 , Where c is the scale factor, d is the wavelet coefficients, and h and g are coefficient of the filter H and G respectively in Fig. 4.

Wavelet base choice
Orthogonal wavelet transforms have fewer coefficients at coarse scales, and therefore the transformed signal is uncorrelated across the scales. Lack of a translation invariant will make filtering by orthogonal wavelet transforms exhibit visual artifacts. The non-orthogonal wavelet offers much better edge detection because the signal is correlated across the scales. Nonorthogonal wavelet is better than orthogonal wavelet for SNR [9]. In this paper we use nonorthogonal wavelets first introduced by Mallat et al. [17]. This base wavelet is shown in Fig. 5. The shift along s in the Radon domain corresponds to the translation of the input image. So we apply a 1D wavelet transform along s .

Threshold-based filtering
Threshold-based filtering is very simple and gives a satisfactory performance. It can be divided into three steps: (1) transform the noisy signal into wavelet coefficient w , (2) employ a hard or soft threshold j t at each scale j , (3) transform back to the original domain, and get the estimated signal. Although Donoho [8] proved the optimality of a soft threshold in theory, a hard threshold has shown better results for certain applications.
As can be seen from filtering theory, during the wavelet filtering process a small threshold value will leave behind all the noisy coefficients, and subsequently, the resultant filtered image may still be noisy. On the other hand, a large threshold value will make more of the coefficients zero, and the resultant image may be blurred and have artifacts. So the threshold is key to the signal filtering effect, how the threshold is determined is critical. We choose jj tc   , where c is a constant. It is well known that for .. .  will suppress 68.26%, 95.44%, and 99.74% of its values. Therefore, imposing c between 3～4, will give good results.
From the above discussion we can see the choice of threshold depends on the noise variance j  in different scales. So the key is how to get the optimum noise variance j  , which is adaptive to different scale characteristics [18].

Noise variance j  in different scales
In order to get the noise variance in different scales we suppose white noise is  , its variance . After applying filter 1 G (see Fig. 4 where j h and j g are the unit pulse response of j H and j G respectively. h and g depend on the scaling function () t  and wavelet function () t  . They are not determined by specific scale.

Noise variance  estimation
Our algorithms require the calculation of the underlying noise variance  . In a real world application, the variance is usually unknown a priori, so it must be estimated from the data. Some papers used the dark (signal free) regions at the boundaries of each image to estimate the noise power at each wavelet scale. But this method requires manually choosing such regions and sometimes the noise in medical images is not uniformly distributed in both signal and free regions [19].
Based on [12] we are using the spatial correlation 2 (, ) Corr j n between the first two scales of wavelet transform to compute the power of ( ) where s n are the pixels through the image as shown in Fig. 1, and J This threshold can vary adaptively, because it depends on the decomposed scales and the pixels through image in the Radon transform adaptively.

Inverse radon transform
After obtaining the filtered sinogram, we perform an inverse radon transform in order to get the original image. It is defined as where R is the filtered projections. Generally, three different inverse Radon transform methods are direct inverse Radon transform (DIRT), filtered back-projection (FBP) and convolution filtered back-projection (CFBP) [20]. DIRT is computationally efficient, but it introduces some artifact. FBP based on linear filtering model often exhibits degradation in recovering from noisy data [21]. Spline-convolution filtered back-projection (SCFBP) offers better approximation performance than the conventional lower-degree formulation (e.g. piecewise constant or piecewise linear models) [22]. For SCFBP the denoised sinogram in the Radon domain is approximated in the B-spline space, while the resulting image in image domain is in the dual-spline space. We used SCFBP to propagate the deniosed sinogram back into the image space along the projection paths.

Experiments and results
Our filtering method has been evaluated by using simulation brain data. We also applied the method to filter real brain MR images. Single scale optimum linear filter -Wiener filter, www.intechopen.com traditional multiscale wavelet filter [8], anisotropic diffusion (AD) filter [2,3] and bilateral filter [5,25] were applied to these dataset in order to compare the performance of these filtering methods. Here we define an average SNR as quality metric. It is given by  (27) with the results averaged over all images and reported as mean decibels ( dB ), where [,] Xmn is the original image and ˆ [ ,] Xmn is the filtered image.
We transformed the intensity of all testing images to uint8 (0~255) before all images were filtered with the five methods. Neighborhoods of size 3 3  were used to estimate the local image mean and standard deviation of Wiener filter. A Db3 wavelet was used in the traditional wavelet filter, and the images were decompounded to four levels. The constant was imposed 3.78 c  in our method, and the images were also decompounded to four levels. For simulated and real brain data the AD filter was used when diffusion constant was 140. And the bilateral filter was performed when spatial function variance 3 s   , range function variance 20 r   .

Simulation brain MR data
We obtained the brain MR images from the McGill phantom brain database for comprehensive validation of the different filtering methods. The MR volume was constructed by subsampling and averaging a high-resolution (1-mm isotropic voxels) dataset consisting of 27 aligned scans from one individual in the stererotaxic space. The volume contains 181  217  181 voxels and covers the entire brain. Based on the realistic phantom, an MR simulator is provided to generate specified MR images [23].
We obtained T1-weighted, T2-weighted and Pd-weighted MR volumes with an isotropic voxel size of 1mm at different noise levels. The noise effect was simulated as zero mean Gaussian noise adding to the MR volume with its standard deviation equal to the noise percentage multiplied by the brightest tissue intensity [24,26]. Fig. 6 illustrates the visual assessment of filtered results on simulated T1-weighted, T2-weighted and Pd-weighted MR images with 9% noise and the comparison of the five different methods. We can easily see that bilateral filter and our method reduce more noise than other filtering methods. But the results of our method are much closer to original images without noise than those of bilateral filter. Fig. 7 shows output SNR comparison of different filters for different image types (T1-weighted, T2-weighted and Pd-weighted MR image) and noise levels (1%, 3%, 5%, 7%, 9%, 11%, 13%, 15% and 17%). The output SNRs of the five methods for three types of MR images always decrease as the noise in MR image increases, and even some output SNRs of AD filter is less than input SNRs. From the three figures it can be seen that the output SNR of our method is always higher than other methods for almost all noised MR images except Pd-weighted images with 3% and 5% noise. Moreover, output SNRs of all denoising methods are always less than input SNRs of original MR images when original images have been added 1% Gaussian noise. We also see that the more noise (over 11%) exists in MR images, the better our method performs than others.    Fig.6. From Fig.8 we can see that the filtered result of our method is closest to original standard image without noise in all denoised results. It demonstrates our method is effective in reducing noise and preserving details in MR image.
(a) Comparison of horizontal profiles between original standard, noised image and the filtered image after using our method.
(b) Comparison of horizontal profiles between original standard, the filtered images using Wiener filter, wavelet filter, AD filter, bilateral filter and our method respectively.  Fig. 9. Qualitative comparison of filtering results after using the different methods for a slice. The first column is the real brain image and the filtered results after using Wiener filter, wavelet filter, AD filter, bilateral filter and our method from the top to the bottom respectively. The second column is the corresponding residuals between the original image and the filtered image. Qualitative comparisons of ROIs are shown in the third and fourth columns. The third column shows the ROIs of real brain image, and the filtered results using Wiener filter, wavelet filter, AD filter, bilateral filter and our method from top to bottom respectively. The fourth column shows the corresponding residuals in ROIs between the original image and the filtered image.

Real MR brain data
Our filtering method also was applied to real T1 weighted MR images of the human brain. Fig. 9 illustrates the visual assessment of the filtered results of real T1-weighted MR brain images using the five different methods. By comparing the filtered results and corresponding residuals, it can easily be seen that Wiener filter makes image blurred and their residual in the whole image is almost the same except edges. Wavelet filter seems to introduce artifacts in the denoised image. AD filter erases small features and transforms image statistics due to its edge enhancement effect resulting in an unnatural image. And it can also be seen that bilateral filter and our method are excellent in reducing the noise, and enhancing the image contrast. But bilateral filter loses many small edges and features while our method preserves more of the details. The enlarged view of the region of interest (ROI) is also shown in Fig.9. It is evident that our method suppresses noise effectively while keeping more of the image details by comparing the two ROIs (50  50 voxels).

Conclusion
A wavelet domain filtering method based on the Radon transform for noise reduction in medical images was presented. We performed the Radon transform for input images, and based on our image model we validated that the Radon transform can improve the SNR and zero mean white noise has no effect on it. Wavelet analysis decomposes a signal into a hierarchy of scales ranging from the coarsest scale to the finest one. Orthogonal wavelet transforms have fewer coefficients at course scales, which prevents the transformed signal from correlated across the scales. Lack of a translation invariant makes filtering by orthogonal wavelet transforms exhibit visual artifacts. Therefore, we apply a 1D nonorthogonal wavelet transform along s , and use a threshold-based method to filter. We use a spatial correlation function to enhance significant structures and dilute noise, and then estimate the original noise variance from the first two scales. Dissimilar to traditional threshold methods that apply the same threshold to the wavelet coefficients in all scales, the proposed method provides an adaptive and robust threshold for every scale. The images in the Radon domain do not have more high-frequency parts or more edges than the corresponding images in the time domain, so our method for noise variance estimation works well. On the other side, inverse radon transform will introduce some artifacts [27]. Currently, our method produced promising results on a small segment of the image in additive Gaussian noise. Further work needs to be performed to shown if this technique is useful in multiplicative noise, on a larger wake image, and the wake detection problem. Simulated brain MR image and real brain MR image images were used to validate our method. Wiener filter, wavelet filter, AD filter and bilateral filter were compared with our method. Our method performed better than the other methods in improving SNR and in preserving the key image details and features. The experimental results showed the superiority of the proposed method as it outperformed the traditional denoising methods. Our method can also be used in other medical imaging.