Tuffy, the thirsty dog A friendly pocket monster
CSC 726, Parallel Algorithms
QR Factorization

Recall that our linear system is:

g = H f

where f is a vector size n*n by 1, H is an matrix size 2*n*n by n*n, and g is a vector size 2*n*n by 1. Note that the system is over determined, so you will need to compute least-squares solution. QR factoriziation is an important method for least squares problems such as this one. A summary of the method and how it applies to least squares problems is linked here.

The vector f contains the pixel values of an n by n, image stacked by columns. Similarly, g contains the stacked pixel values for two images. The vector g represents two observations, and the matrix H represents the blurring operator for both observations. The inverse problem posed here is to find f. For the following files, use your code to compute an estimate of f using g and H. In the example files below, n = 16. For purposes of code validation, the original vector f is also provided. Your solution should be very close to the given f. Later, you will have the opportunity to solve problems for which f is not available.

The following ".dat" files contain little-endian binary data. Note on file format: These are binary files. The first 4 bytes form a little-endian integer M, representing the number of rows in the object. The second 4 bytes form a little-endian integer N, representing the number of columns in the object. The next M * N * sizeof(double) bytes are little-endian double precision values in column-major order. A sample seqential C function to read this format is included here. Sample sequential C++ code is here. Sample C code to write this format is included here. Sample C++ code to write this format is here.

Also, you will want to convert matricies to a graphics format for viewing. A sequential program named dat2ppm to convert to the PPM graphics format is compiled and (currently) installed in /usr/local/bin on gottlieb.cs.wfu.edu. The program dat2ppm takes two options for controlling the colormap: dat2ppm -gray filename.dat produces a gray-scale image; dat2ppm -jet filename.dat produces a pseudo-color image using a colormap similar to Matlab's "jet" colormap.

You can view PPM graphics files on gottlieb with the display command from the ImageMagick software package (already installed on gottlieb).

Sample Input Data for QR Factorization
A small problem for getting things tested and working.

• g_le.dat , a simulated observation (512 x 1)
• h_le.dat , a simulated blurring operator (512 x 256)
• f_le.dat , the unknown. (256 x 1)

What do these images look like ? The vector f is un-stacked into a square 16 x 16 image. The vector g is split into the top and bottom regions g1 and g2, and each is un-stacked into a square 16 x 16 image. For illustration purposes, each image is enlarged (via block replication). Then, we have:

 f g1 g2

The matrix h is illustrated in the 'jet' colormap below. Blues and cool tones indicate low values; oranges, reds and warm tones indicate high values. The top and bottom halves of h are structured: they are block circulant with circulant blocks.

 Matrix h
A somewhat larger problem.
Now that your code is working smoothly, try a slightly larger problem, based on an image size 32 x 32.

For the 2048 x 1024 matrix "h_32.dat" above, what does the QR factorization look like ?

A more interesting problem.
You can also try a larger problem, based on an image size 64 x 64.

An ill-conditioned, spatially varying problem.
Finally, we have arrived at our first real problem. The following files present an ill-conditioned simulated spatially-varying blur on 64x64 images. You will need to include regularization to produce an acceptable restoration.

Sample C code to read the data file and to assist with Tikhonov regularization is found in tikhonov.h and in tikhonov.c. These functions should also compile easily if you are using the C++ compiler, with a minor modification to the syntax for the #include files.

Images:
 g_0 g_1 f