**Resolutıon Enhancement Based Image Compression Technique using Singular Value Decomposition and Wavelet Transforms**

Gholamreza Anbarjafari, Pejman Rasti, Morteza Daneshmand and Cagri Ozcinar

Additional information is available at the end of the chapter

http://dx.doi.org/10.5772/61335

#### **Abstract**

In this chapter, we propose a new lossy image compression technique that uses singular value decomposition (SVD) and wavelet difference reduction (WDR) technique followed by resolution enhancement using discrete wavelet transform (DWT) and stationary wavelet transform (SWT). The input image is decomposed into four different frequency subbands by using DWT. The low-frequency sub‐ band is the being compressed by using DWR and in parallel the high-frequency subbands are being compressed by using SVD which reduces the rank by ignor‐ ing small singular values. The compression ratio is obtained by dividing the total number of bits required to represent the input image over the total bit numbers obtain by WDR and SVD. Reconstruction is carried out by using inverse of WDR to obtained low-frequency subband and reconstructing the high-frequency sub‐ bands by using matrix multiplications. The high-frequency subbands are being enhanced by incorporating the high-frequency subbands obtained by applying SWT on the reconstructed low-frequency subband. The reconstructed low-fre‐ quency subband and enhanced high-frequency subbands are being used to gener‐ ate the reconstructed image by using inverse DWT. The visual and quantitative experimental results of the proposed image compression technique are shown and also compared with those of the WDR with arithmetic coding technique and JPEG2000. From the results of the comparison, the proposed image compression technique outperforms the WDR-AC and JPEG2000 techniques.

**Keywords:** Lossy Image compression, Singular Value Decomposition, Wavelet Differ‐ ence Reduction, Stationary Wavelet Transform, Discrete Wavelet Transform, Image Super Resolution

© 2015 The Author(s). Licensee InTech. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### **1. Introduction**

With the growing demand for multimedia applications, especially high-definition images, efficient storage and transmission of images have been issues of great concern [1, 2, 3, 4]. Image processing deals with the reduction of the amount of bits used to represent an image. Not only that but also resolution of an image plays an important role in many image processing applications, such as video resolution enhancement [5], feature extraction [6], and satellite image resolution enhancement [7]. In general, there are two types of super resolution ap‐ proaches, multi-image super resolution and single image. Multiple-image super-resolution algorithms, like [8], [9], [10] to name a few, receive a couple of low-resolution images of the same scene as input and usually employ a registration algorithm to find the transformation between them. This transformation information is then used along with the estimated blurring parameters of the input low-resolution images, to combine them into a higher-scale framework to produce a super-resolved output image. For multiple-image super-resolution algorithms to work properly there should be subpixels displacements between input low-resolution images. Furthermore, these subpixels displacements should be estimated properly by the registration algorithm, which is usually a challenging task, especially when complicated motion of nonrigid objects, like human body, needs to be modeled. These algorithms are guaranteed to produce proper higher-resolution details; however, their improvement factors are usually limited by factors close to 2 [11].

Single-image super-resolution algorithms, like [12, 13, 14], to name a few, do not have the possibility of utilizing subpixel displacements, because they only have a single input. Instead, they employ a kind of training algorithm to learn the relationship between a set of highresolution images and their low-resolution counterparts. This learned relationship is then used to predict the missing high-resolution details of the input low-resolution images. Depending on the relationship between the training low- and high-resolution images, these algorithms can produce high-resolution images that are far better than their inputs, by improvement factors that are much larger than 2 [15]. Hence, compression of an image and yet reconstructing the image with good resolution is important. Information theory is playing an important role in image compression. Information theory can be used in order to reduce the dimensionality of data such as histogram [16, 17]

There are two categories of image compression techniques, namely lossless and lossy image compression techniques [18, 19]. In lossless image compression, the original image can be perfectly recovered from the compressed image while in lossy compression the original image cannot be perfectively recovered from the compressed image because some information is lost as a result of compression. Lossless compression is used in applications with high requirements such as medical imaging. Lossy compression techniques are very popular because they offer higher compression ratio. The objective of image compression is to achieve as much compres‐ sion as possible with little loss of information [20, 21].

Wavelets are also playing significant role in many image processing applications [12, 22, 23, 24]. The two-dimensional wavelet decomposition of an image is performed by applying the one-dimensional DWT along the rows of the image first, and then the results are decomposed along the columns. This operation results in four decomposed subband images referred to Low-Low (LL), Low-High (LH), High-Low (HL), and High-High (HH). The frequency components of those subbands cover the full frequency spectrum of the original image. Figure 1 shows different subband images of Lena image where the top-left image is the LL subband and the bottom-right image is the HH subband.

**1. Introduction**

36 Wavelet Transform and Some of Its Real-World Applications

limited by factors close to 2 [11].

of data such as histogram [16, 17]

sion as possible with little loss of information [20, 21].

With the growing demand for multimedia applications, especially high-definition images, efficient storage and transmission of images have been issues of great concern [1, 2, 3, 4]. Image processing deals with the reduction of the amount of bits used to represent an image. Not only that but also resolution of an image plays an important role in many image processing applications, such as video resolution enhancement [5], feature extraction [6], and satellite image resolution enhancement [7]. In general, there are two types of super resolution ap‐ proaches, multi-image super resolution and single image. Multiple-image super-resolution algorithms, like [8], [9], [10] to name a few, receive a couple of low-resolution images of the same scene as input and usually employ a registration algorithm to find the transformation between them. This transformation information is then used along with the estimated blurring parameters of the input low-resolution images, to combine them into a higher-scale framework to produce a super-resolved output image. For multiple-image super-resolution algorithms to work properly there should be subpixels displacements between input low-resolution images. Furthermore, these subpixels displacements should be estimated properly by the registration algorithm, which is usually a challenging task, especially when complicated motion of nonrigid objects, like human body, needs to be modeled. These algorithms are guaranteed to produce proper higher-resolution details; however, their improvement factors are usually

Single-image super-resolution algorithms, like [12, 13, 14], to name a few, do not have the possibility of utilizing subpixel displacements, because they only have a single input. Instead, they employ a kind of training algorithm to learn the relationship between a set of highresolution images and their low-resolution counterparts. This learned relationship is then used to predict the missing high-resolution details of the input low-resolution images. Depending on the relationship between the training low- and high-resolution images, these algorithms can produce high-resolution images that are far better than their inputs, by improvement factors that are much larger than 2 [15]. Hence, compression of an image and yet reconstructing the image with good resolution is important. Information theory is playing an important role in image compression. Information theory can be used in order to reduce the dimensionality

There are two categories of image compression techniques, namely lossless and lossy image compression techniques [18, 19]. In lossless image compression, the original image can be perfectly recovered from the compressed image while in lossy compression the original image cannot be perfectively recovered from the compressed image because some information is lost as a result of compression. Lossless compression is used in applications with high requirements such as medical imaging. Lossy compression techniques are very popular because they offer higher compression ratio. The objective of image compression is to achieve as much compres‐

Wavelets are also playing significant role in many image processing applications [12, 22, 23, 24]. The two-dimensional wavelet decomposition of an image is performed by applying the one-dimensional DWT along the rows of the image first, and then the results are decomposed

**Figure 1.** LL, LH, HL, and HH subbands of Lena image obtained by using DWT.

In this research work, a new lossy compression technique which employs singular value decomposition (SVD) and wavelet difference reduction (WDR) is presented. SVD is a lossy image compression technique which can be regarded as a quantization process where it reduces the physcovisual redundancies of the image [25, 26]. In order to enhance the resolution of the decompressed image, stationary wavelet transform (SWT) is used. WDR is one of the state-of-the-art techniques in image compression which uses wavelet transform. It is a lossy image compression technique which achieves compression by first taking the wavelet transform of the input image and then applying the difference reduction method on the transform values [27, 28, 29, 30].

Wavelet transform based techniques also play a significant role in many image processing applications, in particular in resolution enhancement, and recently, many novel resolution enhancement by using wavelet transforms have been proposed. Demirel and Anbarjafari [31] proposed an image resolution enhancement technique based on the input image and interpo‐ lation of the high-frequency subband images obtained by DWT. In their technique, an SWT technique is used in order to enhance the edges. Then, at the same time input image is decomposed into four frequency subbands image by using DWT. After that the input image, as well as the high-frequency subbands are interpolated. The high-frequency subbands of SWT are used to modify the estimated high-frequency subbands. Finally, inverse DWT (IDWT) is applied to combine all frequency subbands in order to generate a high-resolution image. Figure 2 shows the block diagram of the proposed method in [31].

**Interpolation with factor α/2**

**Figure 2.** Block diagram of the proposed method in [31].

The authors in [32] proposed a learning-based super-resolution algorithm. In their proposed algorithm, a multi-resolution wavelet approach was adopted to perform the synthesis of local high-frequency features. Two frequency subbands, LH and HL, were estimated based on wavelet frame in order to get a high-resolution image. The LH and HL frequency subbands were used to prepare their training sets. Then, they used the training set in order to estimate wavelet coefficients for both LH and HL frequency subbands. Finally, the IDWT was used in order to reconstruct a high-resolution image.

In [33], the authors used a complex wavelet-domain image resolution enhancement algorithm based on the estimation of wavelet coefficients. Their method uses a dual-tree complex wavelet transform (DT-CWT) in order to generate a high-resolution image. First, they estimate a set of wavelet coefficients from the DT-CWT decomposition of the rough estimation of the highresolution image. Then, the inverse DT-CWT is used to combine the wavelet coefficients and the low-resolution input image in order to reconstruct a high-resolution image. Figure 3 shows the block diagram of the proposed method in [33].

**Figure 3.** Block diagram of the proposed method in [33].

reduces the physcovisual redundancies of the image [25, 26]. In order to enhance the resolution of the decompressed image, stationary wavelet transform (SWT) is used. WDR is one of the state-of-the-art techniques in image compression which uses wavelet transform. It is a lossy image compression technique which achieves compression by first taking the wavelet transform of the input image and then applying the difference reduction method on the

Wavelet transform based techniques also play a significant role in many image processing applications, in particular in resolution enhancement, and recently, many novel resolution enhancement by using wavelet transforms have been proposed. Demirel and Anbarjafari [31] proposed an image resolution enhancement technique based on the input image and interpo‐ lation of the high-frequency subband images obtained by DWT. In their technique, an SWT technique is used in order to enhance the edges. Then, at the same time input image is decomposed into four frequency subbands image by using DWT. After that the input image, as well as the high-frequency subbands are interpolated. The high-frequency subbands of SWT are used to modify the estimated high-frequency subbands. Finally, inverse DWT (IDWT) is applied to combine all frequency subbands in order to generate a high-resolution image. Figure

+

+

+

**Estimated LH**

**Estimated HL Estimated HH**

**IDWT**

**Interpolation with factor α/2**

**High Resolution Image (αmxαn)**

transform values [27, 28, 29, 30].

38 Wavelet Transform and Some of Its Real-World Applications

**Low resolution image (mxn)**

2 shows the block diagram of the proposed method in [31].

**LL**

**Figure 2.** Block diagram of the proposed method in [31].

**LH**

**Interpolation with factor 2**

> **Interpolation with factor α/2**

The authors in [32] proposed a learning-based super-resolution algorithm. In their proposed algorithm, a multi-resolution wavelet approach was adopted to perform the synthesis of local

**Interpolation with factor 2**

**Interpolation with factor 2**

**HL**

**DWT**

**SWT**

**HH**

**LL**

**LH**

**HL**

**HH**

Patel and Joshi [34] proposed a new learning-based approach for super resolution using DWT. The novelty of their method lies in designing application-specific wavelet basis (filter coeffi‐ cients). First the filter coefficients and learning the high-frequency details in the wavelet domain is used to initial estimate of super-resolution image. Then, they used a sparsely based regularization framework, in which image there was degradation. Finally, the super-resolution image is estimated by the initial super-resolution estimate and the estimated wavelet filter coefficients. Their algorithm has some advantages such as avoiding the use of registered images while learning the initial estimate, use of sparsity prior to preserving neighborhood dependencies in super-resolution image and use of estimated wavelet filter coefficients to represent an optimal point spread function to model image acquisition process. Figure 4 illustrates the block diagram of the proposed method in [34].

In [35], similar to the proposed method in [30], the authors used wavelet domain in order to generate super-resolution image from a single low-resolution image. They proposed an intermediate stage with the aim of estimating high-frequency subbands. The intermediate stage consists of an edge preservation procedure and mutual interpolation between the input low-resolution image and the HF subband images. Sparse mixing weights are calculated over blocks of coefficients in an image, which provides a sparse signal representation in the low-

**Figure 4.** Block diagram of the proposed method in [34].

resolution image. Finally, they used IDWT to combine all frequency subbands in order to reconstruct a high-resolution image. The block diagram of their proposed method is shown in Fig. 5.

In [36], they proposed a learning-based approach for super-resolving an image captured at low spatial resolution. They used a low resolution and a database of high- and low-resolution images as inputs to the proposed method. First, they used DWT in order to obtain highfrequency details of database images. Then, an initial high-resolution image was decimated by using the high-frequency details. In their observation model, they modelled a lowresolution image as an aliased and noisy version of the corresponding high-resolution image and then the initial high-resolution and test image estimated the aliasing matrix entries. After that, the prior model for the super-resolved image was chosen as an Inhomogeneous Gaussian Markov random field (IGMRF) and the model parameters were estimated using the same initial high-resolution estimate. They used a maximum a posteriori (MAP) estimation in order to arrive at the cost function minimized using a simple gradient descent approach. Figure 6 shows the block diagram of the proposed method in [36].

**Figure 5.** Block diagram of the proposed method in [35].

**Figure 6.** Block diagram of the proposed method in [36].

resolution image. Finally, they used IDWT to combine all frequency subbands in order to reconstruct a high-resolution image. The block diagram of their proposed method is shown in

In [36], they proposed a learning-based approach for super-resolving an image captured at low spatial resolution. They used a low resolution and a database of high- and low-resolution images as inputs to the proposed method. First, they used DWT in order to obtain highfrequency details of database images. Then, an initial high-resolution image was decimated by using the high-frequency details. In their observation model, they modelled a lowresolution image as an aliased and noisy version of the corresponding high-resolution image and then the initial high-resolution and test image estimated the aliasing matrix entries. After that, the prior model for the super-resolved image was chosen as an Inhomogeneous Gaussian Markov random field (IGMRF) and the model parameters were estimated using the same initial high-resolution estimate. They used a maximum a posteriori (MAP) estimation in order to arrive at the cost function minimized using a simple gradient descent approach. Figure 6

shows the block diagram of the proposed method in [36].

**Figure 4.** Block diagram of the proposed method in [34].

40 Wavelet Transform and Some of Its Real-World Applications

Fig. 5.

In the proposed compression technique, the input image is firstly decomposed into its different frequency subbands by using 1 level DWT. The LL subband is then being compressed by using DWR and the high-frequency subbands, i.e., LH, HL, and HH, are being compressed by using SVD. The proposed technique has been tested on several well-known images such as, Lena, Peppers, Boat, and Airfield. The results of this technique have been compared with those of JPEG2000 and WDR with arithmetic coding techniques. The quantitative experimental results based on PSNR show that the proposed technique overcomes the aforementioned techniques. The SVD and WDR image compression techniques are discussed in the next section.

#### **2. Review of singular value decomposition and wavelet difference reduction**

#### **2.1. Singular value decomposition**

From a mathematical point of view, an image can be represented by a matrix, which consists of one or three layers in the case the image is grayscale or RGB, respectively. The results of the implementation of SVD on a grayscale image, which is represented by the singlelayer image *A,* are three matrices *U*, *Σ*, and *V*, where *U* and *V* are orthogonal, and *Σ* is a diagonal matrix containing the singular values of *A*. In what follows, the SVD procedure is briefly reviewed. The relation between the matrix *A*, and the decomposed components, *U*, *Σ*, and *V*, can be mathematically presented through the formulation provided in Eqn. (1), where the dimensions of all the matrices are shown, given that the dimensions of the matrix *A* has been *m×n* [37, 38, 39]:

$$A\_{m \times n} = \mathcal{U}\_{m \times m} \Sigma\_{m \times n} \left(V\_{n \times n}\right)^{\mathrm{T}} \tag{1}$$

Eqn. (2) shows how a matrix Σ¯ *<sup>p</sup>*×*<sup>q</sup>* with smaller dimensions *p* ≤*m* and *q* ≤*n* can be used to approximate the diagonal matrix with the dimensions *m* and *n*:

$$
\Sigma\_{m \times n} = \begin{bmatrix} \overline{\Sigma}\_{p \times q} & \mathbf{0} \\ \mathbf{0} & \ddots \end{bmatrix} \qquad \qquad p \le m \text{ and } q \le n \tag{2}
$$

Some columns of *U* and rows of *V* are then reduced in order to reconstruct the compressed image by multiplication. This is shown in Eqn. (3):

$$\boldsymbol{U}\_{m \times n} = \begin{bmatrix} \overline{\boldsymbol{U}}\_{m \times p} & \widetilde{\boldsymbol{U}}\_{m \times \left(m-p\right)} \end{bmatrix} \quad \text{and} \qquad \boldsymbol{V}\_{n \times n} = \begin{bmatrix} \overline{\boldsymbol{V}}\_{m \times q} & \widetilde{\boldsymbol{V}}\_{n \times \left(n-q\right)} \end{bmatrix} \tag{3}$$

The compressed image is then obtained as shown in Eqn. (4):

$$A\_{m \times u} = \overline{\mathcal{U}}\_{m \times p} \overline{\Sigma}\_{p \times q} \left( \overline{\mathcal{V}}\_{u \times q} \right)^{r} \tag{4}$$

Because the singular matrix has sorted singular values (in descending order), by using the physcovisual concept, ignoring low singular value will not significantly reduce the visual quality of the image. Figure 7 shows Lena's picture being reconstructed by using different amount of singular values. This characteristic that an image can be reconstructed by fewer

amounts of singular values makes SVD suitable for compression. Because after reconstruc‐ tion of the image the ignored singular values cannot be recovered, the compression by SVD is lossy [33]. suitable for compression. Because after reconstruction of the image the ignored singular values cannot

Because the singular matrix has sorted singular values (in descending order), by using the

physcovisual concept, ignoring low singular value will not significantly reduce the visual quality of the

This characteristic that an image can be reconstructed by fewer amounts of singular values makes SVD

*T*

*m p nq p q AU V m n* (4)

**Figure 7.** Lena's image of size 256x256 reconstructed by Eq. (4) (a) orıginal Lena image; (b) reconstructed using128 σ; (c) reconstructed using 64 σ; (d) reconstructed using32 σ.

#### **2.2. Wavelet difference reduction**

be recovered, the compression by SVD is lossy [33].

**2. Review of singular value decomposition and wavelet difference**

From a mathematical point of view, an image can be represented by a matrix, which consists of one or three layers in the case the image is grayscale or RGB, respectively. The results of the implementation of SVD on a grayscale image, which is represented by the singlelayer image *A,* are three matrices *U*, *Σ*, and *V*, where *U* and *V* are orthogonal, and *Σ* is a diagonal matrix containing the singular values of *A*. In what follows, the SVD procedure is briefly reviewed. The relation between the matrix *A*, and the decomposed components, *U*, *Σ*, and *V*, can be mathematically presented through the formulation provided in Eqn. (1), where the dimensions of all the matrices are shown, given that the dimensions of the

( )*<sup>T</sup>*

and 0

Some columns of *U* and rows of *V* are then reduced in order to reconstruct the compressed

é ùé ù

ê úê ú ´ ´ ë ûë û = = ´ - ´ - (3)

( )*<sup>T</sup>*

Because the singular matrix has sorted singular values (in descending order), by using the physcovisual concept, ignoring low singular value will not significantly reduce the visual quality of the image. Figure 7 shows Lena's picture being reconstructed by using different amount of singular values. This characteristic that an image can be reconstructed by fewer

£ £ <sup>S</sup> S = <sup>O</sup> (2)

*m n pm qn* ´

<sup>±</sup> ( ) <sup>±</sup> ( ) and *UU U m m*´ ´ *mp mp m n VV V n n nq nq*

*AU V mn mm mn nn* ´ ´´ ´ = S (1)

*<sup>p</sup>*×*<sup>q</sup>* with smaller dimensions *p* ≤*m* and *q* ≤*n* can be used to

*mp nq p q AU V m n* ´ ´ ´ ´ = S (4)

**reduction**

**2.1. Singular value decomposition**

42 Wavelet Transform and Some of Its Real-World Applications

matrix *A* has been *m×n* [37, 38, 39]:

Eqn. (2) shows how a matrix Σ¯

approximate the diagonal matrix with the dimensions *m* and *n*:

0 *p q*

The compressed image is then obtained as shown in Eqn. (4):

é ù ê ú ê ú ë û

´

image by multiplication. This is shown in Eqn. (3):

The WDR is a compression technique, which is based on the difference reduction method. The wavelet transform of the input image is first made; bit plane encoding is then applied to the transform values. The bit plane encoding procedure starts with the initialization stage, where a threshold *To* is chosen such that *To* is greater than all the transform values, and at least, one of the transform values has a magnitude of *To/2*. The next stage is the initialization, where the threshold *T = Tk-1* is updated to *T = Tk*, where *Tk = Tk-1/2*. New significant transform values (*w(i)*) which are satisfying *T ≤ | w(i) | ≤ 2T* are then identified at the significant pass stage. The transform values of these significant transform values are then encoded using the difference reduction method. At the significant pass stage, already quantized values (*wQ*) which satisfy *|wQ| ≥ 2T* are then refined in order to reduce error [27, 29, 30].

#### **3. The proposed lossy image compression technique**

The proposed image compression technique is a lossy compression technique. Firstly, the image is decomposed into its frequency subbands by using DWT. Among these subbands, LL subband is being compressed by using WDR. The high-frequency subband images are being compressed by using SVD. The number of singular values that are being used in order to reconstruct the high-frequency subbands can be reduced into 1, i.e., the highest singular value is enough to reconstruct the high-frequency subbands. If only one singular value is being used in order to reconstruct a matrix, this means that only one column of U and V matrices are being used. The qualitative loss is not psychovisually noticeable up to some point. In order to obtain the compression ratio of the proposed technique, the total number of bits required to represent the original image is divided by the total of number of bits which is obtained by adding the number of bit streams of WDR for LL and that of the SVD compression for LH, HL, and HH.

Decompression is carried out by taking the inverse WDR (IWDR) of the bit streams in order to reconstruct the LL subband and in parallel the matrix multiplications are conducted in order to reconstruct LH, HL, and HH subbands. Due to the losses by ignoring low-valued singular values, high-frequency subbands need to be enhanced. For this purpose, stationary wavelet transform (SWT) is applied to the LL subband image which results in new low- and highfrequency subbands. These high-frequency subbands will have the same direction as the ones obtained by DWT (e.g., horizontal, vertical, and diagonal), so they will be added to the respective ones reconstructed by matrix multiplications. Now, the LL subband image obtained by IWDR and the enhanced LH, HL, and HH subbands are combined by using inverse DWT (IDWT) in order to reconstruct the decompressed image. The enhancement of high-frequency subbands by using SWT results in more sharpened decompressed image. The block diagram of the proposed lossy image compression technique is shown in Fig. 8. The experimental qualitative and quantitative results are represented and discussed in the next section.

#### **4. Experimental results and discussion**

As it was mentioned in the Introduction, for comparison purposes, the proposed lossy image compression was tested on many benchmark images, namely, Lena, Pepper, Boats, Airfield,

Resolutıon Enhancement Based Image Compression Technique… http://dx.doi.org/10.5772/61335 45

**Figure 8.** Block diagram of the proposed blocked based lossy image compression technique.

transform values. The bit plane encoding procedure starts with the initialization stage, where a threshold *To* is chosen such that *To* is greater than all the transform values, and at least, one of the transform values has a magnitude of *To/2*. The next stage is the initialization, where the threshold *T = Tk-1* is updated to *T = Tk*, where *Tk = Tk-1/2*. New significant transform values (*w(i)*) which are satisfying *T ≤ | w(i) | ≤ 2T* are then identified at the significant pass stage. The transform values of these significant transform values are then encoded using the difference reduction method. At the significant pass stage, already quantized values (*wQ*) which satisfy

The proposed image compression technique is a lossy compression technique. Firstly, the image is decomposed into its frequency subbands by using DWT. Among these subbands, LL subband is being compressed by using WDR. The high-frequency subband images are being compressed by using SVD. The number of singular values that are being used in order to reconstruct the high-frequency subbands can be reduced into 1, i.e., the highest singular value is enough to reconstruct the high-frequency subbands. If only one singular value is being used in order to reconstruct a matrix, this means that only one column of U and V matrices are being used. The qualitative loss is not psychovisually noticeable up to some point. In order to obtain the compression ratio of the proposed technique, the total number of bits required to represent the original image is divided by the total of number of bits which is obtained by adding the number of bit streams of WDR for LL and that of the SVD compression for LH, HL, and HH.

Decompression is carried out by taking the inverse WDR (IWDR) of the bit streams in order to reconstruct the LL subband and in parallel the matrix multiplications are conducted in order to reconstruct LH, HL, and HH subbands. Due to the losses by ignoring low-valued singular values, high-frequency subbands need to be enhanced. For this purpose, stationary wavelet transform (SWT) is applied to the LL subband image which results in new low- and highfrequency subbands. These high-frequency subbands will have the same direction as the ones obtained by DWT (e.g., horizontal, vertical, and diagonal), so they will be added to the respective ones reconstructed by matrix multiplications. Now, the LL subband image obtained by IWDR and the enhanced LH, HL, and HH subbands are combined by using inverse DWT (IDWT) in order to reconstruct the decompressed image. The enhancement of high-frequency subbands by using SWT results in more sharpened decompressed image. The block diagram of the proposed lossy image compression technique is shown in Fig. 8. The experimental

qualitative and quantitative results are represented and discussed in the next section.

As it was mentioned in the Introduction, for comparison purposes, the proposed lossy image compression was tested on many benchmark images, namely, Lena, Pepper, Boats, Airfield,

**4. Experimental results and discussion**

*|wQ| ≥ 2T* are then refined in order to reduce error [27, 29, 30].

44 Wavelet Transform and Some of Its Real-World Applications

**3. The proposed lossy image compression technique**

and Goldhill. All the input images were of resolutions *256* x *256* pixels, 8-bit grayscale. Tables 1, 2, and 3 provide a quantitative comparison between the proposed technique, JPEG2000, and WDR [40, 41] based on PSNR values, in dB, for compression ratios 20:1, 40:1, and 80:1, respectively.

The foregoing tables illustrate the superiority of the proposed method in terms of its capability in leading to significantly higher PSNR values compared to the other techniques proposed, previously, in the literature. It is worth noticing that the improvement in the PSNR values brought about by considering the proposed method might better show its impact while keeping in mind the fact that they are calculated in dB, meaning that a logarithmic function determines them, which clarifies how considerable the difference between the actual values has been. To be more clear, if one calculates the difference between the PSNR values obtained using WDR and JPEG2000, and subsequently, that of JPEG2000 and the proposed method, it can be seen that the latter is much higher than the former, although JPEG2000 has always been deemed of significantly better performance than WDR. Thus, it can be concluded that the proposed method makes an enormous enhancement to the PSNR values compared to the ones obtained upon employing WDR or JPEG2000.


**Table 1.** PSNR values in dB for 20:1 compression


**Table 2.** PSNR values in dB for 40:1 compression


**Table 3.** PSNR values in dB for 80:1 compression

In order to ensure the quality of the output of the proposed technique, and for visual illustra‐ tion, the images resulting from the implementation of the foregoing approach were obtained, along with that of JPEG2000 and WDR. Figure 9 shows a part of the magnified Lena image having been compressed using the foregoing approaches, separately, with compression ratio 40:1. As sought from the outset, the proposed method is competent enough to maintain the quality of the image while compressing it, and at the same time, result in better PSNR, which shows its capability in correctly deciding on a reasonable trade-off between the amount of data needed to be transferred, or kept, and the visibility and authenticity of the details in the image

blocks, which is, probably, the most tricky criterion in devising image compression algorithms. As Fig. 9 illustrates, the overall quality of the Lena image being compressed by the proposed method is satisfactory despite possessing much higher PSNR value compared to the JPEG2000 and WDR techniques, and the details are clear and visible, even better than the output of the WDR. the proposed method is satisfactory despite possessing much higher PSNR value compared to the JPEG2000 and WDR techniques, and the details are clear and visible, even better than the output of the

Boats 26.96 26.76 **30.19** Airfield 22.71 22.64 **27.32** Goldhill 27.72 27.69 **32.64**

In order to ensure the quality of the output of the proposed technique, and for visual illustration, the

images resulting from the implementation of the foregoing approach were obtained, along with that of

JPEG2000 and WDR. Figure 9 shows a part of the magnified Lena image having been compressed using

the foregoing approaches, separately, with compression ratio 40:1. As sought from the outset, the

proposed method is competent enough to maintain the quality of the image while compressing it, and at

the same time, result in better PSNR, which shows its capability in correctly deciding on a reasonable

trade-off between the amount of data needed to be transferred, or kept, and the visibility and authenticity

of the details in the image blocks, which is, probably, the most tricky criterion in devising image

**Figure 9.** Zoomed (a) original Lena image, and compressed images by using (b) JPEG2000, (c) WDR, and (d) the pro‐ posed image compression technique at 40:1 compression ratio.

### **5. Conclusion**

WDR.

**Image**

**Image**

**Image**

**Table 1.** PSNR values in dB for 20:1 compression

46 Wavelet Transform and Some of Its Real-World Applications

**Table 2.** PSNR values in dB for 40:1 compression

**Table 3.** PSNR values in dB for 80:1 compression

**Techniques WDR JPEG2000 Proposed**

**Techniques WDR JPEG2000 Proposed**

**Techniques WDR JPEG2000 Proposed**

Lena 35.72 35.99 **39.14** Pepper 34.21 35.07 **40.07** Boats 32.42 33.18 **35.97** Airfield 27.02 27.32 **31.43** Goldhill 31.76 32.18 **38.05**

Lena 32.44 32.75 **35.98** Pepper 31.67 32.40 **36.45** Boats 29.32 29.76 **32.03** Airfield 24.72 24.88 **29.62** Goldhill 29.43 29.72 **34.19**

Lena 29.71 29.62 **32.46** Pepper 28.93 29.54 **33.07** Boats 26.96 26.76 **30.19** Airfield 22.71 22.64 **27.32** Goldhill 27.72 27.69 **32.64**

In order to ensure the quality of the output of the proposed technique, and for visual illustra‐ tion, the images resulting from the implementation of the foregoing approach were obtained, along with that of JPEG2000 and WDR. Figure 9 shows a part of the magnified Lena image having been compressed using the foregoing approaches, separately, with compression ratio 40:1. As sought from the outset, the proposed method is competent enough to maintain the quality of the image while compressing it, and at the same time, result in better PSNR, which shows its capability in correctly deciding on a reasonable trade-off between the amount of data needed to be transferred, or kept, and the visibility and authenticity of the details in the image In this research work, a new lossy image compression technique which uses singular value decomposition and wavelet difference reduction techniques, followed by resolution enhance‐ ment, using discrete wavelet transform and stationary wavelet transform was proposed.

As the first step in the proposed image compression technique, the input image was decom‐ posed into four different frequency subbands using discrete wavelet transform. The lowfrequency subband was compressed using wavelet difference reduction, and in parallel, the high-frequency subbands were compressed using singular value decomposition. The com‐ pression ratio was obtained by dividing the total number of bits required to represent the input image over the total bit numbers obtained by wavelet difference reduction and singular value decomposition.

Reconstruction was carried out using inverse wavelet difference reduction to obtain lowfrequency subband and reconstructing the high-frequency subbands using matrix multipli‐ cations. The high-frequency subbands were enhanced using high frequency obtained by stationary wavelet transform. The reconstructed low-frequency subband and enhanced highfrequency subbands were used to generate the reconstructed image using inverse discrete wavelet transform.

The visual and quantitative experimental results of the proposed image compression technique showed that the proposed image compression technique outperformed the wavelet difference reduction and JPEG2000 techniques.

### **Acknowledgements**

The research was supported by the ERDF program, "Estonian higher education information and communications technology and research and development activities state program 2011-2015 (ICT program)" and Estonian Research Council grant (PUT638).

### **Author details**

Gholamreza Anbarjafari1,3\*, Pejman Rasti1 , Morteza Daneshmand1 and Cagri Ozcinar2

\*Address all correspondence to: sjafari@ut.ee

1 iCV Group, IMS Lab, Institute of Technology, University of Tartu, Tartu, Estonia

2 Centre for Vision, Speech and Signal Processing, University of Surrey, Guildford, United Kingdom

3 Dept. of Electrical and Electronic, Faculty of Engineering, Hasan Kalynocu University, Gaziantep, Turkey

#### **References**

[1] T. Fujii, D. Shirai, Y. Tonomura, M. Kitamura, T. Nakachi, T. Sawabe, M. Ogawara, T. Yamaguchi, M. Nomura and K. Shirakawa, "Digital cinema and super-high-defini‐ tion content distribution on optical high-speed networks," in *Proceedings of the IEEE*, pp. 140-153, 2013.

[2] A. Qian, C. Rongguan, N. a. G. S. Ning, Z. Xiaochen, Z. Jie and K. Deyu, "High-Defi‐ nition Image Processing Algorithm and Digital Platform Design," in *Computer and In‐ formation Technology (CIT), 2012 IEEE 12th International Conference on*, 2012.

image over the total bit numbers obtained by wavelet difference reduction and singular value

Reconstruction was carried out using inverse wavelet difference reduction to obtain lowfrequency subband and reconstructing the high-frequency subbands using matrix multipli‐ cations. The high-frequency subbands were enhanced using high frequency obtained by stationary wavelet transform. The reconstructed low-frequency subband and enhanced highfrequency subbands were used to generate the reconstructed image using inverse discrete

The visual and quantitative experimental results of the proposed image compression technique showed that the proposed image compression technique outperformed the wavelet difference

The research was supported by the ERDF program, "Estonian higher education information and communications technology and research and development activities state program

, Morteza Daneshmand1

and Cagri Ozcinar2

2011-2015 (ICT program)" and Estonian Research Council grant (PUT638).

1 iCV Group, IMS Lab, Institute of Technology, University of Tartu, Tartu, Estonia

2 Centre for Vision, Speech and Signal Processing, University of Surrey, Guildford, United

3 Dept. of Electrical and Electronic, Faculty of Engineering, Hasan Kalynocu University,

[1] T. Fujii, D. Shirai, Y. Tonomura, M. Kitamura, T. Nakachi, T. Sawabe, M. Ogawara, T. Yamaguchi, M. Nomura and K. Shirakawa, "Digital cinema and super-high-defini‐ tion content distribution on optical high-speed networks," in *Proceedings of the IEEE*,

decomposition.

wavelet transform.

**Acknowledgements**

**Author details**

Kingdom

Gaziantep, Turkey

pp. 140-153, 2013.

**References**

reduction and JPEG2000 techniques.

48 Wavelet Transform and Some of Its Real-World Applications

Gholamreza Anbarjafari1,3\*, Pejman Rasti1

\*Address all correspondence to: sjafari@ut.ee


[29] J. Tian and R. O. Wells Jr, "Image data processing in the compressed wavelet do‐ main," in *Signal Processing, 1996., 3rd International Conference on*, pp. 978-981 1996.

[14] Y. Zhu, Y. Zhang and A. Yuille, "Single Image Super-resolution Using Deformable Patches," in *Computer Vision and Pattern Recognition (CVPR), 2014 IEEE Conference on*,

[15] K. Nasrollahi and T. Moeslund, "Super-resolution: a comprehensive survey," *Machine*

[16] H. Demirel and G. Anbarjafari, "Data fusion boosted face recognition based on prob‐ ability distribution functions in different colour channels," *EURASIP Journal on Ad‐*

[17] H. Demirel and G. Anbarjafari, "Iris recognition system using combined colour statis‐

[18] X. Zhang, "Lossy compression and iterative reconstruction for encrypted image," *In‐ formation Forensics and Security, IEEE Transactions on,* vol. 16, no. 1, pp. 53-58, 2011.

[19] F. Hussain and J. Jeong, "Efficient Deep Neural Network for Digital Image Compres‐

[20] M. Boliek, "Beyond compression: a survey of functionality derived from still image coding," in *Signals, Systems and Computers, 2004. Conference Record of the Thirty-Seventh*

[21] R. C. Gonzalez and R. E. Woods, Digital Image Processing, Prentice Hall Upper Sad‐

[22] A. Karami, M. Yazdi and G. Mercier, "Compression of Hyperspectral Images Using Discerete Wavelet Transform and Tucker Decomposition," *Selected Topics in Applied Earth Observations and Remote Sensing, IEEE Journal of,* vol. 5, no. 2, pp. 444-450, 2012.

[23] H. a. A. G. Demirel, C. Ozcinar and S. Izadpanahi, "Video Resolution Enhancement By Using Complex Wavelet Transform," in *Image Processing (ICIP), 2011 18th IEEE In‐*

[24] M. Ghazel, G. H. Freeman and E. R. Vrscay, "Fractal-wavelet image denoising revisit‐ ed," *{Image Processing, IEEE Transactions on,* vol. 15, no. 9, pp. 2669-2675, 2006.

[25] J.-F. Yang and C.-L. Lu, "Combined Techniques of Singular Value Decomposition and Vector Quantization for Image Coding," *IEEE Transactions on Image Processing,*

[27] J. Tian and R. O. Wells Jr, "Embedded image coding using wavelet difference reduc‐

[28] J. Tian and R. O. Wells Jr, "A lossy image codec based on index coding," in *IEEE Data*

[26] M. D. Greenberg, Differential equations & Linear algebra, Prentice Hall, 2001.

tion," in *Wavelet image and video compression*, pp. 289-301, 2002.

tics.," in *Signal Processing and Information Technology, pp. 175-179, 2008.*

sion Employing Rectified Linear Neurons," *Journal of Sensors,* 2015.

*Vision and Applications,* vol. 25, no. 6, pp. 1423-1468, 2014.

*vances in Signal Processing,* vol. 2009, p. 25, 2009.

*Asilomar Conference on*, pp. 1971-1974, 2003.

*ternational Conference on*, pp. 2093-2096, 2011.

vol. 4, no. 8, pp. 1141-1146, 1995.

*Compression Conference*,p. 456, 1996.

pp. 2917-2924 2014.

50 Wavelet Transform and Some of Its Real-World Applications

dle River, NJ, 2002.


## **Adaptive Wavelet Packet Transform**

## Zuofeng Zhou and Jianzhong Cao

Additional information is available at the end of the chapter

http://dx.doi.org/10.5772/61946

#### **Abstract**

Two-dimensional over-complete wavelet packet transform can better represent the tex‐ ture and long oscillatory patterns in natural images.

In this chapter, combining the doubly Wiener filtering algorithm and Wiener cost func‐ tion, a new best wavelet packet decomposition scheme for image denoising applications is proposed. The experiment results for the test image database show the effectiveness of the proposed image denoising algorithm compared to some existing image denoising methods.

**Keywords:** Adaptive wavelet packet transform, Wiener cost function, image denoising

#### **1. Introduction**

Wavelet with vanishing moments are very effective for representing piecewise smooth images [1]. However, two-dimensional separable wavelets are ill-suited to represent long oscillatory patterns in images with abundant textures, partly owing to their poor directional selectivity in frequency domain. These oscillatory variations of intensity can only be represented by smallscale wavelet coefficients. In some image-processing applications, such as image denoising or image compression, those small-scale coefficients are quantized to zero in the low bit rate image compression and are thresholded or shrunken to zero in image denoising, which degrades compression and denoising performance significantly. To overcome this circum‐ stance, one way is to find the more suitable image representation techniques such as curvelet, contourlet [2], and ridgelet [3]. But these methods need the researchers to design the new directional compact filter or multichannel filter banks, which is also challenging in filter design area. Another way is to improve the idea of wavelet design method to accommodate the new requirement where the over-complete wavelet packet decomposition is proposed.

© 2015 The Author(s). Licensee InTech. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Over-complete wavelet packet contains a mass of libraries of waveforms, from which the best wavelet packet base can be selected to efficiently represent long oscillatory patterns given the corresponding criterion. For example, in image compression application, the best wavelet packet base is usually found by pruning the full wavelet packet decomposition given a user predefined cost function. Recently, a large variety of cost function have been proposed, such as Shannon's entropy cost function and vector entropy cost function, which are used in the rate distortion control strategy. On the other hand, different from finding the best tree structure of the wavelet packet decomposition, some researchers devoted themselves for finding the best wavelet packet decomposition base [4] which is equivalent to find the best filter banks of wavelet packet decomposition.

In image denoising application, due to unknown noiseless image as well as a great diversity of filtering methods in the transform domain, selecting the best wavelet packet base is always a difficult problem. Enlightened by the idea of doubly local Wiener filtering method [5] and spatially adaptive wavelet domain shrinkage image denoising algorithms, the best wavelet packet bases selection method using the local Wiener cost function and its application in image denoising is proposed. In this chapter, we will first give the detail of the best wavelet packetbased selection algorithm and then discuss its application in image denoising.

#### **2. Best wavelet packet base selection**

Let *ψ*<sup>0</sup> (*x*), *ψ*<sup>1</sup> (*x*) be the wavelet function and scaling functions, the basic idea of wavelet packet function can be defined as:

$$\begin{aligned} \boldsymbol{\nu}^{\circ 2u}(\mathbf{x}) &= \sqrt{2} \sum\_{k} h(k) \boldsymbol{\nu}^{\circ}(2\mathbf{x} - k), \\ \boldsymbol{\nu}^{\circ 2u + 1}(\mathbf{x}) &= \sqrt{2} \sum\_{k} g(k) \boldsymbol{\nu}^{\circ}(2\mathbf{x} - k), \quad \text{for } n \ge 1 \end{aligned} \tag{1}$$

where *h* (*k*), *g*(*k*) are the orthonormal filters, respectively. Given a signal *s*(*n*) which is satisfied the Nyquist sampling criterion, subspace *V*<sup>0</sup> <sup>≡</sup>span{*ψ*0(*<sup>x</sup>* <sup>−</sup>*l*):*<sup>l</sup>* <sup>∈</sup>**ℤ**} is usually assumed to be the first-level subspace of signal *s*(*n*). The subspace with depth *j* is defined as *Vj <sup>n</sup>* <sup>≡</sup>span{*ψj*,*<sup>l</sup> <sup>n</sup>* (*x*)=2<sup>−</sup> *<sup>j</sup>*/2*ψn*(2<sup>−</sup> *<sup>j</sup> x* −*l*), *l* ∈**ℤ**}. To understand this subspace easily, the dynamic interval 2<sup>−</sup> *<sup>j</sup> n*, 2<sup>−</sup> *<sup>j</sup>* (*n* + 1)) can be associated with the subspace *Vj <sup>n</sup>*. For example, the interval 0, 1) is equivalent to *V*0.

#### **2.1. Two-dimensional wavelet packet bases**

The two-dimensional separable wavelet packet functions are defined as the tensor products of two one-dimensional wavelet packet functions, that are,

$$\begin{aligned} \, \_0V\_0 &= \text{span}\{\psi \, ^0(\text{x} - l)\psi \, ^0(y - p) : (l, \, p) \in \mathbb{Z}^2\} \\ \, \_0\psi \, ^{m, n}(\text{x}, \, \_yy) &= \psi \, ^m(\text{x})\psi \, ^n(y), \quad m, \, n \in \text{0}, \, 1, \, ... \end{aligned} \tag{2}$$

Similar to the one-dimensional case, the subspace

Over-complete wavelet packet contains a mass of libraries of waveforms, from which the best wavelet packet base can be selected to efficiently represent long oscillatory patterns given the corresponding criterion. For example, in image compression application, the best wavelet packet base is usually found by pruning the full wavelet packet decomposition given a user predefined cost function. Recently, a large variety of cost function have been proposed, such as Shannon's entropy cost function and vector entropy cost function, which are used in the rate distortion control strategy. On the other hand, different from finding the best tree structure of the wavelet packet decomposition, some researchers devoted themselves for finding the best wavelet packet decomposition base [4] which is equivalent to find the best filter banks of

In image denoising application, due to unknown noiseless image as well as a great diversity of filtering methods in the transform domain, selecting the best wavelet packet base is always a difficult problem. Enlightened by the idea of doubly local Wiener filtering method [5] and spatially adaptive wavelet domain shrinkage image denoising algorithms, the best wavelet packet bases selection method using the local Wiener cost function and its application in image denoising is proposed. In this chapter, we will first give the detail of the best wavelet packet-

(*x*) be the wavelet function and scaling functions, the basic idea of wavelet packet

<sup>å</sup> (1)

*<sup>n</sup>*. For example, the interval 0, 1)

(2)

*x* −*l*), *l* ∈**ℤ**}. To understand this subspace easily, the dynamic

based selection algorithm and then discuss its application in image denoising.

( ) 2 ( ) (2 ),

= -

(*n* + 1)) can be associated with the subspace *Vj*

(*x* −*l*)*ψ*<sup>0</sup>

*ψ <sup>m</sup>*,*n*(*x*, *y*)=*ψ <sup>m</sup>*(*x*)*ψn*(*y*), *m*, *n* ∈0, 1, …

 y

*x hk x k*

*n n k n n k*

å

<sup>+</sup>

( ) 2 ( ) (2 ), for 1

= -³

where *h* (*k*), *g*(*k*) are the orthonormal filters, respectively. Given a signal *s*(*n*) which is satisfied the Nyquist sampling criterion, subspace *V*<sup>0</sup> <sup>≡</sup>span{*ψ*0(*<sup>x</sup>* <sup>−</sup>*l*):*<sup>l</sup>* <sup>∈</sup>**ℤ**} is usually assumed to be the first-level subspace of signal *s*(*n*). The subspace with depth *j* is defined as

The two-dimensional separable wavelet packet functions are defined as the tensor products

(*y* − *p*):(*l*, *p*)∈**ℤ**2}

*x gk x k n*

 y

wavelet packet decomposition.

54 Wavelet Transform and Some of Its Real-World Applications

Let *ψ*<sup>0</sup>

*Vj*

*<sup>n</sup>* <sup>≡</sup>span{*ψj*,*<sup>l</sup>*

is equivalent to *V*0.

interval 2<sup>−</sup> *<sup>j</sup>*

(*x*), *ψ*<sup>1</sup>

function can be defined as:

**2. Best wavelet packet base selection**

2 2 1

y

y

*<sup>n</sup>* (*x*)=2<sup>−</sup> *<sup>j</sup>*/2*ψn*(2<sup>−</sup> *<sup>j</sup>*

**2.1. Two-dimensional wavelet packet bases**

of two one-dimensional wavelet packet functions, that are,

*<sup>V</sup>*<sup>0</sup> =span{*ψ*<sup>0</sup>

*n*, 2<sup>−</sup> *<sup>j</sup>*

$$V\_j^{m,n} \equiv \text{span}\{\psi\_{j,l,p}^{m,n}(\mathbf{x}) = \mathbf{2}^{-j}\psi^m(\mathbf{2}^{-j}\mathbf{x} - l)\psi^n(\mathbf{2}^{-j}y - p), \ (l, \ p) \in \mathbb{Z}^2\}\tag{3}$$

is referred to as the subspace with depth *j* and a wavelet packet function *ψ <sup>m</sup>*,*n*(*x*, *y*). Figure 1 gives the subspace of 2D wavelet packet base diagram. It can be seen that each node is divided into two branches. Let us associate the dyadic square 2<sup>−</sup> *<sup>j</sup> m*, 2 *<sup>j</sup>* (*m* + 1))× 2<sup>−</sup> *<sup>j</sup> n*, 2<sup>−</sup> *<sup>j</sup>* (*n* + 1)) with the subspace *Vj <sup>m</sup>*,*<sup>n</sup>*, for example, *V*<sup>0</sup> associates with the square 0, 1)<sup>2</sup> , *V*<sup>1</sup> 1,0 associates with the rectangle 1 / 2, 1)× 0, 1 / 2), and *V*<sup>1</sup> 1,1 associates with the square 1 / 2, 1)<sup>2</sup> . It has been proved that: if the squares 2<sup>−</sup> *<sup>j</sup> m*, 2<sup>−</sup> *<sup>j</sup>* (*m* + 1))× 2<sup>−</sup> *<sup>j</sup> n*, 2<sup>−</sup> *<sup>j</sup>* (*n* + 1)), (*m*, *n*, *j*)∈Ω are a partition of the unit square 0, 1)<sup>2</sup> , then the family of functions {<sup>Ψ</sup> *<sup>j</sup>*;*l*, *<sup>p</sup> <sup>m</sup>*,*<sup>n</sup>* (*x*, *<sup>y</sup>*):(*l*, *<sup>p</sup>*)∈**ℤ**<sup>2</sup> , (*m*, *n*, *j*)∈Ω} constitutes an orthogonal wavelet packet base of *V*0.

**Figure 1.** The subspace of 2D wavelet packet base

#### **2.2. Best wavelet packet base selection under wiener cost function**

Similar to the tree structure illustrated in Fig. 1, 2D wavelet packet decomposition can be easily expressed by the quad-tree structure with the root node {0, 0, 0}. Each node (*m*, *n*, *j*) (except the last level node) has four child nodes (2*m*, 2*n*; *j* + 1), (2*m*, 2*n* + 1, *j* + 1), (2*m* + 1, 2*n*, *j*) and (2*m* + 1, 2*n* + 1, *j* + 1). This quad-tree structure facilitates the best wavelet packet decomposition procedure. Given the image and the noise level *σ* (if unknown, the noise level can be estimated by the MAD estimator in the wavelet domain), the wavelet packet coefficient at the node (*m*, *n*, *j*) under J-level wavelet packet decomposition is represented as *yj <sup>m</sup>*,*n*(*p*, *q*), *m*, *n* =0, 1, …, 2 *<sup>j</sup>* −1; *j* =0, 1, …, *J*. For each node, its Wiener cost function is computed by

$$J\left(m,n;j\right) = \sigma^2 \sum\_{\boldsymbol{\nu}} \sum\_{\boldsymbol{\nu}} \frac{\left[\boldsymbol{\mathcal{Y}}\_{j}^{m,n}(\boldsymbol{\mathcal{P}},\boldsymbol{\eta})\right]^2}{\left[\boldsymbol{\mathcal{Y}}\_{j}^{m,n}(\boldsymbol{\mathcal{P}},\boldsymbol{\eta})\right]^2 + \sigma^2}, \quad m,n = 0,1,\ldots,2^j - 1; j = 0,1,\ldots,J. \tag{4}$$

Assuming all the node in the quad-tree decomposition forms a set Class (**A**)={Ω={(*m*, *n*, *j*)}, our task is to find a subset Ω\* from set Ω which have the smallest total cost function (the summation of the Wiener cost function at all wavelet packet decomposition node).

Similar to most of the search algorithms of the best wavelet packet base, the best wavelet packet base under Wiener cost function is obtained by pruning a J-level full quad-tree from bottom to top. The searching algorithm can be described as follows:


$$J\left(m,n;j\right) = \sigma^2 \sum\_{p} \sum\_{q} \frac{\left[\underline{y}^{m,n}(p,q)\right]^2}{\left[\underline{y}^{m,n}\_{/}\left(p,q\right)\right]^2 + \sigma^2} \tag{5}$$

as well as the summation of the Wiener cost functions of its four child nodes by

$$\begin{aligned} \tilde{J}(m, n, j) &= J(S(2m, 2n, j+1)) + J(S(2m, 2n+1, j+1)) + \\ &+ J(S(2m+1, 2n, j+1)) + J(S(2m+1, 2n+1, j+1)) \end{aligned}$$

**iii.** If *J* ˜(*m*, *n*, *j*)< *J*(*m*, *n*, *j*), then

$$\mathcal{S}(m,n,j) = \mathcal{S}(2m, 2n, j+1) \bigcup \mathcal{S}(2m, 2n+1, j+1) \bigcup \mathcal{S}(2m+1, 2n, j+1) \bigcup \mathcal{S}(2m+1, 2n+1, j+1)$$

otherwise, *S*(*m*, *n*, *j*)={(*m*, *n*, *j*)}.

**iv.** Calculate the Wiener cost function of each set *S*(*m*, *n*, *j*) by

$$J\left(\mathcal{S}\left(m,n,j\right)\right) \equiv \sum\_{\{u,v,r\} \in \mathcal{S}\left(u,v,l\right)} J\left(u,v,r\right), \quad m,n = 0,1,...,2^j - 1,\tag{6}$$

**v.** If *j* >0, set *j* ← *j* −1 and return the step (ii); otherwise, output Ω\* =*S*(0, 0, 0).

The set Ω\* is composed of all leaf nodes of the best quad-tree decomposition structure. In this way, given the input noisy image and the noise level, the optimal orthogonal wave‐ let packet decomposition structure and the minimal total Wiener cost function can be obtained. Figure 2 shows a demo of best wavelet packet decomposition where Fig. 2(a) is its square representation and Fig. 2 (b) is the corresponding tree structure of this best wavelet packet decomposition.

**Figure 2.** A demo of best wavelet packet decomposition

Assuming all the node in the quad-tree decomposition forms a set Class (**A**)={Ω={(*m*, *n*, *j*)},

Similar to most of the search algorithms of the best wavelet packet base, the best wavelet packet base under Wiener cost function is obtained by pruning a J-level full quad-tree from bottom

**i.** Let *S*(*m*, *n*, *J*)={(*m*, *n*, *J*)}, *m*, *n* =0, 1, ⋯, 2*<sup>J</sup>* −1 represent 22J leaf nodes in the J-level full

**ii.** For 0≤ *j* < *J* and each node (*m*, *n*, *j*) in the j-th level, calculate the corresponding Wiener

*m n p q <sup>j</sup>*

ë û <sup>=</sup>

( , , ) ( (2 ,2 , 1)) ( (2 ,2 1, 1)) ( (2 1,2 , 1)) ( (2 1,2 1, 1)) *Jmnj JS m nj JS m n j JS m nj JS m n j*

*Smnj S m nj S m n j S m nj S m n j* ( , , ) (2 ,2 , 1) (2 ,2 1, 1) (2 1,2 , 1) (2 1,2 1, 1) = + ++ + + + ++ UUU

, , , , , , 0,1, ,2 1, *<sup>j</sup>*

this way, given the input noisy image and the noise level, the optimal orthogonal wave‐ let packet decomposition structure and the minimal total Wiener cost function can be obtained. Figure 2 shows a demo of best wavelet packet decomposition where Fig. 2(a) is its square representation and Fig. 2 (b) is the corresponding tree structure of this best

is composed of all leaf nodes of the best quad-tree decomposition structure. In

+ + ++ + ++

( )

é ù <sup>+</sup> ë û

= ++ +++

*m n j*

*y pq*

,

é ù

<sup>2</sup> ,

<sup>2</sup> , 2 (,)

s

º = å <sup>K</sup> - (6)

=*S*(0, 0, 0).

åå (5)

summation of the Wiener cost function at all wavelet packet decomposition node).

to top. The searching algorithm can be described as follows:

( )

%

˜(*m*, *n*, *j*)< *J*(*m*, *n*, *j*), then

otherwise, *S*(*m*, *n*, *j*)={(*m*, *n*, *j*)}.

wavelet packet decomposition.

, ;

2

s

**iv.** Calculate the Wiener cost function of each set *S*(*m*, *n*, *j*) by

( ) ( ) ( ) ( , ,) ( , ,)

*uvr Smnj J S mnj J uvr mn* Î

**v.** If *j* >0, set *j* ← *j* −1 and return the step (ii); otherwise, output Ω\*

*y pq J mnj*

as well as the summation of the Wiener cost functions of its four child nodes by

from set Ω which have the smallest total cost function (the

our task is to find a subset Ω\*

56 Wavelet Transform and Some of Its Real-World Applications

quad-tree.

cost function

**iii.** If *J*

The set Ω\*

#### **3. Image denoising algorithms based on best wavelet packet decomposition**

This section gives two image denoising algorithms based on the best wavelet packet decom‐ position. In the above section, the optimal wavelet packet decomposition structure is obtained under the Wiener cost function. When computing the Wiener function, the noiseless image and the noise level are assumed as known. In practice, this assumption is unreasonable. It is well known that the noise level can be accurately estimated by the MAD estimator from the input noisy image. The key problem is how to estimate the noiseless image. Enlightened by the empirical Wiener filtering and the doubly local Wiener filtering image denoising algo‐ rithms, we first filter the noisy image to get the pilot image, and then the pilot image is used to get the best wavelet packet decomposition structure.

It is well known that undecimated wavelet packet decomposition can achieve better image denoising performance than the decimated wavelet packet decomposition. This is owing to the property that they are shift invariance and robustness. So, the undecimated wavelet packet decomposition can ameliorate some unpleasant phenomenon that appears in the maximally decimated wavelet packet decomposition, such as Gibbs-like ringing around edges and specks in smooth region.

In what follows, we will give the detail of the image denoising algorithm using the undeci‐ mated wavelet decomposition. Let *x*(*p*, *q*) be a input noisy image of size 2*<sup>N</sup>* ×2*<sup>N</sup>* , the operator Shift*<sup>k</sup>* ,*<sup>l</sup>* (*x*) denotes circularly shifting the input image *x*(*p*, *q*) by *k* indices in the vertical direction and *l* indices in the horizontal direction, and the operator Unshift*<sup>k</sup>* ,*<sup>l</sup>* (*x*) is a similar operation but in the opposite direction. In the proposed algorithm, we first use the decimated best wavelet packet decomposition algorithm to denoise all possible shift versions of the noisy image to get a set of denoised images and then unshifting these denoised images and average them to get the final denoised image. The new algorithm is referred to as the local Wiener filtering using the best undecimated wavelet packet decomposition (LWF-BUWPD), which can be summarized as follows:


$$E\_{\boldsymbol{\beta}}^{m,n}(\boldsymbol{p},\boldsymbol{q}) = \left\{ \frac{1}{\#\mathcal{W}} \sum\_{(\boldsymbol{s},\boldsymbol{t}) \in \mathcal{W}} \hat{\mathbf{s}}\_{\boldsymbol{\beta},\boldsymbol{w},n}^{2}(\boldsymbol{p}+\boldsymbol{s},\boldsymbol{q}+\boldsymbol{t}) \right\}\_{\boldsymbol{\epsilon}}, \boldsymbol{\epsilon}(\boldsymbol{m},\boldsymbol{n},\boldsymbol{\bar{j}}) \in \Omega^{\*} \tag{7}$$

where *W* and *s* ^ *<sup>j</sup>*,*m*,*n*(*p*, *q*) represent the directional window and the pilot image's best wavelet packet decomposition coefficients, respectively.

**iv.** Using the estimated energy distributions and noise level, the local Wiener filtering is operated on the base wavelet packet decomposition coefficients of the noisy image, that is,

$$\mathfrak{S}\_{j}^{m,n}(p,q) = \frac{E\_{j}^{m,n}(p,q)}{E\_{j}^{m,n}(p,q) + \sigma^{2}} \, y\_{j}^{m,n}(p,q) \,, \text{(\$m,n\$,j\$)} \in \Omega^{\*} \tag{8}$$

where *yj <sup>m</sup>*,*n*(*p*, *q*) are the best wavelet packet decomposition coefficients of the noisy image in the subspace *Vj <sup>m</sup>*,*<sup>n</sup>*.

**v.** Unshift all the shifted denoised images and average them to obtain the final denoised image *s*˜(*p*, *q*)

$$\mathfrak{S}(p\_{\prime}q) = \frac{1}{2^{2^{l}}} \sum\_{k=0}^{2^{l}-1} \sum\_{l=0}^{2^{l}-1} \text{Linshiff}\_{k,l}(\mathfrak{s}\_{k,l})(p\_{\prime}q) \tag{9}$$

#### **4. Experimental results**

We choose the 8-bit 512 × 512 grayscale images "Lena" and "Barbara," and a 256 × 256 texture image, as the test images. In the proposed image denoising algorithm, the best wavelet packet decomposition structure is varied with different shift noisy images and different noise level. For better illustration, the best wavelet packet decomposition structure for different noise level is shown in Fig. 3 for "Barbara" image.

them to get the final denoised image. The new algorithm is referred to as the local Wiener filtering using the best undecimated wavelet packet decomposition (LWF-BUWPD), which

**i.** Using the local Wiener filtering image denoising algorithm to each shift of input noisy

**ii.** Given the pilot image and noise level, the best wavelet packet decomposition structure can be obtained by using the searching algorithm in section 2.

**iii.** Given the best wavelet packet decomposition structure, the empirical energy

ì ü <sup>=</sup> í ý + + ÎW

**iv.** Using the estimated energy distributions and noise level, the local Wiener filtering is

\* = ÎW

**v.** Unshift all the shifted denoised images and average them to obtain the final denoised

We choose the 8-bit 512 × 512 grayscale images "Lena" and "Barbara," and a 256 × 256 texture image, as the test images. In the proposed image denoising algorithm, the best wavelet packet decomposition structure is varied with different shift noisy images and different noise level. For better illustration, the best wavelet packet decomposition structure for different noise level

21 21 2 0 0 , , <sup>1</sup> (,) ( )( , ) <sup>2</sup> *J J <sup>J</sup> k l kl kl spq Unshift s p q* - -

operated on the base wavelet packet decomposition coefficients of the noisy image,

*<sup>m</sup>*,*n*(*p*, *q*) are the best wavelet packet decomposition coefficients of the noisy image in

^(*<sup>p</sup>*, *<sup>q</sup>*) ;

+

<sup>+</sup> % (8)

= = % % <sup>=</sup> å å (9)

*<sup>j</sup>*,*m*,*n*(*p*, *q*) represent the directional window and the pilot image's best wavelet

î þ <sup>å</sup> (7)

\*

can be summarized as follows:

58 Wavelet Transform and Some of Its Real-World Applications

*m n*

packet decomposition coefficients, respectively.

^

that is,

*<sup>m</sup>*,*<sup>n</sup>*.

image *s*˜(*p*, *q*)

**4. Experimental results**

is shown in Fig. 3 for "Barbara" image.

where *W* and *s*

where *yj*

the subspace *Vj*

image shift (*x*(*p*, *q*)) to get a pilot image *s*

, 2

distribution of the pilot image can be estimated by

( )

, , , , 2 (,) (,) ( , ),( , , ) (,)

*E pq s pq y pq mnj E pq* s

*m n m n j m n j j m n j*

Î

, , , <sup>1</sup> (,) <sup>ˆ</sup> ( , ) ,( , , ) #

*<sup>j</sup> st W jmn E pq s p sq t mnj <sup>W</sup>*

**Figure 3.** The best wavelet packet trees for the "Barbara" image with different noise levels: (a) *σ* =10 ; (b) *σ* =15 ; (c) *σ* =20 ; (d) *σ* =25.

In Table 1 and Fig. 4, we give the denoising performance of the image denoising algorithms using the undecimated best wavelet packet decomposition. The experimental results show that for images of structural textures, for example, "Barbara" and texture images, the proposed algorithm greatly improves denoising performance as compared with the existing state-of-theart algorithms.


**Table 1.** The performance comparison of the LWF-BUWPD and several state-of-the-art image denoising algorithms

LWF-BUWPD 35.60 33.87 32.65 31.66 34.35 32.28 30.80 29.63 34.97 32.92 31.72 30.65

8

DFB-GSM algorithm (the left bottom corner, PSNR = 30.93 dB), and the denoised image by the LWF-BUWPD algorithm (the right bottom corner, PSNR = 31.70 dB); (b). Zoomed local regions of the four images in (a). **Figure 4.** (a) The noiseless image (the left top corner), the noisy image (the right top corner, noise level 20), the denoised image by the DFB-GSM algorithm (the left bottom corner, PSNR = 30.93 dB), and the denoised image by the LWF-BUWPD algorithm (the right bottom corner, PSNR = 31.70 dB); (b). Zoomed local regions of the four images in (a).

[1] Shui, P.-L., Zhou, Z.-F., and Li, J.-X., "Image denoising based on adaptive wavelet packets using wiener cost function," IET

[2] Minh, D., and Martin. V., "The contourlet transform: an efficient directional multiresolution image representation," IEEE

Fig. 4 (a) The noiseless image (the left top corner), the noisy image (the right top corner, noise level 20), the denoised image by the

#### Image Processing, 2007, 1(3), pp: 311−318. **Author details**

References

Transactions on Image Processing, Dec. 2005, 14(12), pp: 2091−2106. Zuofeng Zhou\* and Jianzhong Cao

[3] Zhang, B., Jalal M.F., and Starck J.-L., "Wavelets, ridgelets, and curvelets for poisson noise removal," IEEE Transactions on Image Processing, 2008, 17(7), pp: 1093−2009. \*Address all correspondence to: zfzhou@opt.ac.cn

 [4] Mahalakshmi, B.V., and Anand, M.J., "Adaptive wavelet packet decomposition for efficient image denoising by using Xi'an Institute of Optics and Precision Mechanics of Chinese Academy of Sciences, Xi'an, Peoples Republic of China

neighsure shrink method," International Journal of Computer Science & Information Technology 2014, 5(4), pp: 5003−5007. [5] Shui, P.-L., "Image denoising algorithm via doubly local Wiener filtering with directional windows in wavelet domain," IEEE

#### Signal Processing Letters, 2005, 12(10), pp: 681–684. **References**


[3] Zhang, B., Jalal M.F., and Starck J.-L., "Wavelets, ridgelets, and curvelets for poisson noise removal," IEEE Transactions on Image Processing, 2008, 17(7), pp: 1093–2009.

8

LWF-BUWPD 35.60 33.87 32.65 31.66 34.35 32.28 30.80 29.63 34.97 32.92 31.72 30.65

 (a) (b) Fig. 4 (a) The noiseless image (the left top corner), the noisy image (the right top corner, noise level 20), the denoised image by the DFB-GSM algorithm (the left bottom corner, PSNR = 30.93 dB), and the denoised image by the LWF-BUWPD algorithm (the right

[1] Shui, P.-L., Zhou, Z.-F., and Li, J.-X., "Image denoising based on adaptive wavelet packets using wiener cost function," IET

**Figure 4.** (a) The noiseless image (the left top corner), the noisy image (the right top corner, noise level 20), the denoised image by the DFB-GSM algorithm (the left bottom corner, PSNR = 30.93 dB), and the denoised image by the LWF-BUWPD algorithm (the right bottom corner, PSNR = 31.70 dB); (b). Zoomed local regions of the four images in

[2] Minh, D., and Martin. V., "The contourlet transform: an efficient directional multiresolution image representation," IEEE

[3] Zhang, B., Jalal M.F., and Starck J.-L., "Wavelets, ridgelets, and curvelets for poisson noise removal," IEEE Transactions on

 [4] Mahalakshmi, B.V., and Anand, M.J., "Adaptive wavelet packet decomposition for efficient image denoising by using neighsure shrink method," International Journal of Computer Science & Information Technology 2014, 5(4), pp: 5003−5007. [5] Shui, P.-L., "Image denoising algorithm via doubly local Wiener filtering with directional windows in wavelet domain," IEEE

> [1] Shui, P.-L., Zhou, Z.-F., and Li, J.-X., "Image denoising based on adaptive wavelet packets using wiener cost function," IET Image Processing, 2007, 1(3), pp: 311–318.

> [2] Minh, D., and Martin. V., "The contourlet transform: an efficient directional multire‐ solution image representation," IEEE Transactions on Image Processing, Dec. 2005,

Xi'an Institute of Optics and Precision Mechanics of Chinese Academy of Sciences, Xi'an,

bottom corner, PSNR = 31.70 dB); (b). Zoomed local regions of the four images in (a).

Transactions on Image Processing, Dec. 2005, 14(12), pp: 2091−2106.

60 Wavelet Transform and Some of Its Real-World Applications

and Jianzhong Cao

\*Address all correspondence to: zfzhou@opt.ac.cn

References

(a).

Image Processing, 2007, 1(3), pp: 311−318.

**Author details**

Zuofeng Zhou\*

**References**

Image Processing, 2008, 17(7), pp: 1093−2009.

Peoples Republic of China

Signal Processing Letters, 2005, 12(10), pp: 681–684.

14(12), pp: 2091–2106.

