Wavelet Denoising of Multicomponent Images ... - Semantic Scholar
IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 16, NO. 7, JULY 2007
1865
Wavelet Denoising of Multicomponent Images Using Gaussian Scale Mixture Models and a Noise-Free Image as Priors Paul Scheunders and Steve De Backer
Abstract—In this paper, a Bayesian wavelet-based denoising procedure for multicomponent images is proposed. A denoising procedure is constructed that 1) fully accounts for the multicomponent image covariances, 2) makes use of Gaussian scale mixtures as prior models that approximate the marginal distributions of the wavelet coefficients well, and 3) makes use of a noise-free image as extra prior information. It is shown that such prior information is available with specific multicomponent image data of, e.g., remote sensing and biomedical imaging. Experiments are conducted in these two domains, in both simulated and real noisy conditions. Index Terms—Bayesian wavelet-based denoising, Gaussian scale mixture model (GSM), multimodal medical images, multispectral images.
I. INTRODUCTION ITH the evolution of imaging technology, an increasing number of imaging modalities becomes available. In remote sensing, sensors are available that can generate multispectral or hyperspectral data, involving high numbers of bands. In medical imagery, distinct image modalities reveal different features of the internal body. Examples are MRI images, acquired by using different imaging parameters (T1, ), different CT and nuclear T2, proton density, diffusion, medicine imaging modalities. In both cases, multicomponent images are generated. In multicomponent imagery, noise handling may be important. Due to physical constraints, a tradeoff exists between spatial and spectral resolution and SNR. Also, noise filtering and image enhancement can drastically facilitate the processing and analysis of multicomponent imagery. In particular, this holds for the classification and segmentation of multispectral images for identification purposes [1], [2] and for the segmentation of different tissue regions in medical imagery [3], [4]. Multicomponent image noise is usually treated as stochastic Gaussian distributed. Noise in the different bands is not necessarily independent. Decorrelating is traditionally performed using principal component analysis (PCA). The maximum noise fractions transform consists of a transformation in which the signal and noise covariance matrices are diagonalized [5]. Other
W
Manuscript received June 29, 2006; revised February 16, 2007. The associate editor coordinating the review of this manuscript and approving it for publication was Dr. Michael Elad. The authors are with the Vision Lab, Department of Physics, University of Antwerp, 2610 Wilrijk, Belgium (e-mail: [email protected]; steve. [email protected]). Digital Object Identifier 10.1109/TIP.2007.899598
techniques use wavelet transform for spatial and PCA for spectral decorrelation [6]. In this paper, a wavelet-based denoising strategy of multicomponent images is proposed. The paradigm of wavelet shrinking originates from [7]. There, a threshold value was derived using the discrete wavelet transform. The threshold value was universal, i.e., independent of the subband, and soft, i.e., all wavelet coefficients below the threshold were removed, and all others were shrunk by the threshold value. Later on, threshold values were derived that became adaptive, i.e., dependent of the specific subband. For this, Bayesian approaches were applied, where a prior model for the noise-free signal pdf was assumed. Popular priors are generalized Laplacian models (see, e.g., [8] and references therein) and Gaussian scale mixture (GSM) models [9], [10]. Recently, several wavelet-based procedures for multicomponent images were proposed [11]–[16]. Denoising was performed that accounted for the multicomponent image covariances, applying wavelet thresholding, mean squared-error estimation and Bayesian estimation, using different prior models: Gaussian models, GSM models, Laplacian models, and Bernouilli– Gaussian models. From these papers, it is clear that a superior denoising of multicomponent images is obtained when accounting for the full multicomponent image covariance structure, and when a good approximation of the wavelet coefficients marginals is used. In [17], a comparative study of these different techniques showed that all heavy tailed prior models outperformed the multinormal model. When compared with each other, they provided similar results and can, thus, be regarded as state of the art denoising techniques for multicomponent images. In this paper, a Bayesian estimator is constructed that makes use of such a heavy tailed multicomponent model for the noisefree signal. Moreover, other prior information is applied in the form of a “noise-free” single-component image (of course, no image is completely noise-free, with “noise-free” we mean “of high SNR”). In many practical situations involving multicomponent images, a single-component image of higher spatial resolution and/or SNR is available. In remote sensing, e.g., a higher-resolution sensor might be available. In the past, such auxiliary image was applied for fusion with a multispectral image to improve the latter’s spatial resolution while maintaining the spectral information. Several fusion procedures were worked out in the wavelet domain [18], [19]. An example is given by the fusion of Panchromatic data of a SPOT satellite system with Landsat Thematic Mapper multispectral images [20]–[24]. In hyperspectral images, the infrared (IR) part of the
1057-7149/$25.00 © 2007 IEEE
1866
IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 16, NO. 7, JULY 2007
spectrum contains noisy bands near the water vapor absorption band. To denoise these bands, image bands from other parts of the spectrum could be applied as noise-free images. In medical imaging, different images from the same or different modalities with various resolution and SNR can be acquired, e.g., in functional MR or diffusion tensor MR studies, noisy multicomponent images are acquired using fast acquisition procedures. During the same examination, anatomical information is obtained, that can be applied as an extra diagnostic tool. Such auxiliary image information is regularly applied for visualization or registration purposes. In [25], a technique is proposed to enhance the resolution of hyperspectral remote sensing images, with the aid of a high-resolution auxiliary sensor. The goal of our approach is denoising, rather than resolution enhancement. For this, the technique is extended to wavelet domain. In this way, we obtain a multiresolution Bayesian framework that accounts for the correlation between a multicomponent image and an auxiliary noise-free image, in order to improve the SNR of the first. In order to model the wavelet coefficients, multicomponent models should be applied. Since the difference between the three proposed multicomponent heavy-tailed priors in [17] were minor, we conjecture that the choice of the prior in this paper is not too critical as long as it is a heavy tailed multicomponent prior. Since the Laplacian model is not rigorously extendable to the multicomponent case, and the multicomponent Bernouilli–Gaussian model was originally introduced by other authors, in this paper we choose to apply the GSM model. Since the proposed procedure makes use of correlations, it is important that the multicomponent image and the auxiliary image are geometrically coregistered. In some applications, this causes no real problem, since the same sensor is applied for the images. As an example, we will demonstrate the procedure on a hyperspectral remote sensing image to denoise the bands near the water vapor absorbtion band. In some applications the misalignment is small, such as, e.g., in intramodal medical acquisition systems; in others, it should be performed a priori. In order to study the misregistration effects on the propoesed procedure, we include a simulation experiment with misaligned images. In the next section, the multiresolution Bayesian denoising framework is elaborated, and in Section III, the implementation is discussed. Experiments and discussion follow in Section IV. II. DENOISING FRAMEWORK A. Nondecimated Wavelet Transform The wavelet transform reorganizes image content into a low-resolution approximation and a set of details of different orientations and different resolution scales. A fast algorithm for the discrete wavelet transform is an iterative filter bank algorithm of Mallat [26], where a pair of high-pass and low-pass filters followed by downsampling by two is iterated on the low-pass output. The outputs of the low-pass filter are the so-called scaling coefficients and the outputs of the high-pass filter are the wavelet coefficients. At each decomposition level, the filter bank is applied sequentially to the rows and to the columns of the image. Low-pass filtering of both the rows and the columns yields the low-pass LL subband (i.e., approximation subband consisting
of the scaling coefficients). Other combinations of low-pass and high-pass filtering yield the wavelet subbands at different orientations: high-pass filtering of rows and low-pass filtering of columns (HL) yields horizontal edges and the opposite combination (LH) yields vertical edges, while high-pass filtering of both the rows and the columns (HH) yields diagonal edges. The th decomposition level yields the coefficients at the resolution scale . A full signal representation consists of the scaling coefficients at the resolution level and of all the wavelet coefficients at the resolution levels 1 to . The total number of these coefficients, due to downsampling by 2 at each stage, equals the number of the input image samples. This is a critically sampled transform with precisely enough coefficients for the perfect reconstruction. Despite their mathematical elegance and a remarkable signal/ image compression ability, critically sampled representations are less attractive for denoising. These representations lack shift invariance, meaning that the wavelet coefficients of a shifted signal differ from the shifted wavelet coefficients of the unshifted signal. In image denoising, better results are offered by redundant wavelet representations. In this paper, we use a nondecimated wavelet transform implemented with the algorithm zeroes between the àtrous [27]. The algorithm inserts filter coefficients at the resolution level . Downsampling of the filter outputs is excluded, so the size of each wavelet subband equals the size of the input image. B. Wavelet Processing of Multicomponent Data A natural way of exploiting the multicomponent correlations is by vector-based processing, operating on all the components denote the noise-free wavelet coefsimultaneously. Let ficient at spatial position , resolution level , orientation subband , and image component . , and are the corresponding wavelet coefficients of the observed noisy image and noise, respectively. A vector processing approach groups of all the components at a the wavelet coefficients given spatial position, within a subband of a given orientation and resolution level into a -dimensional vector (1) Equivalent processing is typically applied to all the wavelet subbands, and, hence, we shall omit the indices that denote the resolution level and orientation . In each wavelet subband, a multicomponent pixel obeys the additive noise model (2) where the probability density function (i.e., density) of is a multivariate Gaussian of zero mean and covariance matrix . is defined as the zero mean muldenotes tivariate Gaussian distribution with covariance . the estimate the covariance matrix of the noisy vector , and of the covariance matrix of the unknown noise-free vector . The noise covariance in each wavelet subband is, in general, a scaled version of the input image noise covariance, where the scaling factors depend on the wavelet filter coefficients (see, e.g., [28]). With the orthogonal wavelet families [27] that we use in this paper, the noise covariance in all the wavelet subbands is
SCHEUNDERS AND DE BACKER: WAVELET DENOISING OF MULTICOMPONENT IMAGES
equal to the input image noise covariance. Let denote the noise covariance in the input image. When conducting experiments with real noisy data, the noise covariance will be estimated separately. The signal covariance matrix is estimated as (3) Since is a covariance matrix, it needs to be positive semi-positive definite. This is assured by performing an eigenvalue decomposition and clipping possible negative eigenvalues to zero. In the rare event of negative eigenvalues occurring, they were observed to be negligible. C. Estimation Approach and Optimization Criterion Various linear and nonlinear (adaptive) methods can be applied in data denoising. We focus on the Bayesian approach, where a priori knowledge about the distribution of the noisefree data is assumed. In particular, we impose a multicomponent prior distribution (to be called hereafter prior) on the noise-free wavelet coefficients in a given subband. As an optimization criterion, we adopt minimization of the mean squared error, i.e., the Bayesian risk is a quadratic loss function. The minimum mean-squared error (MMSE) estimate is the posterior conditional mean
density for [29], [30]
1867
is required. A widely used prior is Jeffrey’s prior (8)
E. Noise-Free Image A single-component image of the same scene as captured by the multicomponent image is assumed to be available. In this paper, we want to employ the correlation between the multicomponent image and this single-component image to improve the denoising of the first. For this, both images need to be geometrically registered and the images should be correlated. This correlation will be accounted for in a Bayesian denoising procedure, by treating the single-component image as extra prior information. Denote the wavelet detail coefficients of a noise-free singlecomponent image. We assume that and are jointly GSM as the vector, concatenating and distributed. If we write , we can rewrite (6) as (9) where becomes
is zero mean Gaussian distributed; (7) then (10)
where
is a zero-mean Gaussian with covariance , or by taking the expectations over , with . The conditional pdf of given can be written as a mixture of conditional Gaussian distributions (4) (11) Assuming, e.g., a Gaussian prior for the noise-free signal , the above MMSE estimate is the Wiener filter
where
is the Gaussian
with [31]
(5) D. GSM Model
(12)
A GSM model models the noise-free wavelet coefficients as (6) is an indewhere is zero-mean Gaussian distributed and is pendent positive scalar random variable. The probability then a mixture of Gaussians:
Notice that is not dependent on , and we can, therefore, . Also, its value is generally not equal to zero, denote it as making (11) is not strictly a GSM, but a shifted version. F. Bayesian Estimation The Bayes least squares estimate
is given by
(7) is the mixing density, and is a zero mean where or, by taking expectaGaussian with covariance . GSM densities are tions over , with symmetric and zero-mean and heavier tailed than Gaussians. These are known to better model the shape of the wavelet coefficient marginals than Gaussians. In [10], GSMs were applied to model local spatial neighborhoods of wavelet coefficients to denoise greylevel images. In this paper, we apply them to model the pixels of a multicomponent image, where they account for the strong correlation between the different components. A prior
(13) The first term in (13) can be obtained, using Bayes’ rule (14)
1868
with
IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 16, NO. 7, JULY 2007
given by (8) and
leading to (15)
We will now calculate the term rule, one obtains
(22)
. Using Bayes Now define: . Then
and the vectors:
, and
(16) (23) Since conditioned on is given by and since we have assumed that is independent of both , and , it follows that and , conditioned on are independent
, and are independent of , they need to be calSince culated only once. Remark that also for the calculation of the posterior distribution of in (15), their definition can be applied (24)
(17) The second expression appears when using Bayes rule on . Finally, canceling and skipping and because they do not depend on , the last expression is obtained. comes from the observation model The pdf (18) while the second term is given by (11). Then
where , and are the components of and the diagonal elements of respectively. When the covariances are estimated globally over the entire image, is determined globally. Then, , and are calculated for each pixel, after which the integral over is performed. For the calculation of the integral over , we follow [10], where is sampled logarithmically, which requires fewer samples than linear sampling. Optimal results were obtained for an (only interval from 20.5 until 3.5, with a step size of 2 for 13 samples). IV. EXPERIMENTS AND DISCUSSION
(19)
III. IMPLEMENTATION The terms (14) and (19) can be calculated by estimating (12) from the observed image and the noise-free image . In this paper, all estimations are done globally, over the entire image. The expression (19) is dependent on , and should be recalculated for each in the integral of (13). In order to simplify the dependence on of this expression, let us define: , i.e., . Then, (19) can be rewritten as
(20) Let be the square root of the matrix (i.e., ), and let denote the eigenvector/eigenvalue decomposition of . Then the matrix
(21)
To validate and compare the proposed techniques, we perform experiments on different data sets. The data sets contain a multicomponent image and a noise-free single component image of the same scene, that is geometrically registered with the multicomponent image. In the first experiment, MRI data is used. A typical MR imaging system is capable of acquiring image data at different acquisition modes. This leads to series of images with complementary information. Also, modern, fast acquisition procedures allow for performing time-critical experiments, as, e.g., functional MRI and diffusion tensor MRI, where multicomponent images are formed. Of course, spatial resolution and SNR are limited in these cases. During the complete acquisition time of such series of images, a high-quality image may be formed as well. Such image can be applied as a noise-free image in the proposed procedure. In the other two experiments, remote sensing data was used. Multispectral or hyperspectral sensors acquire multicomponent images. A higher quality image can be obtained from, e.g., another sensor, or another part of the reflectance spectrum with higher SNR. We will apply such an image as noise-free image in the proposed procedure. The proposed technique has some advantages over standard denoising of each image component. It accounts for the full covariance structure of the multicomponent image, it makes use of a prior that approximates well the marginal distributions of the wavelet coefficients, and it makes use of a noise-free image as extra prior information. In order to demonstrate the impact of
SCHEUNDERS AND DE BACKER: WAVELET DENOISING OF MULTICOMPONENT IMAGES
Fig. 1.
T
1869
2; P D; and T 1 noisy components of the MRI image and noise-free T 1 image.
these three advantages on the improvement of the performance, the following denoising strategies are applied to compare. • Single component denoising, where each band image is treated independently. A Gaussian prior model is assumed. This is the Wiener filter in wavelet domain. • Multicomponent denoising, where the correlations between the different components are exploited. A Gaussian prior is assumed. • Multicomponent denoising, assuming the GSM model as a prior to better model the marginal distributions of the wavelet coefficients. • Multicomponent denoising, using a noise-free single component image, using a Gaussian prior. • Multicomponent denoising, using a noise-free single component image, using a GSM prior. If the noise is simulated, the performance of the different denoising techniques can be quantitatively described by the PSNR (in dB), defined for 8-bit images as
where is the mean squared error between the original noiseless image and the denoised image . In case of multicomponent images, all bands are separately histogram-stretched to 8 bit, and the mse is averaged over all bands to obtain the PSNR. For the wavelet decomposition, a four-level decomposition of the nondecimated wavelet transform, as described in Section II-A, is used. For the wavelet filters, we use the Daubechies symmlets family [32].
TABLE I PSNR (IN DECIBELS) VALUES AFTER DENOISING OF MRI DATA SET
the image, by using the GSM prior and by using the noise-free image as extra prior information. When including the auxiliary image, sharper edges are obtained, as visible in Fig. 2(d) versus Fig. 2(e) and Fig. 2(f) versus Fig. 2(g). From Fig. 2(f) and (g), one can clearly observe that the noise is better removed using the GSM model. To asses the dependency of the results on the registration accuracy between the auxiliary image and the multicomponent image, we shifted the auxiliary image and repeated the experdB. Fig. 3 shows iment at initial noise level the PSNR after denoising in function of the shift. The dashed line show the results without the use of auxiliary image. After a steep drop in performance at subpixel misregistration errors, the quality slowly reduces until a shift of three pixels, where the misregistration error removes the added value of incorporating the auxiliary image. This demonstrates that for optimal results an almost perfect registration is required, but it is still worthwhile to use the auxiliary image even without perfect registration.
A. Experiment 1: MRI Data
B. Experiment 2: Multispectral Remote Sensing Data
The MRI data set consists of three different simulated MRI , and PD). The data set acquisitions of the same subject ( was obtained from the the Brainweb site [33]. On these images Rician distributed noise of different amplitude is simulated. The noise was independent for each component. A fourth image was , and kept noise-free. The images are shown in generated Fig. 1. Table I shows the obtained PSNR values. Fig. 2 shows a decomponent, for the different denoising tailed window of the procedures. In both Table I and Fig. 1, one can observe a gradual improvement of the results, by accounting for the covariances of
In the second experiment, a seven-band Landsat TM image was taken over the Winnipeg area containing both rural as urban areas. From this image, the thermal band (band 6) was removed and five spectral bands were corrupted with additive Gaussian noise. The noise was chosen independently for each , with varying standard deviation. The red band: band (band 3) served as the noise free image. Table II shows the results of the denoising procedures for this image: a) original Landsat image, b) detail of Landsat image, c) noisy Landsat image, d) single-component denoising with Gaussian prior, e) multicomponent denoising with Gaussian prior, f) multicomponent denoising with GSM prior, g) multicomponent denoising
1870
IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 16, NO. 7, JULY 2007
Fig. 2. (a) Detail of original T 2 MRI image of Fig. 1; (b) T 2 image with simulated Rician noise (PSNR = 34:13); (c) single-component denoising with Gaussian prior (PSNR = 36:21); (d) multicomponent denoising with Gaussian prior (PSNR = 37:49); (e) multicomponent denoising with Gaussian prior and noise-free image (PSNR = 38:71); (f) multicomponent denoising with GSM prior (PSNR = 39:54); (g) multicomponent denoising with GSM prior and noise-free image (PSNR = 40:92).
(Fig. 4). Here, the performance increase due to better description of the prior distribution is not as large as in the MRI case.
Fig. 3. Denoising quality in function of misregistration of auxiliary image. The dashed line indicates the quality achieved with multicomponent denoising without auxiliary image.
TABLE II PSNR (IN DB) VALUES AFTER DENOISING OF MULTISPECTRAL DATA SET
C. Experiment 3: Hyperspectral Remote Sensing Data In the third experiment, an AVIRIS image over Cuprite, Nevada, was taken. Here, we concentrated on the infra-red (IR) part of the spectrum, containing noisy bands near the water vapor absorption band. For this experiment, five bands were uniformly selected in the 2.3–2.5- m region. No additional noise was added. Fig. 5(a) shows one band of the image near 2.5 m. The noise covariance matrix is estimated as follows: let denote the wavelet coefficients of the first resolution level and orientation subband from spectral band . Then the diagonal elements of the noise covariance are estimated by the classical median estimator [34] (25) The off-diagonal elements between spectral bands estimated using [12]
and
are (26)
with Gaussian prior and noise-free image, and h) multicomponent denoising with GSM prior and noise-free image. Again, a gradual improvement can be observed, when accounting for the covariances of the image, using the GSM prior and using the noise-free image as extra prior information
where , and . For the auxiliary image, an image near 2 m was chosen [Fig. 5(b)]. In this region of the spectrum, less noise is present. A detailed window of the denoised images are shown in Fig. 5. Similar results as in the previous experiments are observed. The single component denoising procedure tends to blur the image excessively, removing most of the details. These details are better preserved in the multicomponent approaches.
SCHEUNDERS AND DE BACKER: WAVELET DENOISING OF MULTICOMPONENT IMAGES
1871
Fig. 4. (a) One band of Landsat multispectral image; (b) detail of image; (c) detail image with simulated Gaussian noise ( = 15); (d) single-component denoising with Gaussian prior; (e) multicomponent denoising with Gaussian prior; (f) multicomponent denoising with GSM prior; (g) multicomponent denoising with Gaussian prior and noise-free image; (h) multicomponent denoising with GSM prior and noise-free image.
Fig. 5. (a) Detail of original 2.479-m band of the AVIRIS cuprite image; (b) 1.991-m band of the cuprite image, used as the noise-free image; (c)–(g) results for the 2.479-m band after denoising; (c) single-component denoising with Gaussian prior; (d) multicomponent denoising with Gaussian prior; (e) multicomponent denoising with GSM prior; (f) multicomponent denoising with Gaussian prior and noise-free image; (g) multicomponent denoising with GSM prior and noise-free image.
When adding the noise-free image, the difference between Gaussian and GSM model becomes clear; with the Gaussian prior, the noise is not completely removed, while when using the GSM model, the noise is effectively removed and sharp edges are clearly retained.
V. CONCLUSION In this paper, a Bayesian wavelet-based denoising procedure for multicomponent images was proposed. The procedure makes use of a noise-free single component image as prior
1872
information. The prior model for the wavelet coefficient marginals is a GSM model. Experiments were performed with simulated noise on MRI imagery, multispectral and hyperspectral remote sensing images. The results show a gradual improvement, when applying a technique that 1) fully accounts for the multicomponent image covariances, 2) makes use of Gaussian Scale Mixtures as prior models that approximate well the marginal distributions of the wavelet coefficients, and 3) makes use of a noise-free image as extra prior information. REFERENCES [1] I. L. Thomas, V. M. Benning, and N. P. Ching, Classification of Remotely Sensed Images. Bristol, U.K.: Hilger, 1987. [2] C. Lee and D. A. Landgrebe, “Analyzing high-dimensional multispectral data,” IEEE Trans. Geosci. Remote Sens., vol. 31, no. 4, pp. 388–400, Apr. 1993. [3] T. Taxt and A. Lundervold, “Multispectral analysis of the brain using magnetic resonance imaging,” IEEE Trans. Med. Imag., vol. 13, no. 3, pp. 470–481, Mar. 1994. [4] C. Busch, “Wavelet based texture segmentation of multi-modal tomographic images,” Comput. Graph., vol. 21, no. 3, pp. 347–358, 1997. [5] A. A. Green, M. Berman, P. Switzer, and M. D. Craig, “A transformation for ordering multispectral data in terms of image quality with implications for noise removal,” IEEE Trans. Geosci. Remote Sens., vol. 26, no. 1, pp. 65–74, Jan. 1988. [6] J. L. Starck and P. Querre, “Multispectral data restoration by the wavelet Karhunen-Loève transform,” Signal Process., vol. 81, pp. 2449–2459, 2001. [7] D. L. Donoho and I. M. Johnstone, “Ideal spatial adaptation via wavelet shrinkage,” Biometrika, vol. 81, pp. 425–455, 1994. [8] S. G. Chang, B. Yu, and M. Vetterli, “Adaptive wavelet thresholding for image denoising and compression,” IEEE Trans. Image Process., vol. 9, no. 9, pp. 1532–1546, Sep. 2000. [9] M. K. Mihgak, I. Kozintsev, K. Ramchandran, and P. Moulin, “Lowcomplexity image denoising based on statistical modeling of wavelet coefficients,” IEEE Signal Process. Lett., vol. 6, no. 12, pp. 300–303, Dec. 1999. [10] J. Portilla, V. Strela, M. J. Wainwright, and E. P. Simoncelli, “Image denoising using scale mixtures of Gaussians in the wavelet domain,” IEEE Trans. Image Process., vol. 12, no. 11, pp. 1338–1351, Nov. 2003. [11] A. Benazza-Benyahia and J. C. Pesquet, “An extended sure approach for multicomponent image denoising,” in Proc. IEEE Int. Conf. Acoustics, Speech, Signal Processing, 2004, vol. 2, p. ii-945-8. [12] A. Benazza-Benyahia and J.-C. Pesquet, “Building robust wavelet estimators for multicomponent images using Steins’ principle,” IEEE Trans. Image Process., vol. 14, no. 11, pp. 1814–1830, Nov. 2005. [13] A. Piˇzurica, P. Scheunders, and W. Philips, “Multiresolution multispectral image denoising based on probability of presence of features of interest,” in Proc. Advanced Concepts in Image and Vision Systems, Brussels, Belgium, 2004, pp. 357–364. [14] A. Piˇzurica, B. Huysmans, P. Scheunders, and W. Philips, “Wavelet domain denoising of multispectral remote sensing imagery adapted to the local spatial and spectral context,” presented at the IGGARS, 2005. [15] P. Scheunders, “Wavelet thresholding of multivalued images,” IEEE Trans. Image Process., vol. 13, no. 4, pp. 475–483, Apr. 2004. [16] P. Scheunders and J. Driesen, “Least-squares interband denoising of color and multispectral images,” in Proc. IEEE Int. Conf. Image Processing, 2004, pp. 985–988. [17] S. de Backer, A. Piˇzurica, B. Huysmans, W. Philips, and P. Scheunders, “Denoising of multispectral images using wavelet least-squares estimators,” Image Vis. Comput., 2006. [18] H. Li, B. S. Manjunath, and S. K. Mitra, “Multisensor image fusion using the wavelet transform,” Graph. Models Image Process., vol. 57, no. 3, pp. 235–245, 1995. [19] D. A. Yocky, “Image merging and data fusion by means of the discrete two-dimensional wavelet transform,” J. Opt. Soc. Amer. A, vol. 12, no. 9, pp. 1834–1841, 1995.
IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 16, NO. 7, JULY 2007
[20] D. Yocky, “Multiresolution wavelet decomposition image merger of landsat thematic mapper and spot panchromatic data,” Photogramm. Eng. Remote Sens., vol. 62, pp. 1067–1074, 1996. [21] J. Zhou, D. Civco, and J. Silander, “A wavelet transform method to merge landsat tm and spot panchromatic data,” Int. J. Remote Sens., vol. 19, no. 4, pp. 743–757, 1998. [22] J. Nunez, X. Otazu, O. Fors, A. Prades, V. Pala, and R. Arbiol, “Image fusion with additive multiresolution wavelet decomposition, applications to spot+landsat images,” J. Opt. Soc. Amer. A, vol. 16, pp. 467–474, 1999. [23] J. Nunez, X. Otazu, O. Fors, A. Prades, V. Pala, and R. Arbiol, “Multiresolution-based image fusion with additive wavelet decomposition,” IEEE Trans. Geosci. Remote Sens., vol. 37, no. 3, pp. 1204–1211, May 1999. [24] T. Ranchin and L. Wald, “Fusion of high spatial and spectral resolution images: The arsis concept and its implementation,” Photogramm. Eng. Remote Sens., vol. 66, pp. 49–61, 2000. [25] R. C. Hardie, M. T. Eismann, and G. L. Wilson, “Map estimation for hyperspectral image resolution enhancement using an auxiliary sensor,” IEEE Trans. Image Process., vol. 13, no. 9, pp. 1174–1184, Sep. 2004. [26] S. Mallat, “A theory for multiresolution signal decomposition: The wavelet representation,” IEEE Trans. Pattern Anal. Mach. Intell., vol. 11, no. 7, pp. 674–693, Jul. 1989. [27] S. Mallat, A Wavelet Tour of Signal Processing. London, U.K.: Academic Press, 1998. [28] M. Jansen, Noise Reduction by Wavelet Thresholding. New York: Springer-Verlag, 2001. [29] G. E. P. Bax and C. Tiao, Bayesian Inference in Statistical Analysis. Reading, MA: Addison-Wesley, 1992. [30] M. Figeiredo and R. Nowak, “Wavelet-based image estimation: An empirical bayes approach using Jeffrey’s noninformativeprior,” IEEE Trans. Image Process., vol. 10, no. 9, pp. 1322–1331, Sep. 2001. [31] S. M. Kay, Fundamentals of Statistical Signal Processing: Estimation Theory. Englewoods Cliffs, NJ: Prentice-Hall, 1993. [32] I. Daubechies, Ten Lectures on Wavelets. Philadelphia, PA: SIAM, 1992. [33] C. A. Cocosco, R. K.-S. Kwan, V. Kollokian, and A. C. Evans, Brainweb: Online Interface to a 3-D MRI Simulated Brain Database. [34] D. L. Donoho and I. M. Johnstone, “Adapting to unknown smoothness via wavelet shrinking,” J. Amer. Statist. Assoc., vol. 90, no. 432, pp. 1200–1224, 1995.
Paul Scheunders graduated in physics in 1983 and received the Ph.D. degree in physics, with work in the field of statistical mechanics, in 1990, from the University of Antwerp, Antwerp, Belgium. In 1991, he became a Research Associate at the Vision Lab, Department of Physics of the University of Antwerp, where he became a Lecturer. He has published over 100 papers in international journals and proceedings in the field of image processing and pattern recognition. His research interest includes wavelets and multispectral image processing.
Steve De Backer received the B.S. degree in physics from the University of Hasselt, Belgium, in 1994, and the M.S. and Ph.D. degrees physics from the University of Antwerp, Antwerp, Belgium, in 1996 and 2002, respectively. Since then, he has been a Postdoctorate Fellow at the Visionlab, University of Antwerp. His research focus is on image modeling, pattern recognition, and hyperspectral image analysis.
1865
Wavelet Denoising of Multicomponent Images Using Gaussian Scale Mixture Models and a Noise-Free Image as Priors Paul Scheunders and Steve De Backer
Abstract—In this paper, a Bayesian wavelet-based denoising procedure for multicomponent images is proposed. A denoising procedure is constructed that 1) fully accounts for the multicomponent image covariances, 2) makes use of Gaussian scale mixtures as prior models that approximate the marginal distributions of the wavelet coefficients well, and 3) makes use of a noise-free image as extra prior information. It is shown that such prior information is available with specific multicomponent image data of, e.g., remote sensing and biomedical imaging. Experiments are conducted in these two domains, in both simulated and real noisy conditions. Index Terms—Bayesian wavelet-based denoising, Gaussian scale mixture model (GSM), multimodal medical images, multispectral images.
I. INTRODUCTION ITH the evolution of imaging technology, an increasing number of imaging modalities becomes available. In remote sensing, sensors are available that can generate multispectral or hyperspectral data, involving high numbers of bands. In medical imagery, distinct image modalities reveal different features of the internal body. Examples are MRI images, acquired by using different imaging parameters (T1, ), different CT and nuclear T2, proton density, diffusion, medicine imaging modalities. In both cases, multicomponent images are generated. In multicomponent imagery, noise handling may be important. Due to physical constraints, a tradeoff exists between spatial and spectral resolution and SNR. Also, noise filtering and image enhancement can drastically facilitate the processing and analysis of multicomponent imagery. In particular, this holds for the classification and segmentation of multispectral images for identification purposes [1], [2] and for the segmentation of different tissue regions in medical imagery [3], [4]. Multicomponent image noise is usually treated as stochastic Gaussian distributed. Noise in the different bands is not necessarily independent. Decorrelating is traditionally performed using principal component analysis (PCA). The maximum noise fractions transform consists of a transformation in which the signal and noise covariance matrices are diagonalized [5]. Other
W
Manuscript received June 29, 2006; revised February 16, 2007. The associate editor coordinating the review of this manuscript and approving it for publication was Dr. Michael Elad. The authors are with the Vision Lab, Department of Physics, University of Antwerp, 2610 Wilrijk, Belgium (e-mail: [email protected]; steve. [email protected]). Digital Object Identifier 10.1109/TIP.2007.899598
techniques use wavelet transform for spatial and PCA for spectral decorrelation [6]. In this paper, a wavelet-based denoising strategy of multicomponent images is proposed. The paradigm of wavelet shrinking originates from [7]. There, a threshold value was derived using the discrete wavelet transform. The threshold value was universal, i.e., independent of the subband, and soft, i.e., all wavelet coefficients below the threshold were removed, and all others were shrunk by the threshold value. Later on, threshold values were derived that became adaptive, i.e., dependent of the specific subband. For this, Bayesian approaches were applied, where a prior model for the noise-free signal pdf was assumed. Popular priors are generalized Laplacian models (see, e.g., [8] and references therein) and Gaussian scale mixture (GSM) models [9], [10]. Recently, several wavelet-based procedures for multicomponent images were proposed [11]–[16]. Denoising was performed that accounted for the multicomponent image covariances, applying wavelet thresholding, mean squared-error estimation and Bayesian estimation, using different prior models: Gaussian models, GSM models, Laplacian models, and Bernouilli– Gaussian models. From these papers, it is clear that a superior denoising of multicomponent images is obtained when accounting for the full multicomponent image covariance structure, and when a good approximation of the wavelet coefficients marginals is used. In [17], a comparative study of these different techniques showed that all heavy tailed prior models outperformed the multinormal model. When compared with each other, they provided similar results and can, thus, be regarded as state of the art denoising techniques for multicomponent images. In this paper, a Bayesian estimator is constructed that makes use of such a heavy tailed multicomponent model for the noisefree signal. Moreover, other prior information is applied in the form of a “noise-free” single-component image (of course, no image is completely noise-free, with “noise-free” we mean “of high SNR”). In many practical situations involving multicomponent images, a single-component image of higher spatial resolution and/or SNR is available. In remote sensing, e.g., a higher-resolution sensor might be available. In the past, such auxiliary image was applied for fusion with a multispectral image to improve the latter’s spatial resolution while maintaining the spectral information. Several fusion procedures were worked out in the wavelet domain [18], [19]. An example is given by the fusion of Panchromatic data of a SPOT satellite system with Landsat Thematic Mapper multispectral images [20]–[24]. In hyperspectral images, the infrared (IR) part of the
1057-7149/$25.00 © 2007 IEEE
1866
IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 16, NO. 7, JULY 2007
spectrum contains noisy bands near the water vapor absorption band. To denoise these bands, image bands from other parts of the spectrum could be applied as noise-free images. In medical imaging, different images from the same or different modalities with various resolution and SNR can be acquired, e.g., in functional MR or diffusion tensor MR studies, noisy multicomponent images are acquired using fast acquisition procedures. During the same examination, anatomical information is obtained, that can be applied as an extra diagnostic tool. Such auxiliary image information is regularly applied for visualization or registration purposes. In [25], a technique is proposed to enhance the resolution of hyperspectral remote sensing images, with the aid of a high-resolution auxiliary sensor. The goal of our approach is denoising, rather than resolution enhancement. For this, the technique is extended to wavelet domain. In this way, we obtain a multiresolution Bayesian framework that accounts for the correlation between a multicomponent image and an auxiliary noise-free image, in order to improve the SNR of the first. In order to model the wavelet coefficients, multicomponent models should be applied. Since the difference between the three proposed multicomponent heavy-tailed priors in [17] were minor, we conjecture that the choice of the prior in this paper is not too critical as long as it is a heavy tailed multicomponent prior. Since the Laplacian model is not rigorously extendable to the multicomponent case, and the multicomponent Bernouilli–Gaussian model was originally introduced by other authors, in this paper we choose to apply the GSM model. Since the proposed procedure makes use of correlations, it is important that the multicomponent image and the auxiliary image are geometrically coregistered. In some applications, this causes no real problem, since the same sensor is applied for the images. As an example, we will demonstrate the procedure on a hyperspectral remote sensing image to denoise the bands near the water vapor absorbtion band. In some applications the misalignment is small, such as, e.g., in intramodal medical acquisition systems; in others, it should be performed a priori. In order to study the misregistration effects on the propoesed procedure, we include a simulation experiment with misaligned images. In the next section, the multiresolution Bayesian denoising framework is elaborated, and in Section III, the implementation is discussed. Experiments and discussion follow in Section IV. II. DENOISING FRAMEWORK A. Nondecimated Wavelet Transform The wavelet transform reorganizes image content into a low-resolution approximation and a set of details of different orientations and different resolution scales. A fast algorithm for the discrete wavelet transform is an iterative filter bank algorithm of Mallat [26], where a pair of high-pass and low-pass filters followed by downsampling by two is iterated on the low-pass output. The outputs of the low-pass filter are the so-called scaling coefficients and the outputs of the high-pass filter are the wavelet coefficients. At each decomposition level, the filter bank is applied sequentially to the rows and to the columns of the image. Low-pass filtering of both the rows and the columns yields the low-pass LL subband (i.e., approximation subband consisting
of the scaling coefficients). Other combinations of low-pass and high-pass filtering yield the wavelet subbands at different orientations: high-pass filtering of rows and low-pass filtering of columns (HL) yields horizontal edges and the opposite combination (LH) yields vertical edges, while high-pass filtering of both the rows and the columns (HH) yields diagonal edges. The th decomposition level yields the coefficients at the resolution scale . A full signal representation consists of the scaling coefficients at the resolution level and of all the wavelet coefficients at the resolution levels 1 to . The total number of these coefficients, due to downsampling by 2 at each stage, equals the number of the input image samples. This is a critically sampled transform with precisely enough coefficients for the perfect reconstruction. Despite their mathematical elegance and a remarkable signal/ image compression ability, critically sampled representations are less attractive for denoising. These representations lack shift invariance, meaning that the wavelet coefficients of a shifted signal differ from the shifted wavelet coefficients of the unshifted signal. In image denoising, better results are offered by redundant wavelet representations. In this paper, we use a nondecimated wavelet transform implemented with the algorithm zeroes between the àtrous [27]. The algorithm inserts filter coefficients at the resolution level . Downsampling of the filter outputs is excluded, so the size of each wavelet subband equals the size of the input image. B. Wavelet Processing of Multicomponent Data A natural way of exploiting the multicomponent correlations is by vector-based processing, operating on all the components denote the noise-free wavelet coefsimultaneously. Let ficient at spatial position , resolution level , orientation subband , and image component . , and are the corresponding wavelet coefficients of the observed noisy image and noise, respectively. A vector processing approach groups of all the components at a the wavelet coefficients given spatial position, within a subband of a given orientation and resolution level into a -dimensional vector (1) Equivalent processing is typically applied to all the wavelet subbands, and, hence, we shall omit the indices that denote the resolution level and orientation . In each wavelet subband, a multicomponent pixel obeys the additive noise model (2) where the probability density function (i.e., density) of is a multivariate Gaussian of zero mean and covariance matrix . is defined as the zero mean muldenotes tivariate Gaussian distribution with covariance . the estimate the covariance matrix of the noisy vector , and of the covariance matrix of the unknown noise-free vector . The noise covariance in each wavelet subband is, in general, a scaled version of the input image noise covariance, where the scaling factors depend on the wavelet filter coefficients (see, e.g., [28]). With the orthogonal wavelet families [27] that we use in this paper, the noise covariance in all the wavelet subbands is
SCHEUNDERS AND DE BACKER: WAVELET DENOISING OF MULTICOMPONENT IMAGES
equal to the input image noise covariance. Let denote the noise covariance in the input image. When conducting experiments with real noisy data, the noise covariance will be estimated separately. The signal covariance matrix is estimated as (3) Since is a covariance matrix, it needs to be positive semi-positive definite. This is assured by performing an eigenvalue decomposition and clipping possible negative eigenvalues to zero. In the rare event of negative eigenvalues occurring, they were observed to be negligible. C. Estimation Approach and Optimization Criterion Various linear and nonlinear (adaptive) methods can be applied in data denoising. We focus on the Bayesian approach, where a priori knowledge about the distribution of the noisefree data is assumed. In particular, we impose a multicomponent prior distribution (to be called hereafter prior) on the noise-free wavelet coefficients in a given subband. As an optimization criterion, we adopt minimization of the mean squared error, i.e., the Bayesian risk is a quadratic loss function. The minimum mean-squared error (MMSE) estimate is the posterior conditional mean
density for [29], [30]
1867
is required. A widely used prior is Jeffrey’s prior (8)
E. Noise-Free Image A single-component image of the same scene as captured by the multicomponent image is assumed to be available. In this paper, we want to employ the correlation between the multicomponent image and this single-component image to improve the denoising of the first. For this, both images need to be geometrically registered and the images should be correlated. This correlation will be accounted for in a Bayesian denoising procedure, by treating the single-component image as extra prior information. Denote the wavelet detail coefficients of a noise-free singlecomponent image. We assume that and are jointly GSM as the vector, concatenating and distributed. If we write , we can rewrite (6) as (9) where becomes
is zero mean Gaussian distributed; (7) then (10)
where
is a zero-mean Gaussian with covariance , or by taking the expectations over , with . The conditional pdf of given can be written as a mixture of conditional Gaussian distributions (4) (11) Assuming, e.g., a Gaussian prior for the noise-free signal , the above MMSE estimate is the Wiener filter
where
is the Gaussian
with [31]
(5) D. GSM Model
(12)
A GSM model models the noise-free wavelet coefficients as (6) is an indewhere is zero-mean Gaussian distributed and is pendent positive scalar random variable. The probability then a mixture of Gaussians:
Notice that is not dependent on , and we can, therefore, . Also, its value is generally not equal to zero, denote it as making (11) is not strictly a GSM, but a shifted version. F. Bayesian Estimation The Bayes least squares estimate
is given by
(7) is the mixing density, and is a zero mean where or, by taking expectaGaussian with covariance . GSM densities are tions over , with symmetric and zero-mean and heavier tailed than Gaussians. These are known to better model the shape of the wavelet coefficient marginals than Gaussians. In [10], GSMs were applied to model local spatial neighborhoods of wavelet coefficients to denoise greylevel images. In this paper, we apply them to model the pixels of a multicomponent image, where they account for the strong correlation between the different components. A prior
(13) The first term in (13) can be obtained, using Bayes’ rule (14)
1868
with
IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 16, NO. 7, JULY 2007
given by (8) and
leading to (15)
We will now calculate the term rule, one obtains
(22)
. Using Bayes Now define: . Then
and the vectors:
, and
(16) (23) Since conditioned on is given by and since we have assumed that is independent of both , and , it follows that and , conditioned on are independent
, and are independent of , they need to be calSince culated only once. Remark that also for the calculation of the posterior distribution of in (15), their definition can be applied (24)
(17) The second expression appears when using Bayes rule on . Finally, canceling and skipping and because they do not depend on , the last expression is obtained. comes from the observation model The pdf (18) while the second term is given by (11). Then
where , and are the components of and the diagonal elements of respectively. When the covariances are estimated globally over the entire image, is determined globally. Then, , and are calculated for each pixel, after which the integral over is performed. For the calculation of the integral over , we follow [10], where is sampled logarithmically, which requires fewer samples than linear sampling. Optimal results were obtained for an (only interval from 20.5 until 3.5, with a step size of 2 for 13 samples). IV. EXPERIMENTS AND DISCUSSION
(19)
III. IMPLEMENTATION The terms (14) and (19) can be calculated by estimating (12) from the observed image and the noise-free image . In this paper, all estimations are done globally, over the entire image. The expression (19) is dependent on , and should be recalculated for each in the integral of (13). In order to simplify the dependence on of this expression, let us define: , i.e., . Then, (19) can be rewritten as
(20) Let be the square root of the matrix (i.e., ), and let denote the eigenvector/eigenvalue decomposition of . Then the matrix
(21)
To validate and compare the proposed techniques, we perform experiments on different data sets. The data sets contain a multicomponent image and a noise-free single component image of the same scene, that is geometrically registered with the multicomponent image. In the first experiment, MRI data is used. A typical MR imaging system is capable of acquiring image data at different acquisition modes. This leads to series of images with complementary information. Also, modern, fast acquisition procedures allow for performing time-critical experiments, as, e.g., functional MRI and diffusion tensor MRI, where multicomponent images are formed. Of course, spatial resolution and SNR are limited in these cases. During the complete acquisition time of such series of images, a high-quality image may be formed as well. Such image can be applied as a noise-free image in the proposed procedure. In the other two experiments, remote sensing data was used. Multispectral or hyperspectral sensors acquire multicomponent images. A higher quality image can be obtained from, e.g., another sensor, or another part of the reflectance spectrum with higher SNR. We will apply such an image as noise-free image in the proposed procedure. The proposed technique has some advantages over standard denoising of each image component. It accounts for the full covariance structure of the multicomponent image, it makes use of a prior that approximates well the marginal distributions of the wavelet coefficients, and it makes use of a noise-free image as extra prior information. In order to demonstrate the impact of
SCHEUNDERS AND DE BACKER: WAVELET DENOISING OF MULTICOMPONENT IMAGES
Fig. 1.
T
1869
2; P D; and T 1 noisy components of the MRI image and noise-free T 1 image.
these three advantages on the improvement of the performance, the following denoising strategies are applied to compare. • Single component denoising, where each band image is treated independently. A Gaussian prior model is assumed. This is the Wiener filter in wavelet domain. • Multicomponent denoising, where the correlations between the different components are exploited. A Gaussian prior is assumed. • Multicomponent denoising, assuming the GSM model as a prior to better model the marginal distributions of the wavelet coefficients. • Multicomponent denoising, using a noise-free single component image, using a Gaussian prior. • Multicomponent denoising, using a noise-free single component image, using a GSM prior. If the noise is simulated, the performance of the different denoising techniques can be quantitatively described by the PSNR (in dB), defined for 8-bit images as
where is the mean squared error between the original noiseless image and the denoised image . In case of multicomponent images, all bands are separately histogram-stretched to 8 bit, and the mse is averaged over all bands to obtain the PSNR. For the wavelet decomposition, a four-level decomposition of the nondecimated wavelet transform, as described in Section II-A, is used. For the wavelet filters, we use the Daubechies symmlets family [32].
TABLE I PSNR (IN DECIBELS) VALUES AFTER DENOISING OF MRI DATA SET
the image, by using the GSM prior and by using the noise-free image as extra prior information. When including the auxiliary image, sharper edges are obtained, as visible in Fig. 2(d) versus Fig. 2(e) and Fig. 2(f) versus Fig. 2(g). From Fig. 2(f) and (g), one can clearly observe that the noise is better removed using the GSM model. To asses the dependency of the results on the registration accuracy between the auxiliary image and the multicomponent image, we shifted the auxiliary image and repeated the experdB. Fig. 3 shows iment at initial noise level the PSNR after denoising in function of the shift. The dashed line show the results without the use of auxiliary image. After a steep drop in performance at subpixel misregistration errors, the quality slowly reduces until a shift of three pixels, where the misregistration error removes the added value of incorporating the auxiliary image. This demonstrates that for optimal results an almost perfect registration is required, but it is still worthwhile to use the auxiliary image even without perfect registration.
A. Experiment 1: MRI Data
B. Experiment 2: Multispectral Remote Sensing Data
The MRI data set consists of three different simulated MRI , and PD). The data set acquisitions of the same subject ( was obtained from the the Brainweb site [33]. On these images Rician distributed noise of different amplitude is simulated. The noise was independent for each component. A fourth image was , and kept noise-free. The images are shown in generated Fig. 1. Table I shows the obtained PSNR values. Fig. 2 shows a decomponent, for the different denoising tailed window of the procedures. In both Table I and Fig. 1, one can observe a gradual improvement of the results, by accounting for the covariances of
In the second experiment, a seven-band Landsat TM image was taken over the Winnipeg area containing both rural as urban areas. From this image, the thermal band (band 6) was removed and five spectral bands were corrupted with additive Gaussian noise. The noise was chosen independently for each , with varying standard deviation. The red band: band (band 3) served as the noise free image. Table II shows the results of the denoising procedures for this image: a) original Landsat image, b) detail of Landsat image, c) noisy Landsat image, d) single-component denoising with Gaussian prior, e) multicomponent denoising with Gaussian prior, f) multicomponent denoising with GSM prior, g) multicomponent denoising
1870
IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 16, NO. 7, JULY 2007
Fig. 2. (a) Detail of original T 2 MRI image of Fig. 1; (b) T 2 image with simulated Rician noise (PSNR = 34:13); (c) single-component denoising with Gaussian prior (PSNR = 36:21); (d) multicomponent denoising with Gaussian prior (PSNR = 37:49); (e) multicomponent denoising with Gaussian prior and noise-free image (PSNR = 38:71); (f) multicomponent denoising with GSM prior (PSNR = 39:54); (g) multicomponent denoising with GSM prior and noise-free image (PSNR = 40:92).
(Fig. 4). Here, the performance increase due to better description of the prior distribution is not as large as in the MRI case.
Fig. 3. Denoising quality in function of misregistration of auxiliary image. The dashed line indicates the quality achieved with multicomponent denoising without auxiliary image.
TABLE II PSNR (IN DB) VALUES AFTER DENOISING OF MULTISPECTRAL DATA SET
C. Experiment 3: Hyperspectral Remote Sensing Data In the third experiment, an AVIRIS image over Cuprite, Nevada, was taken. Here, we concentrated on the infra-red (IR) part of the spectrum, containing noisy bands near the water vapor absorption band. For this experiment, five bands were uniformly selected in the 2.3–2.5- m region. No additional noise was added. Fig. 5(a) shows one band of the image near 2.5 m. The noise covariance matrix is estimated as follows: let denote the wavelet coefficients of the first resolution level and orientation subband from spectral band . Then the diagonal elements of the noise covariance are estimated by the classical median estimator [34] (25) The off-diagonal elements between spectral bands estimated using [12]
and
are (26)
with Gaussian prior and noise-free image, and h) multicomponent denoising with GSM prior and noise-free image. Again, a gradual improvement can be observed, when accounting for the covariances of the image, using the GSM prior and using the noise-free image as extra prior information
where , and . For the auxiliary image, an image near 2 m was chosen [Fig. 5(b)]. In this region of the spectrum, less noise is present. A detailed window of the denoised images are shown in Fig. 5. Similar results as in the previous experiments are observed. The single component denoising procedure tends to blur the image excessively, removing most of the details. These details are better preserved in the multicomponent approaches.
SCHEUNDERS AND DE BACKER: WAVELET DENOISING OF MULTICOMPONENT IMAGES
1871
Fig. 4. (a) One band of Landsat multispectral image; (b) detail of image; (c) detail image with simulated Gaussian noise ( = 15); (d) single-component denoising with Gaussian prior; (e) multicomponent denoising with Gaussian prior; (f) multicomponent denoising with GSM prior; (g) multicomponent denoising with Gaussian prior and noise-free image; (h) multicomponent denoising with GSM prior and noise-free image.
Fig. 5. (a) Detail of original 2.479-m band of the AVIRIS cuprite image; (b) 1.991-m band of the cuprite image, used as the noise-free image; (c)–(g) results for the 2.479-m band after denoising; (c) single-component denoising with Gaussian prior; (d) multicomponent denoising with Gaussian prior; (e) multicomponent denoising with GSM prior; (f) multicomponent denoising with Gaussian prior and noise-free image; (g) multicomponent denoising with GSM prior and noise-free image.
When adding the noise-free image, the difference between Gaussian and GSM model becomes clear; with the Gaussian prior, the noise is not completely removed, while when using the GSM model, the noise is effectively removed and sharp edges are clearly retained.
V. CONCLUSION In this paper, a Bayesian wavelet-based denoising procedure for multicomponent images was proposed. The procedure makes use of a noise-free single component image as prior
1872
information. The prior model for the wavelet coefficient marginals is a GSM model. Experiments were performed with simulated noise on MRI imagery, multispectral and hyperspectral remote sensing images. The results show a gradual improvement, when applying a technique that 1) fully accounts for the multicomponent image covariances, 2) makes use of Gaussian Scale Mixtures as prior models that approximate well the marginal distributions of the wavelet coefficients, and 3) makes use of a noise-free image as extra prior information. REFERENCES [1] I. L. Thomas, V. M. Benning, and N. P. Ching, Classification of Remotely Sensed Images. Bristol, U.K.: Hilger, 1987. [2] C. Lee and D. A. Landgrebe, “Analyzing high-dimensional multispectral data,” IEEE Trans. Geosci. Remote Sens., vol. 31, no. 4, pp. 388–400, Apr. 1993. [3] T. Taxt and A. Lundervold, “Multispectral analysis of the brain using magnetic resonance imaging,” IEEE Trans. Med. Imag., vol. 13, no. 3, pp. 470–481, Mar. 1994. [4] C. Busch, “Wavelet based texture segmentation of multi-modal tomographic images,” Comput. Graph., vol. 21, no. 3, pp. 347–358, 1997. [5] A. A. Green, M. Berman, P. Switzer, and M. D. Craig, “A transformation for ordering multispectral data in terms of image quality with implications for noise removal,” IEEE Trans. Geosci. Remote Sens., vol. 26, no. 1, pp. 65–74, Jan. 1988. [6] J. L. Starck and P. Querre, “Multispectral data restoration by the wavelet Karhunen-Loève transform,” Signal Process., vol. 81, pp. 2449–2459, 2001. [7] D. L. Donoho and I. M. Johnstone, “Ideal spatial adaptation via wavelet shrinkage,” Biometrika, vol. 81, pp. 425–455, 1994. [8] S. G. Chang, B. Yu, and M. Vetterli, “Adaptive wavelet thresholding for image denoising and compression,” IEEE Trans. Image Process., vol. 9, no. 9, pp. 1532–1546, Sep. 2000. [9] M. K. Mihgak, I. Kozintsev, K. Ramchandran, and P. Moulin, “Lowcomplexity image denoising based on statistical modeling of wavelet coefficients,” IEEE Signal Process. Lett., vol. 6, no. 12, pp. 300–303, Dec. 1999. [10] J. Portilla, V. Strela, M. J. Wainwright, and E. P. Simoncelli, “Image denoising using scale mixtures of Gaussians in the wavelet domain,” IEEE Trans. Image Process., vol. 12, no. 11, pp. 1338–1351, Nov. 2003. [11] A. Benazza-Benyahia and J. C. Pesquet, “An extended sure approach for multicomponent image denoising,” in Proc. IEEE Int. Conf. Acoustics, Speech, Signal Processing, 2004, vol. 2, p. ii-945-8. [12] A. Benazza-Benyahia and J.-C. Pesquet, “Building robust wavelet estimators for multicomponent images using Steins’ principle,” IEEE Trans. Image Process., vol. 14, no. 11, pp. 1814–1830, Nov. 2005. [13] A. Piˇzurica, P. Scheunders, and W. Philips, “Multiresolution multispectral image denoising based on probability of presence of features of interest,” in Proc. Advanced Concepts in Image and Vision Systems, Brussels, Belgium, 2004, pp. 357–364. [14] A. Piˇzurica, B. Huysmans, P. Scheunders, and W. Philips, “Wavelet domain denoising of multispectral remote sensing imagery adapted to the local spatial and spectral context,” presented at the IGGARS, 2005. [15] P. Scheunders, “Wavelet thresholding of multivalued images,” IEEE Trans. Image Process., vol. 13, no. 4, pp. 475–483, Apr. 2004. [16] P. Scheunders and J. Driesen, “Least-squares interband denoising of color and multispectral images,” in Proc. IEEE Int. Conf. Image Processing, 2004, pp. 985–988. [17] S. de Backer, A. Piˇzurica, B. Huysmans, W. Philips, and P. Scheunders, “Denoising of multispectral images using wavelet least-squares estimators,” Image Vis. Comput., 2006. [18] H. Li, B. S. Manjunath, and S. K. Mitra, “Multisensor image fusion using the wavelet transform,” Graph. Models Image Process., vol. 57, no. 3, pp. 235–245, 1995. [19] D. A. Yocky, “Image merging and data fusion by means of the discrete two-dimensional wavelet transform,” J. Opt. Soc. Amer. A, vol. 12, no. 9, pp. 1834–1841, 1995.
IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 16, NO. 7, JULY 2007
[20] D. Yocky, “Multiresolution wavelet decomposition image merger of landsat thematic mapper and spot panchromatic data,” Photogramm. Eng. Remote Sens., vol. 62, pp. 1067–1074, 1996. [21] J. Zhou, D. Civco, and J. Silander, “A wavelet transform method to merge landsat tm and spot panchromatic data,” Int. J. Remote Sens., vol. 19, no. 4, pp. 743–757, 1998. [22] J. Nunez, X. Otazu, O. Fors, A. Prades, V. Pala, and R. Arbiol, “Image fusion with additive multiresolution wavelet decomposition, applications to spot+landsat images,” J. Opt. Soc. Amer. A, vol. 16, pp. 467–474, 1999. [23] J. Nunez, X. Otazu, O. Fors, A. Prades, V. Pala, and R. Arbiol, “Multiresolution-based image fusion with additive wavelet decomposition,” IEEE Trans. Geosci. Remote Sens., vol. 37, no. 3, pp. 1204–1211, May 1999. [24] T. Ranchin and L. Wald, “Fusion of high spatial and spectral resolution images: The arsis concept and its implementation,” Photogramm. Eng. Remote Sens., vol. 66, pp. 49–61, 2000. [25] R. C. Hardie, M. T. Eismann, and G. L. Wilson, “Map estimation for hyperspectral image resolution enhancement using an auxiliary sensor,” IEEE Trans. Image Process., vol. 13, no. 9, pp. 1174–1184, Sep. 2004. [26] S. Mallat, “A theory for multiresolution signal decomposition: The wavelet representation,” IEEE Trans. Pattern Anal. Mach. Intell., vol. 11, no. 7, pp. 674–693, Jul. 1989. [27] S. Mallat, A Wavelet Tour of Signal Processing. London, U.K.: Academic Press, 1998. [28] M. Jansen, Noise Reduction by Wavelet Thresholding. New York: Springer-Verlag, 2001. [29] G. E. P. Bax and C. Tiao, Bayesian Inference in Statistical Analysis. Reading, MA: Addison-Wesley, 1992. [30] M. Figeiredo and R. Nowak, “Wavelet-based image estimation: An empirical bayes approach using Jeffrey’s noninformativeprior,” IEEE Trans. Image Process., vol. 10, no. 9, pp. 1322–1331, Sep. 2001. [31] S. M. Kay, Fundamentals of Statistical Signal Processing: Estimation Theory. Englewoods Cliffs, NJ: Prentice-Hall, 1993. [32] I. Daubechies, Ten Lectures on Wavelets. Philadelphia, PA: SIAM, 1992. [33] C. A. Cocosco, R. K.-S. Kwan, V. Kollokian, and A. C. Evans, Brainweb: Online Interface to a 3-D MRI Simulated Brain Database. [34] D. L. Donoho and I. M. Johnstone, “Adapting to unknown smoothness via wavelet shrinking,” J. Amer. Statist. Assoc., vol. 90, no. 432, pp. 1200–1224, 1995.
Paul Scheunders graduated in physics in 1983 and received the Ph.D. degree in physics, with work in the field of statistical mechanics, in 1990, from the University of Antwerp, Antwerp, Belgium. In 1991, he became a Research Associate at the Vision Lab, Department of Physics of the University of Antwerp, where he became a Lecturer. He has published over 100 papers in international journals and proceedings in the field of image processing and pattern recognition. His research interest includes wavelets and multispectral image processing.
Steve De Backer received the B.S. degree in physics from the University of Hasselt, Belgium, in 1994, and the M.S. and Ph.D. degrees physics from the University of Antwerp, Antwerp, Belgium, in 1996 and 2002, respectively. Since then, he has been a Postdoctorate Fellow at the Visionlab, University of Antwerp. His research focus is on image modeling, pattern recognition, and hyperspectral image analysis.
