EURASIP Journal on Applied Signal Processing 2005:8, 1205–1211 c ○ 2005 Hindawi Publishing Corporation A Multivariate Thresholding Technique for Image

Multiwavelets, wavelets with several scaling functions, offer simultaneous orthogonality, symmetry, and short support, which is not possible with ordinary (scalar) wavelets. These properties make multiwavelets promising for signal processing applications, such as image denoising. The common approach for image denoising is to get the multiwavelet decomposition of a noisy image and apply a common threshold to each coefficient separately. This approach does not generally give sufficient performance. In this paper, we propose a multivariate thresholding technique for image denoising with multiwavelets. The proposed technique is based on the idea of restoring the spatial dependence of the pixels of the noisy image that has undergone a multiwavelet decomposition. Coefficients with high correlation are regarded as elements of a vector and are subject to a common thresholding operation. Simulations with several multiwavelets illustrate that the proposed technique results in a better performance.


INTRODUCTION
Multiwavelets are a relatively new addition to the wavelet theory, and have received considerable attention since their introduction [1,2,3,4,5,6,7,8,9,10,11,12,13,14]. Contrary to ordinary wavelets, multiwavelets offer simultaneous orthogonality, symmetry, and short support. Similar to performing wavelet decomposition with filters, multiwavelet decomposition can be realized with filterbanks. The filter coefficients in this case are, however, matrices instead of scalars. Therefore, two or more input streams to the multiwavelet filterbank are required to perform the decomposition. Several methods have been developed for obtaining multiple input streams from a given single input stream [2,3,4,15,16,17].
One of the widely used applications of wavelet decomposition is the removal of additive white Gaussian noise from noisy signals [18,19]. The discrimination between the actual signal and noise is achieved by choosing an orthogonal basis, which efficiently approximates the signal with few nonzero coefficients. A signal enhancement can then be obtained by discarding components below a predetermined threshold value. Although the performance of multiwavelets has been evaluated for image compression and coding (see, e.g., [20] and the references therein), less work about denoising applications of multiwavelets exists [21,22,23]. Most of the existing work concentrates on denoising of one-dimensional signals. A detailed discussion of multiwavelets and their applications to signal and image processing can be found in [5,24]. In these works, the noisy image is firstly decomposed with Geronimo, Hardin, and Massopust (GHM) [6] multiwavelets, then each individual coefficient is thresholded with the universal threshold [18]. This is similar to the approach that is widely used for image denoising with scalar wavelet decomposition. The universal threshold is calculated as λ = 2σ 2 log n for a length-n signal that is corrupted with additive white Gaussian noise with zero mean and variance σ 2 . This thresholding technique assumes that the noise on each coefficient is independent. However, this assumption may not be always valid for the coefficients of a multiwavelet decomposition. Therefore, a multivariate thresholding scheme for one-dimensional signals using multiwavelets was introduced in [17]. Instead of thresholding individual multiwavelet coefficients, coefficient vectors are considered and thresholding operation is applied to the whole vector. In this paper, we design a multivariate thresholding technique designed specifically for image denoising. Simulations with various multiwavelets and preprocessing methods reveal that the new technique performs better than term-byterm thresholding, or the univariate method of [5]. This paper is organized as follows. In Section 2, the multiwavelet theory is reviewed and some information on the implementation of multiwavelet decomposition with filterbanks is given. In Section 3, the proposed multivariate thresholding technique for image denoising is introduced. In Section 4, simulation results are presented and finally in Section 5 conclusions are drawn.

BACKGROUND
Multiwavelets are characterized with several scaling functions and associated wavelet functions. Let the scaling functions be denoted in vector form as is called the multiscaling function, T denotes the vector transpose and φ j (t) is the jth scaling function. Likewise, let the wavelets be denoted as Ψ(t) = [ψ 1 (t), ψ 2 (t), . . . , ψ L (t)] T , where ψ j (t) is the jth wavelet function. Then, the dilation and wavelet equations for multiwavelets take the following forms, respectively: The lowpass filter H and the highpass filter G are L × L matrix filters, instead of scalars. In theory, L could be as large as possible, but in practice it is usually chosen to be two. Mallat's pyramid algorithm [25] for single scaling and wavelet functions extends to the matrix version. The resulting twochannel, 2 × 2 matrix filterbank operates on two input data streams, filtering them into four output streams. Each of these streams is then downsampled by a factor of two. This procedure is illustrated in Figure 1. Because a given signal consists of a single stream but the filterbank needs two streams, a method of mapping the data into two streams has to be developed. This mapping process is called preprocessing and is performed by a prefilter [5,7,16]. The postfilter, on the other hand, maps the data from multiple streams into one stream again. In the design of prefilters, it is desirable that properties of multiwavelet bases such as orthogonality, approximation order, short support, and symmetry are preserved as far as possible. For this reason, research is going on for designing multiwavelet bases, called balanced multiwavelets, for which prefiltering can be avoided [26]. The type of the prefilter used in a specific application is important for the performance. Our simulations have revealed that the best image denoising results are obtained with the repeated row prefilter and the approximation prefilter [5].

THRESHOLDING WAVELET AND MULTIWAVELET COEFFICIENTS
Significant wavelet coefficients are extracted by thresholding as proposed by Donoho and Johnstone [18]. The two mostly used methods of thresholding are soft thresholding and hard thresholding. In hard thresholding, the coefficients below a threshold are set to zero; those which are above the threshold remain untouched. In soft thresholding, the coefficients below a threshold value are set to zero as well, but coefficients above the threshold value are shrunk. The amount of shrinking is equal to the threshold value.
The universal threshold is applied to each individual multiwavelet coefficient separately in [5]. However, the individual multiwavelet coefficients are not independent because using any prefilter other than the identity prefilter produces correlated coefficients. Taking this fact into account, Downie and Silverman [17] have proposed a multivariate thresholding method for one-dimensional signal denoising. When multiwavelet decomposition is applied to a one-dimensional signal after prefiltering, each resultant coefficient is represented by a length-L vector. Assuming L = 2, the coefficient is k is the coefficient index, and T is the vector transpose. The thresholding operation is applied to the whole vector coefficients. Using the same approach, we have developed a multivariate thresholding method for multiwavelets that is applicable to image denoising. Due to the special properties of two-dimensional multiwavelet decomposition, however, a different approach should be followed. This idea is going to be elaborated in the subsequent paragraphs.
During a single level of decomposition of an image using a scalar wavelet, the two-dimensional data is replaced with four blocks. These blocks correspond to the subbands that represent either lowpass filtering or highpass filtering in each direction. The procedure for wavelet decomposition consists of consecutive operations on rows and columns of the twodimensional data. The wavelet transform first performs one step of the transform on all rows. This process yields a matrix where the left side contains downsampled lowpass coefficients of each row, and the right side contains the highpass coefficients. Next, one step of decomposition is applied to all columns; this results in four types of coefficients.
(1) HH represents the diagonal features of the image and is formed by highpass filtering in both directions. (2) HL represents the horizontal features of the image and is formed by lowpass filtering the rows and then highpass filtering the columns. (3) LH represents the vertical features of the image and is formed by highpass filtering the rows and then lowpass filtering the columns. (4) LL represents the coefficients that will be further decomposed in the next step. It is formed by lowpass filtering both the rows and the columns.
The subbands corresponding to a single-level scalar wavelet decomposition are illustrated in Figure 2.  Figure 2: Subbands corresponding to a single-level wavelet decomposition. Figure 3: Subbands corresponding to a single-level multiwavelet decomposition.
The multiwavelet decomposition is similar to the wavelet decomposition, but has some differences. The multiwavelet filterbanks have two channels, so the decomposition will have two sets of scaling coefficients and two sets of wavelet coefficients. For multiwavelets, the L and H labels have subscripts denoting the channel to which the data corresponds. For example, the subband L 1 H 2 represents the data from the second channel highpass filtered in the horizontal direction, and the first channel lowpass filtered in the vertical direction.
It has been shown that there exists a spatial dependence between pixels in the different subbands of a wavelet decomposition [27]. This dependence is in the form of a parentchild relationship. Each parent pixel in a deeper decomposition level has four children in the upper level in the form of a 2 × 2 block of adjacent pixels. This idea has been extensively used in image coding because the statistics of the parent-child pixels are similar. If the parent coefficient has a small value, then the children would likely have small values. If the parent coefficient has a large value, the children might also have large values. This observation can also be used for a multivariate image denoising scheme for multiwavelets. If coefficients with high correlations are handled together and thresholding operation is applied to them simultaneously, we may expect to get better results. However, the parent-child relationship does not hold for multiwavelet decomposition [20]. This is due to the fact that, in a single-level decomposition, each of the three subbands LH, HL, and HH are further split into four smaller subbands. This operation destroys the parent-child relationship. The output structure of an image after a single-level multiwavelet decomposition is illustrated in Figure 3. Observations reveal that there is a large amount of similarity in each of the subbands. This suggests that the spatial dependence of the pixels could be restored. As an example, the Lena image after being exposed to a single-level multiwavelet decomposition is illustrated in Figure 4a and the HH subband of the same image consisting of four smaller subbands is illustrated in Figure 4b. The similarity between the illustrated subbands is observable from the figure.
For image denoising applications, the spatial dependence of the pixels should be restored so that coefficients with high correlation can undergo a common thresholding operation. The same observation has been used for designing a novel quantization scheme for image compression in [20]. To illustrate this point more clearly, the HH subband of Figure 3 is shown in Figure 5. The coefficients represented with dots in Figure 5 would have been placed next to each other if scalar wavelet decomposition had been applied. This observation suggests that vectors of length four can be formed by using a coefficient from each of the four subbands. These coefficients are chosen such that they have the same location in their respective subbands. For example, the four coefficients shown in Figure 5 would form a vector. Then, a multivariate thresholding operation is applied to each of these vectors. This procedure is repeated separately for all of the coefficients in all three subbands HH, HL, and LH. Figure 5: HH subband corresponding to a single-level multiwavelet decomposition.
Applying the multiwavelet transform with a prefilter to a noisy image, and then collecting the coefficients in each subband in vectors of length four, we get vector coefficients of the form [2] where υ is the noise-free multiwavelet coefficient vector, ρ is the multiwavelet coefficient vector of the noise, w represents the multiwavelet coefficient vector of the corrupted signal, j is the decomposition level, and k is the coefficient index. ρ j,k has a multivariate normal distribution N(0, Θ j ). Θ j is the covariance matrix for the noise term and depends on the resolution level j. We would like to whiten the noise so that each coefficient within the vector would be distributed independently. We assume that there is no signal component in the vector. Then, the whitening operation could be achieved by multiplying the noise vector with Θ −1/2 j . If y j,k = Θ −1/2 j w j,k , then it can be easily shown that the covariance matrix of y is the identity matrix. The squared length of the vector y can be computed by where the superscript T denotes the transpose. σ j,k is a positive scalar value and has a Chi-squared distribution with four degrees of freedom. The analogous hard thresholding and soft thresholding rules can be applied as in (4) and (5), respectively:ŵ Similar to the procedure developed in [17], this thresholding technique treats the highly correlated multiwavelet coefficients simultaneously and applies a common threshold to the vector of coefficients. The effect of correlation is compensated by a whitening transformation. The covariance matrix used for the transformation depends on the decomposition level and is calculated for the HH, HL, and LH subbands separately at each decomposition level. We also have to calculate the threshold value λ. We assume that there are N identically and independently distributed χ 2 4 random variables and the  maximum of these random variables is denoted by M. The threshold value λ is the infimum of all sequences λ N such that As the number of pixels go to infinity, the probability of the threshold being greater than the maximum of the noise random variables approaches to one. This guarantees that, with high probability, a signal component exists in coefficients that are larger than the threshold value. The threshold value can be found by using the cdf of M in the limit problem (6) and solving for λ. An appropriate sequence satisfying the relation in (6) has been shown to be λ N = 2 log N + 2 log log N [28]. Therefore, this is the value that is used for thresholding the multiwavelet coefficient vectors.

SIMULATION RESULTS
The performance of the multivariate multiwavelet thresholding method that has been proposed in this paper is investigated with simulations. White Gaussian noise with σ = 25 is added to a 256×256 Lena image and denoising by soft thresholding with the proposed method is carried out with GHM [5], CL [8], and SA4 [29] multiwavelets. The same image and multiwavelets are then used for denoising with the univariate scheme. The performance of the scalar wavelet D4 [30] is also studied under the same conditions. The prefiltering methods used for the multiwavelet decompositions are the approximation prefiltering [5] and the repeated row prefiltering [5].
The simulation results can be evaluated objectively and subjectively. For objective evaluation, the signal to noise ratio (SNR) of each denoised image has been calculated. The SNR of the noisy Lena image before exposed to any denoising operation is 6.3793 dB. The SNR values of the denoised images are listed in Tables 1 and 2. The SNRs of the denoised images which are preprocessed with approximation prefiltering and repeated row prefiltering are given in Tables 1 and 2, respectively. Figures 6 to 8 illustrate some of the denoised images for subjective performance comparison. Some important observations can be made from the simulation results. The objective and subjective results prove that the proposed multivariate image denoising technique performs better than the univariate denoising method. The SNR values of the denoised images with the proposed technique are higher and the quality of the images is superior. The highest SNR of 12.6768 dB with approximation prefiltering is attained with the SA4 multiwavelet and the proposed technique. Similarly, the highest SNR of 13.2839 dB with repeated row prefiltering is attained with the CL multiwavelet and the proposed technique. In general, results for the repeated row prefiltering method seem to be better for both multivariate and univariate thresholding methods. It is also observed that the multiwavelets produce more satisfactory results than the ordinary scalar wavelets. The SNR of the denoised image with the scalar D4 wavelet is 10.0017 dB.

CONCLUSIONS
A multivariate vector-based thresholding technique has been introduced for multiwavelet based image denoising applications. The technique is based on the idea of restoring the spatial dependence of pixels of an image that has been subject to a multiwavelet decomposition. Four of such pixels are thought as elements of a vector and a multivariate thresholding scheme is developed for the whole vector. Simulations are carried to evaluate the performance of the proposed technique. Results prove that the new technique is well oriented for image denoising applications and produces promising results both objectively and subjectively.
Ayşin Ertüzün was born in 1959 in Salihli, Turkey. She received the B.S. degree (with honors) from Bogaziçi University, Istanbul, Turkey, the M. Eng. degree from McMaster University, Hamilton, Ontario, Canada and the Ph.D. degree from Bogazici University, Istanbul, Turkey, all in electrical engineering, in 1981, 1984, and 1989, respectively. Since 1988, she is with the Department of Electrical and Electronic Engineering, Bogaziçi University, where she is currently an Associate Professor. Her current research interests are in the areas of blind signal processing, Bayesian methods, adaptive systems, and filtering with applications to communication systems, image processing, independent component analysis and its applications, and pattern recognition. She has authored and coauthored about 50 scientific papers in journals and conference proceedings. She is a Member of IEEE and IAPR.