**7. Algorithms for Wavelet Transform computation**

This section is concerned with a review of variety of algorithms dedicated to implement wavelet transforms. We focus on both 1-Dimensional and 2-Dimensional systems.

The Wavelet Transform for Image Processing Applications 407

The decomposition and the reconstruction processes for a 2-D signal, as in image processing, is achieved through the use of a 2-D filtering process. In this case, only 1/4 of the original signal is obtained at the output of the reducer (the decimation is performed twice). This

This type of decomposition makes this algorithm suitable for a progressive image

Mallat's pyramid is a direct consequence of the multiresolution concept developed by the same author and presented in section 6. Up to date, it is the most widely used approach both in software and hardware - for implementing the wavelet transform (Masud, 1999). Since the one-dimensional decomposition and reconstruction schemes have been already introduced in section 6, we will focus in this section on two-dimensional schemes, which are more suitable for image analysis and synthesis. The two-dimensional decomposition approach is based on the property of separation of the functions into arbitrary x and y directions. The first step is identical to the one-dimensional approach, however, instead of keeping the low-level resolution and processing the high level resolution, both are processed using two identical filter bank after a transposition of the incoming data. Thus, the image is scanned in both horizontal and vertical directions. This result in an average image (or subimage) and three detail images generated by the following 2-D scaling function φ(x,y) φ(x)φ(y) and the vertical, the horizontal and the diagonal wavelet functions: ψ y)(x, φ(x)ψ(y) <sup>1</sup> , ψ y)(x, ψ(x)φ(y) <sup>2</sup> and ψ y)(x, ψ(x)ψ(y) <sup>3</sup> , respectively. To recover the original image, the inverse process is applied. Figure 12 illustrates the analysis

Along the rows

*h* **~**

*g* **~** +

2

+

+

2

Along the Columns

2

2

Im(x,y) **<sup>~</sup>** Im(x,y)

In this case, the frequency band is halved at each stage by a factor of four as represented by

2

2

*h* **~**

*g* **~**

*h* **~**

*g* **~**

scheme can be represented by the pyramidal structure of Figure 11.

and synthesis stages built using three filter banks each.

Along the Columns

2

2

2

2

*h*

*g*

*h*

*g*

Fig. 12. Two-dimensional Mallat's analysis and synthesis tree

transmission scheme.

**7.2 Mallat's Pyramidal algorithm** 

Along the rows

*h*

*g*

Figure 13.

2

2

#### **7.1 Burt's Pyramid**

Dedicated initially to lossless image coding, the pyramid algorithm was first introduced by Burt (Burt & Adelson, 1983). Basically, it decomposes a signal in a low-resolution signal along with some higher resolution signals through a repetition of reduction and expansion processes. At each level, the reduced and expanded signal is compared with the original signal and the difference is stored. In the same time, the reduced signal is repeatedly decomposed by further using the reducer block in the chain. The analysis/synthesis process is shown in Figure 10.

Fig. 10. Pyramidal analysis and synthesis

The reduction block performs the two basic operations of a low pass filtering and decimating by a factor of 2. The expansion block up samples the signal first, then filters it through the use of a synthesis low pass filter. To reconstruct the original signal, the difference signal at each level is added to a previously expanded signal. Repeatedly, the resulting signal is expanded and added to the corresponding difference signal.

Fig. 11. 2-D Pyramidal Structure

The decomposition and the reconstruction processes for a 2-D signal, as in image processing, is achieved through the use of a 2-D filtering process. In this case, only 1/4 of the original signal is obtained at the output of the reducer (the decimation is performed twice). This scheme can be represented by the pyramidal structure of Figure 11.

This type of decomposition makes this algorithm suitable for a progressive image transmission scheme.

#### **7.2 Mallat's Pyramidal algorithm**

406 Advances in Wavelet Theory and Their Applications in Engineering, Physics and Technology

Dedicated initially to lossless image coding, the pyramid algorithm was first introduced by Burt (Burt & Adelson, 1983). Basically, it decomposes a signal in a low-resolution signal along with some higher resolution signals through a repetition of reduction and expansion processes. At each level, the reduced and expanded signal is compared with the original signal and the difference is stored. In the same time, the reduced signal is repeatedly decomposed by further using the reducer block in the chain. The analysis/synthesis process

**2**

**<sup>~</sup>** x(n) c(n)

*h* **~**

+

+ -

resulting signal is expanded and added to the corresponding difference signal. .

Expander

d(n)

The reduction block performs the two basic operations of a low pass filtering and decimating by a factor of 2. The expansion block up samples the signal first, then filters it through the use of a synthesis low pass filter. To reconstruct the original signal, the difference signal at each level is added to a previously expanded signal. Repeatedly, the

+

x(n) ~

Expander

**2** *h*

**7.1 Burt's Pyramid** 

is shown in Figure 10.

*h* **2**

Fig. 10. Pyramidal analysis and synthesis

.

Fig. 11. 2-D Pyramidal Structure

Reducer

Mallat's pyramid is a direct consequence of the multiresolution concept developed by the same author and presented in section 6. Up to date, it is the most widely used approach both in software and hardware - for implementing the wavelet transform (Masud, 1999). Since the one-dimensional decomposition and reconstruction schemes have been already introduced in section 6, we will focus in this section on two-dimensional schemes, which are more suitable for image analysis and synthesis. The two-dimensional decomposition approach is based on the property of separation of the functions into arbitrary x and y directions. The first step is identical to the one-dimensional approach, however, instead of keeping the low-level resolution and processing the high level resolution, both are processed using two identical filter bank after a transposition of the incoming data. Thus, the image is scanned in both horizontal and vertical directions. This result in an average image (or subimage) and three detail images generated by the following 2-D scaling function φ(x,y) φ(x)φ(y) and the vertical, the horizontal and the diagonal wavelet functions: ψ y)(x, φ(x)ψ(y) <sup>1</sup> , ψ y)(x, ψ(x)φ(y) <sup>2</sup> and ψ y)(x, ψ(x)ψ(y) <sup>3</sup> , respectively. To recover the original image, the inverse process is applied. Figure 12 illustrates the analysis and synthesis stages built using three filter banks each.

Fig. 12. Two-dimensional Mallat's analysis and synthesis tree

In this case, the frequency band is halved at each stage by a factor of four as represented by Figure 13.

The Wavelet Transform for Image Processing Applications 409

Unlike the three previous methodologies, the lifting scheme follows another philosophy. The fact is that the Fourier theory is not involved anymore and the construction of any wavelet system lies only in the spatial domain. If the explanation of the theory relies on the works of Sweldens (Sweldens, 1995, 1996 & Valens, 2004) the lifting approach has links with many other schemes (Burrus et al., 1998; Do & Vetterli, 2003, 2005; Cunha et al., 2006; Lu & Do, 2007; Nguyen & Oraintara, 2008 & Brislawn, 2010). The lifting-based wavelet transform can be seen as a succession of three operations: split, predict and update. In the first operation, data is the split into even and odds parts (known also as the lazy wavelet transform). Then, differences or details are calculated through the usage of a predictor. Finally, to compute the average, the even part is updated using the details previously calculated. Figure 16 shows an analysis and synthesis lifting-based wavelet transform.

**Predict Update Merge**

**Even**

**Odd**

The reconstruction operation does exactly the same, but using the reverse process. The data is first predicted, then updated and finally merged. Figure 17 illustrates split and merge

+

**Update Predict**

+

Fig. 15. Frequency bands of Feauveau's Quincux decomposition


**7.4 Swelden's lifting scheme** 

**Split**

**Xo**

**Xe**

**X**


Fig. 16. Lifting-based Wavelet Transform

operations using the polyphase property (Fliege, 1994).

Fig. 13. Frequency Bands of Mallat's 2-D Analysis Algorithm

#### **7.3 Feauveau's non-dyadic structure**

Based on Adelson's work (Adelson et al., 1987), this approach has been introduced by Feauveau (Feauveau, 1990). This decomposition is also known as Quincux. It differs from Mallat's two-dimensional approach by the fact that only the decimated output from the low pass filter is transposed and then processed through a "similar" filter bank. The result is a low resolution average image along with two different detail images from two different resolution levels. The fact is that the decomposition is not dyadic and the initial resolution of a factor of 2 is replaced by a 2 factor leading to an asymmetrical support. Figure 14 shows an analysis and synthesis stage of a Quincux structure.

Fig. 14. Feauveau's analysis and synthesis tree

Due to the removal of the filter bank at the output of the high pass filter - as reported in (Starck et al., 1998) only a wavelet image is involved at each stage. Recently, this approach has been used in an image compression scheme and found to give often better overall performances than other approaches (Stromme, 1999; Ebrahimi et al, 2002; Smith, 2003; Hankerson et al., 2005; Xiong & Ramchandran, 2005; Nai-Xiang et al., 2006; Raviraj & Sanavullah, 2007 & Oppenheim & Schafer, 2010). The frequency bands of a Quincux analysis is shown in Figure 15.

Fig. 15. Frequency bands of Feauveau's Quincux decomposition

### **7.4 Swelden's lifting scheme**

408 Advances in Wavelet Theory and Their Applications in Engineering, Physics and Technology

Based on Adelson's work (Adelson et al., 1987), this approach has been introduced by Feauveau (Feauveau, 1990). This decomposition is also known as Quincux. It differs from Mallat's two-dimensional approach by the fact that only the decimated output from the low pass filter is transposed and then processed through a "similar" filter bank. The result is a low resolution average image along with two different detail images from two different resolution levels. The fact is that the decomposition is not dyadic and the initial resolution of a factor of 2 is replaced by a 2 factor leading to an asymmetrical support. Figure 14 shows

Along the Columns

2

2

*g* **~**

*h* **~** Along the rows

*h*<sup>2</sup> **~**

*g* **~**

+

**~**

2

+

Along the Columns

2

Im(x,y)<sup>2</sup> Im(x,y)

Due to the removal of the filter bank at the output of the high pass filter - as reported in (Starck et al., 1998) only a wavelet image is involved at each stage. Recently, this approach has been used in an image compression scheme and found to give often better overall performances than other approaches (Stromme, 1999; Ebrahimi et al, 2002; Smith, 2003; Hankerson et al., 2005; Xiong & Ramchandran, 2005; Nai-Xiang et al., 2006; Raviraj & Sanavullah, 2007 & Oppenheim & Schafer, 2010). The frequency bands of a Quincux analysis

*g*

*h*

Fig. 13. Frequency Bands of Mallat's 2-D Analysis Algorithm

an analysis and synthesis stage of a Quincux structure.

**7.3 Feauveau's non-dyadic structure** 

Along the rows

*h* 2

*g*

is shown in Figure 15.

2

Fig. 14. Feauveau's analysis and synthesis tree

Unlike the three previous methodologies, the lifting scheme follows another philosophy. The fact is that the Fourier theory is not involved anymore and the construction of any wavelet system lies only in the spatial domain. If the explanation of the theory relies on the works of Sweldens (Sweldens, 1995, 1996 & Valens, 2004) the lifting approach has links with many other schemes (Burrus et al., 1998; Do & Vetterli, 2003, 2005; Cunha et al., 2006; Lu & Do, 2007; Nguyen & Oraintara, 2008 & Brislawn, 2010). The lifting-based wavelet transform can be seen as a succession of three operations: split, predict and update. In the first operation, data is the split into even and odds parts (known also as the lazy wavelet transform). Then, differences or details are calculated through the usage of a predictor. Finally, to compute the average, the even part is updated using the details previously calculated. Figure 16 shows an analysis and synthesis lifting-based wavelet transform.

Fig. 16. Lifting-based Wavelet Transform

The reconstruction operation does exactly the same, but using the reverse process. The data is first predicted, then updated and finally merged. Figure 17 illustrates split and merge operations using the polyphase property (Fliege, 1994).

The Wavelet Transform for Image Processing Applications 411

The orthogonality between the wavelet coefficients and the scaling coefficients is then only a

0g(n)h(n)

The scaling coefficients, which satisfy equation (33), are called Quadrature Mirror Filters

To achieve perfect reconstruction, the analysed signal has to be identical to the synthesised

**2**

**2**

**2** *g*

*cj-1(k)*

*dj-1(k)*

Biorthogonal wavelet basis can be seen as a generalisation of the orthogonal wavelet basis where some imposed restrictions on the latter have been relaxed. Unlike the case of orthogonal basis, the scaling and the wavelet functions need be neither of the same length, nor even numbered. Hence, the quadrature mirror property is not applicable and is replaced with a dual property. For the perfect reconstruction equation to hold, the scaling and the

It is clear that when the analysis and the synthesis filters are similar, the system becomes

<sup>δ</sup>(k)2k)h(n(n)h~

Previously, in orthogonal basis, only the analysis scaling coefficients (or wavelet coefficients) along with their shifted versions were used. In biorthogonal case, the analysing scaling coefficients are kept unchanged, while their shifted versions are replaced by the shifted versions of the synthesis dual filter. In other words, the analysis filter is orthogonal

to its synthesis dual filter. The biorthogonal denomination comes from this feature.

(33)

*~*

*h ~*

n)h(11)((n)g~ <sup>n</sup> (34)

n)(1h~ 1)(g(n) <sup>n</sup> (35)

(36)

j are the input and the output of a

+

*cj (m) ~*

n

one. In other words, (n)c~ (n)c jj , where (n)cj and (n)c~

*g*

Fig. 18. Two-band analysis and synthesis filter bank

wavelet coefficients have to fulfil the following equations:

orthogonal. The "orthogonality" condition in this case is defined by:

n

*h*

two-band filter (or filter bank) as shown in Figure 18, respectively.

**2**

simple implication:

*cj (n)*

**8.2 Biorthogonal basis** 

(QMF).

Fig. 17. Lazy Wavelet Transform: (a) Split, (b) Merge

## **8. The Wavelet Transform revisited**

In many practical problems, both the orthonormal basis (Daubechies, 1988, 1992, 1993) and the biorthogonal basis (Cody, 1994) can be used. The two bases (or families) present similarities and differences. Another scheme, called wavelet packet, which involves either orthonormal basis or biorthogonal basis is also possible (Wickerhauser, 1994). The following briefly describes the main features of orthonormal and biorthogonal bases together with extension to the wavelet packet scheme. It is worth mentioning that other schemes like undecimated wavelet, adaptive wavelets and multiwavelets exist and are beyond the scope of this brief overview.

#### **8.1 Orthonormal basis**

The orthonormal basis emerged from the work initiated by Mallat and Daubechies (Mallat, 1989 & Daubechies, 1988, 1993). The orthonormality property is somewhat seen as the discrete version of the orthogonality property (Masud, 1999). However, the basis functions are further normalised. These concepts have been mentioned when the multiresolution feature and the scaling function have been introduced. The admissibility and the orthogonality conditions ensure the existence and the orthogonality feature of the scaling function, defined by equation (18). This is achieved if:

$$\sum\_{\mathbf{n}} \mathbf{h}(\mathbf{n}) = \sqrt{2} \tag{30}$$

And

$$\sum\_{\mathbf{n}} \mathbf{h}(\mathbf{n}) \mathbf{h}(\mathbf{n} + 2\mathbf{k}) = \mathbf{\hat{S}}(\mathbf{k}) \tag{31}$$

Furthermore, using the two equations above alongside with equation (23), which defines the wavelet function, the orthogonality of the scaling function and the wavelet function at the same scale can be derived. This can be achieved only if the following equality is verified:

$$\mathbf{g(n) = (-1)^n h(1 - n)}\tag{32}$$

.

**Even**

**Odd**

(a) (b)

In many practical problems, both the orthonormal basis (Daubechies, 1988, 1992, 1993) and the biorthogonal basis (Cody, 1994) can be used. The two bases (or families) present similarities and differences. Another scheme, called wavelet packet, which involves either orthonormal basis or biorthogonal basis is also possible (Wickerhauser, 1994). The following briefly describes the main features of orthonormal and biorthogonal bases together with extension to the wavelet packet scheme. It is worth mentioning that other schemes like undecimated wavelet, adaptive wavelets and multiwavelets exist and are beyond the scope

The orthonormal basis emerged from the work initiated by Mallat and Daubechies (Mallat, 1989 & Daubechies, 1988, 1993). The orthonormality property is somewhat seen as the discrete version of the orthogonality property (Masud, 1999). However, the basis functions are further normalised. These concepts have been mentioned when the multiresolution feature and the scaling function have been introduced. The admissibility and the orthogonality conditions ensure the existence and the orthogonality feature of the scaling

2h(n)

δ(k)2k)h(nh(n)

Furthermore, using the two equations above alongside with equation (23), which defines the wavelet function, the orthogonality of the scaling function and the wavelet function at the same scale can be derived. This can be achieved only if the following equality is verified:

n

n

**2**

**2**

**Z-1**

+

(30)

(31)

n)h(11)(g(n) <sup>n</sup> (32)

**2**

**Z 2**

Fig. 17. Lazy Wavelet Transform: (a) Split, (b) Merge

function, defined by equation (18). This is achieved if:

**8. The Wavelet Transform revisited** 

of this brief overview.

**8.1 Orthonormal basis** 

And

The orthogonality between the wavelet coefficients and the scaling coefficients is then only a simple implication:

$$\sum\_{\mathbf{n}} \mathbf{h}(\mathbf{n}) \mathbf{g}(\mathbf{n}) = 0 \tag{33}$$

The scaling coefficients, which satisfy equation (33), are called Quadrature Mirror Filters (QMF).

To achieve perfect reconstruction, the analysed signal has to be identical to the synthesised one. In other words, (n)c~ (n)c jj , where (n)cj and (n)c~ j are the input and the output of a two-band filter (or filter bank) as shown in Figure 18, respectively.

Fig. 18. Two-band analysis and synthesis filter bank

#### **8.2 Biorthogonal basis**

Biorthogonal wavelet basis can be seen as a generalisation of the orthogonal wavelet basis where some imposed restrictions on the latter have been relaxed. Unlike the case of orthogonal basis, the scaling and the wavelet functions need be neither of the same length, nor even numbered. Hence, the quadrature mirror property is not applicable and is replaced with a dual property. For the perfect reconstruction equation to hold, the scaling and the wavelet coefficients have to fulfil the following equations:

$$\widetilde{\mathbf{g}}(\mathbf{n}) = (-1)^{\mathbf{n}} \mathbf{h}(1 - \mathbf{n}) \tag{34}$$

$$\log(\mathbf{n}) = (-1)^{\mathbf{n}} \widetilde{\mathbf{h}} (\mathbf{1} - \mathbf{n}) \tag{35}$$

It is clear that when the analysis and the synthesis filters are similar, the system becomes orthogonal. The "orthogonality" condition in this case is defined by:

$$\sum\_{\mathbf{n}} \widetilde{\mathbf{h}}(\mathbf{n}) \mathbf{h}(\mathbf{n} + 2\mathbf{k}) = \mathbf{\hat{s}}(\mathbf{k}) \tag{36}$$

Previously, in orthogonal basis, only the analysis scaling coefficients (or wavelet coefficients) along with their shifted versions were used. In biorthogonal case, the analysing scaling coefficients are kept unchanged, while their shifted versions are replaced by the shifted versions of the synthesis dual filter. In other words, the analysis filter is orthogonal to its synthesis dual filter. The biorthogonal denomination comes from this feature.

The Wavelet Transform for Image Processing Applications 413

In comparison to classical wavelet approach, the wavelet packet scheme presents the

 Possibility of using different wavelet from a level to another. This strategy has been used in (Masud, 1999) to implement a two-level orthonormal wavelet packet and a

 Possibility of choosing a particular wavelet packet decomposition from the general generic structure of Figure 19. Thus, one can choose either to preserve the orthonormality feature of the decomposition (Wickerhauser, 1994), or highlight the peculiarities of the signal (Masud, 1999). A binary search for the best decomposition

However, there is a cost to be paid. In this case, the computational complexity of a wavelet packet structure is O(nlog(n)) in contrast to the O(n) of the classical wavelet transform.

Recently, The wavelet transform is being increasingly used, not only in the field of image and signal processing applications but also in many other different areas, ranging from mathematics, physics, astronomy to statistics and economics. In image processing based applications, image compression, image denoising and image watermarking are at the cutting edge, and as such, a brief description of these wavelet-based applications is given in the following subsections (Strang & Nguyen, 1996; Burrus et al., 1998; Stromme, 1999; Ebrahimi et al, 2002; Nibouche et al., 2000, 2001a, 2001b, 2001c, 2001d, 2002, 2003; Smith, 2003; Do & Vetterli, 2003, 2005; Hankerson et al., 2005; Nai-Xiang Yap-Peng, 2005; Xiong & Ramchandran, 2005; Chappelier & Guillemot, 2006; Cunha et al., 2006; Nai-Xiang et al., 2006; Raviraj & Sanavullah, 2007; Hernandez-Guzmane et al., 2008; Firoiu et al., 2009; Mallat, 2009; Brislawn,

2010; Oppenheim & Schafer, 2010; Ruikar & Doye, 2010 & Chen & Qian, 2011).

following features (Daubechies, 1992):

three-level biorthogonal wavelet packet.

tree is also possible (Burrus et al., 1998).

**9. Wavelet-based applications** 

Fig. 20. Two-band analysis and synthesis filter bank

At the expense of the energy partitioning property stated by Perseval's equality, which is a direct consequence of the lack of orthogonality, a greater flexibility can be achieved by using the basis and dual basis (Burrus et al., 1998). One of the most "important" features in the biorthogonal basis is the linear phase property, which leads to the filter coefficients (when implementing a wavelet system) being symmetric. In addition, the difference of length between dual filters must be even, leading either to odd or even length of the low pass and the high pass filters. In general, biorthogonal wavelet systems present the following features (Daubechies, 1992):


#### **8.3 Wavelet packets**

In contrast to the "traditional" Mallat's decomposition, which leads to narrow frequency bandwidths (low frequencies) and wide frequency bandwidths (high frequencies), the wavelet packet approach emerged first as a way of adjusting high frequency resolutions. Hence, the Mallat's decomposition scheme is applied to both parts of a filter bank leading to the split of frequencies in progressive finer resolutions. The generic structure of wavelet packet decomposition is shown in Figure 19 and the frequency bandwidths illustrated by

Figure 20. In this scheme, the number of filters increases by a factor of )2(2 ji at each successive subband, where i and j represent two consecutive resolutions and 1ji .

Fig. 19. Three-stage Wavelet Packet Decomposition

At the expense of the energy partitioning property stated by Perseval's equality, which is a direct consequence of the lack of orthogonality, a greater flexibility can be achieved by using the basis and dual basis (Burrus et al., 1998). One of the most "important" features in the biorthogonal basis is the linear phase property, which leads to the filter coefficients (when implementing a wavelet system) being symmetric. In addition, the difference of length between dual filters must be even, leading either to odd or even length of the low pass and the high pass filters. In general, biorthogonal wavelet systems present the following features

The low pass and the high pass filters used in the filter bank have not the same length;

In contrast to the "traditional" Mallat's decomposition, which leads to narrow frequency bandwidths (low frequencies) and wide frequency bandwidths (high frequencies), the wavelet packet approach emerged first as a way of adjusting high frequency resolutions. Hence, the Mallat's decomposition scheme is applied to both parts of a filter bank leading to the split of frequencies in progressive finer resolutions. The generic structure of wavelet packet decomposition is shown in Figure 19 and the frequency bandwidths illustrated by

ji at each

Figure 20. In this scheme, the number of filters increases by a factor of )2(2

successive subband, where i and j represent two consecutive resolutions and 1ji .

*g* **2**

.

*g* **2**

*h* **2**

*g* **2**

*h* **2**

*g* **2**

*h* **2**

*g* **2**

*h* **2**

*h* **2**

*g* **2**

*h* **2**

 The coefficients of the filters are either real numbers or integers; The filters in this family present either even or odd orders;

The high pass filter is either symmetric or antisymmetric.

The low pass filter is always symmetric;

*g*

**2**

**2**

*h*

Fig. 19. Three-stage Wavelet Packet Decomposition

(Daubechies, 1992):

**8.3 Wavelet packets** 

In comparison to classical wavelet approach, the wavelet packet scheme presents the following features (Daubechies, 1992):


However, there is a cost to be paid. In this case, the computational complexity of a wavelet packet structure is O(nlog(n)) in contrast to the O(n) of the classical wavelet transform.

Fig. 20. Two-band analysis and synthesis filter bank
