Image Compression:
How Math Led to the JPEG2000 Standard
General Wavelet Transformations
Introduction
If you read the subsection on the Haar Wavelet Transformation, then you have a basic idea of how a discrete wavelet transformation works and why we are interested in using it for applications such as image compression. You also learned that a drawback to the HWT is the fact that it decouples data. For example, if we apply the HWT to the list (1000, 1002, 10,12), we obtain \sqrt{2} (1001, 11 | 1, 1). Looking at the differences in the output, we might believe that there is not much change in the input data. The problem is the Haar filter is too short - since it only processes two values at a time, it completely missed the large jump between the second and third values.
One of the first approaches to produce better wavelet filters is due to Princeton mathematician/physicist Ingrid Daubechies. In a 1988 paper, she developed a family of lowpass filters (even-length) that can be used to construct an orthogonal discrete wavelet transformation matrix. Moreover, the highpass filters can be constructed from the lowpass filters. These filters were superior to the Haar filter (the length 2 filter in Daubechies family of filters) in many applications, but for applications such as image compression, the filters were not optimal - the filters in Daubechies family are not symmetric.
In 1992 Daubechies, along with collaborators A. Cohen and J.-C. Feauveau, developed a new family of filters called biorthogonal filters. These filters no longer produced orthogonal matrices for discrete wavelet transformations, but the filters were symmetric and could modified to map integers to integers. We will see that both of these properties are desirable for image compression applications.
In this section, we will build the length-4 Daubechies filters, discuss the construction of longer length filters from her family. Finally, we will motivate the construction of biorthogonal filters.
Daubechies Length-4 Filter
Our goal is to construct a length-4 lowpass filter that can be used to form an orthogonal discrete wavelet transformation matrix. Let's call the filter {\bf h} = (h_0, h_1, h_2, h_3) and form the associated Fourier series
H(\omega) = h_0 + h_1 e^{i\omega} + h_2 e^{2i\omega} + h_3 e^{3i\omega}
Recall from the subsection on the HWT, we want the Fourier series of a lowpass filter to satisfy the two properties H(0) = \sqrt{2} and H(\pi) = 0. If we plug \omega = 0 into the Fourier series above we obtain
\sqrt{2} = H(0) = h_0 + h_1 + h_2 + h_3
and if we put \omega = \pi in the Fourier series we get
0 = H(\pi) = h_0 + h_1 e^{i\pi} + h_2 e^{2i\pi} + h_3 e^{3i\pi} = h_0 - h_1 + h_2 - h_3
since for any integer k, e^{ik\pi} = \cos(k\pi) + i \sin(k\pi) = \cos(k\pi) = (-1)^k. Thus we have two linear equations our filter coefficients must satisfy in order to be a lowpass filter. As you may know from linear algebra, linear systems with more variables than equations possess an infinite amount of solutions, so we will build some more equations to better describe our filter.
Suppose we wished to process a vector of length 8. If we were to maintain the same structure as the Haar transformation matrix, then the top half of our DWT using the length-4 filter might look like
H = \left[\matrix{ h_0 & h_1 & h_2 & h_3 & 0 & 0 & 0 & 0 \\ 0 & 0 & h_0 & h_1 & h_2 & h_3 & 0 & 0 \\ 0 & 0 & 0 & 0 & h_0 & h_1 & h_2 & h_3 \\ h_2 & h_3 & 0 & 0 & 0 & 0 & h_0 & h_1} \right]
Since we want an orthogonal matrix for our transform, each row must dot with itself to produce 1 and dot with every other row to produce 0. You can check that these properties lead to the following two equations:
\begin{eqnarray} h_0^2 + h_1^2 + h_2^2 + h_3^2 &=& 1 \\
h_0 h_2 + h_1 h_3 &=& 0
\end{eqnarray}
Thus we have the following system of equations:
\begin{eqnarray}
h_0 + h_1 + h_2 + h_3 &=& \sqrt{2} \\
h_0 - h_1 + h_2 - h_3 &=& 0 \\
h_0^2 + h_1^2 + h_2^2 + h_3^2 &=& 1 \\
h_0 h_2 + h_1 h_3 &=& 0
\end{eqnarray}
Unfortunately, this system is no longer linear since it contains two quadratic equations. But the good news is we can solve it! As a matter of fact, it can be shown that if we can satisfy the last three equations, then we must have h_0 + h_1 + h_2 + h_3 = \pm \sqrt{2}. It can also be shown that there are an infinite number of solutions to the last three equations.
Since our system still has an infinite number of solutions, let's add another equation and see if it leads to a system with a finite number of solutions. Daubechies' idea was to further flatten the Fourier series at \omega = \pi. We can use ideas from Calculus I to perform this task!! Indeed, if we insist that the derivative of the Fourier series at \omega = \pi were zero, then the graph of the absolute value of H(\omega) would approach 0 at \omega = \pi with a horizontal tangent. So we insist that H^\prime(\pi)=0. We can easily compute the derivative of the Fourier series:
H^{\prime}(\omega) = i h_1 e^{i\omega} + 2i h_2 e^{2i\omega} + 3i h_3 e^{3i\omega}
If we plug in \omega = \pi and recall that e^{ik\pi} = (-1)^k, we see that
0 = H^\prime(\pi) = -i h_1 + 2i h_2 - 3i h_3
Dividing both sides by -i gives us h_1 - 2h_2 + 3h_3 = 0. Thus we have the new system
\begin{eqnarray}
h_0 + h_1 + h_2 + h_3 &=& \sqrt{2} \\
h_0 - h_1 + h_2 - h_3 &=& 0 \\
h_1 -2h_2 + 3h_3 &=& 0 \\
h_0^2 + h_1^2 + h_2^2 + h_3^2 &=& 1 \\
h_0 h_2 + h_1 h_3 &=& 0
\end{eqnarray}
It turns out that this system has two solutions (for more information about solving this system, please see Discrete Wavelet Transformations: An Elementary Approach with Applications). The two solutions are reflections of each other. That is, if we find (a, b, c, d) that solves the system, then the other solution is (d, c, b, a). One of the solutions to the system is
\begin{eqnarray}
h_0 &=& (1 + \sqrt{3})/(4\sqrt{2}) \approx 0.482963 \\
h_1 &=& (3 + \sqrt{3})/(4\sqrt{2}) \approx 0.836516\\
h_2 &=& (3 - \sqrt{3})/(4\sqrt{2}) \approx 0.224144\\
h_3 &=& (1 - \sqrt{3})/(4\sqrt{2}) \approx -0.12941
\end{eqnarray}
We can construct the highpass filter from the lowpass filter. Indeed if we take {\bf g} = (g_0, g_1, g_2, g_3) = (h_3, -h_2, h_1, -h_0) or g_k = (-1)^k h_{3-k}, k=0,1,2,3, then {\bf g} is a highpass filter. Moreover, if we construct the bottom half of the wavelet matrix using {\bf g} in the same way we constructed the top half of the transformation matrix, then the wavelet matrix is orthogonal. We have, for processing length 8 vectors, the matrix
W_8 = \left[\matrix{ h_0 & h_1 & h_2 & h_3 & 0 & 0 & 0 & 0 \\ 0 & 0 & h_0 & h_1 & h_2 & h_3 & 0 & 0 \\ 0 & 0 & 0 & 0 & h_0 & h_1 & h_2 & h_3 \\ h_2 & h_3 & 0 & 0 & 0 & 0 & h_0 & h_1 \\ g_0 & g_1 & g_2 & g_3 & 0 & 0 & 0 & 0 \\ 0 & 0 & g_0 & g_1 & g_2 & g_3 & 0 & 0 \\ 0 & 0 & 0 & 0 & g_0 & g_1 & g_2 & g_3 \\ g_2 & g_3 & 0 & 0 & 0 & 0 & g_0 & g_1} \right]
= \left[\matrix{ h_0 & h_1 & h_2 & h_3 & 0 & 0 & 0 & 0 \\ 0 & 0 & h_0 & h_1 & h_2 & h_3 & 0 & 0 \\ 0 & 0 & 0 & 0 & h_0 & h_1 & h_2 & h_3 \\ h_2 & h_3 & 0 & 0 & 0 & 0 & h_0 & h_1 \\ h_3 & -h_2 & h_1 & -h_0 & 0 & 0 & 0 & 0 \\ 0 & 0 & h_3 & -h_2 & h_1 & -h_0 & 0 & 0 \\ 0 & 0 & 0 & 0 & h_3 & -h_2 & h_1 & -h_0 \\ h_1 & -h_0 & 0 & 0 & 0 & 0 & h_3 & -h_2} \right]
The overlap in the rows in each of the highpass and lowpass portions of the transform matrix helps to avoid the decoupling problem encountered by the HWT. The disadvantage to this and other orthogonal wavelet transformation is the fact that the filters are not symmetric (e.g., for length 4 filters, we would like the first and last terms to be the same and the second and third teams to be equal) and there is a "wrapping" problem with the last row in each portion of the transform matrix. Consider the matrix above - to produce the fourth average and fourth difference, we use elements 1, 2, 7, and 8 - these elements are not typically correlated.
Longer Length Daubechies Filters
We can set up and solve similar systems to obtain Daubechies lowpass filters of length 6, 8, 10, ... The quadratic equations in the system are obtained by writing out the top portion of the matrix and then identifying the orthogonality conditions. If the filter length is M, then there are M/2 orthogonality conditions. The other equations come from the Fourier series: H(0) = \sqrt{2} and H^{(k)}(\pi)=0, k=0,\ldots,M/2 -1. For example, when M=6, the system is
\begin{eqnarray}
h_0 + h_1 + h_2 + h_3 + h_4 + h_5 &=& \sqrt{2} \\
h_0 - h_1 + h_2 - h_3 - h_4 + h_5 &=& 0 \\
h_1 -2h_2 + 3h_3 -4h_4 + 5h_5 &=& 0 \\
h_1 -4h_2 + 9h_3 -16h_4 + 25h_5 &=& 0 \\
h_0^2 + h_1^2 + h_2^2 + h_3^2 + h_4^2 + h_5^2 &=& 1 \\
h_0 h_2 + h_1 h_3 + h_2 h_4 + h_3 h_5 &=& 0 \\
h_0 h_4 + h_1 h_5 &=& 0
\end{eqnarray}
A solution to this system (the reflection is the other solution) is
\begin{eqnarray}
h_0 &\approx& 0.332671 \\
h_1 &\approx& 0.806892 \\
h_2 &\approx& 0.459878 \\
h_3 &\approx& -0.135011 \\
h_4 &\approx& -0.0854413 \\
h_5 &\approx& 0.0352263
\end{eqnarray}
Biorthogonal Wavelet Filters
The Daubechies orthogonal family of filters and other orthogonal filters are an improvement over the Haar filter in many ways - the decoupling problem can be avoided and it turns out that the longer filters provide better approximation properties than shorter ones. There are two problems with these filters when we consider them for use in image compression - they are not symmetric and all values are irrational. They do possess some nice properties - they generate an orthogonal transformation, they are finite length, and they possess good approximation properties. But Daubechies was able to show that no orthogonal filter, other than the Haar filter, satisfies these three properties AND is symmetric. In order to construct a symmetric filter, we must give up one of the other properties. Since finite length and good approximation properties are very important, orthogonality is conceded.
Relinquishing orthogonality is not a huge concession - the advantage of an orthogonal transformation is the inverse is the transpose. But we can also think here that the inverse is the transpose of a wavelet matrix - it just happens to be the same! Why not build TWO wavelet matrices \tilde{W} and W where the inverse of one is the transpose of another? So we seek wavelet matrices \tilde{W}, W such that \tilde{W}^{-1} = W^T. By the term wavelet matrix, we mean a matrix that is built from a lowpass filter and a highpass filter and whose row structures behave exactly the same as those W_8 described above for the Daubechies length four filter.
So we seek TWO lowpass filters \tilde{{\bf h}}, {\bf h} and since we are giving up orthogonality, we will insist that both filters be symmetric. It turns out that the highpass filters \tilde{{\bf g}},{\bf g} can be constructed from {\bf h}, \tilde{{\bf h}}, respectively. But how to find these filters? The system for the orthogonal filters was quadratic and increased in complexity as the length of the filter grew. The solution is simple - just pick one of the filters to be a symmetric lowpass filter and then use the fact that the top half of each wavelet matrix \tilde{H}, H, must satisfy \tilde{H} H^T = I, where I is the identity matrix.
Where do we find a symmetric lowpass filter? It turns out one source is very familiar to almost all math students! If you look at the rows of Pascal's triangle, you will see that they are symmetric and the elements in row k sum to 2^k. So to form a lowpass symmetric filter, we take row k from Pascal's triangle, divide each element by 2^k, and then multiply each element by \sqrt{2}! As an example, consider the second row of Pascal's triangle. The values are 1, 2, 1. If we divide each element by 4 and multiply by \sqrt{2}, we obtain the filter
\tilde{{\bf h}} = (\tilde{h}_{-1}, \tilde{h}_0, \tilde{h}_1) = (\sqrt{2}/4, \sqrt{2}/2, \sqrt{2}/4)
These elements sum to \sqrt{2} and they alternate sum to 0 - a lowpass filter! For the other filter, let's try a length 5 filter (there are results that tell us what lengths will work with length 3 filters). We take {\bf h} = (h_2, h_1, h_0, h_1, h_2) to be symmetric and set up the orthogonality conditions - \tilde{{\bf h}} and {\bf h} must dot with each other and give one and two-shifts of {\bf h} must dot with \tilde{{\bf h}} and give 0. We obtain the system
\begin{eqnarray}
\frac{\sqrt{2}}{4} h_1 + \frac{\sqrt{2}}{2} h_0 + \frac{\sqrt{2}}{4} h_1 &=& 1 \\
\frac{\sqrt{2}}{4} h_1 + \frac{\sqrt{2}}{2} h_2 &=& 0
\end{eqnarray}
We add the lowpass condition H(\pi) = 0 or h_2 - h_1 + h_0 - h_1 + h_2 = 0 and simplify to obtain the system
\begin{eqnarray}
h_0 + h_1 &=& \sqrt{2} \\
h_1 + 2h_2 &=& 0 \\
h_0 -2h_1 + 2h_2 &=& 0
\end{eqnarray}
It is easy to solve this linear system and find the solution
{\bf h} = (h_2, h_1, h_0, h_1, h_2) = (-\sqrt{2}/8, \sqrt{2}/4, 3\sqrt{2}/4, \sqrt{2}/4, -\sqrt{2}/8)
You can also check that h_2 + h_1 + h_0 + h_1 + h_2 = \sqrt{2}. For the wavelet matrix \tilde{W}_8, we form {\tilde{\bf g}} from {\bf h} by the formula \tilde{g}_k = (-1)^k h_{1-k}, k=-1,0,1,2,3. For the wavelet matrix W_8, we form {\bf g} from \tilde{{\bf h}} by the formula g_k = (-1)^k \tilde{h}_{1-k}, k=-1,0,1. Here are the wavelet matrices:
\tilde{W}_8 = \left[\matrix{ \tilde{h}_0 & \tilde{h}_1 & 0 & 0 & 0 & 0 & 0 & \tilde{h}_1 \\ 0 & \tilde{h}_1 & \tilde{h}_0 & \tilde{h}_1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & \tilde{h}_1 & \tilde{h}_0 & \tilde{h}_1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & \tilde{h}_1 & \tilde{h}_0 & \tilde{h}_1 \\ \tilde{g}_0 & \tilde{g}_1 & \tilde{g}_2 & \tilde{g}_3 & 0 & 0 & 0 & \tilde{g}_{-1} \\ 0 & \tilde{g}_{-1} & \tilde{g}_0 & \tilde{g}_1 & \tilde{g}_2 & \tilde{g}_3 & 0 & 0 \\ 0 & 0 & 0 & \tilde{g}_{-1} & \tilde{g}_0 & \tilde{g}_1 & \tilde{g}_2 & \tilde{g}_3 \\ \tilde{g}_2 & \tilde{g}_3 & 0 & 0 & 0 & \tilde{g}_{-1} & \tilde{g}_0 & \tilde{g}_1 } \right]
and
W_8 = \left[\matrix{ h_0 & h_1 & h_2 & 0 & 0 & 0 & h_2 & h_1 \\ h_2 & h_1 & h_0 & h_1 & h_2 & 0 & 0 & 0 \\ 0 & 0 h_2 & h_1 & h_0 & h_1 & h_2 & 0 \\ h_2 & 0& 0 & 0 h_2 & h_1 & h_0 & h_1 \\ g_0 & g_1 & 0 & 0 & 0 & 0 & 0 & g_{-1} \\ 0 & g_{-1} & g_0 & g_1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & g_{-1} & g_0 & g_1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & g_{-1} & g_0 & g_1} \right]
Other biorthogonal filter pairs can be constructed in a similar fashion. Daubechies also has a result that gives an explicit formula for H(\omega) once we have selected \tilde{H}(\omega). Please see Discrete Wavelet Transformations: An Elementary Approach with Applications for more details.
At first glance, it doesn't seem that we have solved the "wrapping row" issue that can cause problems for the first and/or last elements computed in each portion of the transformation. The elements of both matrices above are irrational - hardly good candidates for mapping integers to integers! In the subsection Wavelet Transformations Used by JPEG2000, we will discuss how the symmetry of the filters can be exploited to eliminate the wrapping rows and how to modify the above filter pair to produce one that maps integers to integers.
