Two-Dimensional Discrete Wavelet Transform (Biomedical Image Analysis)

Up to this point, we have covered the one-dimensional discrete wavelet transform. In image processing, a multidimensional wavelet transform is desired. Fortunately, the wavelet transform is a linear operation. Therefore, analogous to the Fourier transform, the two-dimensional wavelet transform can be performed by first computing the row-by-row one-dimensional wavelet transform in the horizontal direction, followed by the column-by-column one-dimensional wavelet transform in the vertical direction. This principle can be extended toward any number of dimensions. As a consequence of the decomposition scheme, however, each subband filter generates four image regions with half of the side length. By convention, the image regions are arranged as depicted in Figure 4.7, and regions with reduced detail in either the x or y direction (L/L, L/H, H/L) are subdivided further in the subsequent filter stage (dashed lines).

FIGURE 4.7 Output of a two-dimensional subband decomposition stage. The image is subdivided into four regions: a region with high detail (H/H), a region with high detail in the x-direction (H/L) but low detail in the y-direction, a region with high detail in the y-direction (L/H), and a fourth region with low detail (L/L).

FIGURE 4.8 Example of the pyramidal decomposition of an image. The original image (A) is a noisy MR slice of the abdomen. The first-stage decomposition is shown in (B), where a subdivision into the four regions shown in Figure 4.7 has taken place. In subsequent decomposition steps (C), all regions except region H/H are subdivided further. Note that for visualization purposes, regions H/L, H/H, and L/H have been contrast-enhanced and white divider lines have been added.

Figure 4.8 shows two steps in the decomposition of a medical image, in this example a MRI slice of the abdomen with artificially added Gaussian noise. In the first decomposition step, the image is subdivided into four regions as described in Figure 4.7. In the next subdivision step, regions H/L and L/H are subdivided into two regions, and region L/L is subdivided into four. Further subdivisions of region L/L are possible in subsequent decomposition steps. From this image, the potential of the wavelet decomposition for filtering purposes becomes clear. Region H/H contains predominantly noise, whereas region L/L appears to contain the major parts of the image information. Wavelet-based filter functions are covered in Section 4.3. In addition, the wavelet transform gives rise to lossy image compression: It would appear as if region L/L in Figure 4.8B represents the original image quite well, although it contains only one-fourth of the image area. In fact, wavelet image compression is now being used widely. Wavelet-based image compression is covered in detail in Section 12.3.

Computer Implementation of the Two-Dimensional Wavelet Transform We now focus on implementation of the discrete wavelet transform. For this purpose, an algorithmic representation of Equations (4.16) and (4.17) needs to be formulated. A suitable approach is to implement a one-dimensional subband coding filter and use this filter for higher-level decomposition functions. The convolution described in Equation (4.16) can be viewed as a matrix multiplication of the data vector x with a convolution matrix A, resulting in the output data vector y, where x corresponds to \(m) in Equation (4.16) and y corresponds to the vector containing the next stage’s X(k) and 7 (k):

The matrix A contains, in alternating rows, the parameters of the scaling and wavelet filters:

Note that the matrix coefficients start to wrap around in the last lines, depending on the number of wavelet parameters. A possible implementation is provided in Algorithm 4.1. By using two different indices for input and output data, it becomes possible to perform the convolution step and the rearrangement step (see Figure 4.6) in the same loop, and the output vector will contain the X(k) in the lower half and the 7 (k) in the upper half. From Algorithm 4.1 it is straightforward to derive the subband decoding filter that performs one step of the inverse wavelet transform. Algorithm 4.2 performs one such step according to Equation (4.17). The input data vector must be organized like the output vector of Algorithm 4.1: that is, the X(k) reside in the lower half and the 7 (k) reside in the upper half. The output vector is one contiguous data vector. Algorithms 4.1 and 4.2 build on a table of filter coefficients. Since the parameters g(k) and h(k) are mutually dependent on each other through Equation (4.11), it is generally sufficient to store one set of filter coefficients (usually, the scale filter coefficients) and precede the wavelet transform by a preparatory step, Algorithm 4.3, to populate both filter coefficient arrays.

Algorithm 4.1 Wavelet subband filter. Pseudocode for one stage of the subband decomposition filter described in Equation (4.16). The input data vector is stored in xdata[], and the output will be returned in ydata[]. The length of both data vectors is N. It is assumed that the wavelet filter parameters G[k] and the scale filter parameters H[k] are provided prior to using this algorithm. The number of parameters is L. For this algorithm, N must be an even number.

At this point, these algorithms can be combined to form a pyramidal wavelet decomposition and reconstruction. We begin by providing a one-dimensional decomposition, because it helps us to understand the two-dimensional decomposition. In Algorithm 4.4, a one-dimensional discrete wavelet transform is performed building on Algorithm 4.3, here referred to as dwt.prep, and Algorithm 4.1, referred to as dwt.decompose. With Algorithm 4.4 performing a full one-dimensional pyramidal decomposition, it is now possible to extend the wavelet transform toward two-dimensional data (i.e., images). Since the discrete wavelet transform is a linear operation, a two-dimensional discrete wavelet transform can be implemented as a one-dimensional wavelet transform in the x-direction followed by a one-dimensional wavelet transform in the y-direction. This strategy is reflected in Algorithm 4.5, which builds on the one-dimensional wavelet transform in Algorithm 4.4, here referred to as dwt_1d. When implementing Algorithms 4.4 and 4.5, care should be taken that dwt_prep, the preparation of the filter parameters, is executed only once, that is, in Algorithm 4.5 and not in Algorithm 4.4.

Algorithm 4.2 Wavelet reconstruction subband filter. Pseudocode for one stage of the subband reconstruction filter described in Equation (4.17). The input data vector is stored in xdata[] and must contain the X(k) in the lower half and the 7 (k) in the upper half. The output will be returned in ydata[] .The length ofboth data vectors is N. It is assumed that the wavelet filter parameters G[k] and the scale filter parameters H[k] are provided to this algorithm in a fashion identical to Algorithm 4.1, with the number of parameters being L. For this algorithm, N must be an even number. The calling function is responsible for providing the ydata array.

Algorithm 4.3 Wavelet filter preparation. Pseudocode for the preparation of the filter parameter arrays G[k] and H[k] for one specific example, the Daubechies-4 wavelet (see Table 4.1). This algorithm computes the wavelet filter parameters g(k) from the scale filter parameters h(k) and sets the filter length L that is required by Algorithms 4.1 and 4.2.

Algorithms 4.1 through 4.5 combined perform a two-dimensional wavelet transform on an image with a freely selectable wavelet and a selectable number of decomposition levels. The inverse wavelet transform is implemented in a straightforward manner by substituting the call to the wavelet decomposition (dwt_decompose, Algorithm 4.1) by a call to Algorithm 4.2.