Wavelet Thresholding of Multivalued Images - Semantic Scholar
1
Wavelet Thresholding of Multivalued Images P. Scheunders Vision Lab, Department of Physics, University of Antwerp Groenenborgerlaan 171, 2020 Antwerpen, Belgium Tel.: +32 3 218 04 39; Fax.: +32 3 218 03 18 Email: [email protected]
Abstract— In this paper a denoising technique for multivalued images exploiting interband correlations is proposed. A redundant wavelet transform is applied and denoising is applied by thresholding wavelet coefficients. Specific functions of the wavelet coefficients are defined that exploit interscale and/or interband correlation of the signal. Three functions are studied: the square of the wavelet coefficients, products of coefficients at adjacent scales and products of coefficients from different bands. For these functions, the signal and noise probability density functions (pdf ) become more separated. The high signal correlation between bands is exploited by summing these products over all bands, in this way separating noise and signal pdfs even more. The noise pdf of the proposed quantities is derived analytically and from this, a wavelet threshold is derived. The technique is demonstrated to outperform single band wavelet thresholding on multispectral remote sensing images and on multimodal MRI images.
EDICS: 2-WAVP, 2-COLO I. Introduction With 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, T2, proton density, diffusion, ...), different CT and nuclear medicine imaging modalities. In multivalued imagery, noise handling may be important. Due to physical constraints, a trade-off exists between spatial and spectral resolution and SNR. Also, noise filtering and image enhancement can drastically facilitate the processing and analysis of multivalued imagery. More 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]. Multivalued image noise is usually 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 to a coordinate system in which the noise covariance matrix is diagonalized, followed by a PCA [5]. Other techniques use wavelet transform for spatial and PCA for spectral decorrelation [6]. In this paper, a denoising strategy of multivalued images is proposed by applying the wavelet thresholding paradigm.
The paradigm of wavelet thresholding originates from [7]. There, a threshold value was derived using the discrete wavelet transform [8]. 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 for the signal pdf was assumed (see e.g. [9] and references therein). Also spatially adaptive thresholds were derived [10]. Further on, it was suggested that shift-invariant versions of the wavelet transform were favoured over the discrete wavelet transform. This was also suggested in [11]. There, a hard, scale dependent threshold was applied, with a threshold value t = cσ, with σ the standard deviation of the noise at a particular scale. For Gaussian noise, a value of c between 3 and 4 will suppress most of the noise contribution. In [12], a signal dependent near optimal value for c was derived. Another approach estimates the true signal by a linear minimal mean squared error (LMMSE) technique instead of thresholding. In [13], this technique was applied by Maximum A Posteriori (MAP) estimation of the (spatially local) signal variance, imposing some prior on it. In [14], the correlation between different wavelet scales was exploited by filtering. In [11], this was further improved. In [15], the interscale statistics were analyzed using information theory. In [16], the noise pdf of interscale products was derived analytically and applied for edge detection. In [17], the interscale correlations were exploited using a hybrid approach between LMMSE and hard thresholding, applying shift-invariant wavelet transforms. In this work, redundant wavelet frames were applied [18]. A scale adaptive hard thresholding procedure is derived for multivalued images, based on the modeling of the noise pdf. Multivalued sensor noise is supposed to be decorrelated a priori and is modeled as a Gaussian random vector with independent components of equal variance. Instead of treating each band separately, the high correlation that typically exists between multivalued bands is exploited. First, we make use of the signal interscale and interband correlation to separate noise from signal pdf. For this, three different products are calculated from the wavelet coefficients. The first function is merely the coefficients square. Secondly, the product of coefficients at adjacent scales separates signal from noise due to the high interscale correlation of the signal. Finally the product of coefficients at
2
different bands separates signal from noise due to the high interband correlation of the signal. Then, the proposed functions are summed over the N bands. This even more separates signal from noise, since √ the standard deviation of the noise grows with N , while the standard deviation of the signal grows with N . The pdfs of the noise for the three sums are calculated analytically. From the α-percentile values of these pdfs, threshold values are defined where α is close to one. The wavelet coefficients of the original bands are then thresholded using these threshold values. In an experimental section the proposed quantities are compared with respect to the signal and noise pdfs and to the obtained SNR. Finally, the three obtained thresholding procedures are compared to each other and to single-band thresholding. Experiments are conducted on a 7-band multispectral Landsat image and on a multimodal MRI image. II. Multivalued image wavelet thresholding Suppose an N -band multivalued image I(x, y) = (I1 (x, y), I2 (x, y), ..., IN (x, y)). The additive image noise is modeled as a Gaussian random vector with independent components of equal variance. For i = 1, · · · , N : Ii (x, y) = Iˆi (x, y) + ni (x, y) 1 −x2 fx (x = ni ) = gx (0, σ 2 ) = √ exp ( 2 ) 2σ 2πσ
(1)
A. The dyadic wavelet transform The wavelet transform employed in this work is based on non-orthogonal (redundant) discrete wavelet frames introduced by Mallat[18], [19]. Define a 2-D differentable smoothing function θ(x, y). Define ψ x (x, y) =
∂θ(x, y) ∂θ(x, y) and ψ y (x, y) = ∂x ∂y
(2)
ψ x (x, y) and ψ y (x, y) are quadratic spline wavelets of compact support. The wavelet transform of an image band Ii (x, y) is then defined by: x Di,s (x, y)
=
Ii ∗ ψsx (x, y)
and
y Di,s (x, y)
=
Ii ∗ ψsy (x, y)
(3)
where ∗ denotes the convolution operator and ψsx (x, y) =
1 1 x x y x y ψ ( , ) and ψsy (x, y) = 2 ψ y ( , ) s2 s s s s s x
(4)
y
denote the dilations of the functions ψ and ψ . s is the scale parameter which commonly is set equal to 2j with j = 1, ..., d. This yields the so called dyadic wavelet transform y x of depth d. Di,j and Di,j are referred to as the detail images of Ii , since they contain horizontal and vertical details of ∗ Ii at scale j. Let us denote Di,j where ∗ stands for x or y. B. Sums of detail coefficients over bands After wavelet transformation, the noise of the detail images will be Gaussian [11]: ∗ ∗ ∗ ˆ i,j Di,j (x, y) = D (x, y) + Ni,j (x, y)
fx (x =
∗ ) Ni,j
=
gx (0, σj2 )
σj =
kψj∗ kσ
(5)
Summation over all bands i = 1, · · · , N gives: N X i=1
∗ Di,j (x, y) =
N X
∗ ˆ i,j D (x, y) +
i=1
N X
∗ Ni,j (x, y)
(6)
i=1
Let us follow an intuitive argumentation to show that the sum of detail coefficients separates the signal from the noise, an effect that increases with N . Since the noise between different bands is supposed to be uncorrelated, the second term in eqn. 6 averages the noise. The √obtained noise pdf has a standard deviation that grows ∼ N . The signals between different bands in multivalued images are usually highly correlated, thus the standard deviation of the pdf of the first term grows ∼ N . To show that the latter holds, even for small correlations, consider N random variables drawn from the same pdf with the same standard deviation σ, and a small but equal correlation ρ between them. The standard deviation of the sum of two such variables, up to first p √ order in ρ, is given by: σ[xi + xj ] = σ 2(1 + ρ) ' 2σ(1 + ρ/2). For large N , we find: PN −1 √ ! √ PN i i=1 √ σ[ n=1 xn ] ' N σ + ρσ N √ 2 ' N σ + ρ( N σ) (7) 3 This shows that even for small correlations, the standard deviation of the sum of N random variables from any pdf is proportional to N . C. Products of detail coefficients The argument that summing over bands separates signal from noise also holds for specific functions of the detail coefficients. Three different products are investigated in this paper: the square of the detail coefficients, products of coefficients at adjacent scales and products of coefficients at different bands. Each of these products separates signal from noise in a specific way. C.1 square of detail coefficients Let us square the detail coefficients. Since the noise is not correlated with the signal: ∗ 2 ˆ ∗ )2 (x, y) + (N ∗ )2 (x, y) (8) ) (x, y) ' (D (Di,j i,j i,j The square of the coefficients separates signal from noise since real edge coefficients tend to become larger for coarser scales, while noise coefficients become smaller. By squaring, this effect becomes more prominent. Let us calculate the noise pdf for this square of detail coefficients. In general, the pdf of the square of a random √ variable x with pdf fx is given by: fy (y = x2 ) = √1y fx ( y), The noise pdf of the square of the wavelet detail coefficients is therefore given by ! 1 y 1 ∗ 2 fy (y = (Ni,j ) ) = √ (9) √ exp − 2 2σj 2πσj y ∗ 2 This pdf has an expected value of E[(Ni,j ) ] = σj2 and a ∗ 2 4 variance of V [(Ni,j ) ] = 2σj .
3
of 0:
C.2 interscale products In [14], it was proposed to denoise based on the signal correlation of wavelet coefficients between adjacent scales. Using this, it was suggested to use the product of two ∗ ∗ )(x, y) wavelet coefficients at adjacent scales (Di,j · Di,j+1 to separate signal from noise. Also for this product holds: ∗ ∗ )(x, y) ' · Di,j+1 (Di,j ∗ ∗ ∗ ∗ ˆ ˆ )(x, y) (10) · Ni,j+1 )(x, y) + (Ni,j (Di,j · Di,j+1
Since real edge coefficients are highly correlated, the signal product will be large. For noise coefficients, the correlation between adjacent scales is low, and therefore the product will be small. To calculate the noise pdf of this product, first notice that, due to the non-orthogonality of the proposed wavelet frames, the noise between adjacent scales is correlated. The ∗ ∗ joint pdf of Ni,j and Ni,j+1 is given by: ∗ ∗ , x2 = Ni,j+1 )= fx1 ,x2 (x1 = Ni,j −x2 σ2 −x2 σ2 +2ρ j,j+1 x1 ·x2 σj σj+1 1 j+1 2 j 1 2βσj σj+1 p e (11) 2π σj σj+1 β
with β = σj σj+1 (1 − ρ2j,j+1 ) and the correlation between scales j and j + 1 is given by: RR ∗ ∗ (x, y)dxdy ψj (x, y)ψj+1 ρj,j+1 = qR R (12) R R ∗ (x, y))2 dxdy (ψj∗ (x, y))2 dxdy (ψj+1
∗ ∗ The pdf fy (y = Ni,j · Ni,j+1 ) can be obtained by using the general property that the pdf of the product of 2 correlated random variables x1 and x2R, having joint pdf f (x1 , x2 ), is given by: fy (y = x1 · x2 ) = f (y/x2 , x2 )dx2 /|x2 | [20]. We find: ∗ ∗ )= · Ni,j+1 fy (y = Ni,j |y| 1 yρj,j+1 p K0 exp β β π σj σj+1 β
(13)
C.3 interband products Finally, the correlation of the signal between different bands can be exploited for separating noise from signal. For this, the product of wavelet coefficients at two different bands i and k is calculated. Also for this product, one finds:
·
∗ Nk,j )
1 = K0 πσj2
|y| σj2
!
(15)
∗ ∗ This pdf has an expected value of E[Ni,j · Nk,j ] = 0 and a ∗ ∗ 4 variance of V [Ni,j · Nk,j ] = σj .
D. Sums of products The products that were defined in the previous section separate the noise from the signal pdf. Now, the argumentation of section II-B can be applied: by summing the products over several bands the signal pdf can be further separated from the noise pdf. To calculate the noise pdfs from the sums of products, we use the general property that the pdf of the sum of 2 independent random variables y1 and y2 with pdfs fy1 is given by the convolution: fz (z = and fy2 respectively, R∞ y1 + y2 ) = −∞ fy1 (z − y2 )fy2 (y2 )dy2 [20]. D.1 sum of squared coefficients
For the sum of squared coefficients, one obtains after calculation: ! N N/2−1 X z −z ∗ 2 (Ni,j ) )= exp fz(N ) (z = (16) 2σj2 (2σj2 )N/2 Γ( N2 ) i=1
R∞ where Γ(x) = 0 tx−1 exp (−t)dt is the Gamma function. The pdf of eqn. 16 is called the Γ-pdf Γ(N, σj ). Its exPN ∗ 2 pected value is: E[ i=1 (Ni,j ) )] = N σj2 and its variance PN ∗ 2 is: V [ i=1 (Ni,j ) )] = 2N σj4 . For σj = 1 it is the χ2 -pdf with N degrees of freedom. For large values of N , the pdf becomes the gaussian g(N σj2 , 2N σj4 ). D.2 sum of interscale products For the sum of products of wavelet coefficients at adjacent scales, one obtains after calculation: zρj,j+1 × · = exp = β i=1 z (N −1)/2 |z| (17) K(N −1)/2 √ N N/2 β πβΓ( 2 )(σj σj+1 )
fz(N ) (z
where K0 (·) is the 0-th order modified Bessel function of ∗ the second kind. This pdf has an expected value of E[Ni,j · ∗ ∗ ∗ Ni,j+1 ] = ρj,j+1 σj σj+1 and a variance of V [Ni,j · Ni,j+1 ] = 2 (1 + 2ρ2j,j+1 )σj2 σj+1 .
∗ ∗ )(x, y) ' · Dk,j (Di,j ˆ ∗ )(x, y) + (N ∗ · N ∗ )(x, y) ˆ∗ · D (D k,j i,j k,j i,j
fy (y =
∗ Ni,j
N X
∗ Ni,j
∗ Ni,j+1 )
where Kn (·) is the n-th order modified Bessel function of the second kind. Since it is only defined for integer n, the number of sums N should be chosen The obPN odd. ∗ ∗ tained pdf has an expected value of E[ i=1 Ni,j · Ni,j+1 ]= PN ∗ ∗ N ρj,j+1 σj σj+1 and a variance of V [ i=1 Ni,j · Ni,j+1 ] = 2 N (1 + 2ρ2j,j+1 )σj2 σj+1 . D.3 sum of interband products
(14)
Since the noise between the bands is supposed to be decorrelated, the noise pdf of the product is similar to equation 13, with equal standard deviations and a correlation
Finally for the sum of products of wavelet coefficients at different bands, one has to be careful to sum independent variables only. Each band may appear only once in the sum. Therefore, the number of summations can only be B = N/2 for N even, or B = (N − 1)/2 for N odd. Of
4
course, one has the freedom to choose which bands i and k are to be multiplied. In order to take optimally advantage of the interband correlation, the following procedure is adopted. The correlation of all possible pairs is calculated and ranked. The pair with highest correlation is chosen first. Then each time, the subsequent highest correlation pairs are added, on the condition that each band of the pair has not been included before. In this way, exactly B summations are performed. For the noise pdf of this sum, one obtains: fz(B) (z =
B X i,k
∗ ∗ Ni,j · Nk,j )=
z (B−1)/2 K(B−1)/2 √ πΓ( B2 )(σj )B+1
|z| σj2
!
(18)
PB where the notation i,k denotes the sum over the pairs, as chosen in the proposed approach. Since Kn (·) is only defined for integer n, an additional limitation is that B should odd. The obtained pdf has an expected Pbe PB value B ∗ ∗ ∗ of E[ i,k Ni,j · Nk,j ] = 0 and a variance of V [ i,k Ni,j · ∗ 4 Nk,j ] = N σj . E. Thresholding From the calculated pdfs a thresholding procedure can be developed that removes most of the contribution of the noise. For this, the same argument as in [11] is applied, i.e. most of the noise contribution will be removed when most of the area of the noise pdf is within the threshold value. In [11], for a gaussian noise pdf, 3 to 4 times the noise standard deviation was chosen. In this work, the α(N ) percentile value of the obtained pdfs, f(α )(x) is chosen as a threshold value, with α close to one. In this way, three different thresholds can be defined: τ1
=
(N ) f(α) (
N X
∗ 2 ) ) (Ni,j
(19)
∗ ∗ ) · Ni,j+1 Ni,j
(20)
∗ ∗ Ni,j · Nk,j )
(21)
i=1
τ2
=
(N )
f(α) (
N X i=1
τ3
=
(N )
f(α) (
B X ik
The thresholds are calculated by numerical integration up till α of the pdfs, obtained from the analytical expressions of equations 16, 17 and 18. The thresholds are not applied on the sum of product quantities, but on the original band detail coefficients. For this, the following multivariate thresholding procedure is defined: at all positions where the sum of product is smaller than the obtained threshold value, the detail coefficients of all bands involved are removed. This procedure is expected to outperform single band thresholding, since signal and noise pdfs of the proposed sums of products are more separated than the pdfs of the original coefficients. If at a
particular position a band occurs with a signal coefficient, it is likely that it is accompanied by other signal coefficients at other bands at the same position and it will be retained after thresholding. On the other hand, a noise coefficient at one band is likely to be accompanied by other noise coefficients at different bands at the same position and it will be removed after thresholding. When for a pixel one or a few bands contain noise coefficients while the other bands are signal coefficients, the noise will not be removed, which results in a lower denoising performance. In a rare case, a pixel contains only one or a few signal coefficients. If these are not too large, they will be removed, which results in a signal deterioration. As will be shown in the experimental section, the latter two disadvantages are largely compensated by the advantages, resulting in an overall improved denoising performance. In the case of the sum of squared coefficients, only positive values occur. The thresholding procedure becomes: for all bands k = 1, · · · N : ∗ ∗ y) Dk,j (x, y) = Dk,j (x,P N ∗ 2 (22) if i=1 (Di,j ) (x, y) > τ1 ∗ Dk,j (x, y) = 0 otherwise In the case of the sum of interscale products, the obtained pdf is asymmetric around zero. The signal sum of products can be negative, revealing negative correlations (this can also be seen from figure 2c and 2d in the next section). To include these negative correlations, a second 0 (N ) PN ∗ ∗ · Ni,j+1 ) is necessary. threshold value τ2 = f(α0 ) ( i=1 Ni,j 0
α would then be close to zero. The thresholding procedure becomes: for all bands k = 1 · · · N : ∗ Dk,j (x, y) = 0 PN 0 ∗ ∗ if τ2 < i=1 (Di,j · Di,j+1 )(x, y) < τ2 (23) ∗ ∗ Dk,j (x, y) = Dk,j (x, y) otherwise
In the case of the sum of interband products, the obtained pdf is symmetric around zero. Here, only positive signal correlations are looked for (the largest positive interband correlations are summed). The thresholding procedure becomes: for all bands k = 1, · · · N : ∗ ∗ y) Dk,j (x, y) = Dk,j (x,P B ∗ ∗ (24) if i,k (Di,j · Dk,j )(x, y) > τ3 ∗ Dk,j (x, y) = 0 otherwise The question remains how to determine an optimal value for α. In this work, every effort was made to separate noise and signal pdfs as much as possible. Therefore, it is expected that an optimal value of α will not depend much on the signal, so that we can fix α to a value close to 1. In [12], a signal dependent threshold was derived for a gaussian noise pdf on wavelet coefficients. For this, a function, equivalent to the MSE was derived. If applied to the redundant wavelet representation of this work, the function at scale j is given by: ! Rt 2 x gx (0, σj2 )dx −t 2 Fj (t) = 2 σj − k R t − V [D˜j∗ ] (25) 2 )dx g (0, σ x j −t
5
with k the proportion of coefficients that has been removed image, including signal and noise as ’x’. From the plots, after applying the threshold t, and V [D˜j∗ ] the variance of one can verify that: the pdf of the coefficients after thresholding. This func- • the simulated noise pdfs fits correctly with the analytical tion has a signal dependent minimum that agrees with the pdfs. minimal MSE at a value of t which is then chosen as the • the signal pdfs of the products are broader than the noise optimal threshold value. The calculations to obtain this pdfs. This effect is limited, but seems to be highest for the function can be easily extended towards noise pdfs, other interscale products. than the normal pdf. For the symmetric sum of interband • when summing the products, the signal pdfs become products pdf, we obtain: much broader than the noise pdfs. This holds for all three products, but the effect seems to be highest for the interZ t B X scale products. k 4 2 ∗ ∗ Fj (t) = 2 Bσj − x fx (x = Ni,j · Nk,j )dx α(t) −t In order to quantify to what extent the proposed prodi,k ucts and sums separate signal from noise, the signal-tog B X noise ratio’s can be calculated. A definition of the signal∗ · D∗ ] (26) to-noise ratio of a signal is given by the ratio between the Di,j −V [ k,j i,k signal variance and the noise variance. Because of equations 8, 10 and 14, we find: Rt PB ∗ ∗ with α(t) = −t fx (x = i,k Ni,j · Nk,j )dx. As an example, s in figure 1 this function is plotted for a Landsat image with σ 2 (observed signal) − σ 2 (noise) simulated noise at scale 1, along with the MSE, in function SN R = (27) σ 2 (noise) of α(t). As can be seen, the two minima appear at approximately the same value for α, which is close to 0.98. Similar results were obtained for different scales, different images where the noise standard deviations are obtained from the ∗ and different noise standard deviation values. Values for analytical expressions. In table 1, the SNR’s for D1,j , P N ∗ the optimal value of α were all in the range 0.96-0.99. For D and the three products and their sums are shown i=1 i,j the sum of squares pdf a similar expression for the function in the case of the Landsat image for the first four scales. F can be obtained, and similar results were obtained. For From the table, one can observe that: the sum of interscale products, the asymmetric pdf poses a difficulty. The function F depends on two thresholds, • for all quantities the SNR grows with the scale. The and an optimum needs to be found in 2D. Results did not jump is largest between the first two scales, and from then consistently reveal a minimum in the expected range. In on it is roughly a factor of 2 for the wavelet coefficients fact, one can expect lesser denoising performance in this and a factor of 4 for the products. This effect is due to the case, due to the asymmetry. One of the consequences is observation that the noise standard deviation shrinks with e.g. that negative signal contributions can be corrupted by a factor of 2 between 2 consecutive scales, while the signal variance roughly remains constants between scales. higher-valued positive noise contributions. • For all three products the SNR is improved compared The time complexity of finding the signal dependent α to the wavelet coefficients. The improvement of the SNR is quite high, since the function F needs to be calculated compared to the squared coefficients is a factor of 2 for the for different values of α(t), with each time needs numeriinterscale products and a factor of 1.5 for the interband cal integrations. Since the optimal α seems to be rather products. This demonstrates that exploiting correlations signal independent, a fixed value for α can be chosen, e.g. between scales or between bands separates better the signal α = 0.98. This requires only one numerical integration of from the noise. the specific sum of product function, after which it can be • When summing over the bands, a further improvement applied to any image. of the SNR is obtained for all products, with roughly a factor of 2. The improvement with the sums over interband III. Experiments and discussion products is somewhat smaller, because of the smaller numFor a qualitative appreciation of the proposed tech- ber of summations. This demonstrates that exploiting the niques, the obtained pdfs from the proposed quantities are correlations between the bands by summing, improves the plotted in the case of a 7-band multispectral Landsat im- separation of signal and noise. age. On each band of the image independent additive gausTo demonstrate the denoising properties of the proposed sian noise is added with σ = 20. In total, 6 plots are shown techniques, the three thresholding procedures 22, 23 and 24 at the particularP scale j = 2: in figure 2a and b the plots are compared to each other and to single band thresholding. N ∗ 2 ∗ 2 ) ; in figure 2c and d the plots for (D1,2 ) and i=1 (Di,2 For the latter, 3 approaches are obtained by making use of PN ∗ ∗ ∗ ∗ and i=1 Di,2 · Di,3 and in figure 2e and the 3 proposed products, without summing. For all three · D1,3 for D1,2 PB ∗ ∗ ∗ ∗ f the plots for D1,2 · D2,2 and i,k Di,2 · Dk,2 . On each products, a threshold can be defined by the α-percentile plot, 3 pdfs are shown: the obtained analytical noise pdf value of the product pdfs, obtained from integration up till as a solid line, the noise pdf for a multivalued empty image α of the analytical expressions of equations 9, 13 and 15. 3 with simulated noise as ’+’ and the pdf of the observed different single band thresholding procedures are obtained:
6
For band i: ∗ ∗ Di,j = Di,j ∗ Di,j = 0 or ∗ Di,j = 0
∗ Di,j
or
=
∗ Di,j
∗ ∗ Di,j = Di,j ∗ Di,j = 0
∗ 2 ∗ 2 if (Di,j ) > f(α) (Ni,j ) otherwise
(28)
∗ ∗ if f(1−α) (Ni,j · Ni,j+1 ) ∗ ∗ ∗ ∗ < Di,j · Di,j+1 < f(α) (Ni,j · Ni,j+1 ) (29) otherwise
∗ ∗ ∗ ∗ if Di,j · Dk,j > f(α) (Ni,j · Nk,j ) otherwise
(30)
respectively, where in the last case band k is the one that is largest correlated with band i. Also, for comparison two simple standard wavelet denoising procedures were applied. The LMMSE technique estimates the signal wavelet coefficient of band i and scale j as: ˆ∗ ) σ ˆ 2 (D i,j ∗ ˆ i,j D (x, y) = 2 D∗ (x, y) ˆ ∗ ) i,j σj + σ ˆ 2 (D i,j
a standard technique, i.e by measuring the variance in a local region in the background, and found to be about 13. In figure 5a, one original band is shown. In figure 5b, the denoised band is shown, using the sum of interband product threshold. The same conclusions as for the Landsat image hold. IV. Conclusions In this paper, different denoising procedures for multivalued images were proposed, based on wavelet thresholding. The procedures take advantage of the interscale and interband correlations of the signals, by looking at interband sums of 3 specific products of wavelet coefficients. For these sums, the specific noise pdfs were analytically derived and from these pdfs, threshold values were obtained. From the experiments, we can conclude that all the multi band thresholding procedures outperform the standard single band wavelet thresholding procedure. The sum of interband products leads to the best thresholding procedure. References
(31)
ˆ ∗ ) an estimate for the signal variance. In [13], a with σ ˆ 2 (D i,j Maximum a Posteriori estimation was done for the signal ˆ∗ ) = variance. We will apply a simple estimation: σ ˆ 2 (D i,j P P ∗2 2 ∗ 2 2 ∗ σ (Di,j ) − σj , with σ (Di,j ) = x y Di,j (x, y), the observed signal variance. The wavelet thresholding technique of [11] will also be applied. The coefficients of band i at scale j are hard thresholded with a threshold value of t = 3.5σj . This technique will be referred to as ’PAN’. The proposed denoising techniques are applied to the same Landsat image as before. On the original image, Gaussian noise with different values of σ was added. The α-percentile values of the single band and multiband pdfs are calculated, with α = 0.98 in all cases. For the single band techniques, the results on the first band are given. In table 2, MSE values are given. In fig. 3, the original and noisy band 1 is shown. In figure 4, the denoised results are shown. From the table and the figures, one can conclude that: • the sum of squares, the interband and the sum of interband product thresholding results are superior to LMMSE and single band thresholding. The sum of interband product thresholding leads to the best results. • the interscale product and sum of interscale product thresholding perform poor. This does not agree with the qualitative results on the SNR. The reason can be found in the negative signal correlations and the asymmetric noise 0 pdf. Playing around with the values of α and α to optimize the resulting MSE never lead to improved results. Finally, an experiment is conducted on a noisy multimodal MRI image, diffusion- weighted at 6 independent imaging gradient directions. Noise is mostly thermal and is known to be approximately independent normal distributed. The noise standard deviation is estimated using
[1] [2] [3] [4] [5]
[6] [7] [8]
[9]
[10]
[11] [12] [13]
[14]
[15]
I.L. Thomas, V.M. Benning, and N.P. Ching, Classification of remotely sensed images, Adam Hilger, Bristol, 1987. C Lee and D.A. Landgrebe, “Analyzing high-dimensional multispectral data,” IEEE Trans. Geosci. Remote Sensing, vol. 31, no. 4, pp. 388–400, 1993. T. Taxt and A. Lundervold, “Multispectral analysis of the brain using magnetic resonance imaging,” IEEE Trans. Med. Imaging, vol. 13, no. 3, pp. 470–481, 1994. C. Busch, “Wavelet based texture segmentation of multi-modal tomographic images,” Comput. & Graphics, vol. 21, no. 3, pp. 347–358, 1997. 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 Transactions on Geoscience and Remote Sensing, vol. 26:(1), pp. 65–74, 1988. J.L. Starck and P. Querre, “Multispectral data restoration by the wavelet karhunen-loeve transform,” Signal Processing, vol. 81:(12), pp. 2449–2459, 2001. D.L. Donoho and I.M. Johnstone, “Ideal spatial adaptation via wavelet shrinkage,” Biometrika, vol. 81, pp. 425–455, 1994. S. Mallat, “A theory for multiresolution signal decomposition : the wavelet representation,” IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 11, no. 7, pp. 674–693, 1989. S. Grace Chang, B. Yu, and M. Vetterli, “Adaptive wavelet thresholding for image denoising and compression,” IEEE Transactions on Image Processing, vol. 9:(9), pp. 1532–1546, 2000. S. Grace Chang, B. Yu, and M. Vetterli, “Spatially adaptive wavelet thresholding with context modeling for image denoising,” IEEE Transactions on Image Processing, vol. 9:(9), pp. 1522–1531, 2000. Q. Pan, L. Zhang, and G. Dai aand H. Zhang, “Two denoising methods by wavelet transform,” IEEE Transactions on Signal Processing, vol. 47:(12), pp. 3401–3406, 1999. L. Zhang, P. Bao, and Q. Pan, “Threshold analysis in waveletbased denoising,” Electronic Letters, vol. 37:(24), pp. 1485–1486, 2001. M.K. Mihcak, I. Kozintsev, K. Ramachandran, and P. Moulin, “Low-complexity image denoising based on statistical modeling of wavelet coefficients,” IEEE Signal Processing Letters, vol. 6:(12), pp. 300–303, 1999. Y. Xu, J.B. Weaver, D.M. Healy, and J. Lu, “Wavelet transform domain filtering: a spatially selective noise filtration technique,” IEEE Transactions on Image Processing, vol. 3:(6), pp. 747–758, 1994. J. Liu and P. Moulin, “Information-theoretic analysis of interscale and intrascale dependencies between image wavelet coefficients,” IEEE Transactions on Image Processing, vol. 10:(11), pp. 1647–1658, 2001.
7
1100
1050
1000
950
900 0.8
0.82
0.84
0.86
0.88
0.9
0.92
0.94
0.96
0.98
1
Fig. 1. The functions F (t) and MSE(t) for the first scale of a 7-band Landsat image, in the case of the sum of interband products. On the x-axis, α(t) is plotted, on the y-axis, MSE(t)as ’+’. The function F (t), plotted as ’x’ is rescaled to fit in the plot.
Fig. 2.
∗ )2 and Plots of the pdfs of: a,b: (D1,2
∗ · D ∗ and D1,2 1,3 ∗ . Dk,2
PN
i=1
PN
i=1
Fig. 3. a: orginal and b: noisy band (σ = 20) of multispectral Landsat image.
∗ )2 ; c, d: (Di,2
∗ · D ∗ ; e, f: D ∗ · D ∗ and Di,2 1,2 2,2 i,3
PB
i,k
∗ · Di,2
[16] B.M. Sadler and A. Swami, “Analysis of multiscale products for step detection and estimation,” IEEE Transactions on Information Theory, vol. 45:(4), pp. 1043–1051, 1999. [17] L. Zhang, P.Bao, and X. Wu, “Hybrid inter- and intra- wavelet scale image restoration,” Pattern Recognition, vol. 36:(8), pp. 1737–1746, 2003. [18] S. Mallat and S. Zhong, “Characterization of signals from multiscale edges,” IEEE Trans. Pattern Anal. Machine Intell., vol. 14, pp. 710–732, 1992. [19] S. Mallat, Ed., A wavelet tour of signal processing, vol. 2-nd edition, Academic Press, 1999. [20] A. Papoulis, Probability, random variables and stochastic processes, Mc. Graw-Hill, N.Y., 1991.
Fig. 4. Denoised band of multispectral Landsat image; a: thresholding squared coefficients; b: thresholding sum of squared coefficients; c: thresholding interscale products; d: thresholding sum of interscale products; e: thresholding interband products; f: thresholding sum of interband products.
8
j 1 2 3 4 j 1 2 3 4
PN
∗ D1,j 0.95 3.04 6.54 13.4
i=1
∗ Di,j 1.91 6.65 15.6 30.3
PN
∗ 2 (D1,j ) 1.60 14.0 68.2 243
∗ 2 i=1 (Di,j )
3.02 28.1 138 481
PN
∗ ∗ D1,j · D1,j+1 4.44 31.1 130 428
i=1
∗ ∗ Di,j · Di,j+1 8.40 62.5 264 852
PB
i,k
∗ ∗ D1,j · D2,j 2.06 18.3 84.5 304 ∗ ∗ Di,j · Dk,j 2.80 25.7 127 434
TABLE I SNR values for the products and sums of products in the case of a 7-band Landsat image.
σ 10 15 20 25 30
LMMSE 87 167 250 326 395
PAN 108 200 285 386 471
PN
∗ 2 i=1 (Di,j )
83 149 217 283 348
∗ ∗ D1,j · D1,j+1 106 203 284 396 496
PN
i=1
∗ ∗ Di,j · Di,j+1 100 176 253 330 404
∗ ∗ D1,j · D2,j 85 148 236 292 368
PB
i,k
∗ ∗ Di,j · Dk,j 82 142 205 269 333
TABLE II Obtained MSE values after denoising, applying the proposed thresholding techniques on a 7-band Landsat image.
Fig. 5. a: Original band of a multimodal MRI image; b: thresholding sum of interband products.
Wavelet Thresholding of Multivalued Images P. Scheunders Vision Lab, Department of Physics, University of Antwerp Groenenborgerlaan 171, 2020 Antwerpen, Belgium Tel.: +32 3 218 04 39; Fax.: +32 3 218 03 18 Email: [email protected]
Abstract— In this paper a denoising technique for multivalued images exploiting interband correlations is proposed. A redundant wavelet transform is applied and denoising is applied by thresholding wavelet coefficients. Specific functions of the wavelet coefficients are defined that exploit interscale and/or interband correlation of the signal. Three functions are studied: the square of the wavelet coefficients, products of coefficients at adjacent scales and products of coefficients from different bands. For these functions, the signal and noise probability density functions (pdf ) become more separated. The high signal correlation between bands is exploited by summing these products over all bands, in this way separating noise and signal pdfs even more. The noise pdf of the proposed quantities is derived analytically and from this, a wavelet threshold is derived. The technique is demonstrated to outperform single band wavelet thresholding on multispectral remote sensing images and on multimodal MRI images.
EDICS: 2-WAVP, 2-COLO I. Introduction With 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, T2, proton density, diffusion, ...), different CT and nuclear medicine imaging modalities. In multivalued imagery, noise handling may be important. Due to physical constraints, a trade-off exists between spatial and spectral resolution and SNR. Also, noise filtering and image enhancement can drastically facilitate the processing and analysis of multivalued imagery. More 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]. Multivalued image noise is usually 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 to a coordinate system in which the noise covariance matrix is diagonalized, followed by a PCA [5]. Other techniques use wavelet transform for spatial and PCA for spectral decorrelation [6]. In this paper, a denoising strategy of multivalued images is proposed by applying the wavelet thresholding paradigm.
The paradigm of wavelet thresholding originates from [7]. There, a threshold value was derived using the discrete wavelet transform [8]. 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 for the signal pdf was assumed (see e.g. [9] and references therein). Also spatially adaptive thresholds were derived [10]. Further on, it was suggested that shift-invariant versions of the wavelet transform were favoured over the discrete wavelet transform. This was also suggested in [11]. There, a hard, scale dependent threshold was applied, with a threshold value t = cσ, with σ the standard deviation of the noise at a particular scale. For Gaussian noise, a value of c between 3 and 4 will suppress most of the noise contribution. In [12], a signal dependent near optimal value for c was derived. Another approach estimates the true signal by a linear minimal mean squared error (LMMSE) technique instead of thresholding. In [13], this technique was applied by Maximum A Posteriori (MAP) estimation of the (spatially local) signal variance, imposing some prior on it. In [14], the correlation between different wavelet scales was exploited by filtering. In [11], this was further improved. In [15], the interscale statistics were analyzed using information theory. In [16], the noise pdf of interscale products was derived analytically and applied for edge detection. In [17], the interscale correlations were exploited using a hybrid approach between LMMSE and hard thresholding, applying shift-invariant wavelet transforms. In this work, redundant wavelet frames were applied [18]. A scale adaptive hard thresholding procedure is derived for multivalued images, based on the modeling of the noise pdf. Multivalued sensor noise is supposed to be decorrelated a priori and is modeled as a Gaussian random vector with independent components of equal variance. Instead of treating each band separately, the high correlation that typically exists between multivalued bands is exploited. First, we make use of the signal interscale and interband correlation to separate noise from signal pdf. For this, three different products are calculated from the wavelet coefficients. The first function is merely the coefficients square. Secondly, the product of coefficients at adjacent scales separates signal from noise due to the high interscale correlation of the signal. Finally the product of coefficients at
2
different bands separates signal from noise due to the high interband correlation of the signal. Then, the proposed functions are summed over the N bands. This even more separates signal from noise, since √ the standard deviation of the noise grows with N , while the standard deviation of the signal grows with N . The pdfs of the noise for the three sums are calculated analytically. From the α-percentile values of these pdfs, threshold values are defined where α is close to one. The wavelet coefficients of the original bands are then thresholded using these threshold values. In an experimental section the proposed quantities are compared with respect to the signal and noise pdfs and to the obtained SNR. Finally, the three obtained thresholding procedures are compared to each other and to single-band thresholding. Experiments are conducted on a 7-band multispectral Landsat image and on a multimodal MRI image. II. Multivalued image wavelet thresholding Suppose an N -band multivalued image I(x, y) = (I1 (x, y), I2 (x, y), ..., IN (x, y)). The additive image noise is modeled as a Gaussian random vector with independent components of equal variance. For i = 1, · · · , N : Ii (x, y) = Iˆi (x, y) + ni (x, y) 1 −x2 fx (x = ni ) = gx (0, σ 2 ) = √ exp ( 2 ) 2σ 2πσ
(1)
A. The dyadic wavelet transform The wavelet transform employed in this work is based on non-orthogonal (redundant) discrete wavelet frames introduced by Mallat[18], [19]. Define a 2-D differentable smoothing function θ(x, y). Define ψ x (x, y) =
∂θ(x, y) ∂θ(x, y) and ψ y (x, y) = ∂x ∂y
(2)
ψ x (x, y) and ψ y (x, y) are quadratic spline wavelets of compact support. The wavelet transform of an image band Ii (x, y) is then defined by: x Di,s (x, y)
=
Ii ∗ ψsx (x, y)
and
y Di,s (x, y)
=
Ii ∗ ψsy (x, y)
(3)
where ∗ denotes the convolution operator and ψsx (x, y) =
1 1 x x y x y ψ ( , ) and ψsy (x, y) = 2 ψ y ( , ) s2 s s s s s x
(4)
y
denote the dilations of the functions ψ and ψ . s is the scale parameter which commonly is set equal to 2j with j = 1, ..., d. This yields the so called dyadic wavelet transform y x of depth d. Di,j and Di,j are referred to as the detail images of Ii , since they contain horizontal and vertical details of ∗ Ii at scale j. Let us denote Di,j where ∗ stands for x or y. B. Sums of detail coefficients over bands After wavelet transformation, the noise of the detail images will be Gaussian [11]: ∗ ∗ ∗ ˆ i,j Di,j (x, y) = D (x, y) + Ni,j (x, y)
fx (x =
∗ ) Ni,j
=
gx (0, σj2 )
σj =
kψj∗ kσ
(5)
Summation over all bands i = 1, · · · , N gives: N X i=1
∗ Di,j (x, y) =
N X
∗ ˆ i,j D (x, y) +
i=1
N X
∗ Ni,j (x, y)
(6)
i=1
Let us follow an intuitive argumentation to show that the sum of detail coefficients separates the signal from the noise, an effect that increases with N . Since the noise between different bands is supposed to be uncorrelated, the second term in eqn. 6 averages the noise. The √obtained noise pdf has a standard deviation that grows ∼ N . The signals between different bands in multivalued images are usually highly correlated, thus the standard deviation of the pdf of the first term grows ∼ N . To show that the latter holds, even for small correlations, consider N random variables drawn from the same pdf with the same standard deviation σ, and a small but equal correlation ρ between them. The standard deviation of the sum of two such variables, up to first p √ order in ρ, is given by: σ[xi + xj ] = σ 2(1 + ρ) ' 2σ(1 + ρ/2). For large N , we find: PN −1 √ ! √ PN i i=1 √ σ[ n=1 xn ] ' N σ + ρσ N √ 2 ' N σ + ρ( N σ) (7) 3 This shows that even for small correlations, the standard deviation of the sum of N random variables from any pdf is proportional to N . C. Products of detail coefficients The argument that summing over bands separates signal from noise also holds for specific functions of the detail coefficients. Three different products are investigated in this paper: the square of the detail coefficients, products of coefficients at adjacent scales and products of coefficients at different bands. Each of these products separates signal from noise in a specific way. C.1 square of detail coefficients Let us square the detail coefficients. Since the noise is not correlated with the signal: ∗ 2 ˆ ∗ )2 (x, y) + (N ∗ )2 (x, y) (8) ) (x, y) ' (D (Di,j i,j i,j The square of the coefficients separates signal from noise since real edge coefficients tend to become larger for coarser scales, while noise coefficients become smaller. By squaring, this effect becomes more prominent. Let us calculate the noise pdf for this square of detail coefficients. In general, the pdf of the square of a random √ variable x with pdf fx is given by: fy (y = x2 ) = √1y fx ( y), The noise pdf of the square of the wavelet detail coefficients is therefore given by ! 1 y 1 ∗ 2 fy (y = (Ni,j ) ) = √ (9) √ exp − 2 2σj 2πσj y ∗ 2 This pdf has an expected value of E[(Ni,j ) ] = σj2 and a ∗ 2 4 variance of V [(Ni,j ) ] = 2σj .
3
of 0:
C.2 interscale products In [14], it was proposed to denoise based on the signal correlation of wavelet coefficients between adjacent scales. Using this, it was suggested to use the product of two ∗ ∗ )(x, y) wavelet coefficients at adjacent scales (Di,j · Di,j+1 to separate signal from noise. Also for this product holds: ∗ ∗ )(x, y) ' · Di,j+1 (Di,j ∗ ∗ ∗ ∗ ˆ ˆ )(x, y) (10) · Ni,j+1 )(x, y) + (Ni,j (Di,j · Di,j+1
Since real edge coefficients are highly correlated, the signal product will be large. For noise coefficients, the correlation between adjacent scales is low, and therefore the product will be small. To calculate the noise pdf of this product, first notice that, due to the non-orthogonality of the proposed wavelet frames, the noise between adjacent scales is correlated. The ∗ ∗ joint pdf of Ni,j and Ni,j+1 is given by: ∗ ∗ , x2 = Ni,j+1 )= fx1 ,x2 (x1 = Ni,j −x2 σ2 −x2 σ2 +2ρ j,j+1 x1 ·x2 σj σj+1 1 j+1 2 j 1 2βσj σj+1 p e (11) 2π σj σj+1 β
with β = σj σj+1 (1 − ρ2j,j+1 ) and the correlation between scales j and j + 1 is given by: RR ∗ ∗ (x, y)dxdy ψj (x, y)ψj+1 ρj,j+1 = qR R (12) R R ∗ (x, y))2 dxdy (ψj∗ (x, y))2 dxdy (ψj+1
∗ ∗ The pdf fy (y = Ni,j · Ni,j+1 ) can be obtained by using the general property that the pdf of the product of 2 correlated random variables x1 and x2R, having joint pdf f (x1 , x2 ), is given by: fy (y = x1 · x2 ) = f (y/x2 , x2 )dx2 /|x2 | [20]. We find: ∗ ∗ )= · Ni,j+1 fy (y = Ni,j |y| 1 yρj,j+1 p K0 exp β β π σj σj+1 β
(13)
C.3 interband products Finally, the correlation of the signal between different bands can be exploited for separating noise from signal. For this, the product of wavelet coefficients at two different bands i and k is calculated. Also for this product, one finds:
·
∗ Nk,j )
1 = K0 πσj2
|y| σj2
!
(15)
∗ ∗ This pdf has an expected value of E[Ni,j · Nk,j ] = 0 and a ∗ ∗ 4 variance of V [Ni,j · Nk,j ] = σj .
D. Sums of products The products that were defined in the previous section separate the noise from the signal pdf. Now, the argumentation of section II-B can be applied: by summing the products over several bands the signal pdf can be further separated from the noise pdf. To calculate the noise pdfs from the sums of products, we use the general property that the pdf of the sum of 2 independent random variables y1 and y2 with pdfs fy1 is given by the convolution: fz (z = and fy2 respectively, R∞ y1 + y2 ) = −∞ fy1 (z − y2 )fy2 (y2 )dy2 [20]. D.1 sum of squared coefficients
For the sum of squared coefficients, one obtains after calculation: ! N N/2−1 X z −z ∗ 2 (Ni,j ) )= exp fz(N ) (z = (16) 2σj2 (2σj2 )N/2 Γ( N2 ) i=1
R∞ where Γ(x) = 0 tx−1 exp (−t)dt is the Gamma function. The pdf of eqn. 16 is called the Γ-pdf Γ(N, σj ). Its exPN ∗ 2 pected value is: E[ i=1 (Ni,j ) )] = N σj2 and its variance PN ∗ 2 is: V [ i=1 (Ni,j ) )] = 2N σj4 . For σj = 1 it is the χ2 -pdf with N degrees of freedom. For large values of N , the pdf becomes the gaussian g(N σj2 , 2N σj4 ). D.2 sum of interscale products For the sum of products of wavelet coefficients at adjacent scales, one obtains after calculation: zρj,j+1 × · = exp = β i=1 z (N −1)/2 |z| (17) K(N −1)/2 √ N N/2 β πβΓ( 2 )(σj σj+1 )
fz(N ) (z
where K0 (·) is the 0-th order modified Bessel function of ∗ the second kind. This pdf has an expected value of E[Ni,j · ∗ ∗ ∗ Ni,j+1 ] = ρj,j+1 σj σj+1 and a variance of V [Ni,j · Ni,j+1 ] = 2 (1 + 2ρ2j,j+1 )σj2 σj+1 .
∗ ∗ )(x, y) ' · Dk,j (Di,j ˆ ∗ )(x, y) + (N ∗ · N ∗ )(x, y) ˆ∗ · D (D k,j i,j k,j i,j
fy (y =
∗ Ni,j
N X
∗ Ni,j
∗ Ni,j+1 )
where Kn (·) is the n-th order modified Bessel function of the second kind. Since it is only defined for integer n, the number of sums N should be chosen The obPN odd. ∗ ∗ tained pdf has an expected value of E[ i=1 Ni,j · Ni,j+1 ]= PN ∗ ∗ N ρj,j+1 σj σj+1 and a variance of V [ i=1 Ni,j · Ni,j+1 ] = 2 N (1 + 2ρ2j,j+1 )σj2 σj+1 . D.3 sum of interband products
(14)
Since the noise between the bands is supposed to be decorrelated, the noise pdf of the product is similar to equation 13, with equal standard deviations and a correlation
Finally for the sum of products of wavelet coefficients at different bands, one has to be careful to sum independent variables only. Each band may appear only once in the sum. Therefore, the number of summations can only be B = N/2 for N even, or B = (N − 1)/2 for N odd. Of
4
course, one has the freedom to choose which bands i and k are to be multiplied. In order to take optimally advantage of the interband correlation, the following procedure is adopted. The correlation of all possible pairs is calculated and ranked. The pair with highest correlation is chosen first. Then each time, the subsequent highest correlation pairs are added, on the condition that each band of the pair has not been included before. In this way, exactly B summations are performed. For the noise pdf of this sum, one obtains: fz(B) (z =
B X i,k
∗ ∗ Ni,j · Nk,j )=
z (B−1)/2 K(B−1)/2 √ πΓ( B2 )(σj )B+1
|z| σj2
!
(18)
PB where the notation i,k denotes the sum over the pairs, as chosen in the proposed approach. Since Kn (·) is only defined for integer n, an additional limitation is that B should odd. The obtained pdf has an expected Pbe PB value B ∗ ∗ ∗ of E[ i,k Ni,j · Nk,j ] = 0 and a variance of V [ i,k Ni,j · ∗ 4 Nk,j ] = N σj . E. Thresholding From the calculated pdfs a thresholding procedure can be developed that removes most of the contribution of the noise. For this, the same argument as in [11] is applied, i.e. most of the noise contribution will be removed when most of the area of the noise pdf is within the threshold value. In [11], for a gaussian noise pdf, 3 to 4 times the noise standard deviation was chosen. In this work, the α(N ) percentile value of the obtained pdfs, f(α )(x) is chosen as a threshold value, with α close to one. In this way, three different thresholds can be defined: τ1
=
(N ) f(α) (
N X
∗ 2 ) ) (Ni,j
(19)
∗ ∗ ) · Ni,j+1 Ni,j
(20)
∗ ∗ Ni,j · Nk,j )
(21)
i=1
τ2
=
(N )
f(α) (
N X i=1
τ3
=
(N )
f(α) (
B X ik
The thresholds are calculated by numerical integration up till α of the pdfs, obtained from the analytical expressions of equations 16, 17 and 18. The thresholds are not applied on the sum of product quantities, but on the original band detail coefficients. For this, the following multivariate thresholding procedure is defined: at all positions where the sum of product is smaller than the obtained threshold value, the detail coefficients of all bands involved are removed. This procedure is expected to outperform single band thresholding, since signal and noise pdfs of the proposed sums of products are more separated than the pdfs of the original coefficients. If at a
particular position a band occurs with a signal coefficient, it is likely that it is accompanied by other signal coefficients at other bands at the same position and it will be retained after thresholding. On the other hand, a noise coefficient at one band is likely to be accompanied by other noise coefficients at different bands at the same position and it will be removed after thresholding. When for a pixel one or a few bands contain noise coefficients while the other bands are signal coefficients, the noise will not be removed, which results in a lower denoising performance. In a rare case, a pixel contains only one or a few signal coefficients. If these are not too large, they will be removed, which results in a signal deterioration. As will be shown in the experimental section, the latter two disadvantages are largely compensated by the advantages, resulting in an overall improved denoising performance. In the case of the sum of squared coefficients, only positive values occur. The thresholding procedure becomes: for all bands k = 1, · · · N : ∗ ∗ y) Dk,j (x, y) = Dk,j (x,P N ∗ 2 (22) if i=1 (Di,j ) (x, y) > τ1 ∗ Dk,j (x, y) = 0 otherwise In the case of the sum of interscale products, the obtained pdf is asymmetric around zero. The signal sum of products can be negative, revealing negative correlations (this can also be seen from figure 2c and 2d in the next section). To include these negative correlations, a second 0 (N ) PN ∗ ∗ · Ni,j+1 ) is necessary. threshold value τ2 = f(α0 ) ( i=1 Ni,j 0
α would then be close to zero. The thresholding procedure becomes: for all bands k = 1 · · · N : ∗ Dk,j (x, y) = 0 PN 0 ∗ ∗ if τ2 < i=1 (Di,j · Di,j+1 )(x, y) < τ2 (23) ∗ ∗ Dk,j (x, y) = Dk,j (x, y) otherwise
In the case of the sum of interband products, the obtained pdf is symmetric around zero. Here, only positive signal correlations are looked for (the largest positive interband correlations are summed). The thresholding procedure becomes: for all bands k = 1, · · · N : ∗ ∗ y) Dk,j (x, y) = Dk,j (x,P B ∗ ∗ (24) if i,k (Di,j · Dk,j )(x, y) > τ3 ∗ Dk,j (x, y) = 0 otherwise The question remains how to determine an optimal value for α. In this work, every effort was made to separate noise and signal pdfs as much as possible. Therefore, it is expected that an optimal value of α will not depend much on the signal, so that we can fix α to a value close to 1. In [12], a signal dependent threshold was derived for a gaussian noise pdf on wavelet coefficients. For this, a function, equivalent to the MSE was derived. If applied to the redundant wavelet representation of this work, the function at scale j is given by: ! Rt 2 x gx (0, σj2 )dx −t 2 Fj (t) = 2 σj − k R t − V [D˜j∗ ] (25) 2 )dx g (0, σ x j −t
5
with k the proportion of coefficients that has been removed image, including signal and noise as ’x’. From the plots, after applying the threshold t, and V [D˜j∗ ] the variance of one can verify that: the pdf of the coefficients after thresholding. This func- • the simulated noise pdfs fits correctly with the analytical tion has a signal dependent minimum that agrees with the pdfs. minimal MSE at a value of t which is then chosen as the • the signal pdfs of the products are broader than the noise optimal threshold value. The calculations to obtain this pdfs. This effect is limited, but seems to be highest for the function can be easily extended towards noise pdfs, other interscale products. than the normal pdf. For the symmetric sum of interband • when summing the products, the signal pdfs become products pdf, we obtain: much broader than the noise pdfs. This holds for all three products, but the effect seems to be highest for the interZ t B X scale products. k 4 2 ∗ ∗ Fj (t) = 2 Bσj − x fx (x = Ni,j · Nk,j )dx α(t) −t In order to quantify to what extent the proposed prodi,k ucts and sums separate signal from noise, the signal-tog B X noise ratio’s can be calculated. A definition of the signal∗ · D∗ ] (26) to-noise ratio of a signal is given by the ratio between the Di,j −V [ k,j i,k signal variance and the noise variance. Because of equations 8, 10 and 14, we find: Rt PB ∗ ∗ with α(t) = −t fx (x = i,k Ni,j · Nk,j )dx. As an example, s in figure 1 this function is plotted for a Landsat image with σ 2 (observed signal) − σ 2 (noise) simulated noise at scale 1, along with the MSE, in function SN R = (27) σ 2 (noise) of α(t). As can be seen, the two minima appear at approximately the same value for α, which is close to 0.98. Similar results were obtained for different scales, different images where the noise standard deviations are obtained from the ∗ and different noise standard deviation values. Values for analytical expressions. In table 1, the SNR’s for D1,j , P N ∗ the optimal value of α were all in the range 0.96-0.99. For D and the three products and their sums are shown i=1 i,j the sum of squares pdf a similar expression for the function in the case of the Landsat image for the first four scales. F can be obtained, and similar results were obtained. For From the table, one can observe that: the sum of interscale products, the asymmetric pdf poses a difficulty. The function F depends on two thresholds, • for all quantities the SNR grows with the scale. The and an optimum needs to be found in 2D. Results did not jump is largest between the first two scales, and from then consistently reveal a minimum in the expected range. In on it is roughly a factor of 2 for the wavelet coefficients fact, one can expect lesser denoising performance in this and a factor of 4 for the products. This effect is due to the case, due to the asymmetry. One of the consequences is observation that the noise standard deviation shrinks with e.g. that negative signal contributions can be corrupted by a factor of 2 between 2 consecutive scales, while the signal variance roughly remains constants between scales. higher-valued positive noise contributions. • For all three products the SNR is improved compared The time complexity of finding the signal dependent α to the wavelet coefficients. The improvement of the SNR is quite high, since the function F needs to be calculated compared to the squared coefficients is a factor of 2 for the for different values of α(t), with each time needs numeriinterscale products and a factor of 1.5 for the interband cal integrations. Since the optimal α seems to be rather products. This demonstrates that exploiting correlations signal independent, a fixed value for α can be chosen, e.g. between scales or between bands separates better the signal α = 0.98. This requires only one numerical integration of from the noise. the specific sum of product function, after which it can be • When summing over the bands, a further improvement applied to any image. of the SNR is obtained for all products, with roughly a factor of 2. The improvement with the sums over interband III. Experiments and discussion products is somewhat smaller, because of the smaller numFor a qualitative appreciation of the proposed tech- ber of summations. This demonstrates that exploiting the niques, the obtained pdfs from the proposed quantities are correlations between the bands by summing, improves the plotted in the case of a 7-band multispectral Landsat im- separation of signal and noise. age. On each band of the image independent additive gausTo demonstrate the denoising properties of the proposed sian noise is added with σ = 20. In total, 6 plots are shown techniques, the three thresholding procedures 22, 23 and 24 at the particularP scale j = 2: in figure 2a and b the plots are compared to each other and to single band thresholding. N ∗ 2 ∗ 2 ) ; in figure 2c and d the plots for (D1,2 ) and i=1 (Di,2 For the latter, 3 approaches are obtained by making use of PN ∗ ∗ ∗ ∗ and i=1 Di,2 · Di,3 and in figure 2e and the 3 proposed products, without summing. For all three · D1,3 for D1,2 PB ∗ ∗ ∗ ∗ f the plots for D1,2 · D2,2 and i,k Di,2 · Dk,2 . On each products, a threshold can be defined by the α-percentile plot, 3 pdfs are shown: the obtained analytical noise pdf value of the product pdfs, obtained from integration up till as a solid line, the noise pdf for a multivalued empty image α of the analytical expressions of equations 9, 13 and 15. 3 with simulated noise as ’+’ and the pdf of the observed different single band thresholding procedures are obtained:
6
For band i: ∗ ∗ Di,j = Di,j ∗ Di,j = 0 or ∗ Di,j = 0
∗ Di,j
or
=
∗ Di,j
∗ ∗ Di,j = Di,j ∗ Di,j = 0
∗ 2 ∗ 2 if (Di,j ) > f(α) (Ni,j ) otherwise
(28)
∗ ∗ if f(1−α) (Ni,j · Ni,j+1 ) ∗ ∗ ∗ ∗ < Di,j · Di,j+1 < f(α) (Ni,j · Ni,j+1 ) (29) otherwise
∗ ∗ ∗ ∗ if Di,j · Dk,j > f(α) (Ni,j · Nk,j ) otherwise
(30)
respectively, where in the last case band k is the one that is largest correlated with band i. Also, for comparison two simple standard wavelet denoising procedures were applied. The LMMSE technique estimates the signal wavelet coefficient of band i and scale j as: ˆ∗ ) σ ˆ 2 (D i,j ∗ ˆ i,j D (x, y) = 2 D∗ (x, y) ˆ ∗ ) i,j σj + σ ˆ 2 (D i,j
a standard technique, i.e by measuring the variance in a local region in the background, and found to be about 13. In figure 5a, one original band is shown. In figure 5b, the denoised band is shown, using the sum of interband product threshold. The same conclusions as for the Landsat image hold. IV. Conclusions In this paper, different denoising procedures for multivalued images were proposed, based on wavelet thresholding. The procedures take advantage of the interscale and interband correlations of the signals, by looking at interband sums of 3 specific products of wavelet coefficients. For these sums, the specific noise pdfs were analytically derived and from these pdfs, threshold values were obtained. From the experiments, we can conclude that all the multi band thresholding procedures outperform the standard single band wavelet thresholding procedure. The sum of interband products leads to the best thresholding procedure. References
(31)
ˆ ∗ ) an estimate for the signal variance. In [13], a with σ ˆ 2 (D i,j Maximum a Posteriori estimation was done for the signal ˆ∗ ) = variance. We will apply a simple estimation: σ ˆ 2 (D i,j P P ∗2 2 ∗ 2 2 ∗ σ (Di,j ) − σj , with σ (Di,j ) = x y Di,j (x, y), the observed signal variance. The wavelet thresholding technique of [11] will also be applied. The coefficients of band i at scale j are hard thresholded with a threshold value of t = 3.5σj . This technique will be referred to as ’PAN’. The proposed denoising techniques are applied to the same Landsat image as before. On the original image, Gaussian noise with different values of σ was added. The α-percentile values of the single band and multiband pdfs are calculated, with α = 0.98 in all cases. For the single band techniques, the results on the first band are given. In table 2, MSE values are given. In fig. 3, the original and noisy band 1 is shown. In figure 4, the denoised results are shown. From the table and the figures, one can conclude that: • the sum of squares, the interband and the sum of interband product thresholding results are superior to LMMSE and single band thresholding. The sum of interband product thresholding leads to the best results. • the interscale product and sum of interscale product thresholding perform poor. This does not agree with the qualitative results on the SNR. The reason can be found in the negative signal correlations and the asymmetric noise 0 pdf. Playing around with the values of α and α to optimize the resulting MSE never lead to improved results. Finally, an experiment is conducted on a noisy multimodal MRI image, diffusion- weighted at 6 independent imaging gradient directions. Noise is mostly thermal and is known to be approximately independent normal distributed. The noise standard deviation is estimated using
[1] [2] [3] [4] [5]
[6] [7] [8]
[9]
[10]
[11] [12] [13]
[14]
[15]
I.L. Thomas, V.M. Benning, and N.P. Ching, Classification of remotely sensed images, Adam Hilger, Bristol, 1987. C Lee and D.A. Landgrebe, “Analyzing high-dimensional multispectral data,” IEEE Trans. Geosci. Remote Sensing, vol. 31, no. 4, pp. 388–400, 1993. T. Taxt and A. Lundervold, “Multispectral analysis of the brain using magnetic resonance imaging,” IEEE Trans. Med. Imaging, vol. 13, no. 3, pp. 470–481, 1994. C. Busch, “Wavelet based texture segmentation of multi-modal tomographic images,” Comput. & Graphics, vol. 21, no. 3, pp. 347–358, 1997. 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 Transactions on Geoscience and Remote Sensing, vol. 26:(1), pp. 65–74, 1988. J.L. Starck and P. Querre, “Multispectral data restoration by the wavelet karhunen-loeve transform,” Signal Processing, vol. 81:(12), pp. 2449–2459, 2001. D.L. Donoho and I.M. Johnstone, “Ideal spatial adaptation via wavelet shrinkage,” Biometrika, vol. 81, pp. 425–455, 1994. S. Mallat, “A theory for multiresolution signal decomposition : the wavelet representation,” IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 11, no. 7, pp. 674–693, 1989. S. Grace Chang, B. Yu, and M. Vetterli, “Adaptive wavelet thresholding for image denoising and compression,” IEEE Transactions on Image Processing, vol. 9:(9), pp. 1532–1546, 2000. S. Grace Chang, B. Yu, and M. Vetterli, “Spatially adaptive wavelet thresholding with context modeling for image denoising,” IEEE Transactions on Image Processing, vol. 9:(9), pp. 1522–1531, 2000. Q. Pan, L. Zhang, and G. Dai aand H. Zhang, “Two denoising methods by wavelet transform,” IEEE Transactions on Signal Processing, vol. 47:(12), pp. 3401–3406, 1999. L. Zhang, P. Bao, and Q. Pan, “Threshold analysis in waveletbased denoising,” Electronic Letters, vol. 37:(24), pp. 1485–1486, 2001. M.K. Mihcak, I. Kozintsev, K. Ramachandran, and P. Moulin, “Low-complexity image denoising based on statistical modeling of wavelet coefficients,” IEEE Signal Processing Letters, vol. 6:(12), pp. 300–303, 1999. Y. Xu, J.B. Weaver, D.M. Healy, and J. Lu, “Wavelet transform domain filtering: a spatially selective noise filtration technique,” IEEE Transactions on Image Processing, vol. 3:(6), pp. 747–758, 1994. J. Liu and P. Moulin, “Information-theoretic analysis of interscale and intrascale dependencies between image wavelet coefficients,” IEEE Transactions on Image Processing, vol. 10:(11), pp. 1647–1658, 2001.
7
1100
1050
1000
950
900 0.8
0.82
0.84
0.86
0.88
0.9
0.92
0.94
0.96
0.98
1
Fig. 1. The functions F (t) and MSE(t) for the first scale of a 7-band Landsat image, in the case of the sum of interband products. On the x-axis, α(t) is plotted, on the y-axis, MSE(t)as ’+’. The function F (t), plotted as ’x’ is rescaled to fit in the plot.
Fig. 2.
∗ )2 and Plots of the pdfs of: a,b: (D1,2
∗ · D ∗ and D1,2 1,3 ∗ . Dk,2
PN
i=1
PN
i=1
Fig. 3. a: orginal and b: noisy band (σ = 20) of multispectral Landsat image.
∗ )2 ; c, d: (Di,2
∗ · D ∗ ; e, f: D ∗ · D ∗ and Di,2 1,2 2,2 i,3
PB
i,k
∗ · Di,2
[16] B.M. Sadler and A. Swami, “Analysis of multiscale products for step detection and estimation,” IEEE Transactions on Information Theory, vol. 45:(4), pp. 1043–1051, 1999. [17] L. Zhang, P.Bao, and X. Wu, “Hybrid inter- and intra- wavelet scale image restoration,” Pattern Recognition, vol. 36:(8), pp. 1737–1746, 2003. [18] S. Mallat and S. Zhong, “Characterization of signals from multiscale edges,” IEEE Trans. Pattern Anal. Machine Intell., vol. 14, pp. 710–732, 1992. [19] S. Mallat, Ed., A wavelet tour of signal processing, vol. 2-nd edition, Academic Press, 1999. [20] A. Papoulis, Probability, random variables and stochastic processes, Mc. Graw-Hill, N.Y., 1991.
Fig. 4. Denoised band of multispectral Landsat image; a: thresholding squared coefficients; b: thresholding sum of squared coefficients; c: thresholding interscale products; d: thresholding sum of interscale products; e: thresholding interband products; f: thresholding sum of interband products.
8
j 1 2 3 4 j 1 2 3 4
PN
∗ D1,j 0.95 3.04 6.54 13.4
i=1
∗ Di,j 1.91 6.65 15.6 30.3
PN
∗ 2 (D1,j ) 1.60 14.0 68.2 243
∗ 2 i=1 (Di,j )
3.02 28.1 138 481
PN
∗ ∗ D1,j · D1,j+1 4.44 31.1 130 428
i=1
∗ ∗ Di,j · Di,j+1 8.40 62.5 264 852
PB
i,k
∗ ∗ D1,j · D2,j 2.06 18.3 84.5 304 ∗ ∗ Di,j · Dk,j 2.80 25.7 127 434
TABLE I SNR values for the products and sums of products in the case of a 7-band Landsat image.
σ 10 15 20 25 30
LMMSE 87 167 250 326 395
PAN 108 200 285 386 471
PN
∗ 2 i=1 (Di,j )
83 149 217 283 348
∗ ∗ D1,j · D1,j+1 106 203 284 396 496
PN
i=1
∗ ∗ Di,j · Di,j+1 100 176 253 330 404
∗ ∗ D1,j · D2,j 85 148 236 292 368
PB
i,k
∗ ∗ Di,j · Dk,j 82 142 205 269 333
TABLE II Obtained MSE values after denoising, applying the proposed thresholding techniques on a 7-band Landsat image.
Fig. 5. a: Original band of a multimodal MRI image; b: thresholding sum of interband products.
