Noise Removal for Medical X-ray images in Multiwavelet Domain
Noise Removal for Medical X-ray images in Multiwavelet Domain Ling Wang, Jianming Lu, Yeqiu Li, Takashi Yahagi Yahagi, Lu & Sekiya Lab., Graduate School of Science and Technology,Chiba University, 1-33 Yayoi-cho, Inage-ku, Chiba, 263-8522, Japan Email:[email protected]
Abstract— When a signal is embedded in an additive Gaussian noise, its estimation is often done by finding a wavelet basis that concentrates the signal energy in few coefficients and then thresholding the noisy coefficients. However, in many practical problems such as medical X-ray image, astronomical and lowlight images, the recorded data is not modeled by Gaussian noise but as the realization of a Poisson process. Multiwavelet is a new development to the body of wavelet theory. Multiwavelet simultaneously offers orthogonality, symmetry and short support which are not possible in scalar 2-channel wavelet systems. After reviewing this recently developed theory, a new theory and algorithm for denoising medical X-ray images using multiwavelet multiple resolution analysis (MRA) are presented and investigated in this paper. The proposed covariance shrink (CS) method is used to threshold wavelet coefficients. The form of thresholds is carefully formulated which is the key to more excellent results obtained in the extensive numerical simulations of medical image denoising compared to conventional methods.
I. I NTRODUCTION Noise in medical X-ray images [1] is primarily categorized into quantum mottle, which is related to the number of incident X-rays, and due to artificial noise which is due to the grid, etc. The effect of quantum mottle is manifested as an increase in the graininess of an X-ray image as the dose is reduced. Therefore, noise reduction is of great significance in medical X-ray images. The noise of X-ray images obeys a Poisson law and, hence, is highly dependent on the underlying light intensity pattern being imaged [2]. Over the past decade, wavelet thresholding has received a lot of attention from researchers in image denoising. Scalar orthogonal wavelets with a single mother-wavelet function have played an important role in denoising [4]. The denoising is done on the detail coefficients of wavelet transform, it has been shown that this algorithm offers the advantages of smoothness and adaptation. However, it is known that there is a limitation on the time-frequency localization of scalar wavelet functions. Multiwavelet has been developed recently by using translates and dilates of more than one mother wavelet function [9]. They are known to have several advantages over scalar wavelets. In addition, Strela et al.[9] claimed that multiwavelet hard thresholding offers better results than the tranditional scalar wavelet hard thresholding. In this paper, a new shrink theory and denoising algorithm for medical X-ray images using multiwavelet multiple resolu-
tion analysis (MRA) is presented and investigated. We propose a new shrinkage method—covariance shrink (CS) method. The CS method is subband adaptive—a separate threshold is computed for each detail subbad of multiwavelet transform. We propose near-optimal thresholds and a more suitable image denoising by extensive numerical simulations.We show in this paper that multiwavelet denoising gives better results than the conventional scalar wavelet denoising and conventional multiwavelet denoising. II. S CALAR WAVELET AND MULTIWAVELET DENOISING Scalar wavelet and multiwavelet denoising attempt to remove the noise present in the signal while preserving the signal characteristics, regardless of its frequency content. they involve three steps: 1) A linear forward wavelet transform. 2) Nonlinear shrinkage denoising. 3) A linear inverse wavelet transform. Experimental results show that small wavelet coefficients are dominated by noise, while coefficients with a large absolute value carry more signal information than noise. Replacing noisy coefficients (small coefficients below a certain threshold value) by zeros and inverse wavelet transformation may lead to a reconstruction that has less noise. This method is indeed effective and thresholding is a simple and efficient method for noise reduction. The most well-know thresholding methods include VisuShrink [14] and SureShrink [15]. These threshold choices enjoy asymptotic minimax optimalities over function spaces such as Besov spaces. BayesShrink for denoising determines the threshold for each subband assuming a Generalized Gaussian Distribution (GGD). The results obtained for image denoising using BayesShrink are much better than VisuShrink and SureShrink. In this paper, a new multiwavelet thresholding method based on the two types of multivariate thresholding approach is proposed. We call the proposed method covariance shrink (CS) method. BayesShrink and traditional multiwavelet vector thresholding approach will be the main comparison to the proposed method here, and as will be seen later in this paper, the proposed threshold often yields better results.
594
Fig. 2.
Fig. 1.
Flow chart of the proposed method.
III. T HE PROPOSED METHOD In this section, we are going to introduce the proposed system structure and the CS method. A. Proposed system structure We propose a new denoising method in this paper for medical X-ray images in the multiwavelet domain. Fig. 1 shows the proposed system structure. The algorithm effectively performs the process of “transform,”“denoise,” and “inverse transform”. The proposed algorithm can be characterized as consisting of five main components. 1) Prefilter the shifted image into multiple streams by a specific prefilter. 2) Transform the preprocessed image into the multiwavelet domain using an orthogonal periodic multiwavelet transform. 3) Calculate the thresholds by using the covariance shrink (CS) method and apply the bivariate soft thresholding to the wavelet coefficients (leave the scaling coefficients alone). 4) Postprocess the image. 5) Inverse thransform the result. B. The proposed covariance shrink method Firstly, we analysis the coefficients of multiwavelet transform. Fig. 2 shows the decomposition of the multiwavelet in a multiwavelet pyramid [18] obtained after L=2 steps, in the particular case of r = 2. We can conclude that if we use a different threshold for each subband, we can get better result than using the same threshold for all subbands because of the characteristic of coefficients in the multiwavelet domain. Secondly, we introduce the bivariate thresholding method. The success of wavelet thresholding lies in that, usually,
Multiwavelet image decomposition up to two levels, with r=2.
the signal will be compressed into a few large coefficients, whereas the noise components will be in the small coefficients only. Univariate thresholding is developed by Donoho and Johnstone, and two kinds of thresholding methods (soft and hard) are discussed. Univariate thresholding is successfully used in multiwavelet denoising by Strela et al. [10]. Even though univariate thresholding works in multiwavelet denoising, it does not give enough noise reduction. This is because multiwavelet transform normally produces correlated coefficients. Downie and Silverman developed the bivariate thresholding method by setting the universal threshold to λ = 2 log n [16]. Suppose we apply the DMWT with an appropriate prefilter to a noisy function, then we obtain ∗ M stream coefficients, Dj,m,n = Dj,m,n + Ej,m,n , where ∗ Dj,m,n denotes the signal coefficient, and Ej,m,n indicates a multivariate normal distribution N (0, Vj ). The matrix Vj is the covariance matrix for the error term that depends on the resolution scale j. Suppose we use a threshold of λ, then, the hard thresholding rule in bivariate thresholding can be written as ½ Dj,m,n if θj,m,n ≥ λ ˆ j,m,n = D (1) 0 otherwise q T where θj,m,n = Dj,m,n VjT Dj,m,n . The bivariate sof t thresholding can be formulated as ´ ( ³ λ if θj,m,n ≥ λ Dj,m,n 1 − θj,m,n ˆ (2) Dj,m,n = 0 otherwise To compute θj,k , we need Vj , which can be obtained from Vj = V ar(Dj,m,n ), and we can explicitly obtain the covariance structure of the wavelet coefficients at each level. Here, we can use a robust covariance estimation method to estimate Vj directly from the observed coefficients [17]. Finally, we introduce the proposed Covariance Shrink method. The thresholds of the CS method depend on the amount of each subband coefficients, coefficient’s values of
595
TABLE I
each subband, multiwavelet base, and preprocessor. We calculate the thresholds by using the following equations. p λj,h = βj,h 2cσ 2 (m × n) (3a) λj,v = βj,v
p 2cσ 2 (m × n)
Vjmean
m n 1 XX | Hj,m,n | m × n p=1 q=1
m n 1 XX = | Vj,m,n | m × n p=1 q=1
Djmean =
m n 1 XX | Dj,m,n | m × n p=1 q=1
(6b)
MWV 22.10 24.40 23.61
Proposed 22.89 24.72 25.66
VALUES FOR
Dose 4.17691E-08c/kg 7.43034E-08c/kg 1.38678E-07c/kg 2.98914E-07c/k
Noise 9.66 10.68 11.60 12.15
S PINE 1 X- RAY IMAGE BS 10.33 10.95 11.78 12.27
MWV 10.35 10.98 11.80 12.31
Proposed 11.86 11.97 12.43 12.65
a method of image quality measurement. PSNR we used in medical X-ray images is defined as: P SN R = 20log10 (
1 ) (11) M SE
1 X M SE = (d(x, y)/max(D) − o(x, y)/max(O))2 (12) H × L x,y (d(x, y) ∈ D, o(x, y) ∈ O)
(6c)
Hjmean + Vjmean + Djmean (7) 3 3) Calculate βj,h , βj,v , and βj,d . when σ ≤ 10, calculate β by using the following equations.
βj,h
BS 22.52 24.47 25.62
H,L
Averagej =
Averagej Averagej Djmean , βj,v = , βj,d = Hjmean Vjmean Averagej otherwise, when σ > 10
Noise 20.72 22.81 23.41
S PINE 1 X- RAY IMAGE
TABLE II MSR
(6a)
2) Calculate the average value of coefficient’s magnitudes of each level from 1)
βj,h =
VALUES FOR
Dose 4.17691E-08c/kg 7.43034E-08c/kg 1.38678E-07c/kg
(5b)
p λj,d = βj,d 2cσ 2 (m × n) (5c) where j is the decomposition level; h, v, d are the horizontal direction,the vertical direction, and the diagonal direction, respectively; c is the average variance, which depends on the selection of the multiwavelet base and preprocessor [7]; σ 2 is the variance of Gaussian noise; m × n is the size of the subband of interest. To calculate the thresholds, we need to calculate the following equations. 1) Calculate the average value of coefficient magnitudes of each subband Hjmean =
PSNR
(8)
where the high-dose X-ray image is D (d(x, y) ∈ D) that contains H by L pixels and the reconstructed image is O (o(x, y) ∈ O). We evaluate the proposed method by comparing to the following methods. • The traditional multiwavelet method (MWV). The multiwavelet base is CL, the preprocessing is row repeated preprocessing, the thresholding is scalar hard thresholding [7] and the threshold is calculated by the following equation. p T = 2 × cσ 2 log M × N (13) – c is the average variance.
Averagej Averagej Averagej = , βj,v = , βj,d = , (9) Hjmean Vjmean Djmean
– σ 2 is the variance of Gaussian noise.
IV. E XPERIMENTAL R ESULTS In this section, we show the experimental results and analyze the performance of the proposed algorithm. We test the proposed method on different medical X-ray images with different doses. To evaluate the medical X-ray image quality, we compute the mean-to-standard-deviation ratio (MSR) [13] in the desired region of interest (DROI). DROI is the region that the doctor is most interested in. µd (10) M SR = σd where µd and σd are the mean and the standard deviation computed in the DROI. In order to evaluate the efficiency of the proposed method more precisely, we also compute the Peak Signal to Noise Ratio (PSNR) which itself remains largely unsurpassed as
– M × N is the amount of coefficients in the multiwavelet transform. BayesShrink in wavlet domain (BS). The wavelet base is db5(Daubechies wavelet with window is 5), the thresholding is univariate soft thresholding. The PSNRs are given in Tables I. The MSRs are given in Tables II. Output results are given in Figs. 3∼??. From the previously method tables, we can see that the proposed method gives better values for PSNR, MSR, and CNR than the tranditional methods. From Figs. 3∼??, we can see that we can get the best denoised image by using the proposed method. From the figures and tables shown above, we know that the presented algorithm achieves the highest performance.
596
•
can obtain better result than scalar wavelet when noise the level is high and can preserve edges better than scalar wavelet. All experimental results show that the proposed method gives better results than conventional methods. Our future research will focus on finding more accommodating convergent values for wavelet coefficients. We also want to generalize it to all wavelets and employ it for other kinds of noise. R EFERENCES Noisy PSNR=20.72dB
MWV PSNR=22.10dB
BS PSNR=22.52dB
Proposed PSNR=22.89dB
Fig. 3. The results using different methods applied to Spine1, for dose = 4.17691E-08c/kg.
Noisy
BS
MWV
Proposed
Fig. 4. The amplified images using different methods applied to Spine1, for dose = 4.17691E-08c/kg.
V. C ONCLUSION In this paper, a new method for removing noise from medical X-ray images using CS thresholding in multiwavelet domain is proposed. Experimental results show that multiwavlets generally outperform scalar wavelet in image denoising for six groups of noisy 2D test images, and the results are visually very impressive. Chui-Lina scaling functions and wavelets combined with repeated row preprocessing appears to be a good general method. Multiwavelet denoising gives better results than scalar wavelet. In addition, multiwavelet bivariate denoising works better than univariate denoising. Multiwavelet
[1] T. Okamoto, S. Furui, H. Ichii, S. Yoshino, J. Lu and T. Yahagi: “Noise reduction in digital radiography using wavelet packet based on noise characteristics,” Journal of Signal Processing, Vol. 8, No. 6, pp. 485494, 2004. [2] R.D.Nowak and R.G.Baraniuk: “Wavelet-Domain Filtering for Photon Imaging Systems,” IEEE Trans. Image Processing, Vol. 8, No. 5, pp. 666-678, 1999. ¨ [3] M. K. Ozkan, A. T. Erdern, M. I. Sezan and A. M. Telcalp: “Efficient multiframe wiener restoration of blurred and noise image sequencws,” IEEE Trans. Image Processing, Vol. 1, No. 4, pp. 453–476, 1992. [4] Z. Wang and D. Zhang:“ Progressive Switching Median Filter for the Removal of Impulse Noise from Highly Corrupted Image,” IEEE Trans. Circuits and Systems: Analog and Digital Signal Processing, Vol. 46, No. 1, pp. 78-80, 1999. [5] I. Atkinson, F. Kamalabadi, S. Mohan and D. L. Jones: “Wavelet-based 2-D multichannel signal estimation,” Proceedings of IEEE Int. Conf. Image Processing, 2003. [6] D. L. Donoho, I. M. Johnstone, G. Kerkyacharian and D. Picard: Wavelet Shrinkage: Asymptopia, Journal of the Royal Statistical Society, Ser. B, vol. 57, No. 2, pp. 301–369, 1995. [7] R. Coifman and D. L. Donoho: “Translation invariant de-noising,” Lecture Notes in Statistics, 103, pp. 125-150, 1995. [8] C. Chang, B. Yu, and M. Vetterli:“Adaptive wavelet thresholding for image denoising and compression,” IEEE Trans. Image Processing, Vol. 9, No. 9, pp. 1532–1546, 2000. [9] V.Strela and A.T.Walden, “Signal and Image Denoising via Wavelet Thresholding:Orthogonal Biorthogonal,Scalar and Multiple Wavelet Transforms,” Imperical college,Statistic Section,Technical Report,TR-9801 ,1998. [10] T.D.Bui and G.Chen, “Translation-Invariant Denoising Using Multiwavelets,” IEEE Trans.Signal Processing, vol.46, pp. 3414-3420,1998. [11] X.G.Xia, J.S.Geronimo, D.P.Hardin, and B.W.Suter,“Design of Prefilters for Discrete Multiwavelet Transforms,” IEEE Trans.Signal Processing, vol. 44, pp. 25-35, 1996. [12] V.Strela, P.N.Heller, G.Strang, P.Topiwala, and C.Heil, “The application of multiwavelet filter banks to image processing,” Tech. Rep., Mass. Inst. Technol., Cambirdge, 1995. [13] J. S. Geronimo, D. P. Hardin, and P. R. Massopust, “Fractal functions and wavelet expansions based on several scaling function,” J.Approx.Theory, vol. 78, no. 3, pp. 373-401,1994 [14] D.L.Donoho and I.M.Johnstone, “Ideal spatial adaptation via wavelet shrinkage,” Biometrika, vol. 81, pp. 425-455, 1994. J. Approx.Theory, vol. 78, pp. 373.401, 1994. [15] C. K. Chui and J.-L. Lian, “A Study of Orthonormal Multi-Wavelets,” Applied Numer. vol. 20, pp. 273-298, 1996. [16] D.L.Donoho and I.M.Johnstone, “Adapting to unknown smoothness via wavelet shrinkage,” Journal of the American Statistical Assoc., vol. 90, pp. 1200-1224,1995. [17] E.D.Kolaczyk and D.D.Dizon, “Nonparametric estimation of intensity maps using haar wavelets and poisson noise characteristics,” The Astrophysical Journal, vol. 534, pp.490-505, 2000. [18] M.cotronei, D.Lazzaro, L.B.Montefusco, and L.Puccio, “Image compression through embeded multiwavelet transform coding,” IEEE Trans.Image Processing, vol.9, pp.184-189,2000. [19] T.R.Dowine and B.W.Silverman, “The discrete multiple wavelet transform and thresholding methods,” IEEE Trans. Signal Processing, vol. 46, pp.2558-2561, 1998. [20] P.J.Huber, Robust Statistics. New York: Wiley, 1981.
597
Abstract— When a signal is embedded in an additive Gaussian noise, its estimation is often done by finding a wavelet basis that concentrates the signal energy in few coefficients and then thresholding the noisy coefficients. However, in many practical problems such as medical X-ray image, astronomical and lowlight images, the recorded data is not modeled by Gaussian noise but as the realization of a Poisson process. Multiwavelet is a new development to the body of wavelet theory. Multiwavelet simultaneously offers orthogonality, symmetry and short support which are not possible in scalar 2-channel wavelet systems. After reviewing this recently developed theory, a new theory and algorithm for denoising medical X-ray images using multiwavelet multiple resolution analysis (MRA) are presented and investigated in this paper. The proposed covariance shrink (CS) method is used to threshold wavelet coefficients. The form of thresholds is carefully formulated which is the key to more excellent results obtained in the extensive numerical simulations of medical image denoising compared to conventional methods.
I. I NTRODUCTION Noise in medical X-ray images [1] is primarily categorized into quantum mottle, which is related to the number of incident X-rays, and due to artificial noise which is due to the grid, etc. The effect of quantum mottle is manifested as an increase in the graininess of an X-ray image as the dose is reduced. Therefore, noise reduction is of great significance in medical X-ray images. The noise of X-ray images obeys a Poisson law and, hence, is highly dependent on the underlying light intensity pattern being imaged [2]. Over the past decade, wavelet thresholding has received a lot of attention from researchers in image denoising. Scalar orthogonal wavelets with a single mother-wavelet function have played an important role in denoising [4]. The denoising is done on the detail coefficients of wavelet transform, it has been shown that this algorithm offers the advantages of smoothness and adaptation. However, it is known that there is a limitation on the time-frequency localization of scalar wavelet functions. Multiwavelet has been developed recently by using translates and dilates of more than one mother wavelet function [9]. They are known to have several advantages over scalar wavelets. In addition, Strela et al.[9] claimed that multiwavelet hard thresholding offers better results than the tranditional scalar wavelet hard thresholding. In this paper, a new shrink theory and denoising algorithm for medical X-ray images using multiwavelet multiple resolu-
tion analysis (MRA) is presented and investigated. We propose a new shrinkage method—covariance shrink (CS) method. The CS method is subband adaptive—a separate threshold is computed for each detail subbad of multiwavelet transform. We propose near-optimal thresholds and a more suitable image denoising by extensive numerical simulations.We show in this paper that multiwavelet denoising gives better results than the conventional scalar wavelet denoising and conventional multiwavelet denoising. II. S CALAR WAVELET AND MULTIWAVELET DENOISING Scalar wavelet and multiwavelet denoising attempt to remove the noise present in the signal while preserving the signal characteristics, regardless of its frequency content. they involve three steps: 1) A linear forward wavelet transform. 2) Nonlinear shrinkage denoising. 3) A linear inverse wavelet transform. Experimental results show that small wavelet coefficients are dominated by noise, while coefficients with a large absolute value carry more signal information than noise. Replacing noisy coefficients (small coefficients below a certain threshold value) by zeros and inverse wavelet transformation may lead to a reconstruction that has less noise. This method is indeed effective and thresholding is a simple and efficient method for noise reduction. The most well-know thresholding methods include VisuShrink [14] and SureShrink [15]. These threshold choices enjoy asymptotic minimax optimalities over function spaces such as Besov spaces. BayesShrink for denoising determines the threshold for each subband assuming a Generalized Gaussian Distribution (GGD). The results obtained for image denoising using BayesShrink are much better than VisuShrink and SureShrink. In this paper, a new multiwavelet thresholding method based on the two types of multivariate thresholding approach is proposed. We call the proposed method covariance shrink (CS) method. BayesShrink and traditional multiwavelet vector thresholding approach will be the main comparison to the proposed method here, and as will be seen later in this paper, the proposed threshold often yields better results.
594
Fig. 2.
Fig. 1.
Flow chart of the proposed method.
III. T HE PROPOSED METHOD In this section, we are going to introduce the proposed system structure and the CS method. A. Proposed system structure We propose a new denoising method in this paper for medical X-ray images in the multiwavelet domain. Fig. 1 shows the proposed system structure. The algorithm effectively performs the process of “transform,”“denoise,” and “inverse transform”. The proposed algorithm can be characterized as consisting of five main components. 1) Prefilter the shifted image into multiple streams by a specific prefilter. 2) Transform the preprocessed image into the multiwavelet domain using an orthogonal periodic multiwavelet transform. 3) Calculate the thresholds by using the covariance shrink (CS) method and apply the bivariate soft thresholding to the wavelet coefficients (leave the scaling coefficients alone). 4) Postprocess the image. 5) Inverse thransform the result. B. The proposed covariance shrink method Firstly, we analysis the coefficients of multiwavelet transform. Fig. 2 shows the decomposition of the multiwavelet in a multiwavelet pyramid [18] obtained after L=2 steps, in the particular case of r = 2. We can conclude that if we use a different threshold for each subband, we can get better result than using the same threshold for all subbands because of the characteristic of coefficients in the multiwavelet domain. Secondly, we introduce the bivariate thresholding method. The success of wavelet thresholding lies in that, usually,
Multiwavelet image decomposition up to two levels, with r=2.
the signal will be compressed into a few large coefficients, whereas the noise components will be in the small coefficients only. Univariate thresholding is developed by Donoho and Johnstone, and two kinds of thresholding methods (soft and hard) are discussed. Univariate thresholding is successfully used in multiwavelet denoising by Strela et al. [10]. Even though univariate thresholding works in multiwavelet denoising, it does not give enough noise reduction. This is because multiwavelet transform normally produces correlated coefficients. Downie and Silverman developed the bivariate thresholding method by setting the universal threshold to λ = 2 log n [16]. Suppose we apply the DMWT with an appropriate prefilter to a noisy function, then we obtain ∗ M stream coefficients, Dj,m,n = Dj,m,n + Ej,m,n , where ∗ Dj,m,n denotes the signal coefficient, and Ej,m,n indicates a multivariate normal distribution N (0, Vj ). The matrix Vj is the covariance matrix for the error term that depends on the resolution scale j. Suppose we use a threshold of λ, then, the hard thresholding rule in bivariate thresholding can be written as ½ Dj,m,n if θj,m,n ≥ λ ˆ j,m,n = D (1) 0 otherwise q T where θj,m,n = Dj,m,n VjT Dj,m,n . The bivariate sof t thresholding can be formulated as ´ ( ³ λ if θj,m,n ≥ λ Dj,m,n 1 − θj,m,n ˆ (2) Dj,m,n = 0 otherwise To compute θj,k , we need Vj , which can be obtained from Vj = V ar(Dj,m,n ), and we can explicitly obtain the covariance structure of the wavelet coefficients at each level. Here, we can use a robust covariance estimation method to estimate Vj directly from the observed coefficients [17]. Finally, we introduce the proposed Covariance Shrink method. The thresholds of the CS method depend on the amount of each subband coefficients, coefficient’s values of
595
TABLE I
each subband, multiwavelet base, and preprocessor. We calculate the thresholds by using the following equations. p λj,h = βj,h 2cσ 2 (m × n) (3a) λj,v = βj,v
p 2cσ 2 (m × n)
Vjmean
m n 1 XX | Hj,m,n | m × n p=1 q=1
m n 1 XX = | Vj,m,n | m × n p=1 q=1
Djmean =
m n 1 XX | Dj,m,n | m × n p=1 q=1
(6b)
MWV 22.10 24.40 23.61
Proposed 22.89 24.72 25.66
VALUES FOR
Dose 4.17691E-08c/kg 7.43034E-08c/kg 1.38678E-07c/kg 2.98914E-07c/k
Noise 9.66 10.68 11.60 12.15
S PINE 1 X- RAY IMAGE BS 10.33 10.95 11.78 12.27
MWV 10.35 10.98 11.80 12.31
Proposed 11.86 11.97 12.43 12.65
a method of image quality measurement. PSNR we used in medical X-ray images is defined as: P SN R = 20log10 (
1 ) (11) M SE
1 X M SE = (d(x, y)/max(D) − o(x, y)/max(O))2 (12) H × L x,y (d(x, y) ∈ D, o(x, y) ∈ O)
(6c)
Hjmean + Vjmean + Djmean (7) 3 3) Calculate βj,h , βj,v , and βj,d . when σ ≤ 10, calculate β by using the following equations.
βj,h
BS 22.52 24.47 25.62
H,L
Averagej =
Averagej Averagej Djmean , βj,v = , βj,d = Hjmean Vjmean Averagej otherwise, when σ > 10
Noise 20.72 22.81 23.41
S PINE 1 X- RAY IMAGE
TABLE II MSR
(6a)
2) Calculate the average value of coefficient’s magnitudes of each level from 1)
βj,h =
VALUES FOR
Dose 4.17691E-08c/kg 7.43034E-08c/kg 1.38678E-07c/kg
(5b)
p λj,d = βj,d 2cσ 2 (m × n) (5c) where j is the decomposition level; h, v, d are the horizontal direction,the vertical direction, and the diagonal direction, respectively; c is the average variance, which depends on the selection of the multiwavelet base and preprocessor [7]; σ 2 is the variance of Gaussian noise; m × n is the size of the subband of interest. To calculate the thresholds, we need to calculate the following equations. 1) Calculate the average value of coefficient magnitudes of each subband Hjmean =
PSNR
(8)
where the high-dose X-ray image is D (d(x, y) ∈ D) that contains H by L pixels and the reconstructed image is O (o(x, y) ∈ O). We evaluate the proposed method by comparing to the following methods. • The traditional multiwavelet method (MWV). The multiwavelet base is CL, the preprocessing is row repeated preprocessing, the thresholding is scalar hard thresholding [7] and the threshold is calculated by the following equation. p T = 2 × cσ 2 log M × N (13) – c is the average variance.
Averagej Averagej Averagej = , βj,v = , βj,d = , (9) Hjmean Vjmean Djmean
– σ 2 is the variance of Gaussian noise.
IV. E XPERIMENTAL R ESULTS In this section, we show the experimental results and analyze the performance of the proposed algorithm. We test the proposed method on different medical X-ray images with different doses. To evaluate the medical X-ray image quality, we compute the mean-to-standard-deviation ratio (MSR) [13] in the desired region of interest (DROI). DROI is the region that the doctor is most interested in. µd (10) M SR = σd where µd and σd are the mean and the standard deviation computed in the DROI. In order to evaluate the efficiency of the proposed method more precisely, we also compute the Peak Signal to Noise Ratio (PSNR) which itself remains largely unsurpassed as
– M × N is the amount of coefficients in the multiwavelet transform. BayesShrink in wavlet domain (BS). The wavelet base is db5(Daubechies wavelet with window is 5), the thresholding is univariate soft thresholding. The PSNRs are given in Tables I. The MSRs are given in Tables II. Output results are given in Figs. 3∼??. From the previously method tables, we can see that the proposed method gives better values for PSNR, MSR, and CNR than the tranditional methods. From Figs. 3∼??, we can see that we can get the best denoised image by using the proposed method. From the figures and tables shown above, we know that the presented algorithm achieves the highest performance.
596
•
can obtain better result than scalar wavelet when noise the level is high and can preserve edges better than scalar wavelet. All experimental results show that the proposed method gives better results than conventional methods. Our future research will focus on finding more accommodating convergent values for wavelet coefficients. We also want to generalize it to all wavelets and employ it for other kinds of noise. R EFERENCES Noisy PSNR=20.72dB
MWV PSNR=22.10dB
BS PSNR=22.52dB
Proposed PSNR=22.89dB
Fig. 3. The results using different methods applied to Spine1, for dose = 4.17691E-08c/kg.
Noisy
BS
MWV
Proposed
Fig. 4. The amplified images using different methods applied to Spine1, for dose = 4.17691E-08c/kg.
V. C ONCLUSION In this paper, a new method for removing noise from medical X-ray images using CS thresholding in multiwavelet domain is proposed. Experimental results show that multiwavlets generally outperform scalar wavelet in image denoising for six groups of noisy 2D test images, and the results are visually very impressive. Chui-Lina scaling functions and wavelets combined with repeated row preprocessing appears to be a good general method. Multiwavelet denoising gives better results than scalar wavelet. In addition, multiwavelet bivariate denoising works better than univariate denoising. Multiwavelet
[1] T. Okamoto, S. Furui, H. Ichii, S. Yoshino, J. Lu and T. Yahagi: “Noise reduction in digital radiography using wavelet packet based on noise characteristics,” Journal of Signal Processing, Vol. 8, No. 6, pp. 485494, 2004. [2] R.D.Nowak and R.G.Baraniuk: “Wavelet-Domain Filtering for Photon Imaging Systems,” IEEE Trans. Image Processing, Vol. 8, No. 5, pp. 666-678, 1999. ¨ [3] M. K. Ozkan, A. T. Erdern, M. I. Sezan and A. M. Telcalp: “Efficient multiframe wiener restoration of blurred and noise image sequencws,” IEEE Trans. Image Processing, Vol. 1, No. 4, pp. 453–476, 1992. [4] Z. Wang and D. Zhang:“ Progressive Switching Median Filter for the Removal of Impulse Noise from Highly Corrupted Image,” IEEE Trans. Circuits and Systems: Analog and Digital Signal Processing, Vol. 46, No. 1, pp. 78-80, 1999. [5] I. Atkinson, F. Kamalabadi, S. Mohan and D. L. Jones: “Wavelet-based 2-D multichannel signal estimation,” Proceedings of IEEE Int. Conf. Image Processing, 2003. [6] D. L. Donoho, I. M. Johnstone, G. Kerkyacharian and D. Picard: Wavelet Shrinkage: Asymptopia, Journal of the Royal Statistical Society, Ser. B, vol. 57, No. 2, pp. 301–369, 1995. [7] R. Coifman and D. L. Donoho: “Translation invariant de-noising,” Lecture Notes in Statistics, 103, pp. 125-150, 1995. [8] C. Chang, B. Yu, and M. Vetterli:“Adaptive wavelet thresholding for image denoising and compression,” IEEE Trans. Image Processing, Vol. 9, No. 9, pp. 1532–1546, 2000. [9] V.Strela and A.T.Walden, “Signal and Image Denoising via Wavelet Thresholding:Orthogonal Biorthogonal,Scalar and Multiple Wavelet Transforms,” Imperical college,Statistic Section,Technical Report,TR-9801 ,1998. [10] T.D.Bui and G.Chen, “Translation-Invariant Denoising Using Multiwavelets,” IEEE Trans.Signal Processing, vol.46, pp. 3414-3420,1998. [11] X.G.Xia, J.S.Geronimo, D.P.Hardin, and B.W.Suter,“Design of Prefilters for Discrete Multiwavelet Transforms,” IEEE Trans.Signal Processing, vol. 44, pp. 25-35, 1996. [12] V.Strela, P.N.Heller, G.Strang, P.Topiwala, and C.Heil, “The application of multiwavelet filter banks to image processing,” Tech. Rep., Mass. Inst. Technol., Cambirdge, 1995. [13] J. S. Geronimo, D. P. Hardin, and P. R. Massopust, “Fractal functions and wavelet expansions based on several scaling function,” J.Approx.Theory, vol. 78, no. 3, pp. 373-401,1994 [14] D.L.Donoho and I.M.Johnstone, “Ideal spatial adaptation via wavelet shrinkage,” Biometrika, vol. 81, pp. 425-455, 1994. J. Approx.Theory, vol. 78, pp. 373.401, 1994. [15] C. K. Chui and J.-L. Lian, “A Study of Orthonormal Multi-Wavelets,” Applied Numer. vol. 20, pp. 273-298, 1996. [16] D.L.Donoho and I.M.Johnstone, “Adapting to unknown smoothness via wavelet shrinkage,” Journal of the American Statistical Assoc., vol. 90, pp. 1200-1224,1995. [17] E.D.Kolaczyk and D.D.Dizon, “Nonparametric estimation of intensity maps using haar wavelets and poisson noise characteristics,” The Astrophysical Journal, vol. 534, pp.490-505, 2000. [18] M.cotronei, D.Lazzaro, L.B.Montefusco, and L.Puccio, “Image compression through embeded multiwavelet transform coding,” IEEE Trans.Image Processing, vol.9, pp.184-189,2000. [19] T.R.Dowine and B.W.Silverman, “The discrete multiple wavelet transform and thresholding methods,” IEEE Trans. Signal Processing, vol. 46, pp.2558-2561, 1998. [20] P.J.Huber, Robust Statistics. New York: Wiley, 1981.
597
