Wavelets, the "little waves" of signal processing, came to the fore in the early 1990s as an attractive alternative to classical Fourier Transform based signal and image processing. While the underlying concepts behind wavelets have been known for close to a century (Haar described his multiresolution analysis in 1910), the pioneering work of many applied mathematicians have brought new insights into the field and wavelet theory and applications have found use in diverse areas like geophysical signal processing, medical imaging, and information theory. In this topic, we look at wavelets from the image processing perspective, and develop wavelet-based algorithms for two operations previously studied in this topic. The topic of edge detection was covered in 5.1, and in the first part of this topic we develop a wavelet-based edge detection algorithm, and of course implement it on the DSP. The enhancement of images through various filtering schemes was the focus of next topic, where we utilized linear, non-linear, and adaptive filters to denoise images. In the second part of this topic, the topic of wavelet denoising is introduced and a fixed-point implementation tested on the C6416 is discussed.
Prior to embarking on this abbreviated tour of wavelet-based image processing, we must introduce at least a bare minimum of theory so that the algorithms and code make some sense to the reader. If the following discussion in 6.1 feels too abstract or mathematical in nature, I urge the reader to not get discouraged by its tone – it is not intended to be a treatise on wavelet theory, as this subject has been covered in depth by authors far more qualified than I. A large portion of the mathematical and signal processing background is skipped over in an attempt in an effort to convey the general idea behind the wavelet transform, and introduce enough of the theory to develop wavelet-based models of the two aforementioned image processing algorithms. A full treatment of wavelets easily warrants its own topic, and there are many. More complete and mathematically rigorous discussions of wavelet theory can be found in [1-4].
Transform-based methods are of fundamental importance in signal and image processing. A transform is a mathematical device that takes a function or signal represented in one domain and converts it to a form based in another domain. For example, in 4.1.3, we briefly glossed over the Fourier transform, the most common transform in signal processing applications. The forward Fourier transform can be said to take a signal in the "time-domain" to the "frequency-domain". Put more concretely, a one-dimensional time-domain signal, if plotted, would have its x-axis labeled as something akin to time, and in the two-dimensional (image) case, such a signal’s axes might be labeled as position or pixel location. The Fourier transformed signal allows us to change these labels from time to frequency. Quite often it is the case that a transformed signal highlights certain characteristics we may be interested in analyzing, that are not readily apparent in its original form -such behavior is the essence of transform-based signal/image processing.
A transform is said to be invertible if an inverse transform can be used to recover the original signal from the transformed signal. The Fourier transform is invertible because an inverse transform exists that takes a signal in the frequency-domain back to its original time-domain representation. In general, one-dimensional linear transforms are of the form
Here, instead of viewing y = f(x) in the usual Cartesian sense of a sequence of y values plotted against an x-axis, we express f(x) as a linear decomposition of a weighted sum of basis functions
In the Fourier Transform, these basis functions are complex sinusoids – it has been shown that any signal can be expressed as a (possibly infinite) sum of sinusoids of different frequencies. The transform coefficients au are calculated from an inner product off(x) and
The key point for this discussion is that the basis functions for the classical linear transform y/u (.xj have infinite extent, or in other words they are non-zero for all x. Again using the Fourier transform as an example, this means the Fourier coefficients au depend on all parts of f(x), even as |x| -> oo, Thus, we have essentially lost all information about the time domain. So if we were to, say, analyze an audio signal consisting of a digitized symphonic work using the Fourier transform, we would obtain excellent frequency resolution and would be able to determine which frequencies (notes) were prevalent in that work. What we would not know, however, is when in time those notes occurred. This lack of localization with the Fourier and other related transforms is a major drawback, and is partly what led mathematicians to explore wavelet theory.
They found that in certain respects these classical linear transforms were not the ideal mathematical tool for analyzing signals whose behavior changes with time, or in the case of an image, a two-dimensional signal whose behavior changes with spatial location. Such signals are known as non-stationary, and most interesting signals encountered in real life exhibit this characteristic. While there are some means of circumventing the lack of localization in the Fourier transform (e.g., the so-called Short Time Fourier Transform or STFT, which uses overlapping windows and multiple FFTs), it eventually became clear that wavelets offer an extremely attractive alternative for analyzing non-stationary signals, overcoming many of the drawbacks associated with the Fourier transform.
In one dimension, the wavelet transform uses a two-parameter system that gives us the means to decompose a signal into various frequency bands and cut these bands into slices in time. The wavelet transform can be thought of as a specialized case of subband filtering, and the discrete wavelet transform (DWT) is given by
where k and j are integer indices, and
are the wavelet basis functions.
A wavelet expansion of a function fix) yields location information in both time and scale simultaneously, where scale can be related to frequency. Drawing on the symphony example, if fix) is the digitized audio signal, the DWT of that signal is akin to a musical score, because it tells us what frequencies were played and when those tones occurred. This, in a nutshell, is the overwhelming advantage the wavelet transform offers.
We proceed by first developing an understanding of how to implement the one-dimensional DWT, as the two-dimensional DWT – which is of most significance to us – follows naturally. The DWT is a separable transform, meaning we can implement the two-dimensional transform by performing two one-dimensional transforms in series. The manner in which this is typically done is by first transforming all the rows of an image using the one-dimensional DWT, followed by transforming the columns of the image, using the coefficients obtained from the row transforms.
Calculation of the DWT coefficients is done using a filter bank, which is a series of cascading digital filters. In 4.1.2, we described how digital filtering is accomplished via convolution, and in subband filtering the same concept holds. However, in subband filtering we use multiple filters (hence the term "bank") to decompose the input into its constituent parts. Implementing the DWT using a filter bank entails passing the signal through two digital filters: a low-pass or moving average filter (e.g., the box filter used in Algorithm 5-2 to smooth the image histogram) and another complementary high-pass, or moving difference, filter. The dual-channel subband filter that decomposes the input signal into its wavelet coefficients is the "analysis" filter bank and is one way of expressing the forward DWT. In the analysis phase of the DWT, the input signal is passed through both filters and this filtered output is then decimated, or downsampled, by a factor of two. Decimation is the process by which we take every other sample from an input to form the output, i.e. y[n] = x[2n].
The inverse DWT (IDWT) reconstitutes the decomposed signal by passing the wavelet coefficients through a "synthesis" filter bank. The wavelet coefficients are upsampled by a factor of two (zeroes are inserted between every other sample) and then passed through time-reversed low-pass and high-pass filters. These two outputs are added together to reconstruct the original signal. If we let the symbols *h denote convolution of a signal with a filter h, and i m decimation by a factor of m, then the structure of a generalized one-dimensional two-channel perfect reconstruction orthogonal filter bank is given in Figure 6-1.
In an orthogonal filter bank, the synthesis filters gHp and gu> are time-reversed versions of the analysis filters. That is, if /z[n] is a low-pass or high-pass filter, with N the number of filter coefficients in h, then g[n] = /?[N-l-n] is the corresponding "flipped" synthesis kernel. Moreover, in an orthogonal filter bank, the analysis high-pass filter hHp is related to the low-pass analysis filter hLP in that hHP is a time-reversed alternating negation of the low-pass kernel. In mathematical terms, we form the high-pass filter from the low-pass filter according to the following relation:
where N is the number of taps in the low-pass filter. Thus, for this particular class of filter banks, the entire bank is completely defined by a single filter -the low-pass analysis filter. To summarize, if the length N of an analysis low-pass filter is 4, and
then the rest of the filters in the perfect reconstruction orthogonal filter bank are defined by the following kernels:
The obvious next question is where does the low-pass analysis filter hL? come from? Essentially, there is any number of choices and the exact choice comes down to a filter design problem, with some global constraints imposed on the design space. The analysis filters’ impulse response is directly tied to the mother wavelet, which as we stated must be of limited duration and furthermore have zero mean. Even with these constraints the design space is vast, with various other factors coming into play, such as smoothness, length of the filter’s impulse response, vanishing moments (differentiability, i.e. the more vanishing moments in the wavelet the more differentiable it is), and so on. The derivation of hLP from a mother wavelet can be found in the references, and in the remainder of this topic we will use a few of the more common filter banks, in particular the Haar and Daubechies systems.
Figure 6-1. Single level, one-dimensional perfect reconstruction orthogonal filter bank. Taken together, c: (approximation) and dl (detail) are the complete DWT coefficients au. Each of these vectors is one-half the length of the input signal X.
The left-hand side of Figure 6-1 depicted a single stage of a wavelet decomposition, yielding the wavelet coefficients c, and d, at a single scale. We can cascade this filter in a few different fashions to carry the decomposition one level further, with one quite common method known as the pyramid decomposition structure. The pyramid structure is shown in Figure 6-2, where the approximation wavelet coefficients C\ are fed into the same two-channel filter bank, which in turn emits another set of approximation and detail coefficients. Thus as the cascade proceeds in a dyadic (power-of-two) fashion, the DWT effectively tiles the time-frequency space. This tiling is known as a multi-resolution analysis (MRA), and really gets to the heart of the primary advantage of the wavelet transform – it solves the localization problem of the Fourier and other classical transforms. The MRA decomposes the signal into nested subspaces, thereby giving us the ability to know not only if a particular event or characteristic occurred, but when. Using this formulation of the wavelet decomposition, the output of the analysis cascade is arranged as shown in Figure 6-3, using the example of a three-level wavelet decomposition. The synthesis filter bank reconstitutes the input signal using a similar structure to that of Figure 6-2, but with different filters and the downsampling operation replaced with upsampling of the DWT coefficients (see Figure 6-4).
Another less common decomposition structure, not utilized in the algorithms presented in this topic, is the wavelet packet method. In this formulation, both the approximation and detail coefficients are decomposed at each cascade level, with the end result a binary tree (in one dimension). Wavelet packet analysis does offer a richer signal analysis, although at the cost of some complexity.
Now that we have introduced the mechanics behind the wavelet transform, one can now envision how a wavelet-based processing system might work. In a wavelet-based image processing methodology, one can envision first performing the DWT, and then subsequently manipulating the DWT coefficients. After the modifications have been carried out, the altered DWT coefficients are fed into the synthesis filter bank, thereby yielding a processed image. In 6.3, we use this methodology to denoise images using the DWT.
Quadrature Mirror Filters and Implementing the 2D DWT in MATLAB
The oldest and simplest wavelet is the Haar system, defined by the following analysis filters:
Figure 6-2. Analysis phase of a two-level pyramid wavelet filter bank. This cascading bank of digital filters produces the DWT coefficients as the vectors {c2. d2, dt).
Figure 6-3. Typical arrangement of the output of a one-dimensional three-level DWT implemented via a pyramid cascade of digital filters.
Figure 6-4. Two-level pyramidal filter bank that produces in {c2, d2, dt) the two-level DWT.
The filter defined by hL? outputs a normalized moving average of its input, while htiP outputs a normalized moving difference. The Haar basis functions (the mother wavelet) are discontinuous, and for that reason are not used very often in real applications, except for some specialized video processing. However, due to its simplicity, it is instructive to illustrate its usage in the processing of images. The application of more sophisticated wavelets follows from the general idea behind the implementation of the Haar wavelet transform.
The filter bank tree structures shown in Figures 6-1, 6-2, and 6-4 are known as Quadrature Mirror Filters (QMF). In [1], it is shown that a QMF tree can be efficiently represented in matrix form. For the forward DWT, the convolution and decimation operations can be combined into a single matrix operation. In general, a polyphase matrix H is constructed from the analysis low-pass and analysis high-pass filter coefficients. Remembering that hH? is the time-reversed alternating negation of hLP, if h, are the analysis low-pass filter coefficients and the length of h is 4 (for the Haar wavelet the length is 2), then the format of the block-diagonal matrix H is:
The structure of H is such that the length of the output (transformed) signal is exactly the same as the input. Each successive row in H is translated by 2 samples, so that when we multiply the input signal by this forward transform matrix we achieve the downsampling operation. Moreover, multiplication of the input and the transform matrix yields a circular convolution, hence the wrap-around effect at the bottom rows – the input data is treated as a periodic sequence of samples. The forward transform matrix for the Haar system is given below – note that there is no wrap-around in this case because the Haar system consists of only two filter coefficients and we are dealing with dyadic scales:
These block-diagonal transform matrices are also known as block-Toeplitz matrices, and in MATLAB the Kronecker tensor product can be used to quickly construct such matrices in an elegant fashion. The built-in MATLAB function kron5 can be used to form a 32×32 forward Haar transform matrix with the following code:
If sj are the samples of a one-dimensional signal s, then we compute one level of the DWT by multiplying H by s:
The vector on the left-hand side of the above equation contains in c^j and t/i i the approximation and detail wavelet coefficients respectively, but in order to arrange the output in the form shown in Figure 6-3, we need to shuffle the coefficients around. We do this not because it is aesthetically pleasing to do so, but because it is required in order to implement a multilevel DWT. To rearrange the coefficients appropriately, we define the permutation matrix P as a reordered identity matrix that moves all odd-numbered rows to the top-half of the output matrix, and puts all even-numbered rows in the bottom half of the output matrix (reminiscent of the butterfly structure prevalent in the FFT algorithm). Using the definitions of the matrices S, H, and P, a single iteration through one of the cascade stages of the QMF tree is
Now, to perform one level of the DWT on an image S, we take advantage of the fact that the DWT is separable. We first transform the rows of S using the above equation, which holds for both a one-dimensional Nxl signal and a two-dimensional MxN signal (image). After applying successive one-dimensional DWTs on the rows of the image, we transform the columns of the image by transposing the matrices, as shown below:
A variety of MATLAB M-files illustrating how to implement the DWT and IDWT can be found in the Chap6\wavelet directory. Listing 6-1 are the contents of three of those, which collectively implement a single-level Haar wavelet decomposition.
Listing 6-1: MATLAB functions for performing the Haar 2D DWT on a square image.
Implementing the DWT using matrix operations is trivial in MATLAB, as the programming language inherently supports matrix multiplication and inversion (which as we shall soon see is important when performing the reverse transform). In haar_2d_dwt, the apostrophe character, which in MATLAB denotes matrix transposition, is used to effectively move the transform over to the columns of the image. This implementation of the Haar DWT is limited to square images where the length of all sides is a power-of-two, and a single level of the Haar wavelet decomposition of an image I may be visualized using imshow (haar_dwt (I) , []) as shown in Figure 6.5.
The 2D DWT, with its output arranged in the format illustrated graphically in Figure 6-5, effectively splits the input image into four bands: a low-pass filtered approximation "low-low" (LL) subimage, and three high-pass filtered detail subimages, "low-high" (LH), "high-low" (HL), and "high-high" (HH). These subbands are shown in Figure 6-6.
The wavelet approximation coefficients in the LL(1) subband are nothing more than a subsampled version of the input, while the three high-pass filtered subimages LH1and HH1′ accentuate the horizontal, vertical, and diagonal edges respectively. As one might suspect, we will use this fact to our advantage during implementation of a wavelet-based edge detection scheme in 6.2.2. To continue the pyramidal wavelet decomposition, we recursively perform the same DWT matrix operations on the LL subband, which for a two-level DWT produces an output image like that shown in Figure 6-7. We can continue this process until the LL(n) subimage is nothing more than a 2×2 square that is an average of the entire image. Listing 6-2 is a MATLAB function that performs the 2D Haar pyramid DWT on an image for a given number of decomposition levels. Example output from this function is shown in Figure 6-8.
Figure 6-5. One level of the Haar wavelet decomposition, (a) Original Lenna image, (b) The result from the haar_2d_dwt MATLAB function of Listing 6-1. This matrix, consisting of the 2D DWT coefficients, has been shifted, scaled, and contrast-stretched for optimal display using an 8-bit dynamic range. The upper-left quadrant contains the approximation coefficients, while the other three quadrants are the various detail coefficients.
Figure 6-6. Decomposing an image into four subbands using the 2D DWT.
Figure 6-7. Two-level 2D pyramid DWT.
Figure 6-8. Two-level Haar wavelet decomposition, (a) Original Barbara image, (b) The result from the haar_2d_dwt_n MATLAB function. This 2D DWT coefficient matrix has been contrast-stretched for an 8-bit dynamic range, and artificial gray border lines delineate the seven subbands.
Listing 6-2: MATLAB function for performing a multi-level Haar 2D DWT on a square image.
Implementing the 2D IDWT proceeds along similar lines to that shown in Figure 6.4. Using the definition of outpuf) in terms of the H and P matrices, and if S is now treated as a matrix of wavelet coefficients, the IDWT is calculated as
To compute the IDWT from an N-level DWT, one must proceed up the filter tree by first reconstituting LL(n), then LL(|1"1), and so on all the way up to LL(1). This iterative procedure is illustrated in haar_2d_idwt_n, shown in Listing 6-3.
Listing 6-3: MATLAB function for performing a multi-level Haar 2D IDWT on a square image.
