15-859E Programming Assignment 3:
Wavelet Image Compression

Revision 3: 30 Oct. 98.

Earlier, I said that the assignment would be variational surface modeling, but that appeared too involved for the short time we have, so I settled on this simpler assignment.

• Implement (a) wavelet decomposition (analysis) of images, (b) compression of the wavelet coefficients by quantization and zeroing of small coefficients, and (c) reconstruction (synthesis) of images from the compressed wavelet representation.
• As usual, you can implement in any computer language. I think it wouldn't be too hard to do this assignment in matlab, for example.
• Your program should read and write picture files. At minimum, your program should read pictures with 8 bits of grayscale per pixel.
• The wavelet decomposition should not be computed using integer arithmetic (unless you really know what you're doing), but with floating point. Thus, if your picture is n by n pixels, the input representation will be n*n*8 bits and after decomposition into the wavelet basis, it will use n*n*B bits, where B=32 or 64 depending on whether you use "float" or "double". This is not compression, yet! Let's call this data structure the "raw wavelet representation".
• Next, your program should do compression. The first form of compression your program should have is the zeroing of small wavelet coefficients, as described by Stollnitz. I suggest this be controlled with a parameter frac, so that the fraction frac of the coefficients are kept, and the rest are zeroed. This would compress to about n*n*frac*B bits.
• As a second form of compression, your program should have an option for quantization (discretization of a scalar quantity to fewer bits) of each coefficient. This option should be independent of zeroing. Quantization is not in Stollnitz, but would be vital in any practical compression algorithm, so we'll do it to make the exercise more realistic. When this option is on, it would quantize each coefficient to b bits. When off, we'll define b=B. A simple formula for this is suggested under "Tips", below.
• The "compressed wavelet representation" we have at this point will use n*n*frac*b bits. (If b<8 or frac is small then we're actually compressing, else not!)
• Next in the program, reconstruct the image from the compressed wavelet representation. Do this in floating point, also. If any of the reconstructed pixel values come out<0 then set them to 0 or if >255 then set them to 255. Now we've got an 8-bit image again. Write it out to a file.
• Measure its error. Compute the relative L2 error of the output image "out" relative to the input image "in": sqrt(sum over x,y of (in[x,y]-out[x,y])^2) / sqrt(sum over x,y of in[x,y]^2)). This is the reciprocal of the "signal to noise ratio" that is popular in signal processing. We need to compute the error (whereas Stollnitz had it as a parameter) because we're sometimes doing the extra quantization step. If the fraction was large and (quantization was not turned on or b was large), then the picture should be nearly indistinguishable from the original.
• Run your program on a picture provided by me, src/wavelet/pix/roar512.pgm, and a picture of your own choosing. The "roar" picture and another picture are available in PGM format at 16x16, 128x128, and 512x512 resolutions. You might want to do initial debugging on roar16.pgm, and then final testing on roar512.pgm. This and another test picture are in the src/wavelet/pix directory. You can view the pictures here.
• Make a table listing

  picture-name   npixels   frac   float/b  cf    rel.err   time(sec.)
      

  where "float/b" means whether coefficients are float, double, or quantized, and if quantized, what was the value of b; and "cf" means compression factor (frac*b/8). If cf>1 then we have expansion, not compression. For example, five different compression runs on the same picture might be summarized thus:

  face.pgm       512x512   1      float    32    0         20s
  face.pgm       512x512   .1     float    3.2   .1        20s
  face.pgm       512x512   1      32       32    .05       20s
  face.pgm       512x512   1      8        1     .15       20s
  face.pgm       512x512   .01    4        .005  .2        20s
      

  Note that the two right-hand columns in this example are wild guesses. Run each picture at least four times with results ranging from the most compression (smallest cf) that looked indistinguishable from the original down to a highly compressed version with relative error of .2 or so. See Stollnitz figure 3.5 for an example.
• Print out some of the pictures you generate. For each image you work with, print the input image and two representative output images, say at the useful extremes of low compression and high compression. Label each picture with its table statistics.
• Finish your code by the date listed at top. To submit your program, copy your source code, executable, and input picture files to students/yourlastname/wavelet. Include a README file explaining: what you did, what system your executable runs on (e.g. 900 MHz Pentium 4). Also include the comparative table mentioned earlier.
• During class on 3 Nov., turn in a printout of your README and the printed pictures. We'll discuss the outcome in class, but we probably won't do demos, since most of these programs will be non-interactive.

Tips

• I bet float will give sufficient precision, and that doubles won't be necessary.
• Make your pictures square and powers of two in width & height. 512x512 pixels is a good size. You can use the "xv" program or "Photoshop" to resample an image to any size, or to convert it to grayscale.
• A very simple picture file format I recommend is PGM. I have code for reading/writing PGM files that you can see in pic.c, pic.h, and pnm.c. Or more interesting, inout.cc is a main program that uses it to read a picture into an array of doubles. To convert a picture from TIFF format to PGM, try "tifftopnm". The "convert" program is another good picture file format conversion program. There's also a library called libpgm for reading/writing this format. Try "rtfm libpgm" on Andrew. The PGM file format is documented in "rtfm pgm". For more info on PGM, see the software page.
• You could use "xv", "photoshop", or "pnmtops" to print your pictures. Photoshop doesn't understand PGM file format, however.
• I suggest you put picture file names, frac, and b on the command line, to avoid recompiling for each run. The argument parser discussed in the multigrid assignment would be handy for that.
• There won't be any starter code, per se, beyond that mentioned above.
• Here's a simple quantization formula for converting coefficient c[i,j] into quantized coeff qc[i,j] and then back to recovered coefficient rc[i,j]. It simply maps the interval [-cmax,cmax] to b bits, temporarily stores it into an integer variable, and then converts back.

  double c[i,j]
  #define TINY .001
  // cmax is max of |c[i,j]| for all i,j
  //   (probably the coeff of the level 0 scaling fn, I'm guessing)
  double scalefactor = (pow(2,b-1)-.5-TINY)/cmax
  int qc[i,j] = floor(scalefactor*c[i,j] + .5)
      // qc is an integer in [-2^(b-1)+1,2^(b-1)-1]
  double rc[i,j] = qc[i,j]/scalefactor
      // rc is approximately equal to c
      // the more bits you use (large b) the closer it is
      

• Note that even if you set frac=1, quantization will probably cause some of the smaller coefficients to be zeroed.
• You could also calculate the RMS error in the "wavelet domain", as Stollnitz did. That might work out to be faster.
• Please put some effort into optimization. In your inner loops, instead of dividing by "sqrt(2)", multiply by a constant equal to sqrt(1/2). The header file /usr/include/math.h has one. I hope I don't see any code that looks like:

      c2[pow(2,j)/2+i] = (c[2*i-1]-c[2*i])/sqrt(2);
      

  Use "1<<(j-1)" instead of "pow(2,j)/2", or better yet precompute that outside the loop! The fastest implementation would use pointers.
• For the ambitious, after you get the above done, some extra credit options:
  • Color images
  • Interactive program: put frac, b on sliders and display the output picture, instead of writing it to a file. Optimize the heck out of it.
  • Do compression for real. Write out a compressed image file and try to make it as small as possible without sacrificing image quality. You might want to use Huffman coding or arithmetic coding on the quantized values.
  • Try other wavelet bases. Haar is a primitive wavelet, and human vision is very sensitive to its step discontinuities. Stollnitz page 197 talks about some options.

15-859, Hierarchical Methods for Simulation
Paul Heckbert