Image Compression with Wavelets

Steven is a Ph.D. student in the computer science and operations research department of the University of Montreal in Canada. He can be contacted at pigeon@ iro.umontreal.ca.

"If you don't need it, don't keep it." This mantra is the key not only to clean offices and efficient memory usage, but also to modern graphics and sound compression. Differences that are imperceptible to the average human are simply not worth keeping.

The most popular approach today is to separate your data by "frequency." Generally, high-frequency data is less perceptible, and needn't be stored as accurately. Wavelets present the possibility of separating data according to "feature size," an alternative that many people find much more intuitive.

-- Tim Kientzle

Modern image compression schemes rely heavily on frequency decomposition. The eye isn't really that sensitive to weak, fine-grain, high-frequency detail; it filters most of it out as "noise." This is convenient, because noise is also hard to compress. Simply removing the fine details would make it easier to compress. Because the eye ignores these details anyway, removing them doesn't degrade the picture quality.

One way to proceed is to consider the image (or a small part thereof) as a "time signal." A time signal is a function that changes according to time. So pixels number 0,1,2,...,t come from some function x(t), the time function. In terms of pixels, you might as well call it the space function that varies according to the position of the pixel. Time signals can be decomposed into a series of frequency functions, f(i), that adds together to give x(t) back. Those frequency functions form the frequency domain.

So, to remove all those fine-grain, hard-to-code, imperceptible details, you first transform a part of the picture from the time domain (x(t)) to the frequency domain (f(i)). You can then eliminate the frequencies that correspond to noise. The result is a collection of frequency components ("coefficients") that can be efficiently compressed using some standard technique, such as Huffman, Adaptive Huffman, or arithmetic coding. To recover the image, you decode the compressed image to recover the coefficients f(i), and use those to compute an approximate reconstruction. JPEG basically works this way, and it produces significantly compressed image files.

The Fourier Expansion

The most common frequency transforms are the Fourier transform and the Fourier expansion. The Fourier expansion of a function is a sum of sines and cosines. Each sine and cosine has a distinct frequency and a distinct amplitude. So, for all sufficiently smooth functions, you have {a₀,a₁,a₂,...}, {b₁, b₂, b₃,...}, such that you get the equation in Figure 1.

Every sufficiently smooth function can be expressed by such a sum. If you use a limited number of coefficients, the reconstructed function is only an approximation to f(x). As the number of terms goes to infinity, the sum gives a perfect reconstruction of f(x).

The Fourier transform works well on a wide variety of functions but has serious drawbacks. First, it has problems with discontinuities -- points where the function jumps suddenly from one value to another. These cause the Fourier transform to "ring" as it tries to approximate the discontinuity. This is known as the "Gibbs phenomenon"; see Figure 2.

Second, the Fourier decomposition has trouble with transients. Transients are local features; for example, a sharp peak in an otherwise smooth area. This causes Fourier transforms to have far more as and bs; you will need more terms to get a good approximation than if there were no transient. In terms of compression, you won't be able to throw out as many coefficients, unless you are willing to accept a bad reconstruction of the transient. This can cause you to lose important details if you're doing image compression.

The Discrete Cosine Transform

JPEG uses the Discrete Cosine Transform (DCT) instead of the Fourier transform, but it is limited in the same ways. So, if a patch of the image presents sharp details, it costs a lot of coefficients to represent them with sufficient accuracy. Using fewer coefficients will only make the ringing more visible. This phenomenon (Figure 3) has been referred to as a "mosquito artifact."

Wavelets

In contrast, wavelets handle transients and discontinuities rather well. Instead of being constructed out of sines and cosines that span the axis (or plane, in 2D), wavelets are constructed from compactly supported functions. A compactly supported function is one that takes nonzero values only inside a given interval, and is zero outside this interval. So, if a function exhibits a transient, you need to add more wavelets only around the transient. The overall effect is that smooth regions of the function require fewer wavelets than detailed regions to be satisfyingly approximated, and the wavelets used to reconstruct the transient do not affect the smooth regions.

The thing to remember regarding the difference between trigonometric transforms (based on sines, cosines, or both) and wavelets is that trigonometric transforms are global and wavelets are local.

So, What Does a Wavelet Look Like?

The simplest wavelets are Haar wavelets. Although A. Haar invented wavelets around 1909, the term was coined much later. He intended his "wavelets" to be a tool for progressive reconstruction, as they are today, but in a context of function analysis rather than signal processing, since signal processing hadn't yet been invented. Haar's idea was that a function could be approximated at any point with a sum of wavelets overlapping this point. The strength of his approach is that any function can be perfectly reconstructed given enough coefficients (possibly infinitely many), which is not true of the Fourier series or DCT.

While the Fourier expansion is for continuous functions (as were the original Haar wavelets), I will consider only discrete wavelets, as those are defined only on an integer number of points.

Haar Wavelets

The Haar wavelets are built from the "mother wavelet" in Figure 4. Starting from the mother wavelet, you form a wavelet basis by scaling the mother wavelet to different sizes and translating it to different positions. This gives you one full-size wavelet, two half-size, four quarter-size, and so on.

Mathematically, each of these scaled and translated copies specifies a row in a matrix; see Figure 5. After normalizing (each row becomes a vector with length 1), you obtain the basis matrix in Figure 6. A basis matrix lets you decompose any vector as a linear combination of its rows.

This is the classical form of the Haar wavelet basis. The inverse transform is given by the inverse of Figure 6. Numerical inversion of the matrix gives the inverse transform, but it will turn out to be only the transpose of the forward transform matrix (H^-1==H^T).

Fast Transforms

One way to compute the Haar transform is to generate the n×n matrix and then compute the transform by an ordinary matrix-vector product.

If you compute the wavelet transform as a product of a matrix and a vector, you obtain an O(n²) algorithm, which is clearly not very good. You can also see that there's only about O(n lg n) nonzero terms in the matrix. So you ought to be able to find a transform that is at most O(n lg n). If you are clever, you can even reduce this to exactly O(n) operations for the Haar transform.

To obtain a fast transform, you either use a recurrence relation or a factoring of the transform matrix. Both are easy to find by using a tool such as Mathematica to generate the transform matrix and do some symbolic manipulations. Try to see what H {a,b,c,d,e,f,g,h} looks like for unspecified a,b,c,...h. If it is not too big, you can find a pattern that is useful in devising an efficient algorithm. (Available electronically; see "Resource Center," page 5, is a Mathematica notebook you can experiment with.) Start by noting that the denominators are all powers of --2.

Figure 7 is the final result, which shows a fast Haar transform I found using the factoring approach. I checked by hand that this works for matrices of size 2, 4, and 8. It also works for all powers of 2. The inverse transform is basically the same except that the flow is reversed.

I've also included C++ implementations of the transforms (available electronically). These functions don't require n to be a power of two. In real life, the vector you want to transform can be of any length. Pictures are hardly ever exactly 512×512.

The 2D Wavelet Transform

The matrix H only works with one- dimensional vectors. Images are not vectors but matrices. To Haar transform a 2D matrix, you need a 4D matrix (a tensor). But instead of computing a 4D tensor/2D matrix product, you use a separable transform. You can apply the Haar transform to each row of the matrix, which yields a series of Haar transformed vectors. If you apply the Haar transform again to the columns of this row-transformed matrix, you get a complete 2D (separable) Haar transform.

Wavelet Compression

Now that you are able to transform a piece of the image from the time domain to the wavelet/frequency domain, how do you get any compression at all? The general scheme is as follows. Once you have transformed the patch, you get n×n coefficients. These are quantized and coded with an efficient coder. This is usually implemented by reducing the precision of the coefficients in some way that takes into account the perceptual importance of each of the coefficients: Recall that the weak high frequencies are less important. The greater the reduction in precision, the greater the compression ratio, but, alas, the worse the image gets. The quantization and coding steps are what makes all the difference in such a compression scheme. To recover the image, the codes are decoded to give back the quantized coefficients, which are then unquantized, or bit-augmented, back to the original number of bits. The inverse transform is applied to the recovered coefficients, which will yield the approximation of the original patch.

Wavelet Transforms and Signal Processing

Another clear advantage of wavelets over transforms such as the Fourier or discrete cosine transforms is that it is easier to match coefficients in the wavelet transform to specific picture features. JPEG quantization techniques rely on average power per coefficient (a quantity that takes into account both the frequency and amplitude associated with the coefficient) to quantize the coefficients, rather than the precise contribution of coefficients to the picture. On the other hand, if you look at the wavelet transform of the image, you can see -- or at least guess -- what the original image looked like. Typically, each time a region doubles in size, the corresponding frequency doubles. The largest regions, therefore, usually contain the weakest high-frequency components; see Figure 8.

If you look at the DCT coefficients of an image, you really can't tell anything by eye. Ultimately, the ability to match wavelet coefficients directly to picture features may lead to better image compression schemes, or even schemes that are tailored to only preserve certain features. For example, you might only be interested in preserving edges.

Comparing Wavelets to JPEG

To see how this compares to JPEG, look carefully at Figure 9. Figure 9(a) is the original, 9(b) the JPEG version, and 9(c) is compressed with a Haar wavelet compression system. The compression parameters were chosen so that Figures 9(b) and 9(c) are about the same size after compression.

Shape Is Important

Different applications require different wavelet shapes. A square wavelet introduces square artifacts in the approximated reconstruction. Regions that should be a smooth degradation of color show boxes, see Figure 9(c), when the compression ratio is high. While this isn't too bad in an image, it is disastrous on sound. The sharp edges of the boxes produce a horrifying high-pitched hissing sound. For sound, you would want a basis with smooth sine-like wavelets that imitate the natural characteristics of sound and music.

You must therefore carefully choose a wavelet basis that best fits your needs. There is a lot of active research within the wavelet community to find the optimal wavelet bases for different classes of data.

Conclusion

There are numerous other wavelet transforms. There are other square wavelet transforms, such as the Mallat transform (C++ code that implements this is shown in Listing One). There are also other wavelets that exhibit really curious shapes -- like Daubechies wavelets -- and some are not quite reversible. For more information on wavelets, good places to start are at http://www.wavelet.org/ and http://www .iro.umontreal.ca/~pigeon.

DDJ

Listing One

#include <imath.h>
#include <stdlib.h>
#include <math.h>
/***************************************************************************/
/* Those two routines implement the Mallat decimation algorithm. First     */
/* phase computes the means, second phase the differences. All isO(2n),    */
/* which is rather efficient.                                              */
/* keeps ~20 bits of precision (out of 23), or a precision of about 1e07   */
/***************************************************************************/
void FMT(const float source[], float destination[], int length)
 { 
  float *t;
  int i,d,j;
  t = (float *)malloc(2*length*sizeof(float));
  for (i=0;i<length;i++) t[i]=source[i];
  for (i=0,d=length;i<2*length-2;i+=2,d++) t[d] = t[i]+t[i+1];
  destination[length-1]=t[2*length-2];
  for ( i=2*length-2, j=2*length-3, d=length-2; 
        i>length-1; 
        i--, j-=2, d--) 
    destination[d] = t[j]-t[i]/2.0f;
  free(t);
 }
/*******************************/
void iFMT(const float source[], float destination[], int length)
 {
  float *t; 
  int i,d,p;
  t = (float *)malloc(2*length*sizeof(float));
  t[2*length-2]=source[length-1];
  for ( p=length-2, i=2*length-4 , d=2*length-2; 
       i>=0; i-=2, d--, p--)
   {
    float t0 = t[d]/2;
    float s0 = source[p];
    t[i]   = t0 - s0;
    t[i+1] = t0 + s0;
   }
  for (i=0;i<length;i++) destination[i]=t[i];
  free(t);
 }
/***************************************************************************/
/* These routines compute the fast (normalized) Haar transform             */
/* Works of course with a vector of any lenght > 3.                        */
/* Lonely coefficients are left as is, since we can't pair them correctly, */
/* but are put in the result vector. The precision is also about ~20 bits  */
/* (out of 23)                                                             */
/***************************************************************************/

void FHT(const float source[], float destination[], int length)
 {

  int i,l;
  float *t = (float *)malloc(length*sizeof(float));
  for (i=0;i<length;i++) destination[i]=source[i];
  l=length;
  for (i=Log2(length-1);i;i--)
   {
    int j,k;
    int l2 = l/2;
    for(j=0;j<l;j++) t[j]=destination[j];
    for (k=0,j=0;k<l2;k++,j+=2)
     {
      destination[k]    = (t[j]+t[j+1])*0.7071067812f;
      destination[k+l2] = (t[j]-t[j+1])*0.7071067812f;
     }
    l /= 2 ;  // for the leftovers
   }
  free(t);
 }
/*******************************/
void iFHT(const float source[], float destination[], int length)
 {
  int i,j;
  int lim[32];
  
  int    k  = Log2(length-1);
  float *t  = (float *)malloc( length * sizeof(float));

  if (( 1 << k) != length) k--; 

  for ( j=length, i=k; i; i--, j/=2 ) lim[i]=j;
  lim[0]=1;

  for ( i=0; i<length; i++ ) destination[i]=source[i];

  for(i=0;i<k;i++)
   {
    int d;
    int l=lim[i];
    int l2=lim[i+1];

    for(j=0; j<l2; j++) t[j]=destination[j];

    for( d=0, j=0; j<l; j++, d+=2 )
     {
      destination[d]   =(t[j]+t[j+l])*0.7071067812f;
      destination[d+1] =(t[j]-t[j+l])*0.7071067812f;
     }
   }
  free(t);
 }

Back to Article